mantis-delta — Mathematical background & worked calculations

This document is the mathematical companion to the README. It explains what mantis-delta computes and how, with worked numerical examples you can reproduce line-by-line. Every matrix, steady state, and eigenvalue printed here is generated by the library itself (the exact code snippets are shown alongside).

Contents

Notation and the reaction-network object
The stoichiometry matrix N
Mass-action kinetics: from reactions to ODEs
Worked example I — A ⇌ B end to end
Conservation laws as the left null space of N
Worked example II — Michaelis–Menten
Deficiency and the Feinberg theorems
Worked example III — the CHA cascade (δ = 1)
Chemostatting: folding fixed concentrations
Worked example IV — Brusselator & the Hopf bifurcation
How the steady-state solver works
Stability classification from the Jacobian spectrum
Stochastic kinetics: propensities and the τ-leap
References

1. Notation and the reaction-network object

A chemical reaction network is a triple of finite sets (𝒮, 𝒞, ℛ):

Species 𝒮 = {X₁, …, X_d} — the molecular entities whose concentrations evolve. In code these are rn.species, sorted alphabetically, and d = rn.n_species.
Complexes 𝒞 — the multisets of species standing on either side of a reaction arrow. The complex 2X + Y is the formal object {(X,2),(Y,1)}; the empty complex ∅ is allowed. n = rn.n_complexes is the number of distinct complexes.
Reactions ℛ ⊆ 𝒞 × 𝒞 — directed arrows yᵢ → yⱼ between complexes. A reversible string A <-> B is stored as two directed reactions. r = rn.n_reactions.

Each complex y is a vector in concentration space ℝ^d: y = Σ νₛ eₛ where νₛ is the stoichiometric coefficient of species s in that complex. A reaction yᵢ → yⱼ has an associated reaction vector (yⱼ − yᵢ) ∈ ℝ^d — the net molecular change it produces.

Why complexes are the fundamental node type. mantis-delta builds a directed graph whose nodes are complexes and whose edges are reactions (rn.draw() renders exactly this graph). Deficiency, linkage classes, and weak reversibility are all properties of this graph — not of the species or of the reactions individually. This is the central modelling choice of Chemical Reaction Network Theory (CRNT).

2. The stoichiometry matrix N

The stoichiometry matrix N collects every reaction vector as a column:

N[i, j]  =  (coeff of species i in products of reaction j)
          − (coeff of species i in reactants of reaction j)

so N has shape (d species × r reactions) — rn.stoichiometry_matrix. This is the opposite of some textbooks that use the transpose; mantis-delta keeps species on the rows so the mass-action ODE reads as a clean matrix–vector product (next section).

The stoichiometric subspace is S = span(columns of N) = image of N. Its dimension s = rank(N) is the number of independent directions the state can move in. Everything orthogonal to S is conserved — that is the whole content of §5.

build_stoichiometry_matrix (in stoichiometry.py) is a direct transcription of the formula above: it walks each reaction, subtracts reactant coefficients and adds product coefficients into the relevant row. Species not in the dynamic list (e.g. chemostatted ones, §9) simply have no row.

3. Mass-action kinetics: from reactions to ODEs

Under the law of mass action, the rate (flux) of reaction j is its rate constant times the product of reactant concentrations raised to their stoichiometric coefficients:

v_j(x)  =  k_j · Π_s  x_s^(ν_{s,j})

where ν_{s,j} is the coefficient of species s in the reactant complex of reaction j. Stacking the fluxes into a vector v(x) ∈ ℝ^r, the full dynamics are

dx/dt  =  N · v(x)

This single equation is what odes() assembles symbolically (via SymPy) and what simulate() / the steady-state solver evaluate numerically. Concretely, build_odes computes, for each species i,

d x_i / dt  =  Σ_j  N[i, j] · v_j(x)

and sympy.expands the result. With numeric_rates=False the rate constants stay as symbols k_1, k_2, … (numbered in reaction order); with numeric_rates=True (the default) the supplied float values are substituted in.

The Jacobian is J(x) = ∂(N v)/∂x. Symbolically it is just sympy.diff of each ODE w.r.t. each species (jacobian()); numerically the solver uses the closed-form product rule in _full_jacobian_fn, accumulating J[:, k] += N[:, j] · ∂v_j/∂x_k reaction by reaction. Stability is read off the eigenvalues of J evaluated at a steady state (§12).

