Example code — Zero-divisor methods for absolute concentration robustness: completeness, positive geometry, and quantitative certificates

Absolute concentration robustness means a species has the same concentration at every positive steady state, even when conserved totals vary. Algebraic candidate values may instead describe states with zero or negative concentrations. Even a valid setpoint becomes difficult to infer from nearly stationary measurements when a supporting species is scarce.

The example connects those distinctions to usable models. It recomputes the candidate counterexamples, provides the full loaded-reactor and EnvZ/OmpR kinetics, and checks polynomial identities that turn residual measurements into concentration bounds. Editable rate, load, pool and uncertainty inputs sit at the top of the driver.

Increasing load shrinks the multiplier pool and loses its required floor before the positive branch ends; a full reactor trajectory returns near the nominal setpoint. — The load envelope follows exact equilibrium formulas. The transient is numerical; its multiplier floor is justified separately by the nominal invariant-region certificate and initial-state membership.

Residuals can approach zero while concentration error stays fixed if the multiplier disappears; slower private release increases stored material without changing the steady-state quotient. — A small residual needs a multiplier floor to imply accuracy. Private release preserves steady-state geometry but carries explicit storage and residual costs.

In the loaded reactor, the setpoint stays at 2.2 micromolar while increasing load depletes species B. Its concentration multiplies the setpoint error in the identity $b'=k_1b(a-\alpha)$ , linking concentration error to the rate of change of B. The required accuracy margin disappears before the positive equilibrium does. Exact parameter-box bounds and a freshly checked nominal invariant ellipsoid identify where the floor is justified; numerical trajectories illustrate the full dynamics.

Reusable product-release components show another cost: private intermediates preserve the old steady-state algebra while storing material and adding residual terms. Leakage can even remove the positive operating point. Candidate lists, equilibrium floors and trajectory guarantees remain distinct throughout. Rates are illustrative, and Lean is not rerun.

Python source

"""From ACR candidate algebra to a usable, conditioned concentration guarantee."""
# EDIT THESE INPUTS. Reactor concentrations: micromolar; time: minutes.
REACTOR_RATES=('1/20','1/20','1/20','1/20','3/20')
DILUTION='1/100'
FEED_CONCENTRATION='10'
NOMINAL_LOAD='1/50'
REQUIRED_B_FLOOR='1'
B_DERIVATIVE_ERROR='1/1000'  # a derivative bound, not concentration measurement noise
RATE_RELATIVE_UNCERTAINTY='1/100'
REACTOR_INITIAL_DISPLACEMENT=('1/50','-1/100','1/100')
TRAJECTORY_MINUTES=600
ENVZ_RATES=('1',)*9
ENVZ_TOTAL_X='3/2'
ENVZ_TOTAL_Y='9/2'
RELEASE_RATES=('1','1')
RELEASE_LEAKAGE=('1/100','1/100')

