# Minimal cusp networks and a reusable bistable reactor

Example for manuscript entry 40, **Reaction-minimal cusps in planar bimolecular mass-action networks: an exhaustive certified classification**, SHA-256 `9c35a17e2f132b80f00d057d1ee52f1ba89391f78fd752be2d120f43b0ae086c`.

A cusp organizes a change between one and three nearby equilibria as two parameters vary. Here the parameters are independent reaction rate constants. The example combines an exact catalogue of small reaction mechanisms, fresh checks of all 52 positive cusp witnesses, and a fully usable two-species mass-action reactor. The classification counts network classes, not cusp points or parameter regions.

## Run and editable inputs

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 full run takes roughly a minute, depending on the machine; most work is exact enumeration and checking 52 algebraic templates. Scientific tests take about twelve seconds. The package is self-contained, including `catalogue.json`. It does not need the research workspace or Lean.

At the top of `example.py`, edit `EPSILON`, five optional `CUSTOM_RATES`, positive initial concentrations, time horizon and sampling, unfolding sweep inputs, and `QUERY_REACTIONS`. Rational strings retain exact arithmetic. The default epsilon=1/10 produces three rational equilibria. `CUSTOM_RATES` replaces all five rates of template 02; its equilibrium solver remains valid for arbitrary positive values. The unfolding sweep is a separate study of the specified mu/nu family, even when custom reactor rates are supplied.

Set the two exhaustive-check flags false for repeated kinetic exploration; template 02 is still checked. The census budget raises an error if exceeded, never presenting partial counts as exhaustive. All quantities are nondimensional mathematical examples, not fitted chemical measurements. The long default time horizon resolves slow relaxation close to the cusp; additional early-time samples resolve the initial transient.

## Literal chemistry and the working reactor

Template 02 consists of

```text
0 -> X          k1
X -> 0          k2
Y -> 2X         k3
2X -> X+Y       k4
X+Y -> Y        k5
```

The empty complex represents exchange with an external reservoir. Mass action gives

```text
x' = k1 - k2*x + 2*k3*y - k4*x^2 - k5*x*y
y' = k4*x^2 - k3*y.
```

`PlanarNetwork` owns the literal reactants, products, stoichiometry, structural filter, canonical key and equivalence transport. `MassActionReactor` adds independent positive rates, the full nonlinear field, its analytic Jacobian and positive-domain integration in log concentrations. `RationalCuspReactor` specializes the five reactions above and adds an exact cubic equilibrium solver. These components can also be used without running the catalogue demonstration.

For any positive rates, stationary y=(k4/k3)x^2. Substitution into x' gives

```text
P(x) = k1 - k2*x + k4*x^2 - (k5*k4/k3)*x^3.
```

The solver isolates its real roots with rational intervals and multiplicities. It verifies the derivative sign on each simple-root interval. At a positive equilibrium, `det J=-k3*P'(x)` and `tr J=-k2-2*k4*x-k5*y-k3<0`, so a positive determinant gives local asymptotic stability and a negative determinant a saddle. Repeated roots are labelled nonhyperbolic. Root midpoints are used only for decimal coordinates and explicitly labelled midpoint trace/determinant exports.

This is an **equilibrium reduction**. The nullcline y=(k4/k3)x^2 is generally not an invariant curve: y'=0 there, while its tangent requires y'=(2k4/k3)x*x'. The code integrates both original differential equations, not a scalar surrogate.

## Exact unfolding and bistability

At the cusp, x=y=1 and rates are `(1,3,3,3,1)/11`. Hold k3=k4=3/11 and k5=1/11, and write

```text
k1=(1+mu-nu)/11,    k2=(3-nu)/11,    x=1+z.
```

The stationary equation becomes exactly `11*x'|_(y=x^2)=mu+nu*z-z^3`. The rate map is invertible and rates must remain positive. The fold curve is `(mu,nu)=(-2z^3,3z^2)`, with three roots when `27mu^2<4nu^3`. The plot evaluates equilibrium branches; it is not a time-varying parameter protocol or a measured hysteresis loop.

For mu=0 and nu=epsilon^2, 0<epsilon<1, the three equilibria are

```text
(1-epsilon,(1-epsilon)^2), (1,1), (1+epsilon,(1+epsilon)^2).
```

Their determinants are `(6,-3,6)*epsilon^2/121` and their traces are `(8epsilon-13)/11`, `(epsilon^2-13)/11`, `-(8epsilon+13)/11`. Thus the outer equilibria are stable and the middle is a saddle. The script derives these identities symbolically. At epsilon=1/10 the states are `(0.9,0.81)`, `(1,1)` and `(1.1,1.21)`. Two full nonlinear trajectories approach different outer states; Radau and BDF agree within about 4e-10 in concentration. This illustrates bistability without claiming a complete global basin classification.

## What the catalogue classifies

The universe contains the six complexes `0,X,Y,2X,X+Y,2Y` and their 30 directed non-self reactions. `reactions.csv` maps each zero-based ID to its reaction. A network is a set of five distinct reactions, with rank two and independent positive rates. Both reactants and products have molecularity at most two.

The script freshly enumerates all 142,506 five-reaction sets. The positive-flux surrounding condition and at least four distinct reactants leave 60,036 eligible sets. Quotienting by species exchange gives 30,051 orbits. The complete key includes each reactant exponent and the corresponding primitive **directed** stoichiometric vector; its sorted multiset is compared with its species-exchanged copy. This gives 9,999 mechanism classes, with multiplicities retained.