4. Worked example I — A ⇌ B end to end

The smallest non-trivial network. One reversible reaction → two directed reactions.

from mantis import CRNetwork
rn = CRNetwork.from_string(["A <-> B"], rates={"A -> B": 1.0, "B -> A": 0.5})

Reaction order and rate symbols. rn.rate_keys() → ['A -> B', 'B -> A'], so k₁ = 1.0 (forward) and k₂ = 0.5 (back).

Stoichiometry matrix. Species are ['A', 'B']; reactions are [A→B, B→A]:

N =  [[-1.  1.]      row A:  −1 (consumed in A→B),  +1 (produced in B→A)
      [ 1. -1.]]     row B:  +1 (produced in A→B),  −1 (consumed in B→A)

rank(N) = 1 (the two columns are negatives of each other), so s = 1.

Fluxes and ODEs. v₁ = k₁·A = A, v₂ = k₂·B = 0.5 B. Then

dA/dt = −v₁ + v₂ = −1.0·A + 0.5·B
dB/dt = +v₁ − v₂ = +1.0·A − 0.5·B

exactly what rn.odes() returns: {'A': -1.0*A + 0.5*B, 'B': 1.0*A - 0.5*B}.

Conservation law. Since dA/dt = −dB/dt, the total A + B is constant — this is the one conservation law (rn.conservation_laws → [A + B]), and it is why rank(N) = 1 < 2.

Steady state. Setting dA/dt = 0 gives A = (k₂/k₁)·B, i.e. A/B = 0.5. Combined with the conserved total A + B = A₀ + B₀ = 2 (from the initial condition {A:2, B:0}):

A* = k₂/(k₁+k₂) · (A₀+B₀) = 0.5/1.5 · 2 = 0.6667
B* = k₁/(k₁+k₂) · (A₀+B₀) = 1.0/1.5 · 2 = 1.3333

ss = rn.steady_states({"A": 2.0, "B": 0.0})[0]
ss.concentrations   # {'A': 0.6667, 'B': 1.3333}

Jacobian and stability. For this linear system the Jacobian is constant:

J = [[−k₁,  k₂],     symbolic:  Matrix([[-k_1, k_2], [k_1, -k_2]])
     [ k₁, −k₂]]     numeric :  Matrix([[-1.0, 0.5], [1.0, -0.5]])

Its eigenvalues are −(k₁+k₂) = −1.5 and 0. The zero eigenvalue is not a sign of instability — it is the conservation-law direction (motion along A+B = const is forbidden). mantis-delta filters it out (§12) and classifies the steady state as stable (ss.is_stable == True). This is the Deficiency Zero Theorem at work: δ = 2 − 1 − 1 = 0 and the network is weakly reversible, so a unique asymptotically stable steady state is guaranteed for any positive k₁, k₂.

5. Conservation laws as the left null space of N

A conservation law is a vector c ∈ ℝ^d such that the linear combination c·x is constant along every trajectory:

d(c·x)/dt = c·(N v(x)) = (cᵀN) v(x) = 0  for all x   ⟺   cᵀ N = 0

So conservation laws are exactly the left null space of N (equivalently the null space of Nᵀ). There are d − s independent ones.

For a mass-balanced chemical network these always admit a basis with non-negative integer entries — the physical "moiety totals" (a conserved atom group, a conserved backbone, etc.). conservation_laws_sympy computes this in three steps:

Exact null space. Build N as an integer SymPy matrix and call Nᵀ.nullspace() — exact rational arithmetic, no floating-point fuzz.
Clear denominators. _scale_to_integers multiplies each basis vector by the LCM of its denominators and divides out the GCD, giving the smallest integer vector.
Non-negative, elementary form. _make_nonneg_basis adds/subtracts integer multiples of the basis vectors to (a) eliminate negative entries and (b) minimise the number of non-zeros, so the reported laws are the "most elementary" moiety totals.

The result is returned as SymPy expressions, e.g. E + ES — read as "E + ES = const".