import argparse,csv,hashlib,json,platform
from pathlib import Path
import numpy as np
import scipy
import sympy as s
from certificates import Interval,ResidualCertificate,BlockCandidates
from models import Reaction,MassActionModel,LoadedReactor,EnvZOmpR,envz_box_floors,PrivateRelease
from invariant_region import InvariantEllipsoid,P
import paper_checks
R=s.Rational
MANUSCRIPT_SHA256='91f82500b08f1947437c95b688e681504e38648ef338bb8ae46dfda58439d776'

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,x):(out/name).write_text(json.dumps(x,indent=2,default=lambda x:x.tolist() if isinstance(x,np.ndarray) else str(x))+'\n')
    def table(name,rows):
        with (out/name).open('w',newline='') as f:w=csv.DictWriter(f,fieldnames=list(rows[0]));w.writeheader();w.writerows(rows)
    replay=paper_checks.replay();dump('source_arithmetic_replay.json',replay)
    # Original ideals only: no radical, saturation or conservation equations.
    a,b,c,d=s.symbols('a b c d');fields=[(-a*b+b-a*c+3*c,b*(a-2),c*(a-3)+b),(-a*a+c+d-a*b+b-a,d-b*b,a*a-c,c-d)];algebra=[]
    for label,variables,field,point,positive_branch in [('boundary_adds_candidate',(a,b,c),fields[0],(2,1,1),'(2,tau,tau), tau>0: ACR a=2'),('nonpositive_point_supplies_only_candidate',(a,b,c,d),fields[1],(1,1,1,1),'(tau,tau,tau^2,tau^2), tau>0: no ACR')]:
        candidates=BlockCandidates(variables,a).compute(field);algebra.append(dict(name=label,**candidates,positive_geometry=positive_branch,identity_at_positive_point=[str(s.expand(f).subs(dict(zip(variables,point)))) for f in field]))
    u,v,t=s.symbols('u v t');ideal=[u*t-v,v*t-u];mixed=s.groebner(ideal,u,v,t,order=lambda e:(e[0]+e[1],e[2],e[0]));coeff=[]
    for p in mixed.polys:coeff.append(max(s.Poly(p.as_expr(),u,v).terms(),key=lambda term:(sum(term[0]),term[0][0]))[1])
    mixed_roots=sorted({r for c0 in coeff for r in s.Poly(c0,t).real_roots() if r>0});dump('algebraic_candidates_and_positive_geometry.json',algebra);dump('monomial_order_comparison.json',dict(mixed_basis=[str(p.as_expr()) for p in mixed.polys],mixed_coefficients=list(map(str,coeff)),mixed_candidates=list(map(str,mixed_roots)),block=BlockCandidates((u,v,t),t).compute(ideal),scope='The mixed order eliminates u,v but is not a block order. Its empty candidate list misses the actual ACR value 1.'))
    reactor=LoadedReactor(REACTOR_RATES,DILUTION,FEED_CONCENTRATION,NOMINAL_LOAD);model=reactor.model();eq=reactor.equilibrium(REQUIRED_B_FLOOR);dump('configured_reactor_equilibrium.json',eq)
    if eq['status']!='positive_equilibrium':raise ValueError('Configured reactor has no positive operating point. The exact reason is saved; choose admissible parameters for the transient demo.')
    a,b,c=model.variables;k,D,ci,ell=reactor.parameters();alpha=eq['alpha'];cert=ResidualCertificate(model.variables,0,alpha,k[0]*b,(s.Integer(0),s.Integer(1),s.Integer(0)));box=(Interval(0,ci),Interval(R(REQUIRED_B_FLOOR),ci),Interval(0,ci));floor=k[0]*R(REQUIRED_B_FLOOR);eps=(0,R(B_DERIVATIVE_ERROR),0)
    good=cert.evaluate(model.field,box,eps,floor,(-a/50,0,0));dump('reactor_residual_certificate.json',good)
    bad=ResidualCertificate(model.variables,0,alpha,2*k[0]*b,(0,1,0)).evaluate(model.field,box,eps,floor);nofloor=cert.evaluate(model.field,(box[0],Interval(0,ci),box[2]),eps,floor);dump('rejected_and_unresolved_certificates.json',dict(wrong_identity=bad,unsupported_floor=nofloor))
    direct=cert.evaluate(model.field,box,eps,floor,(0,b/100,0));dump('direct_multiplier_load_certificate.json',direct)
    # Disturbances that cancel after q dot d are simplified before enclosure.
    shared=ResidualCertificate((a,b),0,1,b,(1,1));signed=shared.evaluate((b*(a-1),s.Integer(0)),(Interval(0,3),Interval(1,2)),(R(1,1000),R(1,1000)),1,(a*b,-a*b));dump('signed_load_cancellation.json',signed)
    sweep=[]
    for load in [R(j,10000) for j in range(401)]:
        r=LoadedReactor(REACTOR_RATES,DILUTION,FEED_CONCENTRATION,load);e=r.equilibrium(REQUIRED_B_FLOOR);sweep.append(dict(load=float(load),status=e['status'],setpoint=float(e['alpha']),b_pool=float(e['state'][1]) if e['status']=='positive_equilibrium' else '',floor_satisfied=e.get('floor_satisfied',False),local_stability=e.get('locally_stable',False)))
    table('load_operating_envelope.csv',sweep)
    rb=reactor_box(reactor,R(RATE_RELATIVE_UNCERTAINTY));dump('reactor_parameter_box.json',rb)
    ellipsoid=InvariantEllipsoid(model,eq['state']);proof=ellipsoid.verify();dump('invariant_ellipsoid_certificate.json',proof)
    initial=tuple(x+R(v) for x,v in zip(eq['state'],REACTOR_INITIAL_DISPLACEMENT));inside=ellipsoid.contains(initial) if proof['status']=='certified_forward_invariant' else False
    times,Y=model.integrate(initial,TRAJECTORY_MINUTES);_,Y2=model.integrate(initial,TRAJECTORY_MINUTES,method='BDF');center=np.array(eq['state'],float);U=Y-center[:,None];V=np.einsum('it,ij,jt->t',U,np.array(P,float),U);diagnostics=dict(initial_inside_certified_region=inside,solver_difference=float(np.max(abs(Y-Y2))),minimum_b=float(min(Y[1])),initial_V=float(V[0]),maximum_V=float(max(V)),final_error=float(np.linalg.norm(Y[:,-1]-center)),scope='Numerical transient. The fixed-point ellipsoid theorem applies only if its exact checks and initial membership both pass.')
    table('full_reactor_trajectory.csv',[dict(time_min=t,a=Y[0,i],b=Y[1,i],c=Y[2,i],V=V[i]) for i,t in enumerate(times)]);dump('reactor_trajectory_checks.json',diagnostics)
    extinction=[dict(b=str(R(1,10**j)),a_error='1',b_residual=str(k[0]*R(1,10**j))) for j in range(1,10)];dump('small_residual_extinction_counterexample.json',extinction)
    env=EnvZOmpR(ENVZ_RATES);em=env.model();ee=env.equilibrium(ENVZ_TOTAL_X,ENVZ_TOTAL_Y);dump('envz_equilibrium.json',ee)
    if ee['status']=='positive_equilibrium':
        dump('envz_regularity.json',em.regularity(ee['state']));es=np.array([float(v) for v in ee['state']]);ei=es.copy();shift=min(es[0],es[1])/100;ei[0]+=shift;ei[1]-=shift;et,EY=em.integrate(ei,60);poolX=EY[[0,1,2,5,6]].sum(axis=0);poolY=EY[[3,4,5,6]].sum(axis=0);table('envz_full_trajectory.csv',[dict(time_min=t,**{name:EY[j,i] for j,name in enumerate(em.species)}) for i,t in enumerate(et)]);dump('envz_conservation_checks.json',dict(X_total_error=float(max(abs(poolX-float(R(ENVZ_TOTAL_X))))),Y_total_error=float(max(abs(poolY-float(R(ENVZ_TOTAL_Y))))),scope='Numerical trajectory, no global stability assertion.'))
    eb=envz_box_floors((Interval('99/100','101/100'),)*9,Interval(1,2),Interval(4,5));dump('envz_uniform_equilibrium_floors.json',eb)
    x=em.variables;ek,ealpha,*_=env.constants();q1=(0,0,-(ek[7]+ek[8]),ek[7]+ek[8],0,0,ek[8]);q2=(0,0,-ek[6]*x[4],ek[6]*x[4],0,0,ek[2]);domain=tuple(Interval('1/10','5') for _ in x)
    c1=ResidualCertificate(x,4,ealpha,ek[6]*ek[8]*x[1],q1);c2=ResidualCertificate(x,4,ealpha,ek[6]*ek[8]*x[6],q2)
    dump('envz_two_measured_domain_certificates.json',dict(first=c1.evaluate(em.field,domain,(R(1,1000),)*7,ek[6]*ek[8]/10),second=c2.evaluate(em.field,domain,(R(1,1000),)*7,ek[6]*ek[8]/10),scope='A separately stipulated measured-state box, not the equilibrium-only pool floors.'))
    release=release_demo(out,dump,table)
    plot(out,sweep,eq,times,Y,V,release,extinction)
    print('Source symbolic/rational checks:',replay['count'],replay['status'],flush=True)
    print('Candidate lists:',[a['candidates'] for a in algebra],'; mixed/block order:',mixed_roots,'/ [1]',flush=True)
    print('Loaded reactor:',tuple(map(str,eq['state'])),'; floor limit',eq['floor_load_limit'],'; positive-branch limit',eq['critical_load'],flush=True)
    print('Conditional concentration error:',good['absolute_error_bound'],'micromolar; ellipsoid:',proof['status'],'; initial inside:',inside,flush=True)
    print('Equilibrium floors do not automatically hold along transients. Private release preserves steady-state geometry, not trajectories or storage. Illustrative rates; Lean not rerun.',flush=True)
    here=Path(__file__).parent;digest=lambda p:hashlib.sha256(p.read_bytes()).hexdigest();dump('run_metadata.json',dict(manuscript_sha256=MANUSCRIPT_SHA256,source_sha256=digest(Path(__file__)),python=platform.python_version(),sympy=s.__version__,scipy=scipy.__version__,numpy=np.__version__,module_sha256={p.name:digest(p) for p in sorted(here.glob('*.py')) if p.name not in ['example.py','test_example.py']},output_sha256={p.name:digest(p) for p in sorted(out.iterdir()) if p.is_file() and p.name!='run_metadata.json'}))

