# Equilibrium geometry and stable-state capacity

One kinase adds phosphate groups in a fixed order; one phosphatase removes them. Each enzyme releases substrate after every catalytic step. The full n-site model has 3n+3 species, 6n positive reaction rates and three conserved totals. Those totals define the class in which equilibria and stability are compared.

The paper determines two different capacities: at most **2n-1 positive equilibria**, and at most **n locally asymptotically stable equilibria**, in one class of one source. Its construction attains both together. It first fixes the equilibrium geometry using positive polynomial coefficients, then uses kinetic degrees of freedom that preserve every equilibrium to obtain the desired stability. The stable-state ceiling includes nonhyperbolic attractors but does not count periodic orbits or imply that all trajectories converge to equilibria.

## Run

Use Python 3.11 or newer in a virtual environment:

```sh
python -m pip install -r requirements.txt
python -m unittest -v test_example.py
python example.py
```

The complete ZIP includes all modules, exact benchmark data, scientific checks and saved outputs. Direct dependencies are pinned. The driver has editable inputs at the top and accepts `--output another-directory`. The default replays the full rational source checks, including all 48 equilibria for n=2 through n=7, and then runs the paper-local model and recovery demonstrations. Expect roughly a minute, depending on hardware.

## Build and reuse a model

| Component | Role |
| --- | --- |
| `Geometry` | Prescribed auxiliary roots x>1 and enzyme-total ratio r; the positive recurrence produces A, B, D and explicit equilibria with common totals. |
| `Kinetics` | Independent positive catalytic scales, dissociation constants and binding relaxation; retunes rates without moving the constructed equilibrium set. |
| `Reactor` | Full mass-action kinetics for any number of sites, arbitrary positive rates and totals, exact conservation coordinates, analytic Jacobian and numerical integration. |
| `stability` | Exact characteristic polynomial and regular Routh count; an exceptional/zero-pivot case returns unresolved. |
| `ordered_design` | The paper's successive choices of enzyme ratio, final-site current, binding scale and root separation, with explicit computation budgets. |
| `LocalRecovery` | Fresh rational-interval equilibrium displacement, positivity, Lyapunov, recovery and readout bounds for the three-site benchmark. |
| `phos_sharp`, `phos_capacity` | Exact polynomial, reaction, matrix and loaded-inventory kernels retained from the research source. |

For example, `Geometry(('2','5/2','3','7/2','4'), '4')` specifies five equilibria for a three-site model. `Kinetics.from_coalesced_currents(3, Q(4))` supplies the paper's coalesced-current kinetic choice; call `.rates(geometry)` and construct `Reactor(rates, geometry.construct()['totals'])`. Call `stability` at each state to check attraction. Counting roots alone never labels them stable.

The configured geometry and kinetic arrays at the top of the driver form a separate reusable experiment. Changing one requires matching the number of sites in all arrays. The published benchmark retains its own exact source and explicit certificate. For independent changes to its eighteen rates or three totals, use `Reactor` directly; the analytical certificate is available only within the checked relative parameter radius. No common framework across papers is required.

The ordered design returns a certified finite witness or `unknown_budget`. It does not translate a stopped search into mathematical nonexistence. The existence theorem guarantees termination without an a priori complexity bound; this implementation additionally limits sites and iterations to keep an accidental request manageable. Repeated auxiliary roots are allowed for coalescence calculations, but repeated entries are not counted as distinct equilibria.

## What the default computes

The three-site benchmark has five exact positive equilibria at the same rates and totals `(18, 2, 20)`. Rational characteristic polynomials and Routh arrays give unstable-direction counts **0, 1, 0, 1, 0**. The fully phosphorylated readout orders them, so the three sinks provide three distinct equilibrium labels. The output includes every exact concentration, rate, polynomial, slope and stability result.

The source replay checks the positive recurrence, negative interlacing roots, coefficient identities, full equilibrium equations, determinant formula, loaded matrices, feedback deficit and exact small-site witnesses. The loaded mass matrix retains enzyme-bound substrate; replacing it by free-substrate amounts changes the slow system, as `loaded_inventory_comparison.json` demonstrates. A source-preserving kinetic sweep changes binding speed while retaining the same equilibrium geometry.

`splitting_law.json` combines the exact positive coefficient in the paper's eigenvalue splitting law with fresh 60-digit numerical eigenvalues for four root separations. For this three-site family the critical eigenvalue scales as the fourth power of root separation. The inverse-fourth-power recovery cost is a property of this construction, not a universal speed limit. The limiting feedback-margin table makes another limitation visible: the positive margin decays rapidly with site count.

## Local robustness and recovery

The benchmark admits a jointly independent relative variation of all eighteen rates and three totals by `1/671088640000`, about 1.49e-12. The example freshly checks the rational Lyapunov matrices, reconstructs reactionwise interval bounds, verifies equilibrium displacement and positivity, and checks that the stated initial balls lie inside the invariant ellipsoids. It derives new readout intervals directly from the computed displacement bound and the chosen measurement error.

The return bound uses the exact matrix norm and the freshly evaluated Lyapunov loss. The source replay separately preserves the manuscript's slightly looser stored-loss calculation. Both certify tenfold decay of a local envelope within the reported initial neighborhood. They do not predict a typical recovery time or certify arbitrary preparations. The default readout intervals, including ±0.001 micromolar measurement error, are disjoint. Enlarging the parameter radius to 0.001 fails this certificate and is reported as **not certified**, not as proof that the equilibria disappear.

The nonlinear recovery traces use a one-percent reduction of free S3, transferring the same amount to free S0. This preserves all totals and keeps the initial state positive. These perturbations exceed the tiny certified preparation balls; their recovery is numerical and checked with both LSODA and Radau. Their reported solver agreement does not expand the rigorous basin bound.

The declared units are 0.1 micromolar and 10 seconds, with no fitted biochemical data. They give totals `(1.8, 0.2, 2.0)` micromolar. The source's worst local tenfold-envelope bound is below 32,918,000 seconds, roughly 381 days, whereas the nominal slowest numerical e-folding time is about 22 hours. These quantities answer different questions. The tolerance and preparation radii are far too small to represent practical laboratory robustness.

## Scope, provenance and license

The finite algebraic and rational-interval checks are rerun. The universal equilibrium ceiling and the degree/transversality argument for the stable-state ceiling are imported mathematical results, not consequences of the finite numerical sweep. Theorem 1.1 as a whole is not formalized in Lean; this package does not run Lean or claim otherwise. Catalytic arrows describe maintained driving. No explicit finite-fuel, stochastic retention, switching, inheritance or information-rate model is inferred from a count of deterministic stable labels.

Manuscript: `Sharp_Equilibrium_and_Stable_State_Capacities_Multisite_Phosphorylation.pdf`, SHA-256 `2cda97c29c4bf2773242ede84a0c57667a8aba2e1e90466844a6b5cdce51242f`. Source hashes and adaptations appear in `source_provenance.json`; executable, input and output hashes appear in `run_metadata.json`.

MIT is proposed for original example code, pending owner confirmation. This package does not grant that license or alter rights in imported source material.