"""Low-order mass-action instability invisible to every child selection. Exact structural and spectral checks are separate from numerical trajectories. """ from __future__ import annotations # EDITABLE INPUTS: dimensionless manuscript operating conditions -------------- CIRCULATION_RATIO = '100' # T>1 INVERSE_D_CONCENTRATION = '100' # L>0; equilibrium (1,1,1,1/L) THROUGHPUT = '1' # s>0; rescales all original rates and time DILUTION_RATE = '1' # completion keeps the selected equilibrium PERTURBATION_SIZE = 1e-7 TRAJECTORY_END = 0.8 TRAJECTORY_POINTS = 401 RATE_PERTURBATION_FRACTION = 0.001 PARAMETER_T_RANGE = (1.05, 120.0) PARAMETER_L_RANGE = (0.1, 200.0) MANUSCRIPT_SHA256 = '12479870603ad4aafcc2a7cef8d0c32f22a78f07e210ba83c4f904c249063dbf' # --------------------------------------------------------------------------- import argparse from dataclasses import dataclass import hashlib import itertools import json from pathlib import Path import platform import csv import numpy as np import sympy as sp from scipy.integrate import solve_ivp def Q(x): return sp.Rational(str(x)) def vec(xs): return sp.Matrix([sp.sympify(x) if isinstance(x, sp.Basic) else Q(x) for x in xs]) def encoded(xs): return [str(sp.simplify(x)) for x in xs] @dataclass(frozen=True) class Reaction: name: str reactant: tuple[int, ...] product: tuple[int, ...] class Source: """Literal integer complexes; deterministic homodimer event rate is k*x**2.""" def __init__(self, species, reactions): self.species = tuple(species); self.reactions = tuple(reactions) self.n = len(self.species); self.m = len(self.reactions) if not self.n or not self.m or len(set(self.species)) != self.n: raise ValueError('Distinct species and at least one reaction required.') for r in self.reactions: if len(r.reactant) != self.n or len(r.product) != self.n or any( not isinstance(v, int) or v < 0 for v in (*r.reactant, *r.product)): raise ValueError('Complexes must be matching nonnegative integer vectors.') self.Y = sp.ImmutableMatrix.hstack(*[sp.ImmutableMatrix(r.reactant) for r in self.reactions]) self.P = sp.ImmutableMatrix.hstack(*[sp.ImmutableMatrix(r.product) for r in self.reactions]) self.S = self.P-self.Y def monomials(self, x): return sp.Matrix([sp.prod(x[i]**self.Y[i, j] for i in range(self.n)) for j in range(self.m)]) def reconstruct(self, equilibrium, flux): x = vec(equilibrium); v = vec(flux) if len(x) != self.n or len(v) != self.m or any(a <= 0 for a in (*x, *v)): raise ValueError('Positive matching equilibrium and flux vectors required.') if (self.S*v).applyfunc(sp.simplify) != sp.zeros(self.n, 1): raise ValueError('Flux must be stationary, S*v=0.') mono = self.monomials(x) return MassAction(self, [sp.simplify(v[j]/mono[j]) for j in range(self.m)], x) def molecularity(self): return {'reactants': [int(sum(self.Y[:, j])) for j in range(self.m)], 'products': [int(sum(self.P[:, j])) for j in range(self.m)], 'catalyst_free': all(self.Y[i, j]*self.P[i, j] == 0 for i in range(self.n) for j in range(self.m))} class MassAction: def __init__(self, source, rates, equilibrium=None): self.source = source; self.k = vec(rates) if len(self.k) != source.m or any(k <= 0 for k in self.k): raise ValueError('One positive rate constant per reaction required.') self.equilibrium = None if equilibrium is None else vec(equilibrium) self.x = sp.Matrix(sp.symbols(f'x0:{source.n}')) self.symbolic_flux = sp.matrix_multiply_elementwise(self.k, source.monomials(self.x)) self.symbolic_field = source.S*self.symbolic_flux self.symbolic_jacobian = self.symbolic_field.jacobian(self.x) self._field = sp.lambdify([tuple(self.x)], self.symbolic_field, 'numpy') self._jac = sp.lambdify([tuple(self.x)], self.symbolic_jacobian, 'numpy') if self.equilibrium is not None: if len(self.equilibrium) != source.n or any(v <= 0 for v in self.equilibrium): raise ValueError('Equilibrium must be positive and match species.') if self.exact_field(self.equilibrium).applyfunc(sp.simplify) != sp.zeros(source.n, 1): raise ValueError('Supplied state is not an exact equilibrium.') def exact_field(self, x): return self.symbolic_field.subs(dict(zip(self.x, x))) def jacobian(self, x=None): state = self.equilibrium if x is None else vec(x) if state is None: raise ValueError('Supply a state for the Jacobian.') return self.symbolic_jacobian.subs(dict(zip(self.x, state))).applyfunc(sp.simplify) def field(self, x): return np.asarray(self._field(np.asarray(x, float)), float).reshape(-1) def integrate(self, initial, times, method='Radau'): initial = np.asarray(initial, float); times = np.asarray(times, float) if initial.shape != (self.source.n,) or np.any(initial <= 0) or not np.all(np.isfinite(initial)): raise ValueError('Positive finite initial concentrations required.') if len(times) < 2 or times[0] != 0 or not np.all(np.diff(times) > 0): raise ValueError('Strictly increasing times starting at zero required.') solution = solve_ivp(lambda t, x: self.field(x), (0, times[-1]), initial, t_eval=times, jac=lambda t, x: np.asarray(self._jac(x), float), method=method, rtol=1e-11, atol=1e-14) if not solution.success or np.min(solution.y) <= 0: raise ArithmeticError('Numerical integration failed or left the positive region.') return solution.y.T class OperatingFamily: def __init__(self): self.source = Source(('A', 'B', 'C', 'D'), ( Reaction('A+D -> 2C', (1,0,0,1), (0,0,2,0)), Reaction('B+C -> empty', (0,1,1,0), (0,0,0,0)), Reaction('2C -> D', (0,0,2,0), (0,0,0,1)), Reaction('2D -> A+4B+4C', (0,0,0,2), (1,4,4,0)), Reaction('empty -> A', (0,0,0,0), (1,0,0,0)), Reaction('empty -> D', (0,0,0,0), (0,0,0,1)))) def flux(self, T, s=1): T, s = Q(T), Q(s) if T <= 1 or s <= 0: raise ValueError('Positive stationary cone requires T>1 and s>0.') return s*sp.Matrix([T, 4, T, 1, T-1, 2]) def model(self, T=100, L=100, s=1): L = Q(L) if L <= 0: raise ValueError('L must be positive.') return self.source.reconstruct([1, 1, 1, 1/L], self.flux(T, s)) def from_rates(self, rates): """Closed-form positive equilibrium for any six positive original rates. This follows from the complete stationary flux cone, rather than a nonlinear root solver. Stability still needs a separate calculation. """ k = vec(rates) if len(k) != 6 or any(v <= 0 for v in k): raise ValueError('Six positive rates required.') s = k[5]/2; T = 1+k[4]/s D = sp.sqrt(s/k[3]); C = sp.sqrt(s*T/k[2]) state = sp.Matrix([s*T/(k[0]*D), 4*s/(k[1]*C), C, D]) return MassAction(self.source, k, state) @staticmethod def coefficients(T, L): return ((T+4)*L+5*T+8, (14*T+32)*L+4*T*(T+6), 112*T*L+16*T*T, 64*T*T*L) def exact_family_certificate(self): T, L, u = sp.symbols('T L u', positive=True) a = self.coefficients(T, L); H = sp.expand(Quartic.hurwitz(a)) # Derive the Jacobian from flux and the mass-action derivative identity. v = sp.Matrix([T, 4, T, 1, T-1, 2]) J = self.source.S*sp.diag(*v)*self.source.Y.T*sp.diag(1,1,1,L) if any(sp.expand(x-y) != 0 for x, y in zip(J.charpoly().all_coeffs()[1:], a)): raise AssertionError('Family characteristic polynomial mismatch.') ray = sp.Poly(H.subs({T:100, L:100+u}), u) if not all(v < 0 for v in ray.all_coeffs()): raise AssertionError('Ray sign certificate failed.') reduced = (J[:3,:3]-J[:3,3]*(1/J[3,3])*J[3,:3]).applyfunc(sp.factor) if any(L in x.free_symbols for x in reduced): raise AssertionError('Fast reduction retained L.') cubic = reduced.charpoly().all_coeffs()[1:] threshold = (41+sp.sqrt(2577))/4 return {'family_coefficients': encoded(a), 'H': str(H), 'negative_ray_coefficients_descending': encoded(ray.all_coeffs()), 'large_L_leading_coefficient': str(sp.factor(sp.Poly(H,L).LC())), 'fast_D_matrix': [[str(x) for x in reduced.row(i)] for i in range(3)], 'fast_D_cubic_coefficients': encoded(cubic), 'fast_D_hurwitz': str(sp.factor(cubic[0]*cubic[1]-cubic[2])), 'large_L_threshold': str(threshold), 'large_L_threshold_numeric': float(threshold)} class Quartic: @staticmethod def hurwitz(coefficients): a,b,c,d = coefficients return a*b*c-c*c-a*a*d def __init__(self, matrix): if matrix.shape != (4,4): raise ValueError('Quartic check needs a four-species matrix.') self.z = sp.Symbol('z') self.p = matrix.charpoly(self.z).as_expr() self.a = sp.Poly(self.p, self.z).all_coeffs()[1:] def classification(self): if not all(v > 0 for v in self.a): return {'status': 'unknown', 'reason': 'Positive-coefficient quartic criterion not applicable.'} H = sp.simplify(self.hurwitz(self.a)) status = 'strictly_stable' if H > 0 else ('unstable_complex_pair' if H < 0 else 'imaginary_pair_boundary') return {'status': status, 'coefficients': encoded(self.a), 'H': str(H), 'second_hurwitz': str(sp.simplify(self.a[0]*self.a[1]-self.a[2]))} def shifted_root_certificate(self, lo=2, hi=4): lo, hi = Q(lo), Q(hi) if not 0 < lo < hi: raise ValueError('Positive ordered shift bracket required.') t = sp.Symbol('t') a,b,c,d = sp.Poly(sp.expand(self.p.subs(self.z, self.z+t)), self.z).all_coeffs()[1:] H = sp.expand(self.hurwitz((a,b,c,d))) # Positive polynomial coefficients imply A(t), C(t)>0 throughout t>=lo. positive = all(all(v >= 0 for v in sp.Poly(expr.subs(t,t+lo), t).all_coeffs()) and expr.subs(t,lo)>0 for expr in (a,c)) endpoints = [H.subs(t,lo), H.subs(t,hi)] certified = positive and endpoints[0]*endpoints[1] < 0 return {'certified': bool(certified), 'real_part_bracket': encoded([lo,hi]), 'H_at_endpoints': encoded(endpoints), 'shifted_coefficients': encoded([a,b,c,d]), 'H_shift': str(H), 'root_construction': 'At an intervening zero of H(t), omega=sqrt(C(t)/A(t)); p(t+i*omega)=0.'} class ChildLattice: """Enumerate injective consuming assignments, including the empty selection.""" def __init__(self, source): self.source = source def patterns(self, budget=100000): choices = [[None]+[j for j in range(self.source.m) if self.source.Y[i,j]>0] for i in range(self.source.n)] if np.prod([len(c) for c in choices], dtype=object) > budget: raise RuntimeError('Enumeration budget exceeded; no absence-of-core conclusion.') return [p for p in itertools.product(*choices) if len([x for x in p if x is not None]) == len({x for x in p if x is not None})] def matrix(self, pattern): selected = [i for i,r in enumerate(pattern) if r is not None] if not selected: return sp.zeros(0,0) return self.source.S.extract(selected, [pattern[i] for i in selected]) def padded(self, pattern): if self.source.n != 4: raise ValueError('This quartic certificate route requires four species.') M = -sp.eye(4) for j,r in enumerate(pattern): if r is not None: for i,q in enumerate(pattern): M[i,j] = self.source.S[i,r] if q is not None else 0 return M def certificates(self): scales = sp.symbols('a b c d', positive=True) records = [] def nonnegative(expr, strict=False): coefficients = sp.Poly(sp.expand(expr), *scales).coeffs() return all(v >= 0 for v in coefficients) and (not strict or any(v > 0 for v in coefficients)) for pattern in self.patterns(): padded = self.padded(pattern) a,b,c,d = (padded*sp.diag(*scales)).charpoly().all_coeffs()[1:] H = sp.expand(Quartic.hurwitz((a,b,c,d))) passed = nonnegative(a,True) and nonnegative(b) and nonnegative(c,True) and nonnegative(d) and nonnegative(H) records.append({'assignment': list(pattern), 'species': [i for i,r in enumerate(pattern) if r is not None], 'padded_coefficients': encoded([a,b,c,d]), 'hurwitz_polynomial': str(H), 'hurwitz_monomials': [{'exponents': list(m), 'coefficient': str(v)} for m,v in sp.Poly(H,*scales).terms()], 'certified_D_nonunstable': bool(passed), 'singular': d == 0, 'scope': 'All positive diagonal column scalings, by polynomial nonnegativity and the quartic exclusion lemma.' if passed else 'Criterion inconclusive; no instability conclusion.'}) return records def maximal_patterns(self): patterns = self.patterns() return [p for p in patterns if not any(p!=q and all(r is None or r==q[i] for i,r in enumerate(p)) for q in patterns)] class DilutionCompletion: """Four active reactions, four adjusted feeds, and four first-order losses.""" def __init__(self, original, delta): self.original = original; self.delta = Q(delta) if self.delta <= 0 or original.equilibrium is None: raise ValueError('Positive dilution and an exactly stationary original model required.') s = original.source if s.n != 4 or s.m != 6 or s.Y[:,4:] != sp.zeros(4,2): raise ValueError('Completion expects this six-reaction source with its two terminal inflows.') empty = (0,0,0,0); basis = [tuple(int(i==j) for i in range(4)) for j in range(4)] reactions = list(s.reactions[:4])+[Reaction('feed '+s.species[i], empty, basis[i]) for i in range(4)] reactions += [Reaction('loss '+s.species[i], basis[i], empty) for i in range(4)] feeds = s.S[:,4:]*original.k[4:,:]+self.delta*original.equilibrium self.feed_concentrations = feeds/self.delta self.model = MassAction(Source(s.species,reactions), list(original.k[:4,:])+list(feeds)+[self.delta]*4, original.equilibrium) def transport_certificates(self, original_records): old = {tuple(r['assignment']): r['certified_D_nonunstable'] for r in original_records} lattice = ChildLattice(self.model.source); rows = [] for p in lattice.patterns(): losses = [i for i,r in enumerate(p) if r is not None and r>=8] oldpattern = tuple(None if i in losses else r for i,r in enumerate(p)) if any(p[i] != 8+i for i in losses): raise AssertionError('Non-diagonal loss assignment.') retained = [i for i,r in enumerate(oldpattern) if r is not None] # Verify the exact triangular block: loss columns vanish on every other species. if any(self.model.source.S[j,p[i]] != -int(i==j) for i in losses for j in range(4)): raise AssertionError('Loss transport block failed.') if retained and self.model.source.S.extract(retained,[p[i] for i in retained]) != self.original.source.S.extract(retained,[p[i] for i in retained]): raise AssertionError('Retained old-child block changed.') rows.append({'assignment':list(p),'loss_species':losses,'old_assignment':list(oldpattern), 'certified_D_nonunstable':old.get(oldpattern,False)}) return rows def local_departure(model, size=PERTURBATION_SIZE, end=TRAJECTORY_END): if model.equilibrium is None or size <= 0 or end <= 0: raise ValueError('Positive perturbation, duration and reference equilibrium required.') J = np.array(model.jacobian(),float) eigenvalues, vectors = np.linalg.eig(J) index = max(range(len(eigenvalues)),key=lambda i:(eigenvalues[i].real,eigenvalues[i].imag)) lam = eigenvalues[index]; w = vectors[:,index] pivot = int(np.argmax(np.abs(w))); w = w/w[pivot] # deterministic phase and maximum modulus 1 equilibrium = np.array(model.equilibrium,float).reshape(-1) times = np.linspace(0,end,TRAJECTORY_POINTS) nonlinear = model.integrate(equilibrium+size*w.real,times) linear = size*np.real(np.exp(lam*times[:,None])*w[None,:]) deviations = nonlinear-equilibrium independent = model.integrate(equilibrium+size*w.real,times,method='BDF') scale = float(np.max(np.linalg.norm(linear,axis=1))) return times,deviations,linear,{'eigenvalue_real':float(lam.real),'eigenvalue_imag':float(lam.imag), 'relative_linear_error':float(np.max(np.linalg.norm(deviations-linear,axis=1))/scale), 'independent_solver_error_relative_to_mode':float(np.max(np.linalg.norm(nonlinear-independent,axis=1))/scale), 'scope':'Short-time numerical departure only; no periodic-orbit or long-time claim.'} def write_csv(out,name,headers,rows): with (out/name).open('w',newline='') as stream: writer=csv.writer(stream);writer.writerow(headers);writer.writerows(rows) 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) family=OperatingFamily();model=family.model(CIRCULATION_RATIO,INVERSE_D_CONCENTRATION,THROUGHPUT) children=ChildLattice(family.source);records=children.certificates() if not all(r['certified_D_nonunstable'] for r in records): raise AssertionError('Original child exclusion not established.') completion=DilutionCompletion(model,DILUTION_RATE);transport=completion.transport_certificates(records) spectral=Quartic(model.jacobian());source_witness=family.model() source_root=Quartic(source_witness.jacobian()).shifted_root_certificate() times,deviations,linear,dynamics=local_departure(model) cert=family.exact_family_certificate() result={'inputs':{'T':CIRCULATION_RATIO,'L':INVERSE_D_CONCENTRATION,'s':THROUGHPUT,'dilution':DILUTION_RATE}, 'molecularity':family.source.molecularity(),'rates':encoded(model.k),'equilibrium':encoded(model.equilibrium), 'flux':encoded(family.flux(CIRCULATION_RATIO,THROUGHPUT)), 'jacobian':[[str(v) for v in model.jacobian().row(i)] for i in range(4)], 'operating_classification':spectral.classification(), 'source_100_100_root_certificate':source_root,'family_certificate':cert, 'stable_control':Quartic(family.model(2,1).jacobian()).classification(), 'child_count':len(records),'singular_children':sum(r['singular'] for r in records), 'maximal_assignments':[list(p) for p in children.maximal_patterns()], 'completion':{'reaction_count':12,'feed_concentrations':encoded(completion.feed_concentrations), 'child_count':len(transport),'all_children_transported':all(r['certified_D_nonunstable'] for r in transport), 'spectral_classification':Quartic(completion.model.jacobian()).classification(), 'source_margin_applies':bool(Q(CIRCULATION_RATIO)==100 and Q(INVERSE_D_CONCENTRATION)==100 and Q(DILUTION_RATE)<2*Q(THROUGHPUT))}, 'local_dynamics':dynamics, 'mass_obstruction':'2C->D and 2D->A+4B+4C imply mass_A+4*mass_B=0, impossible with positive masses.', 'scope':'Kinetic abstract-species counterexample; products can have nine molecules. No atom-balanced mechanism, Hopf theorem or Lean rerun.'} (out/'results.json').write_text(json.dumps(result,indent=2)+'\n') (out/'child_certificates.json').write_text(json.dumps(records,indent=2)+'\n') (out/'completion_children.json').write_text(json.dumps(transport,indent=2)+'\n') write_csv(out,'local_trajectory.csv',['time',*[f'nonlinear_d{x}' for x in family.source.species],*[f'linear_d{x}' for x in family.source.species]], [[t,*d,*lin] for t,d,lin in zip(times,deviations,linear)]) perturbations=[] for reaction in range(6): for direction in (-1,1): rates=model.k.copy();rates[reaction]*=1+direction*Q(RATE_PERTURBATION_FRACTION) perturbed=family.from_rates(rates) eig=np.linalg.eigvals(np.array(perturbed.jacobian(),float)) perturbations.append([reaction,direction*RATE_PERTURBATION_FRACTION,*[float(x) for x in perturbed.equilibrium],max(eig.real)]) write_csv(out,'rate_perturbations.csv',['reaction','relative_change','A_star','B_star','C_star','D_star','numeric_spectral_abscissa'],perturbations) Ts=np.linspace(*PARAMETER_T_RANGE,141);Ls=np.geomspace(*PARAMETER_L_RANGE,111) TT,LL=np.meshgrid(Ts,Ls);HH=Quartic.hurwitz(family.coefficients(TT,LL)) write_csv(out,'operating_grid.csv',['T','L','numeric_H'],zip(TT.flat,LL.flat,HH.flat)) lines=[f'Operating rates: {encoded(model.k)}; equilibrium: {encoded(model.equilibrium)}.', f'Operating quartic: {spectral.classification()}.', f'{len(records)} original children certified for all positive column scalings; {result["singular_children"]} singular.', f'Source root bracket [2,4]: H endpoints {source_root["H_at_endpoints"]}.', f'Stable control (T,L)=(2,1): H={result["stable_control"]["H"]}.', f'Dilution completion: {len(transport)} transported children; feed={encoded(completion.feed_concentrations)}.', f'Numerical local mode: {dynamics["eigenvalue_real"]:.9g} + {dynamics["eigenvalue_imag"]:.9g} i.', f'Local nonlinear/linear relative error: {dynamics["relative_linear_error"]:.6g}.', 'Unstable linearization is not a sustained oscillation or atom-resolved mechanism.'] (out/'console.txt').write_text('\n'.join(lines)+'\n');print('\n'.join(lines)) import matplotlib matplotlib.use('Agg') import matplotlib.pyplot as plt fig,ax=plt.subplots(figsize=(8,4.8),layout='constrained') ax.contourf(TT,LL,(HH<0).astype(float),levels=[-.5,.5,1.5],colors=['#edf1f3','#bdd8e7']) ax.contour(TT,LL,HH,levels=[0],colors='#216b91',linewidths=1) ax.axvline(cert['large_L_threshold_numeric'],color='#777',ls='--',label='Large-L threshold (not finite-L boundary)') ax.plot([100,100],[100,200],color='#b95024',lw=3,label='Certified ray: T=100, L>=100') ax.scatter([2],[1],color='black',s=30,label='Certified stable control') ax.scatter([float(Q(CIRCULATION_RATIO))],[float(Q(INVERSE_D_CONCENTRATION))],marker='x',color='#b95024',label='Selected operating point') ax.text(62,35,'H < 0: unstable\n(numerical grid)',ha='center',color='#184a66') ax.text(62,.65,'H > 0: strictly stable\n(numerical grid)',ha='center') ax.set(yscale='log',xlabel='Circulation ratio T',ylabel='Inverse D concentration L',title='Equilibrium stability versus circulation\nand concentration',xlim=PARAMETER_T_RANGE,ylim=PARAMETER_L_RANGE) ax.legend(fontsize=8,loc='lower right') fig.savefig(out/'operating_region.png',dpi=180);fig.savefig(out/'operating_region.svg');plt.close(fig) fig,axes=plt.subplots(1,2,figsize=(10,4),layout='constrained') axes[0].plot(times,deviations[:,0]/PERTURBATION_SIZE,label='Nonlinear mass action',color='#216b91') axes[0].plot(times,linear[:,0]/PERTURBATION_SIZE,'--',label='Linear mode',color='#b95024') axes[0].set(xlabel='Dimensionless time',ylabel='(A - A*) / perturbation size',title='Short-time local departure');axes[0].legend(fontsize=8) axes[1].plot(deviations[:,0]/PERTURBATION_SIZE,deviations[:,2]/PERTURBATION_SIZE,color='#216b91') axes[1].scatter([deviations[0,0]/PERTURBATION_SIZE],[deviations[0,2]/PERTURBATION_SIZE],color='#b95024',label='Initial state') axes[1].set(xlabel='(A - A*) / perturbation size',ylabel='(C - C*) / perturbation size',title='Short-time concentration deviations in the\nA–C plane');axes[1].legend(fontsize=8) for ax in axes:ax.grid(alpha=.2) fig.savefig(out/'local_dynamics.png',dpi=180);fig.savefig(out/'local_dynamics.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()