Skip to content

ODE systems

Your ODE model in 3 steps

1. Install.

pip install "twindex[general,report]"

2. Describe the model and the data. model.json holds plain polynomial expressions. Here is a two-compartment pharmacokinetic model:

{
  "states": ["central", "peripheral"],
  "params": ["k10", "k12", "k21"],
  "rhs": [
    "k21*peripheral - (k10 + k12)*central",
    "k12*central - k21*peripheral"
  ],
  "initial": ["1", "0"],
  "observables": {"central": "central"}
}

data.csv has the columns t,observable,y. Decimals are read exactly, so 0.1 means 1/10, and rational strings such as 1/8 also work. Noise is an absolute bound shared by all observations. To get both example files without cloning the repository, run:

mkdir ode_quickstart && cd ode_quickstart
curl -LO https://raw.githubusercontent.com/Arnav-M10/Twindex/main/examples/ode_quickstart/model.json
curl -LO https://raw.githubusercontent.com/Arnav-M10/Twindex/main/examples/ode_quickstart/data.csv
cd ..

3. Run.

twindex ode ode_quickstart/model.json ode_quickstart/data.csv \
  --eps 1/100 --box k10=0,1 --box k12=0,1 --box k21=0,1 \
  --target k10 --output out/ --starts 1 --max-iter 20

From a clone of the repository, use examples/ode_quickstart/ as the path instead.

The checked-in CSV observes the central compartment of a two-compartment pharmacokinetic model over the first unit of time. --box is required exactly once for each named parameter; order does not matter. An optional --estimate is a comma-separated vector of exact values in model parameter order.

The output contains certificate.json, replay.json, report.md, observations.png, numerical-audit.json when available, and result.json. Replay independently with twindex check out/certificate.json. If integration fails before any candidate is generated, the output contains a result and explanatory report, without a certificate. Reusing an output directory replaces its previous workflow artifacts.

The same workflow is available in Python:

from twindex import identify_ode
from twindex.identify import read_ode_model, read_ode_data

result = identify_ode(
    read_ode_model('examples/ode_quickstart/model.json'),
    read_ode_data('examples/ode_quickstart/data.csv'),
    epsilon='1/100', box=[('0', '1')] * 3, target='k10',
    output='out', starts=1, max_iter=20,
)
print(result.status, result.lower)

result.lower is an exact checked pair gap, or None. Status is certified_pair, witness_only, incompatible, or inconclusive. Incompatible refers only to the attempted witness, never the whole prior. A supplied estimate is tried first; numerical best-fit and twin proposals are then checked independently. numerical_audit is explicitly numerical evidence, not proof. The CLI returns 0 for VALID (including witness-only), 1 for INVALID or malformed input, and 2 for INCONCLUSIVE. A failed search, small lower bound, or a single certified witness proves no identifiability or global upper bound.

Twindex certifies parameter ambiguity for any polynomial ODE system