Conservation laws set the steady state. A closed network has one steady state per stoichiometry class — i.e. per choice of the conserved totals. Those totals are fixed by your initial conditions: starting from {E: 1e-6, ES: 0} pins E + ES = 1e-6 forever, and a different total yields a different steady state. This is why steady_states() takes initial conditions rather than just rate constants.

6. Worked example II — Michaelis–Menten

mm = CRNetwork.from_string(["E + S <-> ES", "ES -> E + P"])

Species ['E', 'ES', 'P', 'S']; rate keys ['E + S -> ES', 'ES -> E + S', 'ES -> E + P'] (3 directed reactions).

Stoichiometry matrix (rows E, ES, P, S; columns = the 3 reactions):

       E+S→ES   ES→E+S   ES→E+P
E   [   −1        +1       +1   ]
ES  [   +1        −1       −1   ]
P   [    0         0       +1   ]
S   [   −1        +1        0   ]

rank(N) = 2, and d = 4, so there are 4 − 2 = 2 conservation laws:

mm.conservation_laws    # [E + ES,  ES + P + S]

E + ES = const — enzyme is recycled, never created or destroyed.
ES + P + S = const — total substrate-derived material is conserved.

Structural verdict. n = 3 complexes (E+S, ES, E+P), l = 1 linkage class, s = 2 ⟹ δ = 3 − 1 − 2 = 0. But the network is not weakly reversible (the arrow ES → E+P has no return path), so the Deficiency Zero Theorem does not apply — a δ = 0 network still needs weak reversibility for the DZT guarantee. mantis-delta reports exactly this in crnt_summary(). (Michaelis–Menten of course still has a well-behaved unique steady state — DZT being inapplicable means the theorem is silent, not that the network misbehaves.)

7. Deficiency and the Feinberg theorems

The deficiency of a network is the non-negative integer

δ = n − l − s

n = number of distinct complexes (n_complexes)
l = number of linkage classes = weakly-connected components of the complex graph (n_linkage_classes)
s = rank(N) = dim of the stoichiometric subspace

compute_deficiency is literally n - l - rank_N. Deficiency measures how far the network's complex-space geometry is from its stoichiometric geometry — the slack between the way complexes are linked in the graph and the dimension actually spanned by the reaction vectors. It is always ≥ 0, and it is a property of the graph plus stoichiometry, independent of the rate constants.

Linkage classes. Two complexes are in the same linkage class if they are connected in the undirected version of the complex graph. get_linkage_classes returns networkx.weakly_connected_components.

Weak reversibility. Every weakly-connected component must also be strongly connected: if you can get from complex y to y′ following arrows, you must be able to get back (not necessarily by a single reaction). check_weak_reversibility tests is_strongly_connected on each weakly-connected subgraph.

Deficiency Zero Theorem (Horn–Jackson 1972; Feinberg 1979)

If δ = 0 and the network is weakly reversible, then for every choice of positive rate constants the mass-action system has exactly one positive steady state in each stoichiometry class, that steady state is asymptotically stable, and there are no periodic orbits. No bistability, no oscillations — structurally ruled out.

deficiency_zero_applies = (δ == 0 and weakly_reversible).

Deficiency One Theorem (Feinberg 1995)

mantis-delta checks the structural hypotheses in check_deficiency_one_theorem:

δ = 1, and
each linkage class has its own deficiency δ_ℓ ∈ {0, 1}, and
at most one linkage class has δ_ℓ = 1.

The per-linkage-class deficiency δ_ℓ = (n_ℓ − 1 − s_ℓ) is computed by _per_lc_deficiency, restricting N to the reactions and species of that linkage class and taking the rank. When the hypotheses hold:

the mass-action system admits at most one steady state per stoichiometry class for any positive rate constants — so bistability is structurally impossible, even though stability and oscillations are not constrained by this theorem alone.

deficiency_one_applies is the boolean returned by that check.

8. Worked example III — the CHA cascade (δ = 1)

The Catalytic Hairpin Assembly biosensor — the library's primary validation case:

cha = CRNetwork.from_string([
    "miR21 + H1 <-> miR21_H1",          # toehold binding
    "miR21_H1 + H2 <-> H1H2 + miR21",   # strand displacement, catalyst release
    "H1H2 + CP <-> H1H2_CP",            # capture-probe / reporter formation
    "H1 + H2 <-> H1H2",                 # leakage dimerisation
])