def reactor_box(reactor,uncertainty):
    k,D,ci,ell=reactor.parameters()
    if not 0<=uncertainty<1:raise ValueError('Relative uncertainty must lie in [0,1).')
    box=lambda value:Interval(value*(1-uncertainty),value*(1+uncertainty));kb=list(map(box,k));cb=box(ci);lb=box(ell);alpha=(kb[1]+kb[2]+D)/kb[0];delta=kb[4]+D-kb[3]*alpha;N=cb-alpha*(1+lb/D)
    if delta.lo<=0 or N.lo<=0:return dict(status='not_certified_uniform_positive_region',alpha=alpha.strings(),delta=delta.strings(),pool=N.strings())
    r=delta/kb[2];floor=r.lo*N.lo/(1+r.lo);return dict(status='certified_equilibrium_region',alpha=alpha.strings(),delta_lower=str(delta.lo),pool_lower=str(N.lo),ratio_lower=str(r.lo),b_floor=str(floor),nominal_setpoint_uncertainty=str(max(abs(alpha.lo-(k[1]+k[2]+D)/k[0]),abs(alpha.hi-(k[1]+k[2]+D)/k[0]))))

def release_demo(out,dump,table):
    # A+B -> 3B, B -> A, B -> empty. Positive steady states: a=1, b>0.
    old=MassActionModel(('a','b'),[Reaction((1,1),(0,3),1),Reaction((0,1),(1,0),1),Reaction((0,1),(0,0),1)]);release=PrivateRelease(0,RELEASE_RATES);new,W=release.apply(old);a,b,z1,z2=new.variables;F=a*b;lam=tuple(map(R,RELEASE_RATES));reconstruction={z1:F/lam[0],z2:F/lam[1]};defects=[s.simplify(f.subs(reconstruction)) for f in new.field];identity=s.expand((a-1)*b-new.field[1]/2-new.field[2]*R(3,2)-new.field[3])
    if any(s.expand(a-b)!=0 for a,b in zip(defects[:2],old.field)) or any(defects[2:]) or identity!=0:raise ArithmeticError('Private-release quotient or transported identity failed.')
    dump('release_equilibrium_and_residual_identity.json',dict(old_field=list(map(str,old.field)),extended_field=list(map(str,new.field)),old_equations_recovered=True,transported_identity='(a-1)b = new_b_derivative/2 + 3*private_1_derivative/2 + private_2_derivative',quotient_scope='Triangular elimination in the original steady-state ideal; no fast-release approximation.',zero_divisor_candidates_old=BlockCandidates(old.variables,old.variables[0]).compute(old.field),zero_divisor_candidates_extended=BlockCandidates(new.variables,new.variables[0]).compute(new.field)))
    costs=[]
    for l in map(R,['1/100','1/10','1','10']):
        chain=PrivateRelease(0,(l,l));c=chain.steady_costs(1,weights=(3,2));costs.append(dict(release_rate=float(l),unit_weight_storage=float(chain.steady_costs(1)['weighted_storage']),cargo_weighted_storage=float(c['weighted_storage'])))
    table('release_storage_sweep.csv',costs);leak=PrivateRelease(0,RELEASE_RATES,RELEASE_LEAKAGE);leak_model,LW=leak.apply(old);lc=leak.steady_costs(1);dump('release_leakage.json',lc)
    # At a positive full steady state a'=b(1-a) forces a=1. Any nonzero
    # leak makes the b-balance strictly negative: there is no positive state.
    yields=lc['product_yields'];effective=R(1,2)*((1-yields[0])+2*(1-yields[1]));dump('release_leakage_existence_check.json',dict(certificate_defect_per_input_flux=effective,positive_equilibrium_exists=bool(effective==0),reason='a balance forces a=1; b balance is -b times a strictly positive release loss when any stage leaks. This example demonstrates why a conditional leakage error bound is not an existence proof.'))
    tt,YY=new.integrate((R(4,5),2,R(1,10),R(1,10)),60);ledger=2*YY[0]+YY[1]+3*YY[2]+2*YY[3];dump('release_trajectory_checks.json',dict(weighted_balance_error=float(max(abs(ledger-ledger[0]))),scope='For these chosen equal B-to-A and B-loss rates only, the weighted field balance cancels. This is not an asserted general conservation law under release.'));table('release_full_trajectory.csv',[dict(time_min=t,**{label:YY[j,i] for j,label in enumerate(new.species)}) for i,t in enumerate(tt)])
    return costs