Equivalence allows reaction permutation, positive stoichiometric column scaling with the reactant fixed, and exchange of X/Y. It does not identify arbitrary dynamically equivalent networks or permit row scaling. `transport_from` supplies the actual matching, species orientation and rate factors. For a target column lambda times its matched source column, the target rate is the source rate divided by lambda, preserving the field after the possible species exchange.

The imported classification partition is checked against every regenerated eligible key:

| Certificate layer | Four sources | Five sources | Total |
| --- | ---: | ---: | ---: |
| Determinant obstruction | 7,205 | 1,858 | 9,063 |
| Quadratic-degeneracy obstruction | 663 | 81 | 744 |
| Cubic-degeneracy obstruction | 132 | 8 | 140 |
| Positive cusp template | 0 | 52 | 52 |

`Catalogue.classify` returns the paper's outcome and an explicit transport to the stored representative. **The 9,947 negative obstruction payloads are not re-proved here.** Their labels use the published classification theorem, while this package freshly checks the finite universe, quotient and source coverage. The negative labels describe which sufficient obstruction decided a class, not a hierarchy of degeneracies that the class can realize.

Every positive class has five distinct reactants. The conclusion that five reactions are minimal also uses the published four-reaction exclusion, not a new four-reaction search in this package. Arbitrary larger networks and trimolecular-product networks are outside the classifier's model.

## Fresh positive certificates and the corrected unfolding condition

`AlgebraicTemplate` reads an exact root polynomial m(t), rational isolating interval, five rate polynomials and vectors p,q,h. It checks that exactly one root lies in the interval, or that a rational singleton is the root of a linear polynomial. It then uses rational polynomial arithmetic modulo m to verify the literal equilibrium and jet equations

```text
Jq=0, pJ=0, pq=1, tau=tr J != 0,
pB(q,q)=0, Jh=-B(q,q), ph=0,
c=pB(q,h)/2 != 0.
```

B is the full second derivative, with no factorial factor. The implementation constructs it directly from the reactant exponents. It freshly checks exact interval signs of every rate, tau and c; all 52 witnesses have positive rates, tau<0 and c<0. Sixteen root polynomials have degree one, fifteen degree two and twenty-one degree three. The degree is not a claim of irreducibility.

For a rate direction u, the correct local unfolding rows are

```text
a(u)=p f_k(u)
b(u)=p J_k(u) q - p B(q,f_k(u))/tau.
```

The subtraction accounts for the transverse equilibrium shift when rates change. The formal predicate discussed in the paper uses only the raw row `p J_k(u)q`; it is not silently substituted for the corrected analytic condition here. The package recomputes the trace-scaled corrected minor D for every selected rate pair and verifies a supplied polynomial identity `U*D+V*m=1`. At the isolated root this proves D is nonzero. Strict signs are checked with exact rational interval arithmetic, not floating-point tolerances.

For template 02, tau=-13/11, c=-75/17303 and D=27/65; the corrected determinant for `(k1,k5)` is -297/845. The raw determinant happens to agree for this particular pair. For template 00 the determinants differ, and a test preserves that distinction. The direct cubic unfolding uses `(k1,k2)`, a different valid pair.

`AlgebraicTemplate.numerical_reactor()` provides an ODE model for any positive template by refining its root interval and substituting its rational midpoint. An optional equivalent target network transports those rates. The algebraic root is the exact cusp witness; a rational numerical approximation should not be called an exactly singular cusp. There is no claimed parameter radius of bistability for all 52 instantiated models.

## Why reactant labels matter

The positive network with IDs `(0,5,12,18,20)` and the negative network `(0,5,12,20,24)` have identical multisets of stoichiometric vectors. Replacing `2X -> X+Y` by `X+Y -> 2Y` preserves the vector (-1,1) but changes the kinetic monomial from x^2 to xy.

The negative network has field

```text
x'=a-bx+2dy-(e+f)xy
y'=((f-e)x-d)y.
```

At a positive equilibrium x=d/(f-e), its determinant is `-d(f-3e)y`. Singularity forces f=3e; then x=d/(2e), and the remaining condition is a=bd/(2e), independent of y. It gives a line of equilibria rather than an isolated cubic cusp. The package checks these identities exactly. Ray geometry without reactant labels cannot decide cusp existence.

## Outputs, provenance and verification scope

Outputs include all reaction IDs, the full trajectory/sweep CSVs, all 52 recomputed positive certificate summaries, exact equilibria, census counts, query transport, figures, console output and reproducibility hashes. `catalogue.json` contains complete positive rate/kernel/center data and the imported representative partition. Its provenance records source hashes for the extracted research data. All 52 recorded Lean witness-file hashes were checked against the supplied research workspace during extraction.

Seven scientific test groups cover the exact census, labelled transport and collision, rational/algebraic jets with rejection of a corrupted Bézout witness, cusp/fold/simple equilibrium multiplicities, arbitrary-rate scalar reduction, numerical template instantiation and independent bistable trajectories. The main run verifies all 52 positive templates; tests independently exercise rational, quadratic and cubic witnesses.

No Lean build is rerun. The paper's formal classification uses Lean 4.30.0, its frozen Mathlib environment and compiled `native_decide` evaluation, not exclusively kernel reduction of every finite computation. The analytic Lyapunov–Schmidt bridge and corrected-minor checks are separate conventional/exact-algebra evidence. This Python package preserves that distinction and does not claim to replace the formal proof or establish cusps in arbitrary supernetworks.

MIT is proposed for the example code, pending owner approval; this README grants no license and does not change manuscript rights.