# Attracting oscillations in distributive multisite phosphorylation

Companion to manuscript 64, SHA-256 `7a0bbd174ad9a01c8765e169ed4dc1541cb7680dd3fd18d8f4529569f52c87d9`. Python 3.11+.

```sh
python -m venv .venv
# Activate the environment using your platform's command.
python -m pip install -r requirements.txt
python -m unittest -v
python example.py
```

## Objects and construction

S0,...,Sn are free substrate molecules with progressively more phosphorylated sites. E and F are free kinase and phosphatase. Ci binds S(i-1) to E; Di binds Si to F. Every binding has a dissociation channel and a catalytic exit. The literal source contains 3n+3 species and 6n elementary reactions. Substrate and each enzyme have conserved totals; complexes contribute to both the appropriate enzyme total and the substrate total.

`EquilibriumDesign` takes positive free and bound concentrations, one productive flux qi per site and dimensionless reverse ratios rho_i,sigma_i. It derives binding, dissociation and catalytic constants. For the kinase arm these are `(1+rho)*q/(S*E)`, `rho*q/C`, `q/C`; the phosphatase arm is analogous. This is a constructive equilibrium parameterization, not fitting measured rates. `with_reverse(r)` changes the D1 ratio together with its binding rate. Changing a single physical rate independently is a different experiment.

The chart has n cumulative substrate coordinates u followed by n C and n D coordinates. ui includes free substrate at levels >=i and complexes whose **input** substrate level is >=i. The reactor integrates deviations y from the designed equilibrium and lifts by `x=xstar+P@y`. This preserves the compatibility class without projecting a drifting trajectory. The independent full-species field and stoichiometric conservation matrix are exposed for checking or coupling to a different solver.

The top of `example.py` contains editable concentrations, fluxes, ratios, sweep values, simulation times, extension loads, driving parameter and illustrative units. Fixed reference certificates remain explicitly separate from the editable exploration. The supplied plots assume a three-site starting design; the model classes support arbitrary n.

## Reusable components

* `reactor.py`: immutable design data and site addition; full mass-action reactor; new rate arrays; chart lift; analytic Jacobian; numerical static-complex solve; binding protocol interface; buffered enzymes; uniform complex relaxation; driven reverse catalysis; integration with per-arm net catalytic ledgers; numerical Newton shooting and Floquet multipliers.
* `phos.py`, `interval_arithmetic.py`: **unmodified, attributed manuscript companion kernels**, with original byte hashes in `provenance.json`. The exact reaction builder, rational interval Hopf certifier, local C2 reduction and resonance constants execute afresh. They do not read saved certificate values.
* `example.py`: explicit experiment composition, source comparisons, exact certificate exports, numerical trajectories, resource accounting, dynamic readouts, figures and metadata.
* `test_example.py`: seven scientific groups covering source conservation, arbitrary site composition, static/dynamic disagreement, freshly replayed certificate signs, buffering/relaxation, reversible-current derivatives and finite-time ledgers.

```python
import numpy as np
from example import design
from reactor import MassActionReactor, SinusoidalBinding

four_sites = design().append_site('1/100')
reactor = MassActionReactor(four_sites.exact())
y0 = np.zeros(reactor.dim)
y0[2] = .001
solution, species = reactor.integrate(y0, 1000, ledger=True)

# Independent rate perturbation preserves totals, not the old equilibrium.
on = reactor.kon.copy()
on[0] *= 1.02
changed = reactor.with_rates(on, reactor.koff, reactor.kcat)

# External selective modulation is a separate declared intervention.
forced = MassActionReactor(design().with_reverse('3/2').exact(),
    SinusoidalBinding(amplitude=.001, frequency=.23508))
```

A changed design can remove the oscillation. The demonstration then refuses failed shooting instead of relabeling an equilibrium as a cycle. The supplied driver locates repeated peaks after a transient to initialize the period; its fallback period guess is not evidence. For a broad sweep, catch a failed orbit search and retain it as unresolved numerical evidence, or write a driver that only evaluates trajectories and spectra. `static_complexes` returns a local physical numerical root when it converges; it neither enumerates all branches nor certifies uniqueness globally.

## What the example establishes

The rational certificate reconstructs the characteristic polynomial affine in r, isolates the unique positive frequency-squared root with exact polynomial arithmetic, encloses r and omega, verifies a positive Routh column for the complementary factor, and encloses the crossing derivative and first Lyapunov coefficient using outward-rounded 160-bit dyadic rational arithmetic. Complex left eigenvectors use the kernel's bilinear normalization `p.T*q=1`, not an accidentally mixed Hermitian convention. The coefficient computation includes the two quadratic resolvent terms.

For witness A, r is about 1.433228306811195, omega about .235080199562195, the crossing is negative and l1 is about -.08915942680578. Together with the classical Hopf theorem this yields attracting small cycles on the r<rA side. Witness H has positive l1: finding a Hopf pair alone does not establish attraction. The kinase-buffered A case is likewise subcritical. The configured added-site source is certified independently; the all-n inheritance theorem is imported, not proved by checking a finite selection of site counts. The paper's n<=2 statement excludes Hopf bifurcations, not every possible periodic orbit.

