Example code — Biochemical readouts that preserve native function Reporter storage, finite-time inference and certified recovery in a shared cofactor pool

A biochemical reporter measures a shared cofactor by binding it, temporarily making that cofactor unavailable to the native reaction. This runnable model retains the bound reporter complex, the regenerated pool, native product and reporter product. It compares reporters with equal effective activity and integrates the native-output loss through acquisition and recovery.

With regeneration rate p, native-use rate k and reporter release rate c, the complete native loss per reporter product is k(1 + p/c)/(p + k). The extra factor 1 + p/c comes from cofactor storage: at p = 1 and release rate c = 0.2 it multiplies the instantaneous-turnover prediction by six. Pulsing or delaying the same reporter changes the amount of signal, but cannot remove this cost at fixed final product.

Slow, fast and delayed reporters with equal effective activity impose different peak suppression and cumulative native-output losses. — Literal source trajectories with the reporter complex retained. Native loss is independently integrated against the matched reporter-free reference, including recovery after acquisition.

Bounded-error reporter records give outer intervals for native function; slower association shutoff requires a longer certified recovery wait. — Left: exact envelope inversion, not a confidence interval or a fitted point estimate. Right: exact certificates at the paper's declared decay rates; connecting lines guide the eye and do not certify intermediate rates.

The example reproduces a design that bounds native-flux suppression below 4.73% while separating the two promised functional classes after eight time units. Exact rational certificates check the parameter box, detector margin and recovery deadline. Numerical trajectories independently check the finite accounting identity; they do not establish coverage of the whole box.

If reporter binding switches off gradually, residual binding causes additional loss. The code bounds that tail and distinguishes a small outstanding loss from the extra loss already incurred. Without the two-class promise, it returns an outer interval for native function, with unresolved and incompatible outcomes retained.

Editable rates, reporter properties, schedules and measurement contracts sit above composable reactor and certificate classes. A nonlinear rate-law extension and an explicit native/background ambiguity show where the model can be extended and where reporter data alone cannot answer the question. These are synthetic design requirements, not experimentally validated assay settings.

Python source

"""Preserving native function while measuring a shared cofactor pool."""
from fractions import Fraction as F
from pathlib import Path
import argparse,csv,hashlib,json,platform
from dataclasses import asdict
import numpy as np
from scipy.optimize import brentq
from reactor import LinearRate,MichaelisMenten,Reporter,AssociationSchedule,CofactorReactor
from certificates import DesignBox,RecoveryCertificate,slope_sandwich

# EDITABLE INPUTS: dimensionless synthetic design requirements, not fitted rates.
REGENERATION=1.
NATIVE_USE=1.
COFACTOR_POOL=1.
REPORTER_TOTAL=1.
REPORTER_ACTIVITY=.08
CATALYTIC_RELEASE=20.
DISSOCIATION=1.
PREPARATION_OFFSET=0.
ACQUISITION='8'
RECOVERY='5'
FADING_RATE='21'
RECOVERY_TARGET='0.000003'
PRESERVATION_ALLOWANCE='0.05'
OBSERVED_RECORDS=('0.30549','0.41151','0.361526225','-0.1','0.6')
FUNCTIONAL_THRESHOLD='0.6'
CONCENTRATION_SCALE_MICROMOLAR='1'
TIME_SCALE_SECONDS='60'
VOLUME_MICROLITRES='100'
INTERIOR_DRAWS_PER_CLASS=8
SEED=17092026
# Primitive association is derived from u, beta, c and reporter total; do not
# vary it independently while asserting that effective activity remains fixed.
UNCERTAINTY=DesignBox()
MANUSCRIPT_SHA256='bd27ad3fee701238aedb2442336de04b53187feb6e7b9cd23072c85c3491adbf'