Counting for the deficiency. The eight distinct complexes are

miR21+H1,  miR21_H1,  miR21_H1+H2,  H1H2+miR21,  H1H2+CP,  H1H2_CP,  H1+H2,  H1H2

so n = 8. All eight are distinct complex nodes — even though H1H2 recurs as a sub-term, the complexes H1H2 + miR21, H1H2 + CP, and H1H2 are three different nodes — so no two reactions share a node and the four reversible reactions form l = 4 independent linkage classes (one two-complex component each). The stoichiometric subspace has rank s = 3. Therefore

δ = n − l − s = 8 − 4 − 3 = 1

The network is weakly reversible (every reaction is written <->). The structural check finds δ = 1 with admissible per-linkage-class deficiencies, so the Deficiency One Theorem applies:

cha.n_complexes              # 8
cha.n_linkage_classes        # 4
cha.deficiency               # 1
cha.is_weakly_reversible     # True
cha._crnt_result.deficiency_one_applies   # True

Four conservation laws (cha.conservation_laws):

CP + H1H2_CP                          # capture probe is conserved
H1H2 + H1H2_CP + H2                   # H2 moiety total
miR21 + miR21_H1                      # the miR-21 catalyst is conserved
H1 + H1H2 + H1H2_CP + miR21_H1        # H1 moiety total

The miR-21 conservation law miR21 + miR21_H1 = const is the mathematical statement that the trigger acts catalytically — it is released by the displacement step and never consumed. The D1T guarantee (at most one steady state per stoichiometry class) is what makes the dose–response curve of examples/04_cha_cascade.py single-valued and monotone.

9. Chemostatting: folding fixed concentrations

An open network holds some species at fixed concentrations — fed by a reservoir, a large buffer, or a flow reactor. mantis-delta treats a chemostatted species not as a dynamic variable but as a parameter folded into the rate constants.

If species A is chemostatted at concentration c_A and appears as a reactant with coefficient ν in reaction j, then its contribution c_A^ν is constant, so it is absorbed into an effective rate constant:

k_j^eff  =  k_j · Π_{A ∈ chemostatted}  c_A^(ν_{A,j})

This is exactly fold_chemostatted_into_rates. The consequences ripple through the whole library consistently:

Species list / N: chemostatted species get no variable and no stoichiometry row — the dynamic subsystem is genuinely smaller.
Deficiency, conservation laws: computed on the reduced complex graph (chemostatted species are stripped from each complex via reduce_complex before the graph is built).
ODEs: only the dynamic species get equations; the folded k^eff carry the reservoir concentrations as numeric prefactors.
Steady states: with chemostatting an open network typically has no conservation laws, which changes the solver strategy entirely (§11) and is what lets it find unstable fixed points.

A chemostatted concentration of 0.0 makes the corresponding reaction a pure sink (zero flux) — the standard way the Brusselator's product-removal steps are written.

10. Worked example IV — Brusselator & the Hopf bifurcation

br = CRNetwork.from_string(
    ["A -> X", "2X + Y -> 3X", "B + X -> Y + D", "X -> E"],
    rates={"A -> X": 1.0, "2X + Y -> 3X": 1.0, "B + X -> Y + D": 1.0, "X -> E": 1.0},
    chemostatted={"A": 1.0, "B": 3.0, "D": 0.0, "E": 0.0},
)

With A, B, D, E chemostatted, only X, Y are dynamic. Folding the fixed concentrations into the rates (k=1 throughout, A=1, B=3) gives the classic reduced system:

dX/dt = 1 − (B+1)·X + X²Y  =  1.0 − 4.0·X + X²·Y
dY/dt =      B·X − X²Y      =        3.0·X − X²·Y

returned verbatim by br.odes().

Fixed point. Setting both to zero: adding the equations gives 1 − X = 0 ⟹ X = 1, then Y = B/X = 3. So (X, Y) = (1, 3)* — independent of where you start.