def plot(out,sweep,eq,times,Y,V,release,extinction):
    import matplotlib
    matplotlib.use('Agg')
    import matplotlib.pyplot as plt
    fig,axs=plt.subplots(1,2,figsize=(11,4),layout='constrained');valid=[r for r in sweep if r['status']=='positive_equilibrium']
    axs[0].plot([r['load'] for r in valid],[r['b_pool'] for r in valid],label='Equilibrium multiplier pool b');axs[0].axhline(float(R(REQUIRED_B_FLOOR)),color='C1',ls='--',label='Required floor');axs[0].axvline(float(eq['floor_load_limit']),color='C1',ls=':');axs[0].axvline(float(eq['critical_load']),color='C3',ls=':');axs[0].set(xlabel='Specific load / per minute',ylabel='Concentration / micromolar',title='Supporting pool versus reactor load');axs[0].legend(fontsize=8)
    for i,name in enumerate(['a','b','c']):axs[1].plot(times,Y[i],label=name)
    axs[1].set(xlabel='Time / minutes',ylabel='Concentration / micromolar',title='Full loaded-reactor trajectory');axs[1].legend(fontsize=8)
    for ax in axs:ax.grid(alpha=.2)
    fig.savefig(out/'operating_region_recovery.png',dpi=180);fig.savefig(out/'operating_region_recovery.svg');plt.close(fig)
    fig,axs=plt.subplots(1,2,figsize=(11,4),layout='constrained');axs[0].loglog([float(R(r['b'])) for r in extinction],[float(R(r['b_residual'])) for r in extinction],'o-',label='Residual magnitude');axs[0].loglog([float(R(r['b'])) for r in extinction],[1]*len(extinction),'--',label='Concentration error stays 1');axs[0].set(xlabel='Multiplier pool b / micromolar',ylabel='Value (respective concentration/time units)',title='Rate residual and concentration error as B\nvanishes');axs[0].legend(fontsize=8)
    axs[1].loglog([r['release_rate'] for r in release],[r['unit_weight_storage'] for r in release],'o-',label='Intermediate concentration');axs[1].loglog([r['release_rate'] for r in release],[r['cargo_weighted_storage'] for r in release],'s-',label='Cargo-weighted inventory');axs[1].set(xlabel='Release rate / per minute',ylabel='Storage at input flux 1 micromolar/min',title='Intermediate storage versus product-release\nrate');axs[1].legend(fontsize=8)
    for ax in axs:ax.grid(alpha=.2)
    fig.savefig(out/'conditioning_and_storage.png',dpi=180);fig.savefig(out/'conditioning_and_storage.svg');plt.close(fig)

if __name__=='__main__':main()

Run output

Source symbolic/rational checks: 74 PASS
Candidate lists: [['2', '3'], ['1']] ; mixed/block order: [] / [1]
Loaded reactor: ('11/5', '17/10', '17/10') ; floor limit 29/1100 ; positive-branch limit 39/1100
Conditional concentration error: 1/50 micromolar; ellipsoid: certified_forward_invariant ; initial inside: True
Equilibrium floors do not automatically hold along transients. Private release preserves steady-state geometry, not trajectories or storage. Illustrative rates; Lean not rerun.