Example code — Certified selective interventions and closed-RAF structure in a pooled prokaryotic network

Finding a self-supporting reaction set does not necessarily identify the reactions needed to make a particular metabolite. In the paper's ten-reaction valine subsystem, two singleton reactions are the complete irreducible RAF catalogue, yet all ten reactions are necessary to produce valine. The example reproduces this distinction and the selective intervention in the supplied pooled prokaryotic reaction dataset.

The package includes the pinned reaction data, its 68-entry food medium, positive catalyst formulas and grouped reversible actions. Reusable classes compute food closure, maxRAF, target capability, small exact catalogues, closed extensions and intervention certificates. Inputs identify the target, cut and food additions at the top of the code.

Constructive activation reaches the same parent and postcut reaction sets as maxRAF pruning: 2148 and 2085 directions respectively. — Constructive lower bounds and RAF upper bounds meet as reaction sets in the actual parent source. The code checks set equality; matching counts alone would not certify uniqueness.

A reaction-membership grid contrasts the two singleton seeds, the five-reaction closed subset, the full ten-reaction target support and the eleven-reaction pooling extension. — Minimal catalytic seeds, target support and ambient closedness carry different information. Indices refer to the literal source directions listed in the downloadable CSV.

Deleting R01209, R01210 and R04441 removes four directed reactions. An 18-species group containing valine becomes inaccessible (a siphon): all 58 surviving directions that produce a member also require a member as a reactant, and none is food. This obstruction does not depend on catalyst annotations.

Restoring the deleted actions individually gives valine-producing RAF witnesses of 9, 13 and 13 directions. Together these certify inclusion-minimality, while the food-supported singleton R09639 forward preserves methionine. This is not a claim that the cut has globally minimum cost. Adding precursor C00141 to the food restores valine through a four-direction witness, showing the environmental boundary explicitly.

The ten-reaction subsystem has 22 RAFs and two closed RAFs: $L=\{0,1,4,7,8\}$ and the full set $H$ . Each already contains every further reaction enabled by its available molecules. Adding one food-enabled pooling direction supplies NADs and leaves one closed RAF containing all eleven directions. Closedness is relative to the ambient reaction set; the formal pooling row represents source catalyst availability.

In the parent source, constructive activation and maxRAF pruning meet at the same set of 2,148 directions, or 2,085 after intervention, with 23 activation layers in each case. Equality is checked by reaction identifiers. These matching lower and upper bounds certify a unique closed RAF without enumerating every parent irreducible RAF.

The implementation keeps two semantics separate: RAF catalysts are checked in the final reactant closure, while constructive activation requires them in the current prefix pool. Seven scientific test groups check these semantics, the source, restoration ranks, subcuts, catalogue counts, precursor rescue and parent set equality.

The results concern structural generation in a specified medium. They do not establish depletion, growth, flux or a unique concentration equilibrium. The package recomputes finite source certificates; it does not rerun Lean.

Python source

"""Selective interventions, catalytic seeds and ambient closedness on literal data.

Run: python example.py --output outputs. Presence/absence structure, not kinetics.
"""
# EDITABLE INPUTS. The source medium has 68 entries; additions change that medium.
TARGET = 'C00183'                 # L-valine
PRESERVATION_TARGET = 'C00073'     # L-methionine
DELETED_ACTIONS = ('R01209','R01210','R04441')  # removes all directions of each row
EXTRA_FOOD = ()
RESCUE_FOOD = 'C00141'             # demonstrated precursor rescue, not any barrier member
SMALL_DIRECTIONS = ('R00014::f','R04672::r','R10916::r','R10985::f','R00585::f',
                    'R01215::r','R04441::f','R00221::r','R00589::r','R_NADs_2::f')
