Example code

"""Full multisite phosphorylation: identical static behavior, different dynamics.

Run python example.py. Inputs below are dimensionless paper witnesses, not fitted
biochemistry. Output JSON distinguishes freshly certified facts from numerics.
"""
# ------------------------- EDITABLE INPUTS -------------------------
SUBSTRATE=('2','12','1/5','2/5')
FREE_KINASE='23/10'
FREE_PHOSPHATASE='2/5'
KINASE_COMPLEXES=('23/50','3/10','17/10')
PHOSPHATASE_COMPLEXES=('23','9/50','4')
PRODUCTIVE_FLUX=('1/10','32/5','6/5')
KINASE_REVERSE=('1/100',)*3
PHOSPHATASE_REVERSE=('1','1/100','1/100')
REVERSE_PARAMETER_SWEEP=(1.,1.3,1.5,2.)
TRANSIENT_TIME=1000.
ORBIT_PARAMETER=1.
ORBIT_PERIOD_GUESS=39.59
ADDED_SITE_LOADS=('1/100',)
REVERSIBLE_DRIVE=20.
CONCENTRATION_UNIT_MICROMOLAR=.1
TIME_UNIT_SECONDS=100.
OUTPUT_DIRECTORY='outputs'
# Certificate reference cases A, H, C2 reduction and site addition are fixed
# paper reproductions; edits above change the separate exploration model.
# ------------------------------------------------------------------
import argparse,csv,hashlib,json,platform
from pathlib import Path
import numpy as np
import sympy as sp
from scipy.signal import find_peaks
import phos
from reactor import EquilibriumDesign,MassActionReactor,ConstantBinding

MANUSCRIPT_SHA256='7a0bbd174ad9a01c8765e169ed4dc1541cb7680dd3fd18d8f4529569f52c87d9'


def design():
    return EquilibriumDesign(SUBSTRATE,FREE_KINASE,FREE_PHOSPHATASE,KINASE_COMPLEXES,PHOSPHATASE_COMPLEXES,PRODUCTIVE_FLUX,KINASE_REVERSE,PHOSPHATASE_REVERSE)


def certificate_record(c):
    keys=('t','r','omega','crossing','l1','eps_derivative')
    return dict(evidence='Fresh exact polynomial isolation and outward-rounded rational arithmetic, using the attributed manuscript kernel; not Lean.',
        **{k:phos.bounds(c[k]) for k in keys if k in c},routh_first_column=[phos.bounds(v) for v in c['routh']],
        p0=[str(v) for v in c['p0']],p1=[str(v) for v in c['p1']],
        right_eigenvector=[phos.bounds(v) for v in c['right']],left_eigenvector=[phos.bounds(v) for v in c['left']],
        constant_coefficient_positive=c['const_coeff_positive'])


def certify_references():
    results={};raw={}
    for name,builder in [('A',lambda r:phos.witness('A',r)),('H',lambda r:phos.witness('H',r)),('A_kinase_clamped',lambda r:phos.witness('A',r,clampE=True)),('A_added_site',lambda r:phos.extended('A',r,list(ADDED_SITE_LOADS)))]:
        m=builder(sp.Rational(1,2));norm=2*m.n+2 if name=='A_added_site' else None
        c=phos.certify_hopf(builder,norm_index=norm);raw[name]=c;results[name]=certificate_record(c)
    reduced=phos.certify_reduced(lambda r:phos.witness('A',r),fast=4,norm_index=8)
    results['A_C2_eliminated']=certificate_record(reduced)
    results['A_C2_eliminated']['scope']='Local algebraic reduction with induced quadratic AND cubic terms; not a trajectory-error guarantee.'
    assert raw['A']['l1'].b<0 and raw['H']['l1'].a>0 and raw['A_added_site']['l1'].b<0
    assert raw['A_kinase_clamped']['l1'].a>0 and reduced['l1'].b<0
    c=raw['A'];m=phos.witness('A',sp.Rational(1,2));L,M,A=m.blocks();Jstatic=L*A.inv()*M
    cp=Jstatic.charpoly().all_coeffs()
    assert all(v>0 for v in cp) and cp[1]*cp[2]>cp[3]
    other=phos.witness('A',sp.Rational(7,3));L2,M2,A2=other.blocks();assert Jstatic==L2*A2.inv()*M2
    sensors={'S3':[0,0,1,0,0,0,0,0,-1],'S3+D3':[0,0,1,0,0,0,0,0,0]}
    gains,identity=phos.gain_constants(lambda r:phos.witness('A',r),c,sensors)
    from interval_arithmetic import I,C,dot
    a=-c['crossing'].re;scale=(a/(-c['omega']*c['l1'])).sqrt();amplitudes={}
    for name,v in sensors.items():
        z=dot([C(x) for x in v],c['right']);amplitudes[name]=phos.bounds(4*(z.re*z.re+z.im*z.im).sqrt()*scale)
    results['mechanism']=dict(static_characteristic=[str(v) for v in cp],static_independence='Exact at two distinct family values; row-factor cancellation proves the full static field independence.',
        resonance_constants={k:phos.bounds(v) for k,v in gains.items()},rank_one_identity=phos.bounds(identity*C(c['r'])),
        amplitude_prefactors=amplitudes,radial_recovery_prefactor=phos.bounds(2*a),
        all_n_scope='The positive-load site-addition theorem is imported; this run freshly certifies the configured finite extension. n<=2 means no Hopf, not no periodic orbit.',
        finite_amplitude_scope='Manuscript interval orbit proofs and validated pulse proof are not rerun here. Simulations, shooting, Floquet spectra and pulse trajectories below are numerical only.')
    return results,raw


