# Coexistence, composition and a consumed precursor

Example for **Permanence and composition control for self-limited consumers coupled to a dynamic resource network**, manuscript entry 31, SHA-256 `8f00f7290e3e9a3f9a9d437d54ae386a3755ec295b37a4ef171d3ea3d449a835`.

Four resident species A, B, z and H produce a resource z. Consumers Xi copy themselves using that resource. Each has its own quadratic loss, so a temporarily abundant consumer limits itself more strongly. System A maintains the copying precursor externally. System B tracks a precursor R which is supplied, washed out and consumed by copying. The paper proves that every consumer eventually stays above a positive concentration from every strictly positive initial state in its specified parameter regime, even with the dynamic reservoir.

This package supplies the literal reactor, independent directed rates, prescribed composition, integrated resource ledgers, numerical stationary continuation, finite-inventory and competitive-exclusion controls, symbolic global floor formulas, and exact rational replays of two conventional certificates. It is deterministic concentration dynamics, not a finite-molecule survival model.

## Run

Python 3.11 or newer:

```sh
python -m venv .venv
# Windows: .venv\Scripts\activate
# macOS/Linux: source .venv/bin/activate
python -m pip install -r requirements.txt
python -m unittest -v test_example.py
python example.py --output outputs
```

The default run takes tens of seconds. It writes five trajectory CSVs, a stationary-branch CSV, two figure pairs, two certificate JSONs, results, console output and provenance hashes. Seven scientific test groups check independently written field equations, exact identities, certificates, numerical balances and limiting cases. The ordinary differential equations are solved with Radau; an independent BDF integration checks the main trajectory. No Lean proofs or historical receipts are rerun.

## Inputs and reusable components

Edit the inputs at the top of `example.py`: resident reversible rate e, reservoir supply d, initial concentrations, desired consumer shares, shared self-limitation normalization kappa, integration horizons and supply sweep. The defaults reproduce the paper's numerical illustration: e=1e-5, d=0.05, residents (3,10,0.05,3), R=0.001 and consumers (1e-8,1). Concentrations and time are nondimensional. These are illustrative source parameters, not measured chemistry.

`Resident` stores thirteen independently editable directed rates, ordered as (a,b,p,q,alpha,beta,ef,er,u,v,h1,h2,dH). `Consumers` stores copying, linear loss and quadratic self-limitation vectors. `Reservoir(feed,washout)` separates precursor feed and washout; `Reservoir.reference(d)` sets both to d. `Reactor(resident,consumers,reservoir=None)` selects system A. Supplying a reservoir selects B. Changing a model does not automatically preserve a theorem's hypotheses.

The thirteen resident reactions are:

| Reaction | Reference rate |
| --- | --- |
| empty -> A; empty -> B | 6; 27 |
| A -> B+z; B+z -> A | 1; 1 |
| A -> empty; B -> empty | 1; 1 |
| B -> 2A; 2A -> B | e; e |
| z -> H; 2z -> H | 16; 2 |
| H -> z; H -> 2z; H -> empty | 1; 1; 1e-4 |

For each consumer, system A adds Xi+z -> 2Xi, Xi -> empty and 2Xi -> Xi, at rates ki, mui and rhoi. System B uses Xi+z+R -> 2Xi and adds precursor feed and washout. These are monomial deterministic rates: the 2Xi loss is rhoi*Xi^2, without a factorial divisor. R washout is not a common reactor dilution applied to all species.

`Reaction` records integer reactant/product multiplicities and implements fluxes and analytic derivatives. `Reactor.field` and `jacobian` accept exact rational or symbolic values; `numerical_field` uses the same reaction list efficiently. `with_rates` replaces the complete ordered vector: thirteen resident rates, then (copy,loss,limit) for each consumer, then reservoir feed/washout if present. This permits independent parameter sweeps without rewriting the equations.

`Simulator` integrates logarithms of consumers, keeping rare populations positive without clipping or artificial immigration. Residents and R are ordinary concentration coordinates. Invalid initial data, solver failure, nonpositive saved values or exponential underflow raise errors. Very long extinction simulations may require a shorter horizon or a different state representation. Zero consumer faces are invariant; the paper's permanence claim requires strictly positive initial values.

`Trajectory` stores concentrations and six integrated rewards: total abundance, copying uptake, reservoir washout, linear loss, quadratic loss and excess quadratic loss due to composition mismatch. The balance check uses integrated rewards rather than a coarse CSV quadrature. Variance loss is defined only when consumers have common growth and positive self-limitation; its column is zero in the deliberately incompatible negative control.

## Composition is a design parameter

With common environmental growth g, consumers obey Xi'=Xi(g-rhoi Xi). Define kappa=1/sum(1/rhoi) and qi=kappa/rhoi. The paper proves Xi/sum(Xj) -> qi. Choosing rhoi=kappa/qi therefore prescribes any strictly positive composition. On that target manifold the total S obeys S'=S(g-kappa*S).

```python
from example import *
consumers = Consumers.target(('1/4', '3/4'), kappa=1)
model = Reactor(Resident.reference(), consumers, Reservoir.reference(.05))
run = Simulator(model).run(model.initial(), np.linspace(0,1200,1201))
```

Using exact fraction strings avoids ambiguous share sums. For more consumers, give matching initial concentrations and target shares. The two-consumer exclusion control remains a separate two-consumer model.

The equal target uses rho=(2,2); the 1:3 target uses (4,4/3). Both have kappa=1 and the same reduced operating equations. Their off-target transients can differ: composition mismatch adds positive quadratic loss. For equal consumers, the composition equation becomes autonomous in the accumulated-abundance clock tau=integral S. Low abundance delays restoration in physical time.