def main():
    parser=argparse.ArgumentParser(description=__doc__);parser.add_argument('--output',type=Path,default=Path('outputs'))
    out=parser.parse_args().output;out.mkdir(parents=True,exist_ok=True)
    def dump(name,data):
        (out/name).write_text(json.dumps(data,indent=2,default=lambda x:{'exact':str(x),'decimal':float(x)} if isinstance(x,F) else str(x))+'\n',encoding='utf-8')
    def table(name,header,rows):
        with (out/name).open('w',newline='',encoding='utf-8') as f:w=csv.writer(f);w.writerow(header);w.writerows(rows)
    T=float(F(ACQUISITION));tau=float(F(RECOVERY))
    def nominal(c=CATALYTIC_RELEASE,p=REGENERATION,k=NATIVE_USE,background=0):
        return CofactorReactor(LinearRate(p),LinearRate(k),Reporter(REPORTER_TOTAL,c,DISSOCIATION,REPORTER_ACTIVITY),COFACTOR_POOL,background)
    comparisons=[];traces={}
    for name,c,delay in [('slow',1.,0.),('fast',20.,0.),('late',20.,T/2)]:
        schedule=AssociationSchedule(T,((delay,T),));t,y,row=nominal(c).run(schedule,tau,PREPARATION_OFFSET)
        _,_,hi=nominal(c,p=2.).run(schedule,tau,PREPARATION_OFFSET)
        row.update(name=name,release=c,delay=delay,nominal_signal_gap_minus_002=hi['acquisition_signal']-row['acquisition_signal']-.02)
        comparisons.append(row);traces[name]=(t,y)
        table(f'trajectory_{name}.csv',['time',*CofactorReactor.columns],zip(t,*y.T))
    configured=nominal().run(AssociationSchedule(T,((0.,T),),float(F(FADING_RATE))),tau,PREPARATION_OFFSET)[2]
    schedules=[]
    for name,windows,gamma in [('constant',((0.,T),),None),('late',((T/2,T),),None),
                               ('pulses',((0.,T/4),(T/2,3*T/4)),None),('fade',((0.,T),),5.),('late-fade',((T/4,T),),1.)]:
        _,_,r=nominal().run(AssociationSchedule(T,windows,gamma),40.,PREPARATION_OFFSET)
        schedules.append(dict(name=name,**r))
    # Equal final product, adjusted exposure; chemistry still determines the cost.
    fixed_product=[]
    for c in [.2,1.,5.,20.]:
        reactor=nominal(c)
        def output(duration):return reactor.run(AssociationSchedule(duration,((0.,duration),)),30.,PREPARATION_OFFSET,samples=31)[2]['complete_product_enclosure_numeric'][0]
        target=.1;right=1.
        while output(right)<target and right<128:right*=2
        if output(right)<target:fixed_product.append(dict(release=c,status='duration-budget-exhausted'));continue
        duration=brentq(lambda d:output(d)-target,1e-5,right,xtol=1e-10)
        result=reactor.run(AssociationSchedule(duration,((0.,duration),)),30.,PREPARATION_OFFSET,samples=101)[2]
        fixed_product.append(dict(release=c,duration=duration,target_product=target,**result))
    box=UNCERTAINTY;general=box.contract(F(ACQUISITION),False,F(PRESERVATION_ALLOWANCE))
    rounded=box.contract(F(ACQUISITION),True,F(PRESERVATION_ALLOWANCE))
    rec=RecoveryCertificate(box);deadlines=[]
    for gamma,n0,n1,wait in [('0.5','0.5','0.5','22.7'),('1','1','1','11.2'),('2','1.85','2','6.8'),('5','1.93','5','5.2'),('10','1.93','10','4.9'),('21','1.93','12','4.9')]:
        gamma,n0,n1,wait=map(F,(gamma,n0,n1,wait));r=rec.fading(gamma,wait,n0,n1);previous=rec.fading(gamma,wait-F(1,10),n0,n1)
        deadlines.append(dict(gamma=gamma,nu0=n0,nu1=n1,wait=wait,**r,certified=r['loss_upper']<F(RECOVERY_TARGET),
                              previous_envelope_above_target=previous['envelope_lower']>F(RECOVERY_TARGET)))
    fading=[];base=nominal().run(AssociationSchedule(T,((0.,T),)),40.,PREPARATION_OFFSET)[2]
    for gamma in [1.,5.,21.]:
        row=nominal().run(AssociationSchedule(T,((0.,T),),gamma),40.,PREPARATION_OFFSET)[2]
        row.update(gamma=gamma,additional_complete_loss_numeric=row['integrated_plus_tail_loss_enclosure_numeric'][0]-base['integrated_plus_tail_loss_enclosure_numeric'][0]);fading.append(row)
    rng=np.random.default_rng(SEED);interior=[]
    for label,prange in [('low',box.low),('high',box.high)]:
        for i in range(INTERIOR_DRAWS_PER_CLASS):
            draw=lambda iv:rng.uniform(float(iv[0]),float(iv[1]))
            p=draw(prange);k=draw(box.native);C=draw(box.pool);u=draw(box.activity);c=draw(box.release);beta=draw(box.dissociation);R=draw(box.reporter);s=rng.uniform(-float(box.preparation),float(box.preparation))
            t,y,row=CofactorReactor(LinearRate(p),LinearRate(k),Reporter(R,c,beta,u),C).run(AssociationSchedule(T,((0.,T),)),tau,s)
            # Both extreme gain/noise records are checked, not a lucky noise draw.
            records=[F(str(row['acquisition_signal']))*box.gain[0]-box.error,F(str(row['acquisition_signal']))*box.gain[1]+box.error]
            interior.append(dict(label=label,parameters=dict(p=p,k=k,C=C,u=u,c=c,beta=beta,R=R,s=s),
                classified_endpoints=[rounded.classify(v,True) for v in records],**row))
    inference=[]
    for y in OBSERVED_RECORDS:
        r=rounded.inverse(F(y));flux=r['flux'];threshold=F(FUNCTIONAL_THRESHOLD)
        decision='incompatible' if flux is None else ('below' if flux[1]<threshold else 'above' if flux[0]>threshold else 'unresolved')
        inference.append(dict(record=F(y),**r,binary_with_promise=rounded.classify(F(y),True),binary_without_promise=rounded.classify(F(y)),functional_decision=decision))
    invrows=[]
    for y in np.linspace(.2,.5,180):
        r=rounded.inverse(F(str(y)))
        if r['p'] is not None:invrows.append([y,*map(float,r['p']),*map(float,r['flux'])])
    table('inference_sweep.csv',['record','p_lower','p_upper','flux_lower','flux_upper'],invrows)
    nonlinear=[]
    for scale in [1,10,100,1000]:
        K=F(10*scale);lo=K*K/(K+1)**2;sandwich=slope_sandwich(lo,1,lo,1,20)
        law=MichaelisMenten(float(K),float(K));model=CofactorReactor(law,law,Reporter(),1.)
        t,y,r=model.run(AssociationSchedule(),40.)
        # This is a long finite horizon. No linear tail formula is applied to MM.
        nonlinear.append(dict(affinity=K,slope_range=(lo,F(1)),cost_sandwich=sandwich,finite_horizon_ratio=float(y[-1,3]/y[-1,2]),final_deficit=float(y[-1,4]-y[-1,0]),evidence='Finite-horizon numerical ratio; infinite-horizon sandwich is a theorem under slope and integrability premises.'))
    alias=[];aliaspaths=[]
    for k,l in [(.5,1.5),(1.5,.5)]:
        t,y,r=nominal(k=k,background=l).run(AssociationSchedule(T,((0.,T),),1.),20.)
        aliaspaths.append(y[:,:3]);alias.append(dict(native=k,background=l,native_stationary_flux=k*REGENERATION*COFACTOR_POOL/(REGENERATION+k+l),signal=r['acquisition_signal']))
    # Native deficits do differ; identical x,b,q alone are the observation alias.
    alias_difference=float(np.max(abs(aliaspaths[0]-aliaspaths[1])))
    def contract_record(c):return dict(rounded=c.rounded,AL=c.AL,AU=c.AU,U=c.U,L=c.L,margin=c.margin,threshold=c.threshold,certified=c.certified)
    summary=dict(constants=box.constants(),box=asdict(box),rounded=contract_record(rounded),unrounded=contract_record(general),
        configured_trajectory=configured,ideal_recovery_upper=rec.ideal(F(RECOVERY)),fast_fading_recovery=rec.fading(21,5,F('1.93'),12),
        comparison=comparisons,schedules=schedules,fixed_product=fixed_product,fading=fading,interior=interior,
        alias=dict(cases=alias,max_observable_difference=alias_difference,note='The source equations depend on native plus background consumption. Native-specific output is not identifiable from these common records.'),
        illustrative_units=dict(acquisition_seconds=F(ACQUISITION)*F(TIME_SCALE_SECONDS),recovery_seconds=F(RECOVERY)*F(TIME_SCALE_SECONDS),
            donor_cap_micromolar=box.high[1]*box.pool[1]*(F(ACQUISITION)+F(RECOVERY))*F(CONCENTRATION_SCALE_MICROMOLAR),
            native_cosubstrate_cap_micromolar=box.native[1]*box.pool[1]*(F(ACQUISITION)+F(RECOVERY))*F(CONCENTRATION_SCALE_MICROMOLAR),
            reporter_cosubstrate_cap_micromolar=box.activity[1]*(1+box.preparation)*box.high[1]*box.pool[1]/(box.high[1]+box.native[0])*(F(ACQUISITION)+F(RECOVERY))*F(CONCENTRATION_SCALE_MICROMOLAR),
            low_complete_loss_picomoles=base['integrated_plus_tail_loss_enclosure_numeric'][0]*float(F(CONCENTRATION_SCALE_MICROMOLAR)*F(VOLUME_MICROLITRES)),
            note='Illustrative scale conversion and finite material caps; not stock recipes or empirical calibration.'))
    dump('summary.json',summary);dump('recovery.json',deadlines);dump('inference.json',inference);dump('nonlinear.json',nonlinear)
    plot(out,traces,invrows,deadlines)
    print(f'Uniform suppression upper: {float(box.constants()["suppression"]):.10%}')
    print(f'Rounded residual detector margin: {float(rounded.margin):.12f}; joint design certified: {rounded.certified}')
    print('Native loss is integrated independently; fading tails include remaining association.',flush=True)
    digest=lambda p:hashlib.sha256(p.read_bytes()).hexdigest();here=Path(__file__).parent
    dump('run_metadata.json',dict(manuscript_sha256=MANUSCRIPT_SHA256,python=platform.python_version(),source_sha256=digest(Path(__file__)),module_sha256={n:digest(here/n) for n in ['reactor.py','certificates.py']},output_sha256={p.name:digest(p) for p in sorted(out.iterdir()) if p.is_file() and p.name!='run_metadata.json'}))