The C2-eliminated reference solves one fast algebraic equation locally. Substitution induces a cubic term even though the original field is quadratic. The kernel retains that term and obtains a nearby attracting Hopf threshold about 1.43082896554. This is a local bifurcation comparison, not a bound on trajectory approximation.

`phos.py` and its interval dependency remain the manuscript implementations, so this is an executable replay rather than an independent formal verification. Correct Python rational arithmetic, SymPy characteristic polynomials and root isolation, and the mathematical Hopf theorem are part of the trust boundary. No Lean compilation, finite-amplitude interval orbit proof, or uncertain pulse-flow proof is performed. All JSON outputs state their evidence level.

## Why static elimination misses it

Multiplying a complex equation by a positive relaxation factor changes its timing but not its zero set. The matched ratio sweep therefore preserves the equilibrium, totals, productive fluxes and the whole local statically eliminated field. Its Jacobian is `L*A^-1*M`, whereas the full dynamic feedback is frequency-dependent: `L*(sI+A)^-1*M`. Removing the complex relaxation loses phase information. Exact reference arithmetic checks the static characteristic polynomial and Hurwitz inequalities; the editable experiment also compares nonlinear local static fields.

At r=1.3 the full equilibrium is unstable; at r=1.5 it is stable. Both share a stable static reduction. Uniformly accelerating all complex equations with relaxation=.5 removes the default oscillation numerically while retaining the static behavior. Arbitrarily changing one reaction is not covered by the uniform-relaxation conclusion. The code rejects relaxation factors that would require nonpositive dissociation rates.

Buffering both **free** enzymes turns their concentrations into external controls. This changes the physical conservation contract: kinase and phosphatase are supplied or removed by buffers, while substrate remains conserved. The resulting substrate/complex model is linear and has no oscillatory attractor by the compartmental argument in the paper. One dynamic enzyme pool can still admit Hopf, as the separate subcritical certificate demonstrates.

## Cycles, costs and readouts

At r=1, numerical shooting gives period 39.590676906, free S3 peak-to-peak about 3.34868 and largest nontrivial Floquet multiplier about .237083. The shooting return residual is around 2e-11. These are numerical results, not an enclosure, global basin claim, minimal-period proof or proof that this orbit is connected to the local Hopf branch. Extremely small computed Floquet multipliers are limited by floating-point precision.

The integrated catalytic ledger checks `ui(T)-ui(0)=integral(c_i*C_i)-integral(gamma_i*D_i)`. On the numerical cycle the per-site differences are below 1e-9 and forward kinase turnover is about 262.065 per period, or 5.924 per conserved substrate. Forward turnover is an ATP-equivalent count in the irreversible source, not automatically thermodynamic dissipation.

For the driven reversible construction, reverse catalytic rates are `kcat*kon/koff*exp(-M)`. Assigned standard potentials make each reaction satisfy local detailed balance, with ATP/ADP/Pi held at unit activity and cycle affinity 2M. The example checks those ratios and integrates **net** catalytic currents. The theorem guarantees persistence only for sufficiently large M and supplies no numerical M0. The default M=20 trajectory is illustrative numerical evidence. This construction varies the designed potentials with M; it is not changing ATP concentration in fixed chemistry. Closed finite fuel is a different system. The saved finite-time net fuel count is not multiplied by affinity and mislabeled entropy production: that cycle identity requires returning to the same state.

Near the stable side, modulation of the D1 binding rate has an exact linear gain law at the Hopf frequency: `(1+rA+delta)*K/delta`, with K about .483638 for free S3. The example recomputes K by interval arithmetic and compares direct matrix response against the rank-one identity. Large linear gain does not validate finite forcing or experimental feasibility. Nominal step trajectories at r=1.3 and 1.5 illustrate the paper's short-pulse discrimination, but this run does not certify preparation uncertainty, actuator variation or measurement-error separation. The paper's stronger interval experiment remains a separately stated result.

Illustrative C0=.1 micromolar and T0=100 seconds map the computed period to about 66 minutes. These are uncalibrated units. At the designed equilibrium almost all phosphatase is bound, catalytic constants span thousands-fold, and a femtoliter contains only about 24 free phosphatase molecules at this scale. A deterministic model may be unsuitable in such a small compartment.

## Outputs and license

JSON files contain rational certificate endpoints, polynomials, eigenvectors, numerical diagnostics and scope. CSVs expose all species and cycle ledgers. PNG/SVG figures are generated from those computations. `run_metadata.json` hashes the manuscript identity, source modules and outputs. Reproducibility is checked under the pinned dependencies; floating trajectories are not portable bitwise guarantees across arbitrary numerical platforms.

MIT is proposed for the new example, pending the owner's decision. No new license grant is asserted; the supplied manuscript kernels and repository/manuscript terms remain separate.