# Overlapping siphons, resident invasion and uniform permanence

A reusable implementation of the fourteen-reaction model in *Overlapping siphons, resident-dependent invasion, and uniform permanence in a two-strain coinfection reaction network* (14 September 2026). Original PDF SHA-256: `3f0fc394b76fb14eae5f5e93292c0feea62e542904517384b4f69d86c36d3df3`.

The compartments are susceptible s, strain-1-only a, strain-2-only b and coinfected c. The model distinguishes three questions: how overlapping extinction sets share linear modes at one boundary state; whether a missing strain can invade a resident population; and when every positive solution eventually keeps every compartment above a common positive floor. A disease-free calculation alone cannot answer the latter two.

## Run and edit

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

Use Python 3.11 or later. The fourteen rates, initial states, simulation duration and independent uncertainty radius are editable at the top. Defaults are dimensionless manuscript examples, not a fit to an epidemic. `Parameters` stores exact rational rates and rejects zero/negative values because the support and the stated theorem require all fourteen rates to be positive.

The default output contains exact resident/covector certificates, the shared-mode polynomial identity, the one-percent parameter-box certificate, normalized growth weights, numerical trajectory and invasion sweeps, an unequal-mortality entropy example, PNG/SVG figures, and hashed source provenance. No numerical universal persistence floor is invented; the JSON records it as null.

```python
from dataclasses import replace
import numpy as np
from example import Parameters, CoinfectionModel, ParameterBox, NormalizedLift

p = Parameters(beta2='1')
model = CoinfectionModel(p)
print(model.permanence())
print(ParameterBox(p, radius='1/100').certificate())
states = model.simulate([0.5, 3.5, 1e-5, 1e-5], np.linspace(0, 180, 721))
low_transmission = CoinfectionModel(replace(p, beta2='1/10'))
print(low_transmission.permanence())
growth = NormalizedLift(model)
print(growth.growth(growth.lift([1, .2, .3, .1])))
```

## Reaction model and rates

`ReactionSource` holds the literal reactant/product complexes; `CoinfectionModel` composes that source with `Parameters`, exact symbolic derivatives and a numerical flow. The reaction order is recruitment; the three alpha infection channels; two eta superinfection channels; two gamma superinfection channels; two beta single-strain transmissions from a coinfected individual; and four removals. Inspect `source.reactions` for the full fourteen rows.

The field is

```text
ds = Lambda - s*(mu0 + alpha1*a + alpha2*b + (alpha3+beta1+beta2)*c)
da = a*(alpha1*s - gamma1*b - eta1*c - mu1) + beta1*s*c
db = b*(alpha2*s - gamma2*a - eta2*c - mu2) + beta2*s*c
dc = (gamma1+gamma2)*a*b + c*(eta1*a + eta2*b + alpha3*s - mu3)
```

Every contact term cancels from the total-population balance: `d(s+a+b+c)=Lambda-mu0*s-mu1*a-mu2*b-mu3*c`. Recruitment and removal make this an open balance, not a conservation law. `bounds` returns the paper's eventual total bound `R=Lambda/min(mu)+1` and susceptible floor `ell=Lambda/(2*Q)`, with Q from the stated loss estimate.

If densities have units D and time has units T, recruitment has units D/T, removal rates 1/T, and contact coefficients 1/(D*T). Removal can aggregate mortality and permanent removal; recovery, immunity and return to susceptibility are absent. Editing the rate vector explores this model; adding compartments or reactions requires rederiving the dynamical theorem.

Numerical integration uses logarithmic concentrations, an analytic transformed Jacobian and a stiff solver. This retains positivity without clipping small missing compartments. Solver failures are explicit. It is still floating-point integration, not a validated trajectory enclosure. The initial states must be strictly positive; exact boundary states are handled by the structural and resident interfaces.

## Shared modes at a common boundary

A siphon is a species set whose simultaneous absence is invariant. Exhaustive source analysis gives exactly the empty set, `{a,c}`, `{b,c}`, and `{a,b,c}`. The shared set `{c}` is not itself a siphon: the two private strains can produce coinfection.

`overlap_certificate` checks both input sets against the source, requires a common nonnegative state at which their union is absent, and derives the normal Jacobian there. A nonzero entry from column j to row i can only decrease siphon membership: `membership(i)` is a subset of `membership(j)`. Grouping by increasing membership cardinality gives upper triangular blocks. Distinct signatures of equal cardinality do not couple. The code retains the factor for the empty signature when supplied normal indices include an uncovered missing compartment.

For the two fully covered siphons, the exact identity is

`chi_union * chi_intersection = chi_S1 * chi_S2`.

At `(s,0,0,0)` the three modes have diagonal values `alpha1*s-mu1`, `alpha2*s-mu2`, `alpha3*s-mu3`. The last is shared and is counted once in the union. The implementation multiplies polynomials and verifies a zero residual. It repeats the check at `s=mu3/alpha3`, where the shared eigenvalue is zero. No determinant is divided out, so a threshold does not break the calculation. Matrices at different resident states cannot be combined by this identity.

