Example Code
A steady state is a set of concentrations at which every species is produced and removed at equal overall rates. The paper proves that, for a minimal Type II autocatalytic core—a branching cycle with no smaller autocatalytic part—under its positive-rate and positive-degradation assumptions, there cannot be two different positive steady states. The examples below show why a larger network with similar reaction connections need not have that property.
Download the complete Python program. Its full source appears at the end of this page. The calculations and figures below were regenerated from that program; you do not need to run it to inspect them.
Two different concentration states, each in exact balance
Within each panel, the two colors refer to the same network and the same fixed reaction and degradation constants. Each species has zero concentration derivative at both states. The seven-species panel uses a logarithmic vertical scale, with exact fractions printed above the orange bars.
| Example | State y | State x | Exact balance check | |---|---|---|---| | Six species, Theorem 8.1 | (1,1,1,1,1,1) | (31/50,57/50,22/25,3/5,14/25,3/5) | All six derivatives are zero at each state. | | Seven species, Theorem 8.3 | (1,1,1,1,1,1,1) | (1,1/16,1/9,1/5,1/3,4/3,12) | All seven derivatives are zero at each state. |
The program builds the product and stoichiometric matrices from the reaction descriptions, computes every forward and reverse rate, and subtracts each degradation loss. It uses exact fractions: a reported zero is the rational number zero, not a small floating-point error. Both examples therefore have at least two positive steady states. These calculations do not establish that either state is stable, or that these are the only positive steady states.
To see the balance directly, a reaction of the form X₀ ⇌ 2X₁ + X₅ has net current
Its contribution is −J₀ to X₀, +2J₀ to X₁, and +J₀ to X₅. The program adds the contributions of all reactions and subtracts dᵢzᵢ for each species. All rate constants remain fixed when it evaluates the second state.
Read every exact species balance, or download all exact parameters and results. Fractions in the JSON file are stored as strings to preserve their full precision.
Why these larger networks are not minimal cores
In the six-species example, retaining X₁ and X₂ gives the forward reactions X₁ → 2X₂ and X₂ → X₁. Using them in amounts 2 and 3 gives a surplus of 1 of each retained species. In the seven-species example, retaining X₂, X₃ and X₄ gives X₂ → X₃, X₃ → X₄ + X₂, and X₄ → X₃. Amounts 2, 3 and 2 give a surplus of 1 of each retained species.
These are smaller autocatalytic parts under the paper's allowed restriction: omitted species are treated as externally maintained. The stated reaction amounts establish structural net production. They are not claimed to be the currents at either of the original network's steady states.
A minimal core can have one positive steady state, or none
This phase diagram concerns a deliberately symmetric example within the paper's family. Its reactions are X₀ ⇌ X₁+X₅, X₁ ⇌ X₂, X₂ ⇌ X₃+X₁, X₃ ⇌ X₄, X₄ ⇌ X₅+X₃, and X₅ ⇌ X₀, together with degradation of every species.
The three fork reactions have forward constant a. The other three have forward constant b=4. All reverse constants are 1, and every degradation constant is d. These are chosen model units, rather than experimental measurements.
Theorem 1.1 guarantees at most one positive steady state. Because rotating this network by two species leaves all reactions and constants unchanged, its unique positive state, when present, must have one common fork concentration f and one common intermediate concentration t. Their exact formulas are
There is one positive steady state precisely when 4a−d(a+5)−d²>0. There is none when this expression is zero or negative. This is a consequence derived here for this family using the paper's theorem, rather than a phase diagram reproduced from the paper.
| Forward constant a | Degradation constant d | Result | |---:|---:|---| | 3 | 1 | All six concentrations are 1. | | 5 | 1 | All six concentrations are 3. | | 3 | 3/2 | No positive steady state. | | 2 | 1 | On the boundary; no positive steady state. |
The all-zero state is not counted because it is not positive. The region labeled “no positive steady state” makes no assertion about where trajectories go over time. The region counts follow from the exact formula and the uniqueness theorem; the 5,800 checked parameter pairs are additional reproducible examples.
Exact parameter grid · Coordinates of the drawn boundary
Inspecting concentrations between the two known states
Here the horizontal coordinate t selects the concentration vector (1−t)y+tx. It is not time. At t=0 and t=1, every balance is exactly zero. The intermediate samples generally have nonzero balances, even though the two endpoints balance. Each panel uses its own vertical scale; all constants in that network have been divided by their largest value, a common rescaling that preserves stationary concentrations.
The displayed lines connect 101 exact rational samples per network. They do not prove that the full network has exactly two steady states or that a trajectory follows this segment. Download these samples.
What else the program checks
The teaching core's default state has currents (2,3,2,3,2,3), producing exactly the degradation loss (1,1,1,1,1,1). The program examines all 664 proper species/reaction restrictions of its unit-weight forward source, checking that none preserves the relevant strongly connected branching structure. This finite check supplies the absence of a smaller autocatalytic core needed to apply the paper's theorem to this example.
It also works through one exact instance of the paper's current-difference calculation. After constructing the matrix K from a base steady state and a chosen candidate second state, it obtains det(K)=1011717/1250>0. The corresponding reduced three-fork determinant is 37471/10000>0. Finally, it checks a unit return-chain elimination with two internal species. The detailed fractions and matrices are retained in results.json.
These worked calculations illustrate parts of the general proof. The script does not reproduce the entire theorem for arbitrary network sizes, and it does not compile Lean. The supplied paper reports Lean verification of the main theorem and the six-species certificate, and exact rational verification of the seven-species certificate.
Run the program yourself
With Python, NumPy, and Matplotlib installed:
python type_ii_demo.py
To install the plotting dependencies:
python -m pip install numpy matplotlib
To run just the exact calculations using the Python standard library:
python type_ii_demo.py --no-figures
The default destination is an outputs directory next to the script. Use --output-dir PATH to choose another location. Re-running replaces files with the same output names. All inputs are included in the source; there is no random seed, account, external service, or private data requirement.
The saved run used Python 3.11.0, NumPy 2.4.1 and Matplotlib 3.10.5; exact versions and the source SHA-256 are recorded in results.json. Both the figure-generating and standard-library-only modes passed.
Python source
#!/usr/bin/env python3
"""Exact examples for Type II core unistationarity (paper dated 7 September 2026).
Run: python type_ii_demo.py
Optional plotting dependencies: numpy, matplotlib. For arithmetic only, use
python type_ii_demo.py --no-figures
All certificates use fractions.Fraction, without a numerical tolerance.
The general uniqueness theorem is supplied by the paper, not proved by sampling.
The phase map is a derived consequence for the symmetric family defined here.
"""
import argparse
import csv
from dataclasses import dataclass
from fractions import Fraction as Q
import hashlib
import itertools
import json
import math
from pathlib import Path
import platform
def require(condition, message):
"""Keep checks active even if Python is run with -O."""
if not condition:
raise AssertionError(message)
def matvec(a, x):
return [sum((v * w for v, w in zip(row, x)), Q(0)) for row in a]
def matmul(a, b):
return [[sum((a[i][k] * b[k][j] for k in range(len(b))), Q(0))
for j in range(len(b[0]))] for i in range(len(a))]
def solve(a, b):
"""Gaussian elimination with exact fractions; singular input raises an error."""
n = len(a)
m = [[Q(v) for v in row] + [Q(b[i])] for i, row in enumerate(a)]
for c in range(n):
pivot = next((r for r in range(c, n) if m[r][c]), None)
require(pivot is not None, "Singular exact linear system")
m[c], m[pivot] = m[pivot], m[c]
p = m[c][c]
m[c] = [v / p for v in m[c]]
for r in range(n):
if r != c:
p = m[r][c]
m[r] = [v - p * w for v, w in zip(m[r], m[c])]
return [row[-1] for row in m]
def determinant(a):
m = [[Q(v) for v in row] for row in a]
result = Q(1)
for c in range(len(m)):
pivot = next((r for r in range(c, len(m)) if m[r][c]), None)
if pivot is None:
return Q(0)
if pivot != c:
m[c], m[pivot] = m[pivot], m[c]
result = -result
p = m[c][c]
result *= p
for r in range(c + 1, len(m)):
t = m[r][c] / p
for j in range(c + 1, len(m)):
m[r][j] -= t * m[c][j]
return result
def inverse(a):
columns = [solve(a, [int(i == j) for i in range(len(a))])
for j in range(len(a))]
return [list(row) for row in zip(*columns)]
@dataclass
class Network:
name: str
weights: tuple
backs: dict
forward: list
reverse: list
degradation: list
def __post_init__(self):
self.n = len(self.weights)
require(self.n >= 3, "At least three species required")
require(all(isinstance(w, int) and w >= 1 for w in self.weights),
"Product coefficients must be positive integers")
require(all(0 <= r < self.n and 0 <= b < self.n
and b not in (r, (r + 1) % self.n)
for r, b in self.backs.items()), "Invalid back target")
for key in ("forward", "reverse", "degradation"):
values = [Q(v) for v in getattr(self, key)]
require(len(values) == self.n and all(v > 0 for v in values),
f"{key} must contain one positive value per species/reaction")
setattr(self, key, values)
# Column r records the full product side of reaction r.
self.P = [[0] * self.n for _ in range(self.n)]
for r, w in enumerate(self.weights):
self.P[(r + 1) % self.n][r] += w
if r in self.backs:
self.P[self.backs[r]][r] += 1
self.N = [[self.P[i][r] - int(i == r) for r in range(self.n)]
for i in range(self.n)]
def products(self, z):
return [math.prod(Q(z[i]) ** self.P[i][r] for i in range(self.n))
for r in range(self.n)]
def balance(self, z):
require(len(z) == self.n and all(v > 0 for v in z),
"This check expects a strictly positive concentration vector")
products = self.products(z)
current = [self.forward[r] * z[r] - self.reverse[r] * products[r]
for r in range(self.n)]
production = matvec(self.N, current)
loss = [d * v for d, v in zip(self.degradation, z)]
return current, production, loss, [p - d for p, d in zip(production, loss)]
def six_species_certificate():
"""Table 1 and equations (8.2)-(8.4), transcribed as integers/fractions."""
net = Network("Six-species example from Theorem 8.1", (2, 2, 1, 3, 1, 3),
{0: 5, 2: 1, 4: 3},
[704180941994460, 619799798320080, 1064148026967810,
1379889114510795, 5001118197194658, 1928841702241720],
[18981846318750, 1952234297034000, 3747998870714400,
18981846318750, 937378238937273, 98424792637500],
[4806051633136950, 18981846318750, 18981846318750,
18981846318750, 18981846318750, 2918522144328875])
x = [Q(31, 50), Q(57, 50), Q(22, 25), Q(3, 5), Q(14, 25), Q(3, 5)]
require(net.products(x) == [Q(9747, 12500), Q(484, 625), Q(171, 250),
Q(2744, 15625), Q(9, 25), Q(29791, 125000)],
"Table 1 product monomials disagree")
return net, x
def seven_species_certificate():
"""Table 2: note that the final two concentrations are 4/3 and 12."""
D = 3546790543
net = Network("Seven-species example from Theorem 8.3", (1,) * 7,
{0: 6, 2: 1, 3: 2, 4: 3, 5: 4, 6: 5},
[Q(79873083, D), Q(4870030864, 11 * D), Q(6049869075, 11 * D),
Q(562292665, 2 * D), Q(560635312, D), Q(76005216, D),
Q(80869 * 43, D)],
[Q(150591284, D), Q(1809648, D), Q(401938000, 11 * D),
Q(1809648, D), Q(486297815, D), Q(1809648, D), Q(1809648, D)],
[Q(72385920, D), Q(1809648, D), Q(413617671, 2 * D),
Q(6785880025, 22 * D), Q(558389511, 2 * D),
Q(1809648, D), Q(1809648, D)])
require(sum(net.forward + net.reverse + net.degradation) == 1,
"Table 2 normalization disagrees")
require(min(net.forward + net.reverse + net.degradation) == Q(1809648, D),
"Table 2 smallest constant disagrees")
return net, [Q(1), Q(1, 16), Q(1, 9), Q(1, 5), Q(1, 3), Q(4, 3), Q(12)]
def symmetric_core(a=Q(3), b=Q(4), d=Q(1)):
"""A teaching example: three identical forks, each separated by one species."""
return Network("Symmetric minimal six-species core", (1,) * 6,
{0: 5, 2: 1, 4: 3}, [a, b] * 3, [1] * 6, [d] * 6)
def symmetric_state(a, b, d):
"""Derived in the plan. Uniqueness + cyclic symmetry makes this exhaustive."""
a, b, d = Q(a), Q(b), Q(d)
require(min(a, b, d) > 0, "Positive kinetic and degradation constants required")
h = a * b - d * (a + b + 1) - d * d
if h <= 0:
return None
tail = h / (1 + 2 * d)
fork = (b - d) * tail / (1 + 2 * d)
require(min(fork, tail) > 0, "Positive branch formula failed")
return [fork, tail] * 3
def strongly_connected(net, species, reactions):
"""Reachability in the forward species graph, using only retained products."""
adjacency = {i: set() for i in species}
for r in reactions:
if r in species:
adjacency[r].update(i for i in species if net.P[i][r])
for start in species:
seen, todo = {start}, [start]
while todo:
for j in adjacency[todo.pop()] - seen:
seen.add(j)
todo.append(j)
if seen != set(species):
return False
return True
def check_small_core_minimality(net):
"""Exhaust all retained species/reaction sets for this small unit-weight example.
A productive strongly connected component in any restriction is itself one
of these retained sets. Thus checking connected sets with a two-product fork
exhausts (Top) witnesses. Reverse currents belong to the kinetic model;
the classified diluted source graph uses the displayed forward reactions.
This is not a general minimal-core recognizer for other source conventions.
"""
require(all(w == 1 for w in net.weights), "This finite check uses unit weights")
full = set(range(net.n))
require(strongly_connected(net, full, full), "Full source is not connected")
require(any(sum(net.P[i][r] for i in full) >= 2 for r in full), "No fork")
checked = 0
for k in range(1, net.n + 1):
for species in itertools.combinations(range(net.n), k):
# Reactions with omitted sources are external feeds, not autonomous
# reactions of the retained diluted source, so do not retain them.
for mask in range(1, 1 << k):
reactions = [species[j] for j in range(k) if mask >> j & 1]
if set(species) == full and set(reactions) == full:
continue
checked += 1
if not any(sum(net.P[i][r] for i in species) >= 2 for r in reactions):
continue
require(not strongly_connected(net, species, reactions),
f"Proper (Top) restriction found: {species}, {reactions}")
return checked
def restriction_witness(net, species, reactions, flux):
s = [[net.N[i][r] for r in reactions] for i in species]
gain = matvec(s, flux)
require(all(v > 0 for v in gain), "Restriction is not strictly productive")
return {"species": species, "reactions": reactions, "matrix": s,
"flux": flux, "production": gain}
def current_difference_check(net, rho):
"""Equations (5.1)-(5.7), on the teaching core at reference state y = 1.
rho is a candidate second concentration vector. This checks the reduction
at that candidate; it does not turn finite tests into a universal proof.
"""
n = net.n
p, q, e = net.forward, net.reverse, net.degradation
require(not any(net.balance([Q(1)] * n)[3]), "Reference state not stationary")
C = [[Q(0)] * n for _ in range(n)]
for r, w in enumerate(net.weights):
successor = (r + 1) % n
secant = sum((rho[successor] ** j for j in range(w)), Q(0))
if r in net.backs:
back = net.backs[r]
C[r][back] = 1
C[r][successor] = rho[back] * secant
else:
C[r][successor] = secant
offset = [v - 1 for v in rho]
require(matvec(C, offset) == [v - 1 for v in net.products(rho)],
"Exact product difference identity failed")
A = [[p[r] * int(r == i) - q[r] * C[r][i] for i in range(n)] for r in range(n)]
scaled_N = [[Q(net.N[i][r]) / e[i] for r in range(n)] for i in range(n)]
AN = matmul(A, scaled_N)
K = [[int(r == s) - AN[r][s] for s in range(n)] for r in range(n)]
base_current = [u - v for u, v in zip(p, q)]
defect = matvec(K, base_current)
require(defect == [q[r] * (sum(C[r]) - 1) for r in range(n)],
"Stationary-current comparison identity failed")
forks, gaps = sorted(net.backs), [i for i in range(n) if i not in net.backs]
block = lambda rows, cols: [[K[i][j] for j in cols] for i in rows]
H = block(gaps, gaps)
H_inverse = inverse(H)
require(all(v >= 0 for row in H_inverse for v in row), "Gap inverse has wrong sign")
correction = matmul(matmul(block(forks, gaps), H_inverse), block(gaps, forks))
R = [[K[i][j] - correction[r][s] for s, j in enumerate(forks)]
for r, i in enumerate(forks)]
m = len(forks)
margins = [R[j][j] - R[j][(j + 1) % m] for j in range(m)]
require(all(R[j][(j - 1) % m] <= 0 and R[j][(j + 1) % m] > 0
and margins[j] >= 1 for j in range(m)), "Fork inequalities failed")
require(all(v > 0 for v in matvec(R, [base_current[i] for i in forks])),
"Positive fork comparison failed")
det_K, det_R = determinant(K), determinant(R)
require(det_K == determinant(H) * det_R and det_K > 0, "Determinant check failed")
delta = matvec(A, offset)
residual = net.balance(rho)[3]
# If rho were stationary, residual would vanish and K delta would be zero.
require(matvec(K, delta) == [-v for v in matvec(A, [r / d for r, d in zip(residual, e)])],
"Current-difference equation failed")
require(any(residual), "Intended nonstationary candidate unexpectedly stationary")
return {"candidate": rho, "C": C, "K": K, "gap_block": H, "fork_matrix": R,
"diagonal_margins": margins, "determinant_K": det_K,
"determinant_fork_matrix": det_R, "comparison": defect,
"candidate_balance": residual}
def passive_chain_check():
"""Lemma 4.1 in a two-internal-species example; coefficients are exact."""
H = [[5, -1], [-3, 7]]
left = solve(H, [2, 0])
right = solve(H, [0, 2])
c, beta = 4 * left[1], right[0]
leak_left, leak_right = 2 - left[0] - c, 2 - 4 * right[1] - beta
require(c > 0 and beta > 0 and min(leak_left, leak_right) >= 0, "Tail is not passive")
X, Y = Q(2), Q(1)
internal = solve(H, [2 * X, 2 * Y])
j_left, j_right = 2 * X - internal[0], 4 * internal[1] - 2 * Y
require(j_left == c * X - beta * Y + leak_left * X, "Left endpoint formula failed")
require(j_right == c * X - beta * Y - leak_right * Y, "Right endpoint formula failed")
require(j_left - j_right == internal[0] + 2 * internal[1], "Internal loss mismatch")
return {"c": c, "beta": beta, "left_loss_coefficient": leak_left,
"right_loss_coefficient": leak_right, "internal_state": internal,
"entering_current": j_left, "leaving_current": j_right}
def save_csv(path, header, rows):
with path.open("w", newline="", encoding="utf-8") as f:
writer = csv.writer(f)
writer.writerow(header)
writer.writerows(rows)
def vector_text(values):
return "(" + ", ".join(str(v) for v in values) + ")"
def make_figures(out, certificates, boundary):
import numpy as np
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
plt.rcParams.update({"font.family": "DejaVu Sans", "font.size": 11,
"axes.spines.top": False, "axes.spines.right": False,
"axes.edgecolor": "#a7b2b5", "axes.labelcolor": "#23363c",
"text.color": "#23363c", "xtick.color": "#44565b",
"ytick.color": "#44565b", "savefig.facecolor": "white"})
teal, orange = "#087e8b", "#cb6632"
versions = {"numpy": np.__version__, "matplotlib": matplotlib.__version__}
fig, axes = plt.subplots(1, 2, figsize=(12, 4.8), layout="constrained")
for ax, (net, x) in zip(axes, certificates):
idx = np.arange(net.n)
ax.bar(idx - .19, np.ones(net.n), .34, color=teal, label="State y")
ax.bar(idx + .19, [float(v) for v in x], .34, color=orange, label="State x")
if net.n == 7:
ax.set_yscale("log", base=2)
ax.set_ylim(1 / 32, 32)
ax.set_yticks([1 / 16, 1 / 4, 1, 4, 16], ["1/16", "1/4", "1", "4", "16"])
else:
ax.set_ylim(0, 1.42)
for i, value in enumerate(x):
ax.annotate(str(value), (i + .19, float(value)), xytext=(0, 5),
textcoords="offset points", ha="center", fontsize=9)
ax.set_xticks(idx, [f"$X_{i}$" for i in idx])
ax.set_ylabel("Concentration / concentration in state y")
if net.n == 7:
ax.set_ylabel("Concentration / concentration in state y\n(logarithmic scale)")
ax.set_title(f"{net.n} species · two exactly verified steady states", loc="left", fontsize=12, pad=14)
ax.grid(axis="y", alpha=.14)
ax.set_axisbelow(True)
ax.legend(frameon=False, loc="upper left" if net.n == 7 else "upper right")
fig.suptitle("The same rate constants can support two different concentration states", fontsize=16)
fig.savefig(out / "two_stationary_states.png", dpi=190)
fig.savefig(out / "two_stationary_states.svg")
plt.close(fig)
fig, ax = plt.subplots(figsize=(9, 6.1), layout="constrained")
a = np.array([row[0] for row in boundary])
critical = np.array([row[1] for row in boundary])
ax.fill_between(a, .05, critical, color="#d6eeee")
ax.fill_between(a, critical, 2.5, color="#f1ece3")
ax.plot(a, critical, color=teal, lw=2.3, label="Exact boundary formula")
ax.scatter([3], [1], c=teal, edgecolors="white", s=100, zorder=4)
ax.annotate("Default: every concentration = 1", (3, 1), xytext=(3.45, .75),
fontsize=10, arrowprops={"arrowstyle": "-", "color": teal})
ax.scatter([3], [1.5], marker="x", c=orange, s=75, zorder=4)
ax.annotate("Comparison: no positive steady state", (3, 1.5), xytext=(.5, 1.8),
fontsize=10, arrowprops={"arrowstyle": "-", "color": orange})
ax.text(4.65, .32, "ONE positive\nsteady state", color=teal, fontsize=14, ha="center", weight="bold")
ax.text(2.6, 2.15, "NO positive steady state", color="#806249", fontsize=14, ha="center", weight="bold")
ax.set(xlim=(.25, 6), ylim=(.05, 2.5), xlabel="Forward rate constant at each fork, a",
ylabel="Degradation constant for every species, d")
ax.set_title("A minimal core: where does a positive steady state exist?\n"
"Six species · b = 4 · all reverse rate constants = 1", loc="left", pad=18)
ax.text(.02, .98, "Boundary: 4a − d(a + 5) − d² = 0", transform=ax.transAxes,
va="top", fontsize=11)
fig.savefig(out / "positive_state_phase_map.png", dpi=190)
fig.savefig(out / "positive_state_phase_map.svg")
plt.close(fig)
# The segment connects two exact states; it is not a time trajectory.
fig, axes = plt.subplots(1, 2, figsize=(12, 4.8), layout="constrained")
line_rows = []
for ax, (net, x), color in zip(axes, certificates, [teal, orange]):
scale = max(net.forward + net.reverse + net.degradation)
normalized = Network(net.name, net.weights, net.backs,
[p / scale for p in net.forward],
[q / scale for q in net.reverse],
[d / scale for d in net.degradation])
errors = []
for j in range(101):
t = Q(j, 100)
z = [1 + t * (v - 1) for v in x]
residual = normalized.balance(z)[3]
value = max(abs(v) for v in residual)
errors.append(float(value))
line_rows.append([net.n, str(t), str(value), float(value)])
ax.plot(np.linspace(0, 1, 101), errors, color=color, lw=2.4,
label=f"{net.n}-species example")
ax.scatter([0, 1], [0, 0], color=color, s=45, clip_on=False)
ax.set(xlabel="Position t on the segment from y to x",
ylabel="Largest absolute concentration derivative\n(after a common rescaling of all rates)")
ax.set_title(f"{net.n}-species example", loc="left", fontsize=12, pad=12)
ax.grid(alpha=.13)
fig.suptitle("Both endpoints balance exactly; intermediate concentrations need not balance", fontsize=15)
fig.savefig(out / "between_stationary_states.png", dpi=190)
fig.savefig(out / "between_stationary_states.svg")
plt.close(fig)
save_csv(out / "between_states.csv", ["species_count", "t_exact", "max_derivative_exact", "max_derivative_decimal"], line_rows)
return versions
def main():
parser = argparse.ArgumentParser(description=__doc__)
parser.add_argument("--output-dir", type=Path, default=Path(__file__).resolve().parent / "outputs")
parser.add_argument("--no-figures", action="store_true", help="Run exact checks without NumPy or Matplotlib")
args = parser.parse_args()
out = args.output_dir
out.mkdir(parents=True, exist_ok=True)
report, data, balances = [], {}, []
emit = report.append
emit("TYPE II CORES: EXACT CHECKS AND WORKED EXAMPLES")
emit("Paper: type_II_l_unistationarity_arxiv.pdf, 7 September 2026")
emit("All reported balance zeros use exact rational arithmetic.")
certificates = [six_species_certificate(), seven_species_certificate()]
for net, x in certificates:
states = {"y": [Q(1)] * net.n, "x": x}
require(x != states["y"] and all(v > 0 for v in x), "States must be positive and distinct")
emit(f"\n{net.name}")
emit(" x = " + vector_text(x))
emit(" All forward, reverse and degradation constants: strictly positive")
state_data = {}
for name, z in states.items():
current, production, loss, residual = net.balance(z)
require(residual == [0] * net.n, f"{net.name}, state {name}: nonzero balance")
emit(f" State {name}: exact balance residual = {vector_text(residual)}")
state_data[name] = {"concentrations": z, "current": current,
"net_reaction_production": production, "degradation_loss": loss,
"balance_residual": residual}
for i in range(net.n):
balances.append([net.n, name, i, z[i], production[i], loss[i], residual[i]])
data[f"certificate_{net.n}"] = {"name": net.name, "P": net.P, "N": net.N,
"forward": net.forward, "reverse": net.reverse,
"degradation": net.degradation, "states": state_data}
six, seven = certificates[0][0], certificates[1][0]
witnesses = [restriction_witness(six, [1, 2], [1, 2], [2, 3]),
restriction_witness(seven, [2, 3, 4], [2, 3, 4], [2, 3, 2])]
data["proper_autocatalytic_restrictions"] = witnesses
emit("\nWhy these examples are not minimal cores:")
for witness in witnesses:
emit(" Retain species " + ", ".join(f"X{i}" for i in witness['species'])
+ ": flux " + vector_text(witness['flux'])
+ " gives production " + vector_text(witness['production']))
# Lemma 3.2's formula at M=2 is (3,5), a different valid witness.
require(matvec(witnesses[0]["matrix"], [3, 5]) == [2, 1], "Lemma 3.2 witness failed")
data["lemma_3_2_length_two_witness"] = {"flux": [3, 5], "production": [2, 1]}
core = symmetric_core()
count = check_small_core_minimality(core)
require(core.balance([Q(1)] * 6)[3] == [0] * 6, "Teaching core default balance failed")
core_gain = matvec(core.N, [2, 3] * 3)
require(core_gain == [1] * 6, "Full core production witness failed")
data["teaching_core"] = {"P": core.P, "N": core.N, "forward": core.forward,
"reverse": core.reverse, "degradation": core.degradation,
"proper_restrictions_checked": count, "productive_flux": [2, 3] * 3,
"production": core_gain}
emit(f"\nTeaching core: {count} proper species/reaction restrictions checked; none retain (Top).")
emit(" At a=3, b=4, d=1, the state (1,1,1,1,1,1) balances exactly.")
emit(" Its uniqueness follows from Theorem 1.1 after checking the core assumptions.")
presets = [(3, 4, 1), (5, 4, 1), (3, 4, Q(3, 2)), (2, 4, 1)]
data["phase_presets"] = []
for a, b, d in presets:
z = symmetric_state(a, b, d)
if z is not None:
require(symmetric_core(a, b, d).balance(z)[3] == [0] * 6, "Symmetric formula failed")
data["phase_presets"].append({"a": a, "b": b, "d": d, "positive_state": z})
emit(f" a={a}, b={b}, d={d}: " + ("positive state " + vector_text(z) if z else "no positive steady state"))
kernel = current_difference_check(core, certificates[0][1])
data["current_difference_example"] = kernel
emit("\nExact current-difference calculation on the teaching core:")
emit(f" det(K) = {kernel['determinant_K']} > 0")
emit(f" det(reduced fork matrix) = {kernel['determinant_fork_matrix']} > 0")
emit(" This rejects the selected candidate; the paper proves nonsingularity generally.")
tail = passive_chain_check()
data["passive_chain_example"] = tail
emit(f"\nUnit-chain elimination: c={tail['c']}, beta={tail['beta']}, "
f"left loss={tail['left_loss_coefficient']}, right loss={tail['right_loss_coefficient']}.")
# Classify a rational grid exactly. Only drawing the boundary uses floats.
grid = []
for ia in range(25, 601, 5):
for id_ in range(5, 251, 5):
a, d = Q(ia, 100), Q(id_, 100)
z = symmetric_state(a, 4, d)
if z is not None:
require(symmetric_core(a, 4, d).balance(z)[3] == [0] * 6, "Grid balance failed")
grid.append([a, d, int(z is not None), z[0] if z else "", z[1] if z else ""])
from math import sqrt
boundary = []
for j in range(576):
a = (25 + j) / 100
s = a + 5
boundary.append([a, 8 * a / (sqrt(s * s + 16 * a) + s)])
save_csv(out / "exact_balances.csv", ["species_count", "state", "species_index", "concentration", "net_reaction_production", "degradation_loss", "balance_residual"], balances)
save_csv(out / "phase_grid.csv", ["a_exact", "d_exact", "positive_state_count", "fork_concentration_exact", "intermediate_concentration_exact"], grid)
save_csv(out / "phase_boundary.csv", ["a_decimal", "critical_d_decimal"], boundary)
versions = {} if args.no_figures else make_figures(out, certificates, boundary)
emit(f"\nPhase map: {len(grid)} parameter pairs classified exactly; every displayed positive state checked.")
emit(" Entire regions follow from the derived formula + the paper's uniqueness theorem.")
emit("\nScope: two counterexamples verified, not a count of all their steady states.")
emit("No stability, convergence, or thermodynamic-consistency conclusion is asserted.")
emit("The Lean sources were not supplied or compiled by this program.")
emit("ALL CHECKS PASSED")
text_report = "\n".join(report) + "\n"
(out / "run_output.txt").write_text(text_report, encoding="utf-8")
data["provenance"] = {"paper": "type_II_l_unistationarity_arxiv.pdf", "paper_date": "2026-09-07",
"python": platform.python_version(), "plotting_versions": versions,
"script_sha256": hashlib.sha256(Path(__file__).read_bytes()).hexdigest(),
"arithmetic": "exact rational for all checks; floating point for drawing only",
"phase_map_scope": "symmetric unit-weight six-species core; theorem + derived formula"}
(out / "results.json").write_text(json.dumps(data, indent=2, default=str) + "\n", encoding="utf-8")
print(text_report, end="")
if __name__ == "__main__":
main()
Run output
TYPE II CORES: EXACT CHECKS AND WORKED EXAMPLES Paper: type_II_l_unistationarity_arxiv.pdf, 7 September 2026 All reported balance zeros use exact rational arithmetic. Six-species example from Theorem 8.1 x = (31/50, 57/50, 22/25, 3/5, 14/25, 3/5) All forward, reverse and degradation constants: strictly positive State y: exact balance residual = (0, 0, 0, 0, 0, 0) State x: exact balance residual = (0, 0, 0, 0, 0, 0) Seven-species example from Theorem 8.3 x = (1, 1/16, 1/9, 1/5, 1/3, 4/3, 12) All forward, reverse and degradation constants: strictly positive State y: exact balance residual = (0, 0, 0, 0, 0, 0, 0) State x: exact balance residual = (0, 0, 0, 0, 0, 0, 0) Why these examples are not minimal cores: Retain species X1, X2: flux (2, 3) gives production (1, 1) Retain species X2, X3, X4: flux (2, 3, 2) gives production (1, 1, 1) Teaching core: 664 proper species/reaction restrictions checked; none retain (Top). At a=3, b=4, d=1, the state (1,1,1,1,1,1) balances exactly. Its uniqueness follows from Theorem 1.1 after checking the core assumptions. a=3, b=4, d=1: positive state (1, 1, 1, 1, 1, 1) a=5, b=4, d=1: positive state (3, 3, 3, 3, 3, 3) a=3, b=4, d=3/2: no positive steady state a=2, b=4, d=1: no positive steady state Exact current-difference calculation on the teaching core: det(K) = 1011717/1250 > 0 det(reduced fork matrix) = 37471/10000 > 0 This rejects the selected candidate; the paper proves nonsingularity generally. Unit-chain elimination: c=3/4, beta=1/16, left loss=13/16, right loss=11/16. Phase map: 5800 parameter pairs classified exactly; every displayed positive state checked. Entire regions follow from the derived formula + the paper's uniqueness theorem. Scope: two counterexamples verified, not a count of all their steady states. No stability, convergence, or thermodynamic-consistency conclusion is asserted. The Lean sources were not supplied or compiled by this program. ALL CHECKS PASSED