# Cycle structure, capacity and resource costs of phosphorylation memory

Example for manuscript 72, SHA-256 `7970bdb2f32090091d3a3a5dc2bc73e835db7b7d23165410ff053eb5aa2175ea`.

The paper distinguishes equilibrium count, stable-state capacity, operational labels and information capacity. This example reproduces its seven-equilibrium/four-sink source, then exposes the cost of turning those sinks into locally recoverable, continuously retained, readable labels. It also supplies a full reaction model suitable for independent numerical and stochastic exploration.

## Run

Python 3.11 is the tested runtime. In a fresh environment:

```sh
python -m venv .venv
# Activate .venv for your shell.
python -m pip install -r requirements.txt
python example.py
python -m unittest -v
python example.py --output replay
```

Allow a few minutes. The default driver replays 776 source checks, builds the exact equilibrium census, runs a moderate-rate full reactor and one finite-count path, and writes figures and tables. No saved output is accepted as a computed result. The supplied rational matrices, rate lists and count vectors are witness inputs and are checked afresh.

## Two distinct sources

The fixed published source has three sites, eight phosphoforms, twelve edges, 34 species, 72 irreversible reaction channels and a 31-dimensional conserved class. Its totals are `(E_T,F_T,S_T)=(10,1,264490005424)`. Its scale separation is extreme. Double-precision eigenvalues are not used to classify its equilibria, and simulating a short trajectory would not establish its very slow recovery.

The editable source is a **separate symmetric lift** with moderate illustrative rates, totals `(2,2,20)`, and initially unmodified substrate. These are transparent placeholders in arbitrary concentration/time units, not measured biochemical parameters. It provides a practical full-network starting point; it does not claim four stable labels. Edit `CHAIN_RATES`, inventories, starting phosphoform, clock, reverse activity, horizon or system size at the top of `example.py`. Supply one six-rate tuple per site. Decimal strings become exact rational constants. The explicit network builder supports up to four sites; exact cofactor elimination is limited to eight phosphoforms.

The certificate size and clock multipliers operate on the **published** source. The loading-kernel shift is a separate fixed-target family exploration. None of these edits silently transfers a stability or probability certificate to a new source.

## Reusable objects

`EdgeKinetics` holds the six elementary rates on one edge: kinase association, dissociation and catalysis, followed by the phosphatase counterparts. Its loadings are `p=a/(b+c)` and `q=alpha/(beta+gamma)`; effective rates are `c*p` and `gamma*q`. Elementary rates, loadings and effective rates have different roles.

`PhosphorylationNetwork` constructs a chain or cube with **one C and one Y complex per edge**, one shared free kinase and one shared free phosphatase. Species order is all free phosphoforms, `E,F`, all `C`, then all `Y`. Each binding step consumes one free substrate and one free enzyme; each release or catalytic step returns enzyme. Its inventory matrix retains substrate in every complex. It exposes stoichiometry, reaction rates, the analytic Jacobian, a conserved-class chart, the full ODE and a literal direct SSA. Binding reactants are distinct species, so propensity is `k*N_i*N_j/Omega`, with no factor one half. Unimolecular propensity is `k*N_i`. Event-budget exhaustion raises an error instead of returning an apparently completed path.

`symmetric_lift` divides kinase association at level `k` by `n-k` and phosphatase association by `k+1`. At zero reverse activity, summed propensities for every aggregated jump equal those of the original chain, even for nonuniform microstates. This is tested with exact fractions, independently of the numerical trajectory comparison. The aggregation does not cover arbitrary site-specific perturbations or observations. At positive reverse activity the extension remains a full model, but the induced reverse aggregation factors differ; the driver reports no same-rate chain comparison in that case.

`EquilibriumGeometry` constructs matrix-tree polynomials, the loading polynomials and the eliminant. It isolates all real roots exactly and checks positivity separately. An algebraic root is not automatically a chemical equilibrium. The removable reconstruction case `u=E/F=E_T/F_T` is treated by its remaining substrate-total equation rather than discarded for dividing by zero. High-precision state reconstruction supplies useful model inputs; it is labelled numerical, not an interval certificate. The irreversible eliminant explicitly rejects a reversible model.

`PublishedSource` exposes the full exact elementary rates and a loading-kernel shift. `QuadraticMemoryContract` turns **already verified** local constants into sufficient size and error accounts. `PublishedOperatingContract` connects these reusable formulas to the source's checked constants. Passing arbitrary positive numbers to the formula class does not verify a Lyapunov hypothesis.

Example of an independent full-network model:

```python
from fractions import Fraction as F
import numpy as np
from reactor import EdgeKinetics, PhosphorylationNetwork

levels = [EdgeKinetics(F(1), F(1), F(2), F(1), F(1), F(1)) for _ in range(2)]
model = PhosphorylationNetwork.symmetric_lift(levels)
x = np.zeros(model.size)
x[0], x[model.E], x[model.F] = 20, 2, 2
path = model.simulate(x, np.linspace(0, 100, 501))
print(model.inventory @ path[-1])
```

To use the exact published rate list, instantiate `PublishedSource().network`. Its enormous inventories and very slow modes demand appropriate scaling and precision; the exact stability replay is the evidence for its four sinks.

## What is rebuilt

`paper_checks.py` is an attributed adaptation of the supplied read-only checker. The data path is flattened, and `main()` returns its report instead of writing beside the witnesses. All arithmetic is retained. `ExactPaperReplay` selects the full mode, which checks both original and final source sink certificates. Scientific output removes wall-clock timing fields.