Jacobian and the Hopf condition. The symbolic Jacobian (with k's symbolic) is

J = [[2XY·k₂ − 3k₃ − k₄ ,  X²·k₂ ],
     [ −2XY·k₂ + 3k₃     ,  −X²·k₂]]

Evaluated at (1, 3) with all k = 1, this is [[2, 1], [−3, −1]], whose eigenvalues are

λ = 0.5 ± 0.866 j

ss = br.steady_states({"X": 1.0, "Y": 1.0})[0]
ss.concentrations    # {'X': 1.0, 'Y': 3.0}
ss.eigenvalues       # [0.5+0.866j, 0.5-0.866j]
ss.is_stable         # False
ss.is_oscillatory    # True   (complex eigenvalues)

Positive real part + non-zero imaginary part = an unstable focus: trajectories spiral outward and settle onto a limit cycle. The general Brusselator Hopf condition is B > 1 + A² (here 3 > 1 + 1² = 2 ✓). At B = 1.5 the same construction gives eigenvalues with negative real part — a stable spiral, is_stable == True, is_oscillatory == True.

Why an algebraic solver matters here. Forward ODE integration can never land on an unstable fixed point — it spirals away to the limit cycle. The fixed point (1, 3) is found only because the solver solves N v(x) = 0 algebraically (§11). This is the headline capability of the open-network path: locating Hopf centres that simulation cannot reach.

11. How the steady-state solver works

A steady state solves N v(x) = 0 with x ≥ 0. The strategy (find_steady_states) branches on whether the network has conservation laws.

Closed / semi-open networks (have conservation laws)

The steady state is unique only within a stoichiometry class, so the conserved totals must be pinned. The solver:

ODE integration (_integrate_to_ss, SciPy stiff Radau) from the supplied IC. Integration automatically stays on the conservation manifold and relaxes to the stable attractor. The horizon t_end defaults to 1e4 s — deliberately short so oscillatory systems don't hang; raise it for very slow (e.g. leakage-dominated) kinetics.
Algebraic least_squares from the same IC. This can converge onto unstable fixed points that integration diverges away from. The algebraic system (_build_augmented_system) stacks:
rank(N) independent ODE rows — chosen by QR factorisation with column pivoting of Nᵀ so the rows are genuinely independent (no redundant equations), and
one equation per conservation law, c·x − total = 0, pinning the stoichiometry class.

That makes a square, well-posed root-find solved with bound-constrained (x ≥ 0) trust-region least squares. 3. Multi-start: random initial points generated on the constraint manifold (_make_feasible_initial rescales a random vector to match each conserved total), then ODE-integrated, falling back to algebraic solve. This is how multiple steady states (in genuinely multistable networks) are discovered.

Open / chemostatted networks (no conservation laws)

ODE integration is skipped entirely: a limit-cycle attractor would make it orbit forever. The solver uses only algebraic least_squares, first from the supplied IC, then from scale-aware random restarts (sampled around the magnitude of the initial conditions). This is the path that finds the Brusselator's unstable focus (§10).

Filtering the results

Every candidate is post-processed:

Non-negativity: any x_i < −tol is rejected.
Deduplication: a candidate within 10 % relative L2 distance of an existing one is dropped.
Residual ranking: solutions are sorted by the true ODE residual ‖N v(x)‖₂, and candidates whose residual is more than 1e4× worse than the best are discarded. Using a relative* residual threshold (not a fixed absolute one) is essential for stiff systems where all fluxes are tiny in absolute terms — it correctly separates a genuinely converged state from one "stuck at the initial condition".

The surviving SteadyState objects carry .concentrations, .eigenvalues, .is_stable, .is_oscillatory, and .residual (rule of thumb: < 1e-6 good, < 1e-10 excellent).

12. Stability classification from the Jacobian spectrum

For each steady state the solver forms the numeric Jacobian J(x*) (_full_jacobian_fn, the analytic product-rule derivative of N v) and takes its eigenvalues with scipy.linalg.eigvals. classify_steady_state then decides:

Filtering conservation-law zeros. A closed network's Jacobian has structurally-zero eigenvalues, one per conservation law (motion off the stoichiometry class is forbidden). These would muddy the stability test, so eigenvalues with

|λ| < 1e-4 · max_k |λ_k|

are treated as zero and dropped before classification. The remaining "significant" eigenvalues decide stability.

Condition on significant eigenvalues	Classification	`is_stable`	`is_oscillatory`
all Re(λ) < 0, all Im = 0	stable node	✓	✗
all Re(λ) < 0, some Im ≠ 0	stable spiral	✓	✓
some Re(λ) > 0, some Im ≠ 0	unstable focus (Hopf)	✗	✓
some Re(λ) > 0, all Im = 0	unstable node	✗	✗
some Re(λ) = 0	non-hyperbolic (borderline)	✗	—

is_stable = all significant eigenvalues have Re(λ) < 0.
is_oscillatory = any significant eigenvalue has |Im λ| > 1e-10·|λ|. Note this is true for both stable spirals (Re < 0, damped oscillation) and unstable foci (Re > 0, the Hopf case) — it flags "the system rotates in phase space near this point", which is the right signal for identifying Hopf-bifurcation candidates.

13. Stochastic kinetics: propensities and the τ-leap

When molecule counts are small (single-cell volumes, ≲10³ molecules), the deterministic ODE is the wrong model and mantis-delta offers two exact/approximate stochastic simulators that work in molecule counts, not concentrations.

From rate constant to stochastic propensity

A deterministic rate constant k (with concentration units) must be converted to a stochastic rate constant c (per-molecule). For a reaction of total reactant order m, with reaction volume V and Avogadro's number N_A,

c  =  k / (N_A · V)^(m − 1)

(first-order reactions, m = 1, are unchanged). The propensity — the probability per unit time that the reaction fires given the current counts — is a(x) = c · h(x), where h(x) is the number of distinct reactant combinations:

h(x)  =  Π_s  C(n_s, ν_s)  =  Π_s  n_s·(n_s−1)···(n_s−ν_s+1) / ν_s!

This combinatorial h (not the mass-action power n_s^{ν_s}) is what makes the SSA exact. Both simulators compute it directly in their propensity helper.

Exact SSA — `stochastic_simulate`

Gillespie's direct method: at each step, with total propensity a₀ = Σ aⱼ, draw the time to the next event as an exponential, τ = −ln(r₁)/a₀, and select which reaction fires with probability aⱼ/a₀ (linear search on a second uniform r₂). Apply the reaction's integer change vector, advance time, repeat. Exact for any population size, but slow when reactions fire frequently.

Approximate τ-leap — `tau_leap_simulate`

Instead of one reaction per step, fire many at once over a leap τ by sampling each reaction's firing count from Poisson(aⱼ · τ). Roughly N× faster than direct SSA, where N is the mean number of firings per leap.

Adaptive τ (default, Cao–Gillespie–Petzold 2006): τ is chosen to bound the fractional change of every propensity per leap by epsilon (default 0.03), using the per-species drift μ_i = Σ_j ν_{ji} a_j and variance σ²_i = Σ_j ν_{ji}² a_j. Pass tau=… to override with a fixed leap.
Leap rejection: if a tentative leap would drive any species count negative, that leap is discarded and replaced by a single exact-SSA step. This keeps counts non-negative without the heavier binomial Tian/Burrage machinery — adequate for kinetic-design work.

Both return a StochasticResult with .times, .counts, .concentrations, plus .at(t) (step-function lookup) and .final(). As populations grow, a τ-leap trajectory's mean converges to the deterministic ODE solution — a sanity check exercised in the test suite.

14. References

Horn F, Jackson R (1972). General mass action kinetics. Arch. Rational Mech. Anal. 47, 81–116. (Deficiency Zero Theorem.)
Feinberg M (1979). Lectures on Chemical Reaction Networks. University of Wisconsin. Available at crnt.osu.edu.
Feinberg M (1995). The existence and uniqueness of steady states for a class of chemical reaction networks. Arch. Rational Mech. Anal. 132, 311–370. (Deficiency One Theorem.)
Feinberg M (2019). Foundations of Chemical Reaction Network Theory. Springer.
Gillespie D T (1977). Exact stochastic simulation of coupled chemical reactions. J. Phys. Chem. 81, 2340–2361. (Direct-method SSA.)
Cao Y, Gillespie D T, Petzold L R (2006). Efficient step size selection for the tau-leaping simulation method. J. Chem. Phys. 124, 044109. (Adaptive τ.)

Every numerical value in this document is produced by mantis-delta itself; reproduce them by pasting the snippets into a Python session with pip install -e . from the package root. See the README for the user-facing API and the examples/ directory for complete runnable scripts.