## Resident-dependent invasion and exact scope decisions

The disease-free resident is `E0=(Lambda/mu0,0,0,0)`. Single-strain residents have `si=mui/alphai`, `ui=(Lambda-mu0*si)/mui`; a resident exists only if ui>0. `invasion_blocks` differentiates the full source field at E1 and E2 and extracts the missing `(b,c)` and `(a,c)` blocks respectively.

`InvasionBlock` uses the exact two-dimensional Metzler criterion. For `M=[[A,B],[C,D]]` with B,C>0, strict invasion holds when `A>=0` or `D>=0` or `A*D<B*C`. Strict decay holds when both diagonal entries are negative and `B*C<A*D`. The remaining threshold has dominant eigenvalue zero. The code constructs and checks rational **left** covectors `(1,r)` and their componentwise relative margins. Numerical eigenvalues are supplementary readouts.

With all rates positive and both residents present, strict mutual invasion invokes the paper's uniform permanence theorem. A strictly decaying missing block invokes its converse: that resident is a sink attracting some positive states, so uniform permanence fails. A threshold or absent resident is reported as undecided/outside these hypotheses, not rounded into a success/failure classification. Permanence does not mean global convergence to an interior equilibrium, and the time needed to reach a floor can depend on the initial state.

For the fixed published beta2 family, all disease-free reproduction numbers remain `(8,4,4)` while invasion of E1 switches at `beta2=51/140`. At beta2=0.1 its dominant missing-block eigenvalue is negative; at beta2=1 it is about 0.187386. E2 remains strictly invadable. Changing beta2 changes off-diagonal disease-free Jacobian entries, so that sweep alone does not keep the full Jacobian fixed.

The stronger saved example holds beta2=0.1 and changes eta1 from 0.1 to 0.2. This leaves the **entire disease-free Jacobian** and the resident coordinates unchanged, yet changes invasion at E1 from decay to growth. Superinfection is quadratic in missing compartments at E0 but becomes linear in the missing strain when the other strain is resident.

## A continuous independent parameter box

`ParameterBox` uses exact rational interval arithmetic on all fourteen independently varying rates. It encloses the resident coordinates and expands the weighted matrix columns before bounding, preserving cancellation such as the gamma1 term in `(1,1)*M2`. It checks fixed positive covectors `(1,2)` at E1 and `(1,1)` at E2, as well as resident existence throughout the box.

At the reference vector and radius 1%, the reproduced relative margins are exactly `5601/101000` and `8701/25250`. This is a certificate for the entire continuous box, not a corner sample. Other centers/radii can be explored; a failed enclosure is inconclusive because these fixed weights and interval bounds are only sufficient. The conclusion is permanence at each parameter vector, with its own existential floor. A single common floor across the box is not supplied by this paper or this program.

## The normalized growth and recovery components

`NormalizedLift` constructs resident weights `ru,rv`, a compensating disease-free weight rj and integer k. The masses are `U=a+ru*c`, `V=b+rv*c`, `J=a+b+rj*c`, with product `P=U*V*J**k`. The lift carries `(a/U,b/V,a/J,b/J)` as four independent coordinates and exposes the polynomial eight-dimensional field and continuous relative growth functions `(HU,HV,HJ,G)`. It verifies positive lower bounds for G over the equilibrium fibers. Tests include a case where `HU+HV` is negative in the pure-coinfected disease-free direction and the factor `J**k` is essential.

These objects make the proof's mechanism available for further modelling, but they do not compute its compactness-derived common growth window or product floor q. Arbitrary extended boundary directions need not lie in the physical compact closure K. Evaluating the polynomial field at such a point does not establish physical realizability or invariance of an arbitrary direction box.

`conditional_recovery(q)` implements the exact algebra that converts an **externally justified** eventual product floor into weighted-mass floors and then individual compartment floors. It labels q as an assumption and its outputs as conditional. Supplying a simulated minimum does not justify that assumption. The default run never calls this routine with a guessed floor.

`ResidentEntropy` constructs the corrected relative entropy on either single-strain face and its exact derivative formula. The additional epsilon times the product of the two composition deviations makes the eventual drift negative definite. The saved face example uses unequal losses `mu0=6/5`, `mu1=13/10`. This explains resident convergence without assuming a floor for the infected coordinate in advance.

## Validation

Seven test groups check the literal source and total balance, exact resident stationarity, every siphon, overlap at a zero shared mode, an added uncovered normal mode, covector signs and threshold scope, the fourteen-dimensional interval certificate, quotient-field and product-growth identities, the disease-free compensation case, the entropy and recovery algebra, and numerical trajectories against a second stiff solver and the exact total-population solution when all losses equal one.

The figures compare finite trajectories from several initial states and a near-sink extinction example. Their minima are observations, not universal lower bounds. The general permanence and overlap results remain the paper's theorems; this package rechecks their finite hypotheses and algebra and supplies reusable modelling components. It does not rerun Lean.

Licensing: MIT is proposed, pending the owner's decision. No new license grant is made by this example.