POOLING_DIRECTION = 'R_NADs_1::f'
MAX_SUBSETS = 65536                # exact enumeration is only for small ambient sets
MANUSCRIPT_SHA256 = 'd271244c721f70339c8e26094d51d691ef1915e556e742e3dc0648a64c0aa50c'
PINNED_SOURCE_SHA256 = 'f6834863b2d272d074af5befcb6149c11edf5a64b1aa0ba716ebeb8f10b3700b'

import argparse
import csv
from dataclasses import dataclass, replace
import hashlib
import itertools
import json
from pathlib import Path
import platform
import re


@dataclass(frozen=True)
class CatalystFormula:
    """Positive formula; conjunction is never flattened into alternative catalysts."""
    kind: str
    children: tuple = ()
    atom: str = ''

    def __post_init__(self):
        if self.kind not in ('true','atom','and','or'): raise ValueError('Positive formulas only.')
        if self.kind=='atom' and (not self.atom or self.children): raise ValueError('Malformed atom.')
        if self.kind in ('and','or') and len(self.children)!=2: raise ValueError('Binary connective required.')
        if self.kind=='true' and (self.children or self.atom): raise ValueError('Malformed true constant.')
        if self.kind in ('and','or') and self.atom: raise ValueError('Connectives have no atom label.')

    def holds(self,pool):
        if self.kind=='true': return True
        if self.kind=='atom': return self.atom in pool
        if self.kind=='and': return all(c.holds(pool) for c in self.children)
        return any(c.holds(pool) for c in self.children)

    def mandatory_atoms(self):
        if self.kind=='true': return frozenset()
        if self.kind=='atom': return frozenset((self.atom,))
        a,b=(c.mandatory_atoms() for c in self.children)
        return a|b if self.kind=='and' else a&b

    def to_data(self):
        return [self.kind,self.atom] if self.kind=='atom' else [self.kind,*[c.to_data() for c in self.children]]

    @classmethod
    def parse(cls,text):
        text=text.strip().replace('|',',').replace('*','&')
        text=re.sub(r'\s*([(),&])\s*',r'\1',text)
        text=re.sub(r'\s+',',',text)
        if not text: return cls('true')
        tokens=re.findall(r'[(),&]|[^(),&]+',text); pos=0
        def term():
            nonlocal pos
            if pos>=len(tokens): raise ValueError('Incomplete catalyst formula.')
            token=tokens[pos]; pos+=1
            if token=='(':
                result=disjunction()
                if pos>=len(tokens) or tokens[pos]!=')': raise ValueError('Missing closing parenthesis.')
                pos+=1; return result
            if token in ',&)': raise ValueError('Expected catalyst atom.')
            return cls('atom',atom=token)
        def conjunction():
            nonlocal pos
            result=term()
            while pos<len(tokens) and tokens[pos]=='&':
                pos+=1; result=cls('and',(result,term()))
            return result
        def disjunction():
            nonlocal pos
            result=conjunction()
            while pos<len(tokens) and tokens[pos]==',':
                pos+=1; result=cls('or',(result,conjunction()))
            return result
        result=disjunction()
        if pos!=len(tokens): raise ValueError('Unconsumed catalyst tokens.')
        return result


@dataclass(frozen=True)
class Reaction:
    id: str
    action: str
    reactants: tuple             # (species, integer coefficient), retained for reuse
    products: tuple
    catalyst: CatalystFormula
    source_line: int = 0

    def __post_init__(self):
        if not self.id or not self.action: raise ValueError('Direction and action identifiers required.')
        if any(not s or not isinstance(n,int) or n<=0 for s,n in self.reactants+self.products):
            raise ValueError('Positive integer stoichiometric coefficients required.')

    @property
    def inputs(self): return frozenset(s for s,n in self.reactants)
    @property
    def outputs(self): return frozenset(s for s,n in self.products)