\[ u'(t) = F(t, u; \theta), \qquad u(0) = u_0(\theta), \]

with observations that are affine combinations of the states. Polynomial vector fields cover a wide range of models:

  • mass-action chemical kinetics;
  • compartment models in epidemiology (SIR, SEIR) and pharmacokinetics;
  • Lotka–Volterra ecology;
  • FitzHugh–Nagumo neurons;
  • the Lorenz system.

Signs are unrestricted, and an unknown initial value is simply another parameter.

Derived targets such as R₀

Often the quantity that matters is a combination of parameters. A pair claim can name any rational expression, and Twindex certifies the gap in that quantity directly:

$ twindex check examples/study/ode/sir_r0_pair.json
examples/study/ode/sir_r0_pair.json: VALID: twins certified; beta/gamma differs by 3245/4178 across data-compatible parameters

In Python, use check_pair_ode(problem, w1, w2, gap, target_expr="beta/gamma") or make_certificate_ode(..., target_expr="beta/gamma"). The expression is evaluated exactly, and if it divides by zero the certificate is INVALID.

Beyond polynomials

Right-hand sides may use division, negative powers, exp, sin and cos:

from twindex.ode import ODEModel

enzyme = ODEModel(["S"], ["Vmax", "Km"], ["-Vmax*S/(Km + S)"], ["10"])      # Michaelis–Menten
reactor = ODEModel(["T"], ["k", "E"], ["-k*exp(-E/T)*(T - 1)"], ["2"])      # Arrhenius
gene = ODEModel(["p"], ["b", "K"], ["b/(1 + (p/K)**2) - p"], ["0"])         # Hill repression
forced = ODEModel(["x", "v"], ["w", "a"], ["v", "-w**2*x + a*sin(t)"], ["1", "0"])
print(enzyme.states)   # ('S', '__inv1'): the checker adds z = 1/(Km + S) exactly

Twindex rewrites each model as an exactly equivalent polynomial system by adding auxiliary states, and proves that the rewrite is faithful (THEORY Theorem 12). Every tool, including twindex ode, reports and the browser checker, then works unchanged. The solver picks state scales automatically, because auxiliary states such as \(1/(K_m + S)\) live on very different scales from the original ones.

Ready-made templates with units, priors, schedules and stated assumptions live in examples/ode_templates/. Each is tested end to end, from synthetic data to a certificate that replays VALID.

template field non-polynomial terms
sir epidemiology none
pk2 pharmacokinetics none
enzyme biochemistry (mass action) none
michaelis_menten enzyme kinetics, the classic \(V_{max}/K_m\) ridge division
hill_repression gene regulation division, square
predator_prey_holling ecology (Rosenzweig–MacArthur) division
seasonal_sir seasonal epidemics sine forcing

Early-epidemic data cannot determine R₀

from fractions import Fraction as F

from twindex.checker import Verdict
from twindex.ode import ODEModel, ODEProblem, check_pair_ode, make_certificate_ode
from twindex.solver_ode import candidate_ode, fit_witness_ode

sir = ODEModel(states=["S", "I"], params=["beta", "gamma"],
               rhs=["-beta*S*I", "beta*S*I - gamma*I"],
               initial=["0.99", "0.01"], observables={"I": "I"})

# Synthetic data: infections at t = 0.5, 1, ..., 4 from beta = 0.5, gamma = 0.2.
truth = [F(1, 2), F(1, 5)]
p = candidate_ode(sir, truth, horizon=4, segments=4)
times = [F(k, 2) for k in range(1, 9)]
problem = ODEProblem(sir, tuple((t, "I", p(t)[1]) for t in times),
                     epsilon=F(1, 1000), box=((0, 2), (0, 1)), target=0)

# A twin with the same early growth rate beta*S0 - gamma but a different R0.
twin = [F(18, 25), F(18, 25) * F(99, 100) - (truth[0] * F(99, 100) - truth[1])]
q, _ = fit_witness_ode(problem, twin, segments=4)

report = check_pair_ode(problem, (truth, p), (twin, q), F(11, 50))
assert report.verdict is Verdict.VALID
print(report.reasons[0])   # beta differs by 11/50 across data-compatible parameters
certificate = make_certificate_ode(problem, [(truth, p), (twin, q)], F(11, 50))

The two parameter sets have reproduction numbers \(R_0 = \beta/\gamma = 2.5\) and about \(1.72\). Both reproduce every measurement within \(10^{-3}\), so no method can recover \(R_0\) from these data. The certificate proves this as an exact statement under the bounded-noise model; it does not rely on a profile likelihood or a linearisation. The saved certificate, examples/study/ode/sir_pair.json (9 KB), replays with twindex check and needs no numpy.

How the proof works

The candidate is a piecewise polynomial with a tube of radius \(\delta\) around it. On each segment the checker bounds, in exact rational arithmetic:

  • the residual \(\rho \ge \sup|p' - F(t,p)|\);
  • the logarithmic norm \(\mu\) of \(\partial F/\partial u\) over the tube.

From these it derives the error bound

\[ |u - p| \le \frac{E}{1 - \mu s} + \frac{\rho s}{1 - \mu^+ s}. \]

If this bound stays inside the tube, the true solution exists and obeys it (THEORY Theorem 11). For states with very different units, pass scales, which switches to the weighted norm \(\max_i |e_i|/w_i\).

Prior art covers related ground. Guaranteed bounded-error estimation with validated integration and set inversion (T. Raïssi, N. Ramdani and Y. Candau, Automatica, 2004) and structural identifiability tools (DAISY, SIAN, StructuralIdentifiability.jl) both address this area. Twindex's contribution is a portable, exact certificate of ambiguity together with a small, independent checker that anyone can replay.

Suggesting an additional measurement time

from fractions import Fraction as Q
from twindex.ode import parse_certificate_ode
from twindex.certificate import read_certificate_json
from twindex.measurement import rank_measurement_times

_, problem, pair, _ = parse_certificate_ode(
    read_certificate_json('examples/study/ode/sir_pair.json'))
ranked = rank_measurement_times(problem, pair, [Q(k, 10) for k in range(41)])

This is a heuristic ranking of time/observable pairs by the absolute candidate-polynomial difference divided by epsilon. It replays both witnesses first and accepts exact times only within both candidate horizons; it does not extrapolate or optimize over the compatible set. By default all model observables are ranked; observables=['I'] restricts the list. Rows contain exact Fraction time, score, difference, the two error_bounds, and separates_twins.

The optional exact statement is deliberately narrow. If separates_twins is true, the difference exceeds 2*epsilon + error1 + error2 strictly, using exact observable bounds from fresh VALID witness replays. The two closed intervals of possible noisy observations are disjoint. Every future datum at that time, for that observable and the same absolute noise bound, is incompatible with at least one of these two fixed parameter vectors. Equality does not suffice. There is no conclusion about diameters, uniqueness, other parameters, or the best experimental design. A high heuristic score alone proves no separation.