The replay reconstructs tree weights by Laplacian cofactors and verifies the loading target `B-10D=product(u-j,j=1..6)`. It recomputes the degree-15 eliminant, finds eight positive roots and rejects the inadmissible one above 10. Four rational Lyapunov matrices and rational whiteners prove strict positivity and strict full-class Hurwitz inequalities using outward interval arithmetic. The separate geometry census also excludes zeros of its reconstruction denominators throughout each root interval. Simple roots and the paper's degree theorem then identify the remaining three equilibria as unstable with index -1; a scalar crossing alone is not used as a stability test.

The checker reconstructs the local quadratic remainder, positive-class radii, readout spacing, local and global jump-noise bounds, sufficient system size and error bounds. It refines root intervals to verify all four explicit 34-species integer preparations at 220-digit interval resolution, including exact inventories and membership in their preparation ellipsoids. Physical conversions and robustness/driving bounds are replayed. Floating-point comparisons in the source checker are confined to display rounding of the resource table; exact inequalities determine its scientific certificates.

The source has four normalized readouts approximately `3.31e-21, .316355, .558771, .670451`, with minimum gap above `.1116`. Effective square affinities are nonzero, but their value is not an ATP free energy. Square-balanced effective rates have capacity exactly `n`; unrestricted capacity is bounded, not known exactly. In particular, the example shows `4 <= C_3 <= 8`, not `C_3=4`.

## Operating accounts and their limits

The local theorem uses an ellipsoid `V=h^T P h < v`, drift decay `lambda` and jump-noise bound `C`. Continuous retention failure is bounded by `1/4096 + C*T/(Omega*v)` from the preparation set. Recovery from the declared larger disturbance set has its own bound, and their sum bounds recovery followed by storage. An excursion out of the safe set counts as failure even if the path returns. These are sufficient local bounds, not a global basin theorem or a writing protocol.

For the final source the sufficient size is about `1.31e100`; total substrate molecules are about `3.5e111`. At total substrate 1 micromolar and after uniformly slowing all rates to the stated upper ceilings, certified recovery takes about `1.4e31` seconds and volume about `5.8e93` litres. Retention and recovery bounds are respectively below `.00105` and `.008987`, and their composite is below `.011`. These numbers diagnose the construction and its conservative certificate. They are not minimum resource requirements, even for this source, and are not impossibility results for phosphorylation memory.

The resource comparison uses each source's **own** matched horizon `T=100*t_recovery`. Multiplying all rates by a common factor scales both drift and noise and inversely scales those times, leaving the molecule budget unchanged. The driver verifies this identity exactly. It would not hold as the same statement with a fixed absolute storage horizon. When an edited size is too small, failure bounds above one remain visible as uninformative bounds, rather than being converted into a claim of failure.

The fixed-target obstruction is separate: with the specified effective rates, target, enzyme totals and ratio interval `[4.1,4.2]`, free substrate exceeds `1.49e11`, forcing kinase/substrate below `6.8e-11`. A loading-kernel shift cannot remove it. Such a shift preserves free substrate **at a fixed enzyme ratio**, not the equilibrium locations for a fixed substrate total. Changing effective rates, target, retained branch or enzyme totals falls outside the obstruction.

`crowding_scaling_reference.csv` illustrates only the exponents for the paper's different three-site sequential split family: barrier proportional to `delta^6`, recovery to `delta^-4` and readout gap to `delta`. Prefactors and finite-size errors are not calculated. The theorem concerns a specified local tube and first takes size to infinity at fixed separation, then separation to zero. It is not a uniform joint-limit or global exit theorem. The four-label source here is well separated, so its cost is not attributed to crowding.

## Driving and observation

Setting positive `REVERSE_ACTIVITY=theta` adds the two reverse catalytic channels on every edge. Their constants are `a*c/b` and `alpha*gamma/beta`, multiplied by theta. At unit free activities the closed pairs balance with `C=a/b` and `Y=alpha/beta`; tests verify this exactly. The shared chemostat potential is `-2*log(theta)` in thermal units. The source certificate guarantees four deterministic sinks for `theta <= 1e-102`, a very conservative sufficient threshold, not a minimum energetic requirement. Its finite-molecule guarantee is **not** transferred: reverse channels change equilibria, drift and noise.

The scalar packing expression is an upper bound for bounded observation error, not an attained information capacity. Ideal-center separation is checked separately from the safe-set variation required by the source's actual decoder. Passive observation does not improve chemical retention. A binding reporter would need extra reactions and retained inventory.

## Reproduction, files and attribution

The package contains source witness JSON, four paper-local modules, the driver, tests, pinned dependencies and computed outputs. Outputs include the exact replay report, independent root census, reconstructed full states, pathway ratios, operational/resource accounts, reaction list, full ODE path, sampled SSA path, capacity bounds and two PNG/SVG figures. Hashes identify the manuscript, modules and witness inputs.

Seven tests check full inventory/Jacobian identities, exact deterministic and stochastic aggregation, the exceptional equilibrium reconstruction and admissible root census, target invariance and detailed-balance extension, common-clock/composite accounting, actual ODE and SSA conservation, and invalid inputs. The driver separately executes the long exact stability/operating replay.

Original files and SHA-256 values are recorded in `source_provenance.json`; witness bytes are unchanged. The adaptation is local file layout/reporting only. The analytic capacity, degree, perturbation and probability theorems remain conventional mathematics. The supplied Lean algebra is not rerun, and this is not an end-to-end formal proof.

License: MIT is proposed for new example code, pending the owner's choice. No license grant is made here; attribution does not resolve the supplied research code's licensing.