@dataclass(frozen=True)
class ReactionSystem:
    reactions: tuple
    food: frozenset
    source_sha256: str = ''

    def __post_init__(self):
        if len({r.id for r in self.reactions})!=len(self.reactions): raise ValueError('Duplicate directed identifiers.')

    @property
    def by_id(self): return {r.id:r for r in self.reactions}

    @classmethod
    def read(cls,path,expected_hash=None):
        data=Path(path).read_bytes(); digest=hashlib.sha256(data).hexdigest()
        if expected_hash is not None and digest!=expected_hash: raise ValueError('Pinned source hash mismatch.')
        food=set(); reactions=[]; actions=set()
        def side(text):
            result=[]
            for item in text.strip().split('+'):
                parts=item.split()
                if not parts: continue
                if len(parts)==1: result.append((parts[0],1))
                elif len(parts)==2 and parts[0].isdigit() and int(parts[0])>0: result.append((parts[1],int(parts[0])))
                else: raise ValueError(f'Invalid stoichiometry: {item}')
            return tuple(result)
        for line_no,line in enumerate(data.decode('utf-8').splitlines(),1):
            if not line.strip() or line.lstrip().startswith('#'): continue
            if line.startswith('Food:'): food.update(line[5:].split()); continue
            match=re.fullmatch(r'([^:]+):\s*(.*?)\s*\[(.*?)\]\s*(<=>|=>|<=)\s*(.*)',line)
            if not match: raise ValueError(f'Unrecognized source row {line_no}.')
            action,left,cat,arrow,right=match.groups(); action=action.strip()
            if action in actions: raise ValueError('Duplicate source action.')
            actions.add(action); a,b=side(left),side(right); formula=CatalystFormula.parse(cat)
            for direction in (('f','r') if arrow=='<=>' else ('r',) if arrow=='<=' else ('f',)):
                inputs,outputs=(a,b) if direction=='f' else (b,a)
                reactions.append(Reaction(action+'::'+direction,action,inputs,outputs,formula,line_no))
        return cls(tuple(reactions),frozenset(food),digest)

    def subset(self,ids):
        ids=frozenset(ids)
        if not ids<=self.by_id.keys(): raise ValueError('Unknown direction in subset.')
        return replace(self,reactions=tuple(r for r in self.reactions if r.id in ids))

    def delete_actions(self,actions):
        actions=frozenset(actions)
        if not actions<={r.action for r in self.reactions}: raise ValueError('Unknown action in deletion.')
        return replace(self,reactions=tuple(r for r in self.reactions if r.action not in actions))

    def add_food(self,species): return replace(self,food=self.food|frozenset(species))