The default reservoir equilibrium has total S approximately **0.0406741019**. Equal consumers each have approximately **0.0203370509**; the 1:3 target has approximately **0.0101685255** and **0.0305055764**. `OperatingAnalysis.spectral_split` checks the environmental/total block and n-1 composition eigenvalues -kappa*S. The paper prints the kappa=1 case -S; the general expression follows by differentiating the same target equations. Stable composition modes do not prove stability or global convergence of the environmental block.

## Resource accounting and supply limitations

For reference system B, with J=R*z*S,

    R' = d(1-R)-J
    (R+S)' = d-dR-S/2-sum(rhoi*Xi^2).

The default numerical reservoir balance error is below 3e-11; Radau and BDF concentrations differ by less than 2e-8. At the operating state about **43.98%** of precursor inflow enters copying; the remaining **56.02%** leaves as unused R washout. These percentages concern only the reservoir balance, not all material or energy flows in the resident network.

The weighted variance identity is

    sum(rhoi*Xi^2) = kappa*S^2 + sum(rhoi*(Xi-qi*S)^2).

Integrated resource accounting yields the asymptotic time-average ceiling (sqrt(1+16*kappa*d)-1)/(4*kappa). `OperatingAnalysis.ceiling` evaluates its cancellation-resistant equivalent. The default ceiling is **0.0854101966**, greater than the simulated stationary total. It is an asymptotic time-average bound, not an instantaneous ceiling or a sufficient design condition. Finite-time averages include initial inventory divided by the elapsed time.

At d=0 the consumers tend to zero. The finite-inventory control checks integrated uptake <=R(0) and integrated total abundance <=2[R(0)+S(0)]. A positive floor uniform over all d>0 is impossible. Removing quadratic losses and using linear losses (1/2,3/4) instead gives the exact ratio X2/X1=(X2(0)/X1(0))*exp(-t/4), demonstrating exclusion outside the theorem's rate cube.

`OperatingAnalysis.equilibrium` solves positive concentrations in log coordinates and checks per-capita consumer residuals, which avoids accepting a zero-consumer root. The supply sweep is numerical continuation, not a certified global branch. `donor_certificate` recomputes the exact scalar enclosure near z0=0.9957940123, its derivative sign and positivity. This supports the paper's conventional implicit-function argument: locally S*(d)=2(1-1/(2z0))*d+O(d^2). It does not give an all-trajectory floor proportional to d.

## Global and local certificates have different scopes

`PermanenceBounds` evaluates the paper's reference and robust aggregate/recovery formulas as **exact rational prefactor times exp(-rational exponent)**. The exponents are around 3e11, so ordinary floating-point floors would underflow to zero. These formulas invoke the paper's resident-potential and permanence results; the package does not reprove those global results. Reference e must lie in [1/200000,1/50000], n>=2, and B needs a compact positive supply interval. Robust independent-rate radii also include aggregate_floor/1000: taking only 1e-18 is insufficient. The displayed consumer recovery term is not the minimum floor across resident and reservoir species; the manuscript supplies those additional coordinate bounds.

`local_certificate_inputs.json` contains the paper's exact seven-coordinate center and rational P matrix, box radius 1e-5, rate tolerance 1e-14 and rational square-root lower bound. These are fixed certificate inputs, independent of editable simulation defaults. `certificates.py` derives the full field from the reaction list and recomputes Sylvester minors, norm bounds, interval Jacobian/field errors, inward boundary margin and contraction inequalities using exact rational arithmetic. It also verifies that the equilibrium lies in the ellipsoid. A larger attempted box reports `accepted=False`; that means the certificate fails, not that the model is unstable.

The certified ellipsoid is E={(x-c)^T P (x-c) <= p_min*(1e-5)^2}, where p_min=1/||P^-1||_infinity. For every independent perturbation of all 21 rates within 1e-14, trajectories starting in E stay there and converge to its unique equilibrium, with each consumer at least **0.0203270509477**. This is an exact arithmetic replay supporting a conventional analytic proof, not a Lean certificate. It does not certify when a plotted trajectory enters E and does not assert global equilibrium convergence.

## Physical interpretation and extension

`UnitScales` converts nondimensional concentrations, times, rates and fluxes consistently: concentration C*x, time T*t, an order-m rate k/(C^(m-1)*T), and flux C*J/T. The illustrative C=1 micromolar and T=1 hour give copying flux about 0.0219914 micromolar/hour. These chosen units do not calibrate the reaction model.

The resident z<->H and 2z<->H channels cannot both conserve positive molecular masses without omitted species. System B makes one precursor dynamic; it is still an externally driven effective model, not an atom-balanced autonomous apparatus. Its copying step is effectively trimolecular. Adding intermediates requires a justified new model. Loss flux is not harvested product unless a collection channel is specified.

Reuse the reaction and reactor objects for new rate vectors, the simulation/ledger layer for resource budgets, or the rational interval machinery for other finite polynomial certificates. New mechanisms and heterogeneous copying/loss rates need their own composition and permanence analysis.

## Provenance and license

The manuscript PDF is unchanged. Exact certificate inputs transcribe Appendix B and agree with the supplied research certificate; no private proof workspace is needed. Saved metadata hashes the executable source, certificate implementation/input and outputs. Replaying scientific JSON/CSV is checked in a clean extracted package.

MIT is proposed for the original companion code, pending the owner's license decision. This README does not grant a license or change rights to the manuscript or dependencies.