def plot(out,traces,inverse,deadlines):
    import matplotlib
    matplotlib.use('Agg')
    import matplotlib.pyplot as plt
    fig,axs=plt.subplots(1,2,figsize=(10,4),layout='constrained')
    for name,(t,y) in traces.items():
        axs[0].plot(t,100*(y[:,4]-y[:,0])/y[:,4],label=name)
        axs[1].plot(y[:,2],y[:,3],label=name)
    axs[0].axhline(5,color='black',ls=':');axs[0].set(xlabel='Time including recovery',ylabel='Native-flux suppression (%)',title='Native-flux suppression by reporter design')
    axs[1].set(xlabel='Accumulated reporter product',ylabel='Independently integrated native-output loss',title='Native-output loss versus reporter product')
    for ax in axs:ax.grid(alpha=.2);ax.legend(fontsize=8)
    fig.savefig(out/'storage.png',dpi=180);fig.savefig(out/'storage.svg');plt.close(fig)
    fig,axs=plt.subplots(1,2,figsize=(10,4),layout='constrained');iv=np.array(inverse)
    axs[0].fill_between(iv[:,0],iv[:,3],iv[:,4],alpha=.25);axs[0].plot(iv[:,0],iv[:,3]);axs[0].plot(iv[:,0],iv[:,4]);axs[0].axhline(float(F(FUNCTIONAL_THRESHOLD)),color='black',ls=':')
    axs[0].set(xlabel='Bounded-error reporter record',ylabel='Outer interval for stationary native flux',title='Native-flux bounds from reporter\nmeasurements')
    axs[1].plot([float(r['gamma']) for r in deadlines],[float(r['wait']) for r in deadlines],'o-');axs[1].axhline(5,color='black',ls=':')
    axs[1].set(xscale='log',xlabel='Association fade rate',ylabel='Certified recovery wait',title='Required recovery time versus reporter\nshutoff rate')
    for ax in axs:ax.grid(alpha=.2)
    fig.savefig(out/'inference-recovery.png',dpi=180);fig.savefig(out/'inference-recovery.svg');plt.close(fig)


if __name__=='__main__':main()

Run output

Uniform suppression upper: 4.7205487518%
Rounded residual detector margin: 0.005557315359; joint design certified: True
Native loss is integrated independently; fading tails include remaining association.