class RAFAnalysis:
    """Catalyst-free food closure, followed by final-pool RAF support tests."""
    def __init__(self,system):
        self.system=system; self.rows=system.by_id; self.universe=frozenset(self.rows)
        self.inputs={k:r.inputs for k,r in self.rows.items()}; self.outputs={k:r.outputs for k,r in self.rows.items()}

    def validate(self,ids):
        ids=frozenset(ids)
        if not ids<=self.universe: raise ValueError('Direction outside this ambient set.')
        return ids

    def closure(self,ids):
        remaining=set(self.validate(ids)); pool=set(self.system.food); ranks={}; layer=0
        while True:
            ready=sorted(k for k in remaining if self.inputs[k]<=pool)
            if not ready: break
            for k in ready: ranks[k]=layer; pool.update(self.outputs[k])
            remaining.difference_update(ready); layer+=1
        return frozenset(pool),ranks

    def supported(self,k,pool): return self.inputs[k]<=pool and self.rows[k].catalyst.holds(pool)

    def is_raf(self,ids):
        ids=self.validate(ids); pool,_=self.closure(ids)
        return bool(ids) and all(self.supported(k,pool) for k in ids)

    def maximum(self,ids=None):
        current=self.universe if ids is None else self.validate(ids); counts=[len(current)]
        while True:
            pool,_=self.closure(current); kept=frozenset(k for k in current if self.supported(k,pool)); counts.append(len(kept))
            if kept==current: return current,counts
            current=kept

    def capability(self,target,ids=None):
        maximum,trace=self.maximum(ids); pool,_=self.closure(maximum)
        return bool(maximum) and target in pool

    def constructive(self):
        """Stricter than RAF: catalysts must already hold in the prefix pool."""
        pool=set(self.system.food); remaining=set(self.universe); layers=[]
        while True:
            ready=sorted(k for k in remaining if self.supported(k,pool))
            if not ready: break
            layers.append(ready)
            for k in ready: pool.update(self.outputs[k])
            remaining.difference_update(ready)
        return frozenset(self.universe-remaining),layers

    def is_closed(self,ids):
        ids=self.validate(ids)
        if not self.is_raf(ids): return False
        pool,_=self.closure(ids)
        return all(not self.supported(k,pool) for k in self.universe-ids)

    def closed_extension(self,seed):
        current=self.validate(seed)
        if current and not self.is_raf(current): raise ValueError('Closed extension requires an initial RAF or empty seed.')
        while True:
            pool,_=self.closure(current); enlarged=current|frozenset(k for k in self.universe if self.supported(k,pool))
            if enlarged==current: return current
            current=enlarged

    def closed_uniqueness_certificate(self):
        maximum,trace=self.maximum(); active,layers=self.constructive()
        equal=active==maximum
        return {'status':'unique_closed_raf' if equal and maximum else 'no_raf' if not maximum else 'unresolved_bounds',
            'set_equality':equal,'maximum':sorted(maximum),'activated':sorted(active),
            'pruning_counts':trace,'activation_layers':layers}

    def catalogue(self,target,budget=MAX_SUBSETS):
        ids=sorted(self.universe); required=2**len(ids)
        if required>budget: return {'status':'unresolved_budget','required_subsets':required}
        rafs=[]; closed=[]; producing=[]
        for mask in range(1,required):
            candidate=frozenset(ids[i] for i in range(len(ids)) if mask>>i&1)
            if not self.is_raf(candidate): continue
            rafs.append(candidate)
            if self.is_closed(candidate): closed.append(candidate)
            if target in self.closure(candidate)[0]: producing.append(candidate)
        minimal=lambda sets:[s for s in sets if not any(t<s for t in sets)]
        encode=lambda sets:[sorted(s) for s in sets]
        return {'status':'complete','subsets_checked':required-1,'raf_count':len(rafs),
            'irreducible':encode(minimal(rafs)),'closed':encode(closed),'target_rafs':encode(producing),
            'minimal_target_rafs':encode(minimal(producing))}

    def necessary_producers(self,target):
        """Sound unique-producer propagation inside this ambient set; may be incomplete."""
        forced=set(); pending=[target]; seen=set(); reasons=[]; impossible=[]
        while pending:
            species=pending.pop(0)
            if species in seen or species in self.system.food: continue
            seen.add(species); producers=sorted(k for k in self.universe if species in self.outputs[k])
            if not producers: impossible.append(species)
            if len(producers)!=1: continue
            k=producers[0]; reasons.append({'species':species,'only_producer':k})
            if k in forced: continue
            forced.add(k); pending.extend(sorted(self.inputs[k]|self.rows[k].catalyst.mandatory_atoms()))
        return {'forced_reactions':sorted(forced),'reasons':reasons,'impossible_species':impossible,
            'scope':'Necessary reactions for this target in this ambient set; ambiguous producers and optional catalyst branches are not resolved.'}


@dataclass(frozen=True)
class SiphonCertificate:
    species: frozenset

    def check(self,system):
        intersections=[]; violations=[]
        for r in system.reactions:
            if r.outputs&self.species:
                row={'reaction':r.id,'inputs_in_barrier':sorted(r.inputs&self.species),'outputs_in_barrier':sorted(r.outputs&self.species)}
                intersections.append(row)
                if not r.inputs&self.species: violations.append(r.id)
        food_overlap=system.food&self.species
        return {'valid_siphon':not violations,'food_disjoint':not food_overlap,
            'absence_certified':not violations and not food_overlap,'food_overlap':sorted(food_overlap),
            'violations':violations,'surviving_producers':intersections}


