Example code
A red-cell metabolic model can balance peroxide consumption on paper without establishing how fast the cell can actually remove peroxide. This example includes the paper's complete 19,620-column red-cell model, exact certificates, a reusable optimization model, and separate finite-carrier and paired-signal modules. Each answers a different part of that question.
Service counts turnover of the Na/K-ATPase ion pump. The peroxide index counts two specified routes: twice the catalase extent plus the glutathione-peroxidase extent. Chemical rows track inventories; enzyme-allocation rows impose algebraic constraints. The code retains both, with every source decimal interpreted exactly for the certificate audit.
The package freshly verifies a 1,580-entry rational feasible witness and every column of three upper certificates. In the restricted one-hour diagnostic medium, one witness at high service is feasible at every lower service requirement. Combined with a global upper bound, it bounds the maximum peroxide-turnover index at required pump service : for every . This is a uniform enclosure, not a grid of solver outputs. The exact witness gives a variation bound below ; the residual service interval is explicitly left unresolved.
Fresh numerical optimizations then show why that objective is insufficient. Maximizing the two-route index can leave almost all imported peroxide in the terminal pool. Minimizing terminal peroxide can divert it into uncounted reactions instead. A further optimization minimizes unnecessary total reaction activity and removes that extra handling. These are different solutions of the declared relaxation, with complete peroxide accounts; none controls peak burden or damage.
The finite-carrier module supplies a missing kinetic restriction. A carrier must be present to cycle between reduced and oxidized states. Its sharp finite-time turnover envelope depends on total amount, initial oxidation state and both rate limits. Saturated controls attain the bound, a fully oxidized pool starts quadratically, and zero carrier forces zero turnover. The plotted coefficients are illustrative, not fitted red-cell parameters.
A separate paired-signal module cancels shared oxidative input before inferring integrated effective reduction. It preserves a common initial peak and treatment fraction, carries observation and intervention errors, and can reject the paper's synthetic inconsistent pair. Passing that test is only a necessary consistency check.
Inputs at the top support new media, inventories, capacity perturbations, carrier controls and signal scenarios. The download includes the full source representation, exact sparse certificates, seven scientific test groups, computed outputs and data provenance. The published plateau remains tied to its frozen premises when numerical inputs change. These results concern a conditional model; they do not establish an external peroxide dose, a physiological tolerance threshold or clinical safety, and Lean is not rerun.
Python source
"""Energy service, peroxide fates, finite carrier kinetics, and paired signals.
Run python example.py --output outputs. All S7 certificates are checked afresh.
Illustrative kinetic coefficients are NOT fitted red-cell parameters.
"""
from __future__ import annotations
import argparse
import csv
from dataclasses import dataclass
from fractions import Fraction as Q
import hashlib
import json
import math
from pathlib import Path
import platform
import numpy as np
from scipy.integrate import solve_ivp
import sympy as sp
import mpmath as mp
from metabolic_model import SourceModel,RationalAudit,LinearRelaxation,OBJECTIVES
# EDITABLE INPUTS -----------------------------------------------------------
HORIZON_HOURS = 1
REQUIRED_SERVICE = 1.0 # mmol/gDW over horizon; source Na/K ATPase
MEDIUM = 'restricted' # original, sulfur_closed, restricted, internal_blocked
INITIAL_CHEMICAL_AMOUNTS = {} # zero is the paper's algebraic diagnostic, not a viable cell
TERMINAL_CHEMICAL_FLOORS = {}
OXIDATION_CAP = 2.0 # a, per hour; illustrative carrier module
REGENERATION_CAP = 1.0 # b, per hour; illustrative
TOTAL_CARRIER = 1.0 # arbitrary unit; results per unit, no invented physiological pool
INITIAL_REDUCED_FRACTIONS = (0., 1/3, 1.)
KINETIC_HORIZON = 3.0
PAIRED_PEAK = 1.0 # synthetic raw signal units
RESIDUAL_INPUT = .03 # one constant r at every time
SIGNAL_HORIZON = 4.0
RAW_OBSERVATIONS = ((Q('0.8'),Q('0.5')),(Q('0.7'),Q('0.6'))) # (untreated, treated)
RAW_ERROR = Q('0.01')
RESIDUAL_INPUT_INTERVAL = (Q(0),Q('0.05'))
LP_TOLERANCE = 1e-8 # numerical objective-fixing tolerance, not certification
RTOL, ATOL = 1e-10, 1e-12
MANUSCRIPT_SHA256 = 'a8cd8efce97007780e7b88fc7f19dc94a31c3e7da20f22b66dbc3f1a0db97774'
# --------------------------------------------------------------------------
class PlateauCertificate:
"""One exact feasible witness plus two global bounds covers an entire interval."""
def __init__(self,primal,turnover_upper,service_upper):
self.M=primal['M'];self.L=primal['Q'];self.U=turnover_upper['bound_without_reverse'];self.UM=service_upper['bound_without_reverse']
if not 0<=self.M<=self.UM or not 0<=self.L<=self.U:raise ValueError('Inconsistent enclosure inputs.')
def query(self,required):
required=Q(required)
if required<0:raise ValueError('Nonnegative required service expected.')
if required>self.UM:return dict(status='infeasible_by_exact_service_bound')
if required<=self.M:return dict(status='uniformly_enclosed',lower=self.L,upper=self.U)
return dict(status='unresolved_feasibility',global_turnover_upper=self.U)
@dataclass(frozen=True)
class CarrierPool:
oxidation: float = OXIDATION_CAP
regeneration: float = REGENERATION_CAP
total: float = TOTAL_CARRIER
def __post_init__(self):
if not all(np.isfinite(v) for v in (self.oxidation,self.regeneration,self.total)) or min(self.oxidation,self.regeneration)<=0 or self.total<0:raise ValueError('Positive finite a,b and nonnegative finite total required.')
def validate(self,g0,time):
if not np.isfinite(g0) or not 0<=g0<=self.total or not np.isfinite(time) or time<0:raise ValueError('Initial reduced pool and nonnegative horizon out of range.')
def bound(self,time,g0):
self.validate(g0,time);a,b,C=self.oxidation,self.regeneration,self.total;k=a+b;s=k*time
# T - (1-exp(-kT))/k loses its leading quadratic term at short times.
if abs(s)<1e-3:
startup=sum((-1)**n*s**n/math.factorial(n) for n in range(2,13))/k
else:startup=time+math.expm1(-s)/k
return a/k*(-math.expm1(-s))*g0+a*b*C/k*startup
def saturated_state(self,time,g0):
self.validate(g0,time);a,b,C=self.oxidation,self.regeneration,self.total
return b*C/(a+b)+(g0-b*C/(a+b))*math.exp(-(a+b)*time)
@property
def harmonic_rate(self):return self.total/(1/self.oxidation+1/self.regeneration)
def required_total(self,turnover,time,theta):
if not np.isfinite(turnover) or turnover<0 or not np.isfinite(theta) or not 0<=theta<=1:raise ValueError('Nonnegative turnover and reduced fraction in [0,1] required.')
per_unit=CarrierPool(self.oxidation,self.regeneration,1).bound(time,theta)
if per_unit==0:return 0. if turnover==0 else None
return turnover/per_unit
def evolve(self,g0,duration,control=None):
self.validate(g0,duration)
if duration<=0:raise ValueError('Use bound() for the zero-time endpoint.')
def rhs(t,v):
oxidation,regeneration=(1.,1.) if control is None else control(t,float(v[0]))
if not 0<=oxidation<=1 or not 0<=regeneration<=1:raise ValueError('Control fractions outside rate envelopes.')
g=v[0];vout=oxidation*self.oxidation*g;uin=regeneration*self.regeneration*(self.total-g)
return [uin-vout,vout,uin]
t=np.linspace(0,duration,301);sol=solve_ivp(rhs,(0,duration),[g0,0,0],t_eval=t,rtol=RTOL,atol=ATOL,method='DOP853')
if not sol.success:raise RuntimeError(sol.message)
return t,sol.y.T
class PairedRecovery:
"""Same initial peak, one residual-input fraction, common effective reduction."""
def __init__(self,peak=PAIRED_PEAK,residual=RESIDUAL_INPUT):
if not np.isfinite(peak) or peak<=0 or not np.isfinite(residual) or not 0<=residual<1:raise ValueError('Positive common peak and r in [0,1) required.')
self.peak=peak;self.r=residual
def transform(self,untreated,treated):return (np.asarray(treated)-self.r*np.asarray(untreated))/((1-self.r)*self.peak)
def evolve(self,duration,input_rate,reduction,discrepancy=None):
def rhs(t,v):
q=input_rate(t);k=reduction(t);extra=0 if discrepancy is None else discrepancy(t,v[1])
if q<0 or k<0:raise ValueError('Synthetic input and reduction rates must be nonnegative.')
return [q-k*v[0],self.r*q-k*v[1]+extra,k,abs(extra)/((1-self.r)*self.peak)]
t=np.linspace(0,duration,401);sol=solve_ivp(rhs,(0,duration),[self.peak,self.peak,0,0],t_eval=t,rtol=RTOL,atol=ATOL,method='DOP853')
if not sol.success:raise RuntimeError(sol.message)
return t,sol.y.T
@staticmethod
def transformed_error(r,peak,error_C,error_D,error_peak,Zmax=1):
r,peak,ec,ed,ep,z=map(Q,(r,peak,error_C,error_D,error_peak,Zmax))
if not 0<=r<1 or peak<=ep or min(ec,ed,ep,z)<0:raise ValueError('Invalid raw-error assumptions.')
return (ec+r*ed)/((1-r)*(peak-ep))+z*ep/(peak-ep)
@staticmethod
def infer_integrated(z_observed,total_error):
z,e=Q(z_observed),Q(total_error)
if e<0:raise ValueError('Nonnegative discrepancy plus observation allowance required.')
lo=z-e;hi=min(Q(1),z+e)
if hi<=0 or lo>hi:return dict(status='incompatible')
mp.iv.dps=45
iv=lambda v:mp.iv.mpf(v.numerator)/v.denominator
lower=math.nextafter(float((-mp.iv.log(iv(hi))).a),-math.inf)
if lo<=0:return dict(status='no_finite_upper_identification',K_lower=max(0.,lower),K_upper=None)
upper=math.nextafter(float((-mp.iv.log(iv(lo))).b),math.inf)
return dict(status='bounded',K_lower=max(0.,lower),K_upper=upper,z_lower=str(lo),z_upper=str(hi))
@staticmethod
def reject_increase(first,second,error,r_interval,peak=1,discrepancy_sum=0):
"""A sufficient exact rejection with ONE common r and peak.
discrepancy_sum bounds the sum of the two normalized intervention errors.
It is subtracted before declaring rejection; passing is only consistency.
"""
d1,c1=map(Q,first);d2,c2=map(Q,second);e=Q(error);rlo,rhi=map(Q,r_interval);P=Q(peak);gamma=Q(discrepancy_sum)
if e<0 or not 0<=rlo<=rhi<1 or P<=0 or gamma<0:raise ValueError('Invalid common-r observation assumptions.')
margins=[c2-c1-2*e+r*(d1-d2-2*e)-(1-r)*P*gamma for r in (rlo,rhi)]
return dict(rejected=min(margins)>0,minimum_numerator_margin=min(margins),common_r=list(map(str,(rlo,rhi))))
def symbolic_checks():
a,b,C,T,t,g0=sp.symbols('a b C T t g0',positive=True);k=a+b
weight=a/k*(1-sp.exp(-k*(T-t)));g=b*C/k+(g0-b*C/k)*sp.exp(-k*t)
bound=a*b*C*T/k+a/k*(g0-b*C/k)*(1-sp.exp(-k*T))
assert sp.simplify(sp.diff(weight,t)-(k*weight-a))==0 and weight.subs(t,T)==0
assert sp.simplify(sp.diff(g,t)-b*(C-g)+a*g)==0
assert sp.simplify(sp.integrate(a*g,(t,0,T))-bound)==0
assert sp.simplify(sp.limit(bound/T,T,sp.oo)-a*b*C/k)==0
quadratic=sp.expand(sp.series(bound.subs(g0,0),T,0,3).removeO())
assert sp.simplify(quadratic-a*b*C*T*T/2)==0
D,Cs,q,r,kk=sp.symbols('D Cs q r kk')
assert sp.expand((r*q-kk*Cs)-r*(q-kk*D)+kk*(Cs-r*D))==0
return dict(adjoint=str(weight),attaining_state=str(g),finite_time_bound=str(bound),oxidized_start=str(quadratic),paired_cancellation='(Cs-rD)prime = -k(Cs-rD) exactly')
def run_lp_diagnostics(model):
lp=LinearRelaxation(model,MEDIUM,HORIZON_HOURS,INITIAL_CHEMICAL_AMOUNTS,TERMINAL_CHEMICAL_FLOORS)
baseline=lp.solve(required_service=REQUIRED_SERVICE)
if not baseline['finished']:return dict(baseline=baseline,scope='No successful optimum; later cleanup omitted.')
fixed={name:(baseline[key]-LP_TOLERANCE,baseline[key]+LP_TOLERANCE) for name,key in [('service','M'),('turnover','Q')]}
cleaned=lp.solve('peroxide',REQUIRED_SERVICE,fixed)
parsimonious=None
if cleaned['finished']:
fixed['peroxide']=(cleaned['terminal_peroxide']-lp.initial[model.ri['M_h2o2_c']]-LP_TOLERANCE,cleaned['terminal_peroxide']-lp.initial[model.ri['M_h2o2_c']]+LP_TOLERANCE)
parsimonious=lp.solve('turnover',REQUIRED_SERVICE,fixed,parsimony=True)
sensitivity=[dict(fraction=f,result=lp.solve(required_service=REQUIRED_SERVICE,capacity_fractions={'R_GTHOx':f})) for f in (0.,.5,1.)]
return dict(baseline=baseline,minimum_terminal=cleaned,parsimonious=parsimonious,GTHOx_sensitivity=sensitivity,
scope='Fresh floating-point LPs on the configured medium; intermediate-burden and damage constraints absent. Absolute-flux parsimony mixes source column scales and is a declared selection rule.')
def main():
ap=argparse.ArgumentParser(description=__doc__);ap.add_argument('--output',type=Path,default=Path('outputs'));ap.add_argument('--skip-lp',action='store_true');args=ap.parse_args();out=args.output;out.mkdir(parents=True,exist_ok=True)
source=Path(__file__).parent;model,certs,manifest=SourceModel.load(source);audit=RationalAudit(model)
primal=audit.primal(certs['primal']);duals={name:audit.upper(certs[name],OBJECTIVES[name],medium='original' if name=='joint' else 'restricted') for name in ('joint','turnover','service')}
assert duals['joint']['interval_cost']==Q(1283461673726143619667945455169487587,50452437500000000000000000000000000)<Q('25.439042')
plateau=PlateauCertificate(primal,duals['turnover'],duals['service']);endpoint=Q('1.0345238')
assert plateau.M>=endpoint and plateau.L>=Q('6.50048030') and plateau.U<=Q('6.50048050')
assert plateau.U-plateau.L<Q('0.000000193')
exact_plateau=dict(M_witness=plateau.M,Q_witness=plateau.L,turnover_upper=plateau.U,service_upper=plateau.UM,
rounded_service_endpoint=endpoint,variation_bound=plateau.U-plateau.L,rounded_endpoint_gap=plateau.UM-endpoint,witness_endpoint_gap=plateau.UM-plateau.M,
configured_query=plateau.query(Q(str(REQUIRED_SERVICE))),scope='Frozen restricted medium at one hour, zero chemical stocks/floors; independent of editable numerical scenarios.')
lp=None if args.skip_lp else run_lp_diagnostics(model)
pool=CarrierPool();kinetic_rows=[];kinetic_summaries=[]
for theta in INITIAL_REDUCED_FRACTIONS:
g0=theta*pool.total;t,y=pool.evolve(g0,KINETIC_HORIZON)
exact=np.array([pool.bound(float(ti),g0) for ti in t]);diff=float(max(abs(y[:,1]-exact)))
assert diff<1e-8 and np.max(abs(g0+y[:,2]-y[:,1]-y[:,0]))<1e-9
kinetic_rows.extend((theta,ti,*row,bnd) for ti,row,bnd in zip(t,y,exact))
kinetic_summaries.append(dict(theta=theta,turnover_at_one_hour=pool.bound(1.,g0),solver_error=diff,
total_needed_for_witness_GTHP=pool.required_total(float(primal['extent_GTHP']),1.,theta)))
signal=PairedRecovery();input_rate=lambda t:.3+.2*math.sin(2*t)**2;reduction=lambda t:.4+.2*math.cos(t)**2
t,y=signal.evolve(SIGNAL_HORIZON,input_rate,reduction);z=signal.transform(y[:,0],y[:,1]);error=float(max(abs(z-np.exp(-y[:,2]))));assert error<1e-8
tt,yy=signal.evolve(SIGNAL_HORIZON,input_rate,reduction,lambda t,cs:.01*math.sin(t));zz=signal.transform(yy[:,0],yy[:,1]);assert np.all(abs(zz-np.exp(-yy[:,2]))<=yy[:,3]+1e-9)
paired_rows=[(ti,*row,zi,damaged,allowance) for ti,row,zi,damaged,allowance in zip(t,y,z,zz,yy[:,3])]
rejection=PairedRecovery.reject_increase(*RAW_OBSERVATIONS,RAW_ERROR,RESIDUAL_INPUT_INTERVAL)
result=dict(manuscript_sha256=MANUSCRIPT_SHA256,source_hashes=manifest,source_dimensions=dict(rows=len(model.rows),columns=len(model.columns),chemical=len(model.chem),auxiliary=len(model.aux)),
exact_primal=primal,exact_duals=duals,plateau=exact_plateau,bypass=audit.bypass(),symbolic=symbolic_checks(),numerical_lp=lp,
carrier=kinetic_summaries,paired=dict(identity_solver_error=error,synthetic_rejection=rejection,
inference=PairedRecovery.infer_integrated(Q('0.5'),Q('0.02')),raw_error_bound=PairedRecovery.transformed_error(Q('0.03'),1,Q('.01'),Q('.01'),Q('.01'))),
scope='Conditional mathematical model. Exact source verification, numerical LPs, illustrative kinetic envelopes and synthetic signal tests are distinct. No clinical threshold, donor fit or Lean replay.')
def dump(name,data):(out/name).write_text(json.dumps(data,indent=2,default=str)+'\n')
def table(name,header,rows):
with (out/name).open('w',newline='') as f:w=csv.writer(f);w.writerow(header);w.writerows(rows)
dump('results.json',result);table('carrier.csv',['theta','time_h','reduced','oxidation_turnover','regeneration_turnover','sharp_upper'],kinetic_rows)
table('paired_signals.csv',['time','untreated','treated','integrated_k','zero_discrepancy','z','z_with_discrepancy','Gamma'],paired_rows)
table('dual_costs.csv',['certificate','reaction','residual','interval_cost'],[(name,v['reaction'],v['residual'],v['interval_cost']) for name,d in duals.items() for v in d['costs']])
lines=[f'All {len(model.columns)} source columns audited exactly: primal feasible; joint, turnover and service bounds recomputed.',
f'Plateau: exact witness Q={float(plateau.L):.12f}; upper={float(plateau.U):.12f}; service certified through {endpoint}.',
f'Variation <= {float(plateau.U-plateau.L):.9g}; rounded-endpoint service gap {float(plateau.UM-endpoint):.9g}.',
'The sulfur bypass cancels all chemical carriers but retains eight auxiliary terms.',
f'Paired synthetic contradiction rejected: {rejection["rejected"]}; identity solver error {error:.3g}.',
'Finite carrier kinetics are illustrative; source feasibility is not a kinetic or physiological tolerance claim.']
(out/'console.txt').write_text('\n'.join(lines)+'\n');print('\n'.join(lines))
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
fig,axs=plt.subplots(1,2,figsize=(10,4),layout='constrained');x=np.linspace(0,float(endpoint),100)
offset=Q('6.50048030');axs[0].fill_between(x,float((plateau.L-offset)*10**7),float((plateau.U-offset)*10**7),alpha=.25,color='#246d91')
axs[0].axhline(float((plateau.L-offset)*10**7),color='#246d91',label='Exact feasible witness');axs[0].axhline(float((plateau.U-offset)*10**7),color='#bd5a24',label='Exact upper certificate')
axs[0].set(xlabel='Required service (mmol/gDW)',ylabel='10^7 × (turnover − 6.50048030)',title='Certified peroxide-turnover band across\npump service');axs[0].legend(fontsize=8)
if lp and lp.get('parsimonious') and lp['parsimonious']['finished'] and lp['minimum_terminal']['finished']:
stages=[lp['baseline'],lp['minimum_terminal'],lp['parsimonious']];xx=np.arange(3);bottom=np.zeros(3)
for key,label in [('Q','Counted turnover'),('peroxide_other_consumption','Other consumption'),('terminal_peroxide','Terminal pool')]:
vals=np.array([s[key] for s in stages]);axs[1].bar(xx,vals,bottom=bottom,label=label);bottom+=vals
axs[1].set_xticks(xx,['Max turnover','Min terminal','Then parsimonious']);axs[1].legend(fontsize=7)
axs[1].set(ylabel='Peroxide amount (mmol/gDW)',title='Peroxide allocation under three\noptimization objectives')
for ax in axs:ax.grid(axis='y',alpha=.2)
fig.savefig(out/'envelope.png',dpi=180);fig.savefig(out/'envelope.svg');plt.close(fig)
fig,axs=plt.subplots(1,2,figsize=(10,4),layout='constrained')
display_scale=pool.total if pool.total else 1.
for theta in INITIAL_REDUCED_FRACTIONS:
d=np.array([r[1:] for r in kinetic_rows if r[0]==theta]);axs[0].plot(d[:,0],d[:,-1]/display_scale,label=f'Initial reduced fraction {theta:.3g}')
axs[0].plot(t:=np.linspace(0,KINETIC_HORIZON,100),pool.harmonic_rate*t/display_scale,'k:',label='Harmonic rate × time');axs[0].set(xlabel='Time (hours)',ylabel='Oxidation turnover / total carrier' if pool.total else 'Oxidation turnover (zero carrier)',title='Carrier turnover by initial oxidation state');axs[0].legend(fontsize=7)
data=np.array(paired_rows);axs[1].plot(data[:,0],data[:,5],label='Transformed paired signal z');axs[1].plot(data[:,0],np.exp(-data[:,3]),'--',label='exp(−integrated k)');axs[1].fill_between(data[:,0],data[:,6]-data[:,7],data[:,6]+data[:,7],alpha=.2,label='Discrepant z ± Gamma');axs[1].set(xlabel='Synthetic time (hours)',ylabel='Normalized transformed signal',title='Paired-signal estimate of integrated\nreduction');axs[1].legend(fontsize=7)
for ax in axs:ax.grid(alpha=.2)
fig.savefig(out/'kinetics.png',dpi=180);fig.savefig(out/'kinetics.svg');plt.close(fig)
digest=lambda p:hashlib.sha256(p.read_bytes()).hexdigest()
dump('run_metadata.json',dict(manuscript_sha256=MANUSCRIPT_SHA256,source_sha256=digest(Path(__file__)),supporting_source_sha256={n:digest(source/n) for n in ['metabolic_model.py','freeze_source.py','s7_model.json','certificates.json','source_manifest.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'}))
if __name__=='__main__':main()
Run output
All 19620 source columns audited exactly: primal feasible; joint, turnover and service bounds recomputed. Plateau: exact witness Q=6.500480304680; upper=6.500480497416; service certified through 5172619/5000000. Variation <= 1.9273609e-07; rounded-endpoint service gap 9.44791832e-08. The sulfur bypass cancels all chemical carriers but retains eight auxiliary terms. Paired synthetic contradiction rejected: True; identity solver error 7.8e-12. Finite carrier kinetics are illustrative; source feasibility is not a kinetic or physiological tolerance claim.