# Feedback-induced bistability and global state selection

This standalone example implements the four-species reactor in *Feedback-induced bistability and global state selection in coupled autocatalytic cores* (September 2026). The original PDF has SHA-256 `83d9a0ecb93ca62c09b473fd55615b34712b44e5d6077598706b156d4d01cb44`.

Two modules exchange material through a fork reaction. The AB module contains `A <-> B + z` and `B + G <-> 2A`; the ZH module contains `z + F <-> H` and `H <-> 2z`. Reservoirs supply A and B, and H leaves through a waste channel. F, G and the reservoirs have maintained unit activities. These are dimensionless theoretical parameters from the manuscript, not experimentally calibrated chemistry.

The feedback intervention splits the fork into dynamic `A <-> B + z` with coefficients `(g,g)` and buffered `A <-> B + Zb` with `(1-g,1-g)`, where `Zb=1`. Forward and reverse coefficient sums remain one. Actual currents, throughput and dissipation need not stay fixed. At weak feedback the stationary response is unique; sufficiently strong feedback supports two attracting compositions. In the fully coupled flagship family, the same core-activity labels occur at both sinks and therefore cannot identify which composition was selected.

## Run and reproduce

Use Python 3.11 or later in a virtual environment:

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

The default run saves equilibrium states, rational root enclosures, a numerical feedback sweep, two trajectories, activity margins, the potential relative to the saddle, a recomputed exact fold certificate, PNG/SVG figures, console output and hashed run metadata. It needs no manuscript files or network connection after dependencies are installed. There is no random sampling in the main run.

Editable inputs are at the top of `example.py`. Concentrations are ordered `(A,B,z,H)`; time uses the reference reverse-rate unit. The default reverse-core rate is `1/100000`, output loss `1/10000`, and feedback `g=1`. Two initial states are constructed as the saddle plus/minus `0.02` times its unit unstable eigenvector, oriented toward increasing B. `CUSTOM_INITIAL_STATES` adds positive states to this experiment. The 1200-unit horizon is an illustration, not a uniform convergence-time guarantee.

## Components and reuse

`ReversiblePair` owns the mass-action current and stoichiometric change of one reduced reaction. `Parameters` validates rates and tracks whether they lie in the flagship theorem's scope. `CoupledReactor` supplies the complete field, analytic Jacobian, actual-current production margins, and stiff integration. The field is expanded for integration speed; tests compare it with the reaction sum.

`StationaryResponse` reconstructs the other concentrations from a candidate z and isolates every real root of the cleared polynomial using exact rational arithmetic. It rejects nonpositive reconstructions: a positive algebraic root alone is insufficient. It also exposes the independently clamped AB and ZH responses while retaining reservoirs and loads. The stationary residual is **not a scalar ODE for z**. The plotting sweep uses floating-point polynomial roots and full Jacobian eigenvalues; it is explicitly an illustration and can miss extremely close roots near a fold.

`ResponsePotential` implements the definite integrals and residuals in equations (29)-(33). It supplies the potential, its derivative along the full four-dimensional field, and the nonnegative dissipation lower bound. `NumericalSelector` requires membership in the absorbing region, a strict potential gap below the saddle, and a strict B-side test. It reports `unresolved` or `outside_absorbing_region` otherwise. A tolerance excludes tiny numerical gaps, but does **not** turn floating-point integration and quadrature into a rigorous enclosure. No nearest-equilibrium assignment is used.

`Interval` implements rational interval arithmetic. `FoldCertificate` independently recomputes the fixed Appendix B contraction for `(Q,Q_z)`, the derivative signs, positive reconstruction, the determinant identity, and the cubic Hurwitz inequalities. Its implementation adapts the manuscript's `scripts/certify_fold.py`; all acceptance tests use fractions. It proves one local fold in the specified box, not uniqueness of all folds or a complete bifurcation diagram. The certificate remains fixed at the paper's parameters when exploration inputs change.

For example, explore a changed output load without silently applying the flagship potential theorem:

```python
from dataclasses import replace
from fractions import Fraction as Q
from example import Parameters, CoupledReactor, StationaryResponse

base = Parameters()
for loss in (Q(1,10000), Q(1,10), Q(1)):
    reactor = CoupledReactor(replace(base, d=loss))
    for root in StationaryResponse(reactor).exact_roots():
        state = root['state']
        print(loss, state, reactor.activity(state))
    path = reactor.integrate((10, 20, 1, 9), end=200)
    print(path.y[:, -1])
```

For feedback sweeps, replace `g`; for a reversible waste channel, replace `eta` with a positive rational. Rate edits preserve this network's structure. To introduce additional reactions, extend the reaction collection **and** update the expanded field and Jacobian; the tests are intended to catch inconsistencies. Stationary elimination, potential and certificates are specialized to this architecture and must be rederived after structural changes. Composability here separates kinetics, stationary analysis and observation; it does not imply arbitrary networks inherit these formulas.

## What the default results establish

Exact root isolation finds three positive equilibria with z approximately `0.9957940123`, `2.0248399792`, and `2.9763672438`. Numerical Jacobian spectra classify the outer states as sinks and the middle state as a saddle with one growing mode. Tests repeat root isolation at both endpoints of the flagship reverse-rate interval, at weak/strong routed endpoints, and in the `d=1` uniqueness regime. They check the literal reaction field, clamped balances, activity balance, potential derivative, sampled dissipation, opposite trajectory destinations, parameter-scope guards and exact local fold inequalities.

The fold lies in `|z-1.617456101591| <= 1e-9`, `|g-0.970187992295| <= 1e-10`. The exact contraction bound is below `1.10e-7`; `Q_g>0` and `Q_zz<0` orient the locally created sink/saddle pair toward increasing g.

At both sinks the AB B-production is negative; both ZH productions are positive. Their concentrations and flux magnitudes differ. The identity `p_AB,A + 2 p_AB,B + p_ZH,z = dz/dt` explains why both distinguished cores cannot be strictly productive simultaneously at a stationary state.

The global convergence, analytic basin boundary, measure-zero separator and eventual completeness of the strict selection test are manuscript theorems for **g=1, eta=0, (a,b,u,v,d)=(6,27,16,2,1/10000), e in [1/200000,1/50000]**. This package does not rerun the Lean proofs, certify finite-time flow enclosures, compute the four-dimensional separator, or infer inheritance under molecular noise. The global theorem is not transferred to routed or reversible-sink perturbations. The potential class rejects parameters outside its stated scope.

The reaction model, numerical parameter sweep and initial-state experiments remain usable outside that theorem scope; conclusions there require their own analysis. Output metadata distinguishes exact certificate data from numerical diagnostics.

## Reuse terms

MIT licensing is proposed pending the owner's decision. This example does not independently grant a new license; confirm repository terms before redistribution.