class InterventionStudy:
    def __init__(self,system,actions,target,preserve,barrier):
        self.system=system; self.actions=frozenset(actions); self.target=target; self.preserve=preserve; self.barrier=barrier

    def evaluate(self,restoration):
        cut=self.system.delete_actions(self.actions); analysis=RAFAnalysis(cut); barrier=self.barrier.check(cut)
        negative=barrier['absence_certified'] and self.target in self.barrier.species
        witnesses={}
        for action in sorted(self.actions):
            available=self.system.delete_actions(self.actions-{action}); evaluator=RAFAnalysis(available)
            ids=frozenset(restoration.get(action,()))
            valid=ids<=evaluator.universe and evaluator.is_raf(ids)
            pool,ranks=evaluator.closure(ids) if ids<=evaluator.universe else (frozenset(),{})
            witnesses[action]={'valid_target_raf':valid and self.target in pool,'size':len(ids),
                'reactions':sorted(ids),'reactant_ranks':ranks,'final_closure':sorted(pool)}
        return {'deleted_actions':sorted(self.actions),'deleted_directions':sorted(set(self.system.by_id)-set(cut.by_id)),
            'target_absence_by_siphon':negative,'target_capability_by_maxraf':analysis.capability(self.target),
            'preservation_capability':analysis.capability(self.preserve),'barrier':barrier,'restoration':witnesses,
            'inclusion_minimal_cut_certified':negative and all(w['valid_target_raf'] for w in witnesses.values()),
            'global_minimum_cost_claimed':False}


