Example code
Can selective degradation create two positive steady states in a minimal autocatalytic core? For the paper's Type II and Type V families, there is at most one, even when some or all loss coefficients vanish. The example implements their literal reversible mass-action models, including nonlinear reverse reactions and intermediate species along one-to-one reaction chains.
Editable inputs define the topology, losses and a constructive operating state. Reusable components expose exact path reduction, full state reconstruction, stationary-flow certificates and numerical continuation at fixed reversible rates. Inputs are dimensionless examples, without chemical calibration.
Two mechanisms explain the result. In Type V, eliminating intermediate species along passive reversible chains gives three equations that rule out two distinct positive solutions by comparing their concentration ratios. Reverse-fork kinetics supplies the strict inequality, so the proof does not need positive degradation. The code supports all eight linked/collapsed pair patterns and arbitrary finite path lengths.
Path elimination must retain leakage. For with all four edge constants and equal to one, both transmissions and both endpoint leakage coefficients equal . Ignoring those leaks changes the steady-state problem. This reduction is exact at stationarity; the expanded model remains available for transient dynamics.
For a separated Type II skeleton, the example checks the stationary factorization of the scaled Jacobian and the admissibility condition . Its exact six-species example gives
All 64 principal minors in the determinant expansion are checked with rational arithmetic. Zero losses are allowed throughout, and backward net currents are handled by shifting the admissible one-way flows. The determinant bound is specific to the retained separated skeleton; it is not a bound on a smallest singular value or a claim about the expanded model.
The plotted concentration branch and eigenvalues are numerical diagnostics. Uniqueness does not establish existence or stability at arbitrary loss. The paper proves existence and local stability near zero degradation, while global continuation and possible Hopf bifurcations remain outside that conclusion.
Seven scientific test groups check the source equations, all eight Type V patterns, internally degraded paths, exact certificates and independent numerical root agreement. The package does not rerun Lean; its README separates the paper's universal theorem from the finite checks performed here.
Python source
"""Exact stationary reductions and reusable mass-action models for mixed degradation.
Run: python example.py --output outputs
Inputs are dimensionless constructive examples, not fitted chemical constants.
"""
from __future__ import annotations
# EDITABLE INPUTS: zero is a valid loss; all reversible rates must stay positive.
TYPE_II_ACTIVE_WEIGHTS = ((2,), (1,), (3,)) # one or more active edges per gap
TYPE_II_LOSS_BY_NAME = {"F0": "1/5", "T1": "0", "F1": "0", "T2": "1/3", "F2": "0", "T0": "0"}
TARGET_CONCENTRATION = "1" # construct rates with this exact positive state
MIN_ONE_WAY_FLOW = "1" # positive margin above the admissibility boundary
LOSS_SWEEP = tuple(i / 40 for i in range(41)) # fixed rates, d = scale * target d
TYPE_V_PATH_EDGES = (3, 0, 2) # 0 identifies U_i=S_i; positive values give unit paths
TYPE_V_INTERNAL_LOSS = "1/4"
TYPE_V_BASE_LOSS = ("0", "1/5", "0")
TYPE_V_FORK_LOSS = ("0", "0", "1/3")
MAX_PRINCIPAL_MINORS = 1024
MANUSCRIPT_SHA256 = "6d260d5d6596b099430287d32563d91ebf3f22ccbee5c5b9a889b54470ccf3ed"
import argparse
import csv
from dataclasses import dataclass, replace
import hashlib
import itertools
import json
from pathlib import Path
import platform
import numpy as np
import sympy as sp
from scipy.optimize import least_squares
Q = sp.Rational
def vector(values):
return sp.Matrix([Q(v) for v in values])
def positive(values, name, allow_zero=False):
if any(v < 0 if allow_zero else v <= 0 for v in values):
raise ValueError(f"{name} must be {'nonnegative' if allow_zero else 'positive'}.")
@dataclass(frozen=True)
class UnitPath:
"""Exact stationary two-port; only one-to-one reactions may be put here."""
forward: tuple
reverse: tuple
internal_loss: tuple
def __post_init__(self):
for field in ('forward', 'reverse', 'internal_loss'):
object.__setattr__(self, field, tuple(Q(v) for v in getattr(self, field)))
if not self.forward or len(self.forward) != len(self.reverse) or len(self.internal_loss) != len(self.forward)-1:
raise ValueError('A path needs m+1 reversible edges and m internal losses.')
positive(self.forward + self.reverse, 'Path rates')
positive(self.internal_loss, 'Internal losses', True)
def compress(self):
c, b, left, right = self.forward[0], self.reverse[0], Q(0), Q(0)
for a, beta, d in zip(self.forward[1:], self.reverse[1:], self.internal_loss):
pivot = b + right + a + d
c, b, left, right = c*a/pivot, b*beta/pivot, left+c*(right+d)/pivot, beta*(right+d)/pivot
return c, b, left, right
def reconstruct(self, X, Y):
"""Solve the independent tridiagonal balance, not the compression recursion."""
X, Y = Q(X), Q(Y)
positive((X,Y), 'Endpoints')
m = len(self.internal_loss)
if not m:
return sp.Matrix([X,Y])
matrix, rhs = sp.zeros(m), sp.zeros(m,1)
for k in range(m):
matrix[k,k] = self.reverse[k] + self.forward[k+1] + self.internal_loss[k]
if k: matrix[k,k-1] = -self.forward[k]
else: rhs[k] += self.forward[0]*X
if k < m-1: matrix[k,k+1] = -self.reverse[k+1]
else: rhs[k] += self.reverse[-1]*Y
return sp.Matrix([X, *matrix.inv()*rhs, Y])
@dataclass(frozen=True)
class SourceNetwork:
"""Column r consumes one X_r and produces column r of P."""
names: tuple
P: sp.ImmutableMatrix
family: str = 'unclassified'
def __post_init__(self):
object.__setattr__(self, 'P', sp.ImmutableMatrix(self.P))
if len(set(self.names)) != len(self.names) or self.P.shape != (len(self.names),)*2:
raise ValueError('Unique species and square product matrix required.')
if any(v.is_integer is not True or v < 0 for v in self.P):
raise ValueError('Product multiplicities must be nonnegative integers.')
if self.N.det() == 0: raise ValueError('This stationary-flow interface requires invertible N.')
@property
def N(self): return self.P-sp.eye(len(self.names))
def monomials(self, x):
return sp.Matrix([sp.prod(x[i]**self.P[i,r] for i in range(len(x))) for r in range(len(x))])
def realize(self, x, d, margin=MIN_ONE_WAY_FLOW):
"""Construct exact rates at a chosen root; permits backward net currents."""
x, d = vector(x), vector(d)
if len(x) != len(self.names) or len(d) != len(x): raise ValueError('Wrong vector length.')
positive(x, 'Concentrations'); positive(d, 'Losses', True); positive((Q(margin),), 'Flow margin')
currents = self.N.inv()*sp.diag(*d)*x
q = vector([max(0,-j)+Q(margin) for j in currents]); p = q+currents
return MassActionModel(self, tuple(p[i]/x[i] for i in range(len(x))),
tuple(q[i]/m for i,m in enumerate(self.monomials(x))), tuple(d))
def separated_type_ii(weights=TYPE_II_ACTIVE_WEIGHTS):
"""Arbitrary nonempty active gaps; each terminal return is a unit edge."""
if len(weights) < 3 or any(not gap or any(int(m)!=m or m<1 for m in gap) for gap in weights):
raise ValueError('At least three forks, nonempty gaps, positive integer active weights.')
names=[]; sectors=[]
for j,gap in enumerate(weights):
nodes=[f'F{j}', *[f'G{j}_{k}' for k in range(1,len(gap))], f'T{(j+1)%len(weights)}']
names.extend(nodes); sectors.append(nodes)
index={s:i for i,s in enumerate(names)}; P=sp.zeros(len(names))
for j,(gap,nodes) in enumerate(zip(weights,sectors)):
for k,m in enumerate(gap): P[index[nodes[k+1]],index[nodes[k]]] = m
P[index[f'T{j}'],index[f'F{j}']] += 1
P[index[f'F{(j+1)%len(weights)}'],index[nodes[-1]]] = 1
return SourceNetwork(tuple(names), P, 'separated_type_ii_skeleton')
def coincident_type_ii(weights=(2,1,3)):
if len(weights)!=3 or any(int(m)!=m or m<1 for m in weights): raise ValueError('Three positive integer weights required.')
P=sp.zeros(3)
for i,m in enumerate(weights): P[(i+1)%3,i]=m; P[(i+2)%3,i]=1
return SourceNetwork(('X0','X1','X2'), P, 'coincident_type_ii')
@dataclass(frozen=True)
class MassActionModel:
source: SourceNetwork
forward: tuple
reverse: tuple
loss: tuple
def __post_init__(self):
for name in ('forward','reverse','loss'):
object.__setattr__(self,name,tuple(Q(v) for v in getattr(self,name)))
if len(getattr(self,name)) != len(self.source.names): raise ValueError('Wrong rate vector length.')
positive(self.forward+self.reverse, 'Reversible rates'); positive(self.loss,'Losses',True)
def flows(self,x):
return sp.diag(*self.forward)*sp.Matrix(x), sp.diag(*self.reverse)*self.source.monomials(x), sp.diag(*self.loss)*sp.Matrix(x)
def field(self,x):
p,q,e=self.flows(x)
return self.source.N*(p-q)-e
def scaled_jacobian(self,x):
p,q,e=self.flows(x); N=self.source.N
return N*sp.diag(*p)-N*sp.diag(*q)*self.source.P.T-sp.diag(*e)
def detailed_balance(self):
"""Positive zero-loss root for these reversible constants; not the mixed-loss root."""
N=np.array(self.source.N,float)
return np.exp(np.linalg.solve(N.T,np.log(np.array(self.forward,float)/np.array(self.reverse,float))))
def numeric(self,x):
x=np.asarray(x,float)
if np.any(x<=0) or not np.all(np.isfinite(x)): raise ValueError('Numerical state must be finite and positive.')
P=np.array(self.source.P,float); N=P-np.eye(len(x))
p=np.array(self.forward,float)*x; q=np.array(self.reverse,float)*np.exp([email protected](x)); e=np.array(self.loss,float)*x
[email protected](p)[email protected](q)@P.T-np.diag(e)
return N@(p-q)-e, H
def continuation(self,scales,initial=None):
"""Log-coordinate local solves. Failure is unresolved, never nonexistence."""
guess=self.detailed_balance() if initial is None else np.array(initial,float)
rows=[]
for scale in scales:
if scale<0: raise ValueError('Physical loss scales must be nonnegative.')
model=replace(self,loss=tuple(Q(str(scale))*d for d in self.loss))
# Dividing f by x avoids accepting the origin as a tiny-residual root.
def residual(z):
x=np.exp(z); f,_=model.numeric(x); return f/x
def jac(z):
x=np.exp(z); f,H=model.numeric(x); return H/x[:,None]-np.diag(f/x)
try:
fit=least_squares(residual,np.log(guess),jac=jac,bounds=(-25,25),xtol=1e-13,ftol=1e-13,gtol=1e-13,max_nfev=1000)
except (ValueError, FloatingPointError, np.linalg.LinAlgError) as exc:
rows.append({'scale':float(scale),'status':'unresolved','relative_residual':None,'reason':str(exc)}); break
x=np.exp(fit.x); error=float(np.max(np.abs(residual(fit.x))))
if not fit.success or error>1e-8:
rows.append({'scale':float(scale),'status':'unresolved','relative_residual':error}); break
f,H=model.numeric(x); physical=H/x[None,:]
rows.append({'scale':float(scale),'status':'numerical_root','x':x.tolist(),
'relative_residual':error,'spectral_abscissa':float(max(np.linalg.eigvals(physical).real)),
'det_minus_H':float(np.linalg.det(-H))})
guess=x
return rows
def expand_unit_edge(self,start,end,path):
"""Replace a literal unit reaction; active m>1 or fork edges are rejected."""
i=self.source.names.index(start); j=self.source.names.index(end)
if self.source.P[:,i] != sp.eye(len(self.loss))[:,j]: raise ValueError('Only a standalone unit edge may be expanded.')
extra=tuple(f'{start}_path{k}' for k in range(1,len(path.forward)))
if set(extra)&set(self.source.names): raise ValueError('Path names collide.')
names=self.source.names+extra; n=len(names); P=sp.zeros(n)
P[:len(self.loss),:len(self.loss)]=self.source.P; P[:,i]=sp.zeros(n,1)
route=[i,*range(len(self.loss),n),j]; a=list(self.forward)+[Q(1)]*len(extra); b=list(self.reverse)+[Q(1)]*len(extra)
for k,(u,v) in enumerate(zip(route,route[1:])):
P[v,u]=1; a[u]=path.forward[k]; b[u]=path.reverse[k]
model=MassActionModel(SourceNetwork(names,P,'expanded_unit_path'),tuple(a),tuple(b),self.loss+path.internal_loss)
c,beta,left,right=path.compress(); ca=list(self.forward); cb=list(self.reverse); cd=list(self.loss)
ca[i]=c; cb[i]=beta; cd[i]+=left; cd[j]+=right
return model,replace(self,forward=tuple(ca),reverse=tuple(cb),loss=tuple(cd))
class StationaryCertificate:
"""Exact identities on a stationary flow triple; no division by a loss."""
def __init__(self,model,x):
self.model=model; self.x=vector(x); positive(self.x,'State')
if model.field(self.x)!=sp.zeros(len(x),1): raise ValueError('Exact stationarity required; an approximate root is not a certificate.')
self.p,self.q,self.e=model.flows(self.x); self.N=model.source.N; V=self.N.inv(); self.J=V*self.e
self.shift=vector([max(0,-j) for j in self.J]); self.t=self.q-self.shift
self.B=(sp.diag(*self.J)-V*sp.diag(*self.e))*V.T
self.A=sp.diag(*self.shift)-self.B; self.H=model.scaled_jacobian(self.x)
if -self.H != self.N*(sp.diag(*self.q)-self.B)*self.N.T: raise AssertionError('Factorization failed.')
def principal_audit(self,budget=MAX_PRINCIPAL_MINORS):
n=len(self.x)
if 2**n>budget: return {'status':'unresolved_budget','required_minors':2**n}
minors=[]; expansion=Q(0)
for mask in range(2**n):
idx=[i for i in range(n) if mask>>i&1]
value=self.A.extract(idx,idx).det() if idx else Q(1)
minors.append(value); expansion+=value*sp.prod(self.t[i] for i in range(n) if i not in idx)
if expansion!=(self.A+sp.diag(*self.t)).det(): raise AssertionError('Multiaffine expansion failed.')
lower=self.N.det()**2*sp.prod(self.t); actual=(-self.H).det()
return {'status':'all_principal_minors_nonnegative' if min(minors)>=0 else 'negative_principal_minor',
'count':len(minors),'zero_minors':sum(m==0 for m in minors),'minimum':str(min(minors)),
'det_minus_H':str(actual),'lower_bound_from_nonnegative_minors':str(lower) if min(minors)>=0 else None,
'ratio_to_lower_bound':str(actual/lower) if min(minors)>=0 else None,
'scope':'Exact certificate for this triple. The paper proves the universal bound for separated Type II skeletons only.'}
@dataclass(frozen=True)
class TypeVTopology:
"""All eight identifications and arbitrary finite linked paths."""
edges: tuple
def __post_init__(self):
if len(self.edges)!=3 or any(int(n)!=n or n<0 for n in self.edges): raise ValueError('Three nonnegative integer edge counts required.')
def names_for_pair(self,i):
return (f'U{i}',) if self.edges[i]==0 else (f'U{i}',*[f'Z{i}_{k}' for k in range(1,self.edges[i])],f'S{i}')
def source(self):
names=tuple(s for i in range(3) for s in self.names_for_pair(i)); index={s:i for i,s in enumerate(names)}; P=sp.zeros(len(names))
for i in range(3):
route=self.names_for_pair(i)
for a,b in zip(route,route[1:]): P[index[b],index[a]]=1
for j in range(3):
if j!=i: P[index[f'U{j}'],index[route[-1]]]=1
return SourceNetwork(names,P,'type_v')
def reduce(self,model):
if model.source.names!=self.source().names or model.source.P!=self.source().P: raise ValueError('Model does not match this Type V topology.')
return ReducedTypeV(self,model)
class ReducedTypeV:
def __init__(self,topology,model):
self.topology=topology; self.model=model; self.sectors=[]
self.alpha=[]; self.beta=[]; self.lam=[]; self.mu=[]
for i in range(3):
route=topology.names_for_pair(i); ids=[model.source.names.index(s) for s in route]; last=ids[-1]
f,g=model.forward[last],model.reverse[last]; d=model.loss[ids[0]]
if len(ids)==1:
self.sectors.append(None); alpha,beta,lam,mu=f,g,d,Q(0)
else:
path=UnitPath(tuple(model.forward[k] for k in ids[:-1]),tuple(model.reverse[k] for k in ids[:-1]),tuple(model.loss[k] for k in ids[1:-1]))
a,b,left,right=path.compress(); d+=left; h=model.loss[last]+right; D=b+f+h
self.sectors.append((path,a,g,D)); alpha,beta,lam,mu=f*a/D,g*(b+h)/D,d+h*a/D,h*g/D
self.alpha.append(alpha); self.beta.append(beta); self.lam.append(lam); self.mu.append(mu)
def G(self,u):
return sp.Matrix([2*self.beta[i]*u[(i+1)%3]*u[(i+2)%3]/u[i]+self.lam[(i+1)%3]*u[(i+1)%3]/u[i]+self.lam[(i+2)%3]*u[(i+2)%3]/u[i]+self.mu[(i+1)%3]*u[(i+2)%3]+self.mu[(i+2)%3]*u[(i+1)%3] for i in range(3)])
def residual(self,u): return sp.diag(*u)*(2*sp.Matrix(self.alpha)-self.G(u))
def reconstruct(self,u):
u=vector(u); positive(u,'Base concentrations'); values=[]
for i,sector in enumerate(self.sectors):
if sector is None: values.append(u[i])
else:
path,a,g,D=sector; S=(a*u[i]+g*u[(i+1)%3]*u[(i+2)%3])/D
values.extend(path.reconstruct(u[i],S))
return sp.Matrix(values)
def tangent(self,u):
z=sp.symbols('u0:3',positive=True)
return self.G(z).jacobian(z).subs(dict(zip(z,u)))*sp.diag(*u)
def ratio_witness(self,u,v):
"""A strict row defect shows unequal positive u,v cannot both be roots."""
u,v=vector(u),vector(v); positive(u,'u');positive(v,'v'); r=[u[i]/v[i] for i in range(3)]
if all(a==1 for a in r): return {'status':'identical'}
upward=sorted(r)[1]>=1; i=(min if upward else max)(range(3),key=lambda k:r[k]); j,k=(i+1)%3,(i+2)%3
terms=[2*self.beta[i]*v[j]*v[k]*(r[j]*r[k]-r[i]),self.lam[j]*v[j]*(r[j]-r[i]),self.lam[k]*v[k]*(r[k]-r[i]),self.mu[j]*v[i]*v[k]*r[i]*(r[k]-1),self.mu[k]*v[i]*v[j]*r[i]*(r[j]-1)]
direction=1 if upward else -1
if direction*terms[0]<=0 or any(direction*t<0 for t in terms): raise AssertionError('Ratio sign structure failed.')
return {'status':'cannot_both_be_stationary','row':i,'terms':list(map(str,terms)),'strict_defect':str(sum(terms))}
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)
source=separated_type_ii(); x=vector([TARGET_CONCENTRATION]*len(source.names)); d=vector([TYPE_II_LOSS_BY_NAME.get(n,'0') for n in source.names])
model=source.realize(x,d); cert=StationaryCertificate(model,x); audit=cert.principal_audit(); sweep=model.continuation(LOSS_SWEEP)
topology=TypeVTopology(TYPE_V_PATH_EDGES); vs=topology.source(); vd=[]
for i in range(3):
route=topology.names_for_pair(i)
if len(route)==1: vd.append(Q(TYPE_V_BASE_LOSS[i])+Q(TYPE_V_FORK_LOSS[i]))
else: vd.extend([Q(TYPE_V_BASE_LOSS[i]),*[Q(TYPE_V_INTERNAL_LOSS)]*(len(route)-2),Q(TYPE_V_FORK_LOSS[i])])
vx=vector([TARGET_CONCENTRATION]*len(vs.names)); vm=vs.realize(vx,vd); reduced=topology.reduce(vm); u=vector([TARGET_CONCENTRATION]*3)
if reduced.reconstruct(u)!=vx or reduced.residual(u)!=sp.zeros(3,1): raise AssertionError('Type V exact reconstruction failed.')
# Internal-only loss survives as endpoint leakage; no loss division is used.
path=UnitPath((1,1),(1,1),(1,)); path_values=path.reconstruct(2,3)
coincident=coincident_type_ii(); cm=coincident.realize([1,1,1],[0,Q(1,3),0])
result={'scope':'Dimensionless exact stationary examples and numerical local continuation. Universal uniqueness is the paper theorem, not a solver inference; Lean is not rerun.',
'type_ii':{'species':source.names,'product_matrix':[[int(v) for v in row] for row in source.P.tolist()],
'forward':list(map(str,model.forward)),'reverse':list(map(str,model.reverse)),'loss':list(map(str,d)),
'exact_state':list(map(str,x)),'currents':list(map(str,cert.J)),'admissibility_shift':list(map(str,cert.shift)),
'min_one_way_flows':list(map(str,cert.t)),'certificate':audit,'numerical_continuation':sweep},
'type_v':{'edges':TYPE_V_PATH_EDGES,'species':vs.names,'exact_state':list(map(str,vx)),
'forward':list(map(str,vm.forward)),'reverse':list(map(str,vm.reverse)),'loss':list(map(str,vd)),
'alpha':list(map(str,reduced.alpha)),'beta':list(map(str,reduced.beta)),'lambda':list(map(str,reduced.lam)),'mu':list(map(str,reduced.mu)),
'tangent_determinant':str(reduced.tangent(u).det()),'ratio_witness':reduced.ratio_witness([Q(3,2),1,Q(1,2)],u)},
'internal_only_loss':{'compression':list(map(str,path.compress())),'endpoints_2_3_reconstruction':list(map(str,path_values)),
'left_current':str(2-path_values[1]),'right_current':str(path_values[1]-3),'leakage':str(path_values[1])},
'coincident_branch':{'det_N':str(coincident.N.det()),'inverse':[[str(a) for a in row] for row in coincident.N.inv().tolist()],
'exact_root_verified':cm.field([1,1,1])==sp.zeros(3,1),'separated_determinant_bound_claimed':False}}
def csv_file(name,headers,rows):
with (out/name).open('w',newline='') as f:
w=csv.writer(f);w.writerow(headers);w.writerows(rows)
csv_file('loss_continuation.csv',['loss_scale','status',*source.names,'relative_residual','spectral_abscissa','det_minus_H'],
[[r['scale'],r['status'],*r.get('x',['']*len(x)),r['relative_residual'],r.get('spectral_abscissa',''),r.get('det_minus_H','')] for r in sweep])
csv_file('source_reactions.csv',['model','source','products','forward','reverse','degradation'],
[[name,s,' + '.join(f'{network.P[i,r]} {network.names[i]}' for i in range(len(network.names)) if network.P[i,r]),str(m.forward[r]),str(m.reverse[r]),str(m.loss[r])]
for name,m in [('Type II',model),('Type V',vm)] for network in [m.source] for r,s in enumerate(network.names)])
path_rows=[]
for loss in [Q(k,10) for k in range(31)]:
p=UnitPath((1,1),(1,1),(loss,)); c,b,l,r=p.compress(); vals=p.reconstruct(2,3)
path_rows.append([str(loss),str(c),str(b),str(l),str(r),str(vals[1]),str((c+l)*2-b*3),str(c*2-(b+r)*3)])
csv_file('passive_path.csv',['internal_loss','c','beta','lambda_left','lambda_right','internal_concentration','left_current','right_current'],path_rows)
(out/'results.json').write_text(json.dumps(result,indent=2)+'\n')
lines=[f'Type II: exact positive stationary state; {audit["count"]} principal minors checked.',
f'det(-H)={audit["det_minus_H"]}; lower bound={audit["lower_bound_from_nonnegative_minors"]}.',
f'Type V: {len(vs.names)} literal species reduce to three base variables with exact reconstruction.',
'An internal-only loss gives two endpoint leaks of 1/3, even when both direct endpoint losses vanish.',
f'Fixed-rate continuation: {sum(r["status"]=="numerical_root" for r in sweep)}/{len(LOSS_SWEEP)} numerical roots.',
'At most one positive root is the theorem; existence, stability at arbitrary loss and global continuation are not asserted.']
(out/'console.txt').write_text('\n'.join(lines)+'\n'); print('\n'.join(lines))
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
fig,axes=plt.subplots(1,2,figsize=(10,4.2),layout='constrained'); roots=[r for r in sweep if r['status']=='numerical_root']
for i,n in enumerate(source.names): axes[0].plot([r['scale'] for r in roots],[r['x'][i] for r in roots],label=n)
axes[0].set(xlabel='Scale of selective degradation (fixed reversible rates)',ylabel='Positive stationary concentration',title='Numerical local continuation'); axes[0].legend(ncol=2,fontsize=8)
axes[1].plot([r['scale'] for r in roots],[r['spectral_abscissa'] for r in roots],color='#b95024'); axes[1].axhline(0,color='black',ls='--',lw=.8)
axes[1].set(xlabel='Scale of selective degradation',ylabel='Largest real part of physical Jacobian eigenvalues',title='Stability diagnostic for these computed states')
for ax in axes: ax.grid(alpha=.2)
fig.savefig(out/'continuation.png',dpi=180); fig.savefig(out/'continuation.svg'); plt.close(fig)
rows=np.array([[float(Q(v)) for v in row] for row in path_rows]); fig,axes=plt.subplots(1,2,figsize=(10,4.2),layout='constrained')
axes[0].plot(rows[:,0],rows[:,1],label='Transmission c = beta'); axes[0].plot(rows[:,0],rows[:,3],label='Leakage lambda_left = lambda_right')
axes[0].set(xlabel='Loss at internal species Z',ylabel='Effective coefficient',title='X <-> Z <-> Y: all four edge rates = 1'); axes[0].legend(fontsize=8)
axes[1].plot(rows[:,0],rows[:,6],label='Current leaving X'); axes[1].plot(rows[:,0],rows[:,7],label='Current entering Y')
axes[1].plot(rows[:,0],rows[:,0]*rows[:,5],'--',label='Their difference: internal loss')
axes[1].set(xlabel='Loss at internal species Z',ylabel='Stationary current (X=2, Y=3)',title='Endpoint currents after internal-species\nelimination'); axes[1].legend(fontsize=8)
for ax in axes: ax.grid(alpha=.2)
fig.savefig(out/'path_reduction.png',dpi=180); fig.savefig(out/'path_reduction.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__)),
'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
Type II: exact positive stationary state; 64 principal minors checked. det(-H)=90368/1125; lower bound=36. Type V: 8 literal species reduce to three base variables with exact reconstruction. An internal-only loss gives two endpoint leaks of 1/3, even when both direct endpoint losses vanish. Fixed-rate continuation: 41/41 numerical roots. At most one positive root is the theorem; existence, stability at arbitrary loss and global continuation are not asserted.