def encoded(v):
    if isinstance(v,np.ndarray):return v.tolist()
    if isinstance(v,np.generic):return v.item()
    if isinstance(v,complex):return dict(real=v.real,imag=v.imag)
    raise TypeError(type(v).__name__)


def main():
    parser=argparse.ArgumentParser(description=__doc__);parser.add_argument('--output',default=OUTPUT_DIRECTORY)
    out=Path(parser.parse_args().output);out.mkdir(exist_ok=True)
    def dump(name,data):(out/name).write_text(json.dumps(data,indent=2,default=encoded)+'\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)
    reference,raw=certify_references();dump('hopf_certificates.json',reference)
    print('Fresh rational Hopf certificates: A attracting, H subcritical, kinase clamp subcritical, site addition attracting.',flush=True)
    d=design();base=MassActionReactor(d.exact());n=base.n
    if n!=3:raise ValueError('The teaching plots use a three-site starting design; the reactor classes accept arbitrary site count.')
    rows=[];trajectories={};static_fields=[]
    y0=np.zeros(base.dim);y0[2]=.01
    for r in REVERSE_PARAMETER_SWEEP:
        model=MassActionReactor(d.with_reverse(r).exact());J=model.jacobian(0,np.zeros(model.dim));eig=np.linalg.eigvals(J)
        sol,X=model.integrate(y0,TRANSIENT_TIME)
        z,uf=model.static_complexes(np.array([.001,-.001,.001]));static_fields.append(uf)
        rows.append(dict(r=r,equilibrium_residual=float(max(abs(model.field(0,np.zeros(model.dim))))),leading_real=float(max(eig.real)),
            static_leading_real=float(max(np.linalg.eigvals(model.static_jacobian()).real)),static_local_field=uf,
            minimum_concentration=float(X.min()),conservation_drift=float(np.max(abs(model.conservation@X-([email protected])[:,None])))))
        trajectories[str(r)]=(sol,X)
        table(f'trajectory_r_{r}.csv',['time',*model.labels],zip(sol.t,*X))
    assert np.max(np.ptp(np.array(static_fields),axis=0))<1e-10
    dump('static_and_dynamic.json',dict(evidence='Exact reference static comparison above; editable sweep and nonlinear implicit solves are numerical.',sweep=rows))
    # Locate a repeated maximum after a long transient, then solve the return map.
    orbitmodel=MassActionReactor(d.with_reverse(ORBIT_PARAMETER).exact())
    settled,X=orbitmodel.integrate(y0,TRANSIENT_TIME,6001)
    peaks,_=find_peaks(X[3]);peaks=peaks[settled.t[peaks]>.6*TRANSIENT_TIME]
    if len(peaks)<2:raise RuntimeError('No repeated oscillation detected. Adjust the parameter, transient or disable this experiment in your own driver.')
    seed=settled.y[:,peaks[-1]];periodguess=settled.t[peaks[-1]]-settled.t[peaks[-2]]
    anchor,T,multipliers,residual=orbitmodel.periodic_orbit(seed,periodguess if periodguess>0 else ORBIT_PERIOD_GUESS)
    cyc,CX=orbitmodel.integrate(anchor,T,1601,ledger=True);led=cyc.y[orbitmodel.dim:,-1]
    trivial=int(np.argmin(abs(multipliers-1)));nontrivial=np.delete(multipliers,trivial)
    orbit=dict(evidence='Numerical shooting and variational integration, NOT a validated orbit or minimal-period proof.',r=ORBIT_PARAMETER,period=T,
        shooting_residual=residual,multipliers=multipliers,largest_nontrivial_modulus=float(max(abs(nontrivial))),
        minimum_concentration=float(CX.min()),S3_peak_to_peak=float(np.ptp(CX[3])),
        kinase_turnovers=led[:n],phosphatase_turnovers=led[n:],paired_flux_difference=led[:n]-led[n:],
        total_forward_ATP_turnover=float(sum(led[:n])),turnover_per_substrate=float(sum(led[:n])/([email protected])[2]))
    table('periodic_orbit.csv',['time',*orbitmodel.labels,*[f'integral_{s}' for s in ['C1','C2','C3','D1','D2','D3']]],zip(cyc.t,*CX,*cyc.y[orbitmodel.dim:]))
    dump('periodic_orbit.json',orbit)
    # Uniform slowing preserves equilibrium and static behavior; clamping changes
    # the physical conservation contract by replacing enzyme pools with buffers.
    interventions=[]
    for label,model in [('both_free_enzymes_buffered',MassActionReactor(d.with_reverse('1').exact(True,True))),
        ('complexes_faster',MassActionReactor(d.with_reverse('1').exact(),relaxation=.5)),
        ('driven_reversible',MassActionReactor(d.with_reverse('1').exact(),drive=REVERSIBLE_DRIVE))]:
        sol,X=model.integrate(y0,TRANSIENT_TIME)
        interventions.append(dict(name=label,terminal_peak_to_peak=float(np.ptp(X[3,-301:])),minimum_concentration=float(X.min()),
            equilibrium_residual_at_original_design=float(max(abs(model.field(0,np.zeros(model.dim)))))))
        trajectories[label]=(sol,X)
    # Detailed-balance construction: potentials vary with M; not fixed chemistry
    # under changing ATP concentration. Reverse currents keep the three totals.
    rm=MassActionReactor(d.with_reverse('1').exact(),drive=REVERSIBLE_DRIVE);mm=rm.model
    mu=np.zeros(len(rm.xstar));mu[:n+1]=np.arange(n+1)*REVERSIBLE_DRIVE
    mu[rm.cpx]=mu[rm.lev]-np.log(rm.kon/rm.koff)
    binding_error=np.log(rm.kon/rm.koff)-(mu[rm.lev]+mu[rm.enz]-mu[rm.cpx])
    cat_potential=mu[rm.cpx]-mu[rm.out]-mu[rm.enz];cat_potential[:n]+=2*REVERSIBLE_DRIVE
    catalysis_error=np.log(rm.kcat/rm.reverse)-cat_potential
    revsol,revX=rm.integrate(anchor,TRANSIENT_TIME,1201,ledger=True)
    # No period proof here: accumulated net fuel is reported over a finite window.
    net=float(sum(revsol.y[rm.dim:rm.dim+n,-1]))
    dump('interventions.json',dict(experiments=interventions,thermodynamics=dict(driving_M=REVERSIBLE_DRIVE,
        local_detailed_balance_residual=float(max(np.max(abs(binding_error)),np.max(abs(catalysis_error)))),
        finite_window_net_ATP=net,cycle_affinity=2*REVERSIBLE_DRIVE,
        scope='Finite-time net consumption is not an entropy-production identity away from a closed periodic orbit. The theorem gives no numerical M0. Potentials are redesigned together with M; ATP/ADP/Pi remain chemostatted.')))
    c=raw['A'];rc=phos.mid(c['r']);omega=phos.mid(c['omega']);sensor=np.array([0,0,1,0,0,0,0,0,-1]);resonance=[]
    for delta in [.01,.03,.1,.3,1.]:
        m=MassActionReactor(phos.witness('A',str(rc+delta)));J=m.jacobian(0,np.zeros(9));b=np.zeros(9);b[6]=(1+rc+delta)*.1
        gain=abs([email protected](1j*omega*np.eye(9)-J,b))
        bounds=reference['mechanism']['resonance_constants']['S3'];K=float((sp.Rational(bounds[0])+sp.Rational(bounds[1]))/2)
        exactlaw=(1+rc+delta)*K/delta
        resonance.append(dict(delta=delta,direct_gain=gain,rank_one_law=exactlaw,relative_residual=abs(gain/exactlaw-1)))
    pulses=[]
    for r in ['13/10','3/2']:
        pm=MassActionReactor(phos.witness('A',r),ConstantBinding(1.5));ps,PX=pm.integrate(np.zeros(9),.05,51)
        pulses.append(dict(r=r,D1_terminal=float(PX[9,-1]),scope='Nominal numerical pulse only; no preparation/actuator uncertainty enclosure.'))
    dump('dynamic_readouts.json',dict(resonance=resonance,nominal_pulses=pulses,scope='The resonance law is exact for the linearized response at the certified frequency. Finite forcing, measurement error and experiment feasibility require separate analysis.'))
    dump('units.json',dict(calibrated=False,concentration_unit_micromolar=CONCENTRATION_UNIT_MICROMOLAR,time_unit_seconds=TIME_UNIT_SECONDS,
        totals_micromolar=([email protected])*CONCENTRATION_UNIT_MICROMOLAR,period_minutes=T*TIME_UNIT_SECONDS/60,
        free_phosphatase_molecules_per_femtoliter=float(base.xstar[base.model.iF]*CONCENTRATION_UNIT_MICROMOLAR*602.214076),
        caveat='Illustrative scale only; heavy enzyme loading, broad catalytic rates and low molecule numbers can invalidate a biological or deterministic interpretation.'))
    plot(out,trajectories,cyc,CX,rows,resonance)
    print(f'Numerical r={ORBIT_PARAMETER:g} orbit: period {T:.9f}; S3 swing {orbit["S3_peak_to_peak"]:.6f}; largest nontrivial multiplier {max(abs(nontrivial)):.6f}.',flush=True)
    print('Static fields agree while the full stability changes. Reverse catalysis, buffering, fuel and dynamic readouts evaluated.',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(),
        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 plot(out,trajectories,cyc,X,rows,resonance):
    import matplotlib
    matplotlib.use('Agg')
    import matplotlib.pyplot as plt
    fig,axs=plt.subplots(1,2,figsize=(10,4),layout='constrained')
    for key in [str(REVERSE_PARAMETER_SWEEP[1]),str(REVERSE_PARAMETER_SWEEP[2])]:
        sol,Z=trajectories[key];axs[0].plot(sol.t,Z[3],label='r = '+key)
    axs[0].set(xlabel='Time (model units)',ylabel='Free fully phosphorylated S3',title='Phosphorylation trajectories at matched\nequilibria');axs[0].legend()
    axs[1].plot([r['r'] for r in rows],[r['leading_real'] for r in rows],'o-',label='Full mass action')
    axs[1].plot([r['r'] for r in rows],[r['static_leading_real'] for r in rows],'s--',label='Static elimination')
    axs[1].axhline(0,color='black',lw=.7);axs[1].set(xlabel='D1 dissociation ratio r',ylabel='Largest real eigenvalue',title='Stability of full and complex-eliminated\nmodels');axs[1].legend(fontsize=8)
    for ax in axs:ax.grid(alpha=.2)
    fig.savefig(out/'mechanism.png',dpi=180);fig.savefig(out/'mechanism.svg');plt.close(fig)
    fig,axs=plt.subplots(1,2,figsize=(10,4),layout='constrained')
    for i in [3,6,8,9,11]:axs[0].semilogy(cyc.t/cyc.t[-1],X[i],label=['S0','S1','S2','S3','E','F','C1','C2','C3','D1','D2','D3'][i])
    axs[0].set(xlabel='Fraction of numerical period',ylabel='Concentration',title='An attracting numerical cycle');axs[0].legend(fontsize=8,ncol=2)
    axs[1].loglog([r['delta'] for r in resonance],[r['direct_gain'] for r in resonance],'o',label='Matrix response')
    axs[1].loglog([r['delta'] for r in resonance],[r['rank_one_law'] for r in resonance],label='Exact resonance law')
    axs[1].set(xlabel='Distance above Hopf threshold',ylabel='Linear S3 gain at Hopf frequency',title='Linear substrate response near the Hopf\nthreshold');axs[1].legend(fontsize=8)
    for ax in axs:ax.grid(alpha=.2)
    fig.savefig(out/'cycle_and_readout.png',dpi=180);fig.savefig(out/'cycle_and_readout.svg');plt.close(fig)


if __name__=='__main__':main()