def main():
    parser=argparse.ArgumentParser(description=__doc__);parser.add_argument('--output',type=Path,default=Path('outputs'));args=parser.parse_args()
    out=args.output;out.mkdir(parents=True,exist_ok=True);base=Path(__file__).resolve().parent
    system=ReactionSystem.read(base/'prokaryotic-network.txt',PINNED_SOURCE_SHA256).add_food(EXTRA_FOOD)
    inputs=json.loads((base/'certificate_inputs.json').read_text())
    if inputs['source_sha256']!=system.source_sha256: raise ValueError('Certificates use a different source.')
    barrier=SiphonCertificate(frozenset(inputs['barrier']))
    study=InterventionStudy(system,DELETED_ACTIONS,TARGET,PRESERVATION_TARGET,barrier); result=study.evaluate(inputs['restoration'])
    parent=RAFAnalysis(system); cut_system=system.delete_actions(DELETED_ACTIONS); cut=RAFAnalysis(cut_system)
    small=RAFAnalysis(system.subset(SMALL_DIRECTIONS)); augmented=RAFAnalysis(system.subset((*SMALL_DIRECTIONS,POOLING_DIRECTION)))
    small_catalogue=small.catalogue(TARGET); augmented_catalogue=augmented.catalogue(TARGET)
    parent_certificate=parent.closed_uniqueness_certificate(); cut_certificate=cut.closed_uniqueness_certificate()
    L=frozenset(SMALL_DIRECTIONS[i] for i in (0,1,4,7,8)); H=frozenset(SMALL_DIRECTIONS)
    rescue_ids=tuple(SMALL_DIRECTIONS[i] for i in (7,8,4,5)); rescue=RAFAnalysis(cut_system.add_food((RESCUE_FOOD,)))
    preserved=RAFAnalysis(cut_system.subset(('R09639::f',)))
    necessity=small.necessary_producers(TARGET)
    result.update({'manuscript_sha256':MANUSCRIPT_SHA256,'source_sha256':system.source_sha256,
        'food':sorted(system.food),'source_actions':len({r.action for r in system.reactions}),'source_directions':len(system.reactions),
        'small_catalogue':small_catalogue,'augmented_catalogue':augmented_catalogue,'target_necessary_producers':necessity,
        'same_parent_closed_extensions':parent.closed_extension(L)==parent.closed_extension(H)==frozenset(parent_certificate['maximum']),
        'methionine_singleton_valid':preserved.is_raf(preserved.universe) and preserved.capability(PRESERVATION_TARGET),
        'precursor_rescue':{'food_added':RESCUE_FOOD,'reactions':list(rescue_ids),'valid_target_raf':rescue.is_raf(rescue_ids) and TARGET in rescue.closure(rescue_ids)[0]},
        'parent_bounds':{k:v for k,v in parent_certificate.items() if k not in ('maximum','activated','activation_layers')},
        'cut_bounds':{k:v for k,v in cut_certificate.items() if k not in ('maximum','activated','activation_layers')},
        'scope':'Exact finite source computations, not a Lean rerun. Structural capability is not flux, growth or concentration equilibria.'})
    for name,certificate in [('parent',parent_certificate),('postcut',cut_certificate)]:
        (out/f'{name}_closed_certificate.json').write_text(json.dumps(certificate,indent=2)+'\n')
    (out/'results.json').write_text(json.dumps(result,indent=2)+'\n')
    def csv_file(name,headers,rows):
        with (out/name).open('w',newline='') as f:
            writer=csv.writer(f);writer.writerow(headers);writer.writerows(rows)
    csv_file('restoration_witnesses.csv',['restored_action','direction','reactant_rank','catalyst_in_final_closure'],
        [[action,k,w['reactant_ranks'].get(k,''),system.by_id[k].catalyst.holds(set(w['final_closure']))] for action,w in result['restoration'].items() for k in w['reactions']])
    csv_file('siphon_producers.csv',['direction','reactants_in_barrier','products_in_barrier'],
        [[r['reaction'],';'.join(r['inputs_in_barrier']),';'.join(r['outputs_in_barrier'])] for r in result['barrier']['surviving_producers']])
    csv_file('activation_layers.csv',['model','synchronous_layer','new_directions','cumulative_directions'],
        [[name,i+1,len(layer),sum(map(len,c['activation_layers'][:i+1]))] for name,c in [('parent',parent_certificate),('postcut',cut_certificate)] for i,layer in enumerate(c['activation_layers'])])
    csv_file('small_reactions.csv',['index','direction','action','reactants','products','catalyst_formula_tree'],
        [[i,k,system.by_id[k].action,json.dumps(system.by_id[k].reactants),json.dumps(system.by_id[k].products),json.dumps(system.by_id[k].catalyst.to_data())] for i,k in enumerate(SMALL_DIRECTIONS)])
    lines=[f'Source: {len(system.reactions)} directions, {len(system.food)} food species/markers.',
        f'Cut removes {len(result["deleted_directions"])} directions; selective inclusion-minimal certificate: {result["inclusion_minimal_cut_certified"]}.',
        f'Restoration witness sizes: {[w["size"] for w in result["restoration"].values()]}; methionine singleton: {result["methionine_singleton_valid"]}.',
        f'H: {small_catalogue["raf_count"]} RAFs, {len(small_catalogue["irreducible"])} irrRAFs, {len(small_catalogue["closed"])} closed RAFs, {len(small_catalogue["target_rafs"])} target RAF.',
        f'Adding the food-enabled pooling row: {len(augmented_catalogue["closed"])} closed RAF.',
        f'Parent bounds: {parent_certificate["pruning_counts"]}; A=M: {parent_certificate["set_equality"]}.',
        f'Postcut bounds: {cut_certificate["pruning_counts"]}; A=M: {cut_certificate["set_equality"]}.',
        'No complete parent irrRAF catalogue, global cut optimum, depletion or kinetic uniqueness is claimed.']
    (out/'console.txt').write_text('\n'.join(lines)+'\n');print('\n'.join(lines))
    import matplotlib
    matplotlib.use('Agg')
    import matplotlib.pyplot as plt
    import numpy as np
    fig,axes=plt.subplots(1,2,figsize=(10,4.2),layout='constrained')
    for name,c,color in [('Parent',parent_certificate,'#187f91'),('After cut',cut_certificate,'#b95024')]:
        cumulative=np.cumsum([len(layer) for layer in c['activation_layers']])
        axes[0].step(range(1,len(cumulative)+1),cumulative,where='post',label=name,color=color)
        axes[0].axhline(len(c['maximum']),color=color,ls='--',lw=.8)
        axes[1].plot(range(len(c['pruning_counts'])),c['pruning_counts'],'o-',label=name,color=color)
    axes[0].set(xlabel='Constructive activation layer',ylabel='Activated directions',title='Constructive lower bound reaches maxRAF')
    axes[1].set(xlabel='Pruning step (last equality retained)',ylabel='Directions in RAF upper bound',title='Catalysts checked in final reactant closure')
    for ax in axes:ax.grid(alpha=.2);ax.legend(fontsize=8)
    fig.savefig(out/'parent_bounds.png',dpi=180);fig.savefig(out/'parent_bounds.svg');plt.close(fig)
    fig,ax=plt.subplots(figsize=(10,4.2),layout='constrained')
    rows=[set([SMALL_DIRECTIONS[0]]),set([SMALL_DIRECTIONS[7]]),set(L),set(H),set((*H,POOLING_DIRECTION))]
    labels=['irrRAF seed 0','irrRAF seed 7','Closed L in H','Only valine RAF in H','Only closed RAF in H + q']
    matrix=np.array([[k in row for k in (*SMALL_DIRECTIONS,POOLING_DIRECTION)] for row in rows])
    ax.imshow(matrix,cmap=matplotlib.colors.ListedColormap(['#f0f0ed','#187f91']),vmin=0,vmax=1,aspect='auto')
    ax.set_yticks(range(5),labels);ax.set_xticks(range(11),[str(i) for i in range(10)]+['q']);ax.set(xlabel='Literal subsystem reaction index (source mapping in CSV)',title='Reaction membership in catalytic seeds and\nvaline-producing sets')
    ax.set_xticks(np.arange(-.5,11),minor=True);ax.set_yticks(np.arange(-.5,5),minor=True);ax.grid(which='minor',color='white',linewidth=2);ax.tick_params(which='minor',bottom=False,left=False)
    fig.savefig(out/'subsystem.png',dpi=180);fig.savefig(out/'subsystem.svg');plt.close(fig)
    digest=lambda p:hashlib.sha256(p.read_bytes()).hexdigest()
    (out/'run_metadata.json').write_text(json.dumps({'manuscript_sha256':MANUSCRIPT_SHA256,'source_sha256':digest(Path(__file__)),
        'input_sha256':{name:digest(base/name) for name in ('prokaryotic-network.txt','certificate_inputs.json')},'python':platform.python_version(),
        'output_sha256':{p.name:digest(p) for p in sorted(out.iterdir()) if p.is_file() and p.name!='run_metadata.json'}},indent=2)+'\n')


if __name__=='__main__':main()

Run output

Source: 9231 directions, 68 food species/markers.
Cut removes 4 directions; selective inclusion-minimal certificate: True.
Restoration witness sizes: [9, 13, 13]; methionine singleton: True.
H: 22 RAFs, 2 irrRAFs, 2 closed RAFs, 1 target RAF.
Adding the food-enabled pooling row: 1 closed RAF.
Parent bounds: [9231, 2823, 2148, 2148]; A=M: True.
Postcut bounds: [9227, 2584, 2085, 2085]; A=M: True.
No complete parent irrRAF catalogue, global cut optimum, depletion or kinetic uniqueness is claimed.