End-to-end system identification: a symbolic ODE model

This notebook identifies the parameters of an epidemic model written as a system of differential equations with SymPy, and asks the same question its user-defined companion asks: once the model fits, can the parameters actually be recovered?

Here the answer depends on when you stopped collecting data, which makes it a concrete demonstration that identifiability is a property of the experiment and not only of the model.

A MATLAB Live Script version (Examples/sir_symbolic.m) runs the same computation with the same section numbering. Throughout, # MATLAB: comments give the equivalent toolbox call.

The model

The classical SIR system for a closed population of $N=1000$:

$$\frac{dS}{dt}=-\beta\frac{SI}{N},\qquad \frac{dI}{dt}=\beta\frac{SI}{N}-\gamma I,\qquad \frac{dR}{dt}=\gamma I$$

Two factors are to be identified: $\beta$ (transmission rate, 1/day) and $\gamma$ (recovery rate, 1/day). Their ratio is the basic reproduction number $R_0=\beta/\gamma$, the quantity an epidemiologist actually wants.

Only the infected compartment $I(t)$ is observed, and the initial conditions are known: one infected individual in an otherwise susceptible population.

The three expressions below are the whole model.

import warnings
warnings.filterwarnings("ignore")

import numpy as np
import pandas as pd
import sympy as sp
import matplotlib.pyplot as plt

from gsua_csb import (SymbolicODEModel, design_matrix, uncertainty_analysis,
                      coverage_metric, parameter_estimation,
                      identifiability_analysis, profile_likelihood)

N = 1000.0                                  # closed population
t = sp.Symbol("t")                          # time symbol
S, I, R = sp.symbols("S I R")               # state variables
beta, gamma = sp.symbols("beta gamma")      # parameters to identify

# One right-hand side per state, in the same order as the state list below.
# MATLAB: odes = [diff(S) == -beta*S*I/N; diff(I) == beta*S*I/N - gamma*I; diff(R) == gamma*I]
odes = [-beta * S * I / N,                  # dS/dt  susceptible
         beta * S * I / N - gamma * I,      # dI/dt  infected   (this is what we observe)
         gamma * I]                         # dR/dt  recovered

odes

[-0.001*I*S*beta, 0.001*I*S*beta - I*gamma, I*gamma]

1. Preparing the environment

For a symbolic model the factor vector is the initial conditions of the states first, in the order the states are listed, followed by the model parameters — the same convention MATLAB uses. Giving a factor a degenerate range fixes it, so the five rows below leave only $\beta$ and $\gamma$ free.

The model integrates all three states, but only $I(t)$ is measured. output=1 selects it once, at construction, and every routine that follows honours it — the counterpart of MATLAB's 'output',2 argument to gsua_dataprep, or of editing T.Properties.CustomProperties.output afterwards. It can be changed at any time with model.set_output(...).

NAMES  = ["S0", "I0", "R0i", "beta", "gamma"]
truth  = np.array([999.0, 1.0, 0.0, 0.35, 0.10])   # S0, I0, R0, then beta, gamma
ranges = np.array([[999, 999],      # S0    susceptible at t=0  (known -> degenerate)
                   [1,   1],        # I0    infected at t=0     (known -> degenerate)
                   [0,   0],        # R0    recovered at t=0    (known -> degenerate)
                   [0.15, 0.9],     # beta  to estimate
                   [0.02, 0.4]])    # gamma to estimate
FREE = [3, 4]                                       # indices of beta and gamma


def build_model(domain):
    """Build the SIR model over a given time domain, observing I(t) only.

    MATLAB: [T,~] = gsua_dataprep(odes, vars, [0 tmax], 'sirEpidemicModel',
                                  'range', ranges, 'output', 2)
    """
    return SymbolicODEModel(odes, [S, I, R], t, [beta, gamma], domain=domain,
                            names=NAMES, range=ranges, nominal=truth,
                            output=1)                # <-- observe state 1, I(t)


model_full = build_model(np.linspace(1, 80, 20))
print("free factors:", [NAMES[i] for i in FREE], " | true R0 =", truth[3] / truth[4])
print("evaluate returns", model_full.evaluate(truth).shape, "-- one row, because output=1")

free factors: ['beta', 'gamma']  | true R0 = 3.4999999999999996
evaluate returns (1, 20) -- one row, because output=1

2. Synthetic data

The true epidemic uses $\beta=0.35$ and $\gamma=0.10$, so $R_0=3.5$. Case counts are noisier when they are larger, so the measurements carry Poisson-like noise whose spread grows with the count.

xfull = np.linspace(1, 80, 20)              # 80 days of surveillance, 20 reports

rng = np.random.default_rng(3)              # fix the noise draw so the page reproduces
# MATLAB: cleanFull = gsua_eval(truth, T, xfull, [], false, false);
clean_full = model_full.evaluate(truth, xfull)              # already just I(t)
yfull = clean_full + np.sqrt(np.maximum(clean_full, 1)) * rng.standard_normal(clean_full.shape)

tdense = np.linspace(1, 80, 300)            # dense grid, for drawing curves only

plt.figure(figsize=(6.4, 3.8))
plt.plot(tdense, model_full.evaluate(truth, tdense)[0], lw=1.6, label="true epidemic")
plt.plot(xfull, yfull[0], "ko", ms=5, label="surveillance data")
plt.axvline(25, ls="--", c="grey")
plt.text(26, 300, "end of early phase", fontsize=9, color="grey")
plt.xlabel("time (days)"); plt.ylabel("infected individuals I(t)")
plt.title("Simulated epidemic, $R_0$ = 3.5")
plt.legend(); plt.grid(alpha=.3); plt.tight_layout(); plt.show()

3. Can these parameters be estimated at all?

Before spending any optimizer budget it is worth asking whether the data is even reachable: does it fall inside the range of epidemics the model can produce over the $\beta$ and $\gamma$ bounds declared in section 1? If it does not, no amount of optimization will help. This is the reachability check that opens the toolbox's semi-automated identification cycle.

# MATLAB: M0 = gsua_dmatrix(T,300);  Y0 = gsua_ua(M0,T,'xdata',xfull,'ynom',yfull);
M0 = design_matrix(model_full, 300, seed=0)     # 300 samples of the (beta, gamma) box
ua0 = uncertainty_analysis(model_full, M0, xdata=xfull, y_exp=yfull)

cost_data, cost_band, p5, _, p95 = coverage_metric(ua0.Y, yfull, margin=0.1)
contained = float(np.mean((yfull >= p5) & (yfull <= p95)))

pd.DataFrame({"contained": [contained],
              "median_band_width": [float(np.median(p95 - p5))],
              "cost_data": [cost_data],
              "cost_band": [cost_band]}).round(4)

	contained	median_band_width	cost_data	cost_band
0	0.9	313.9139	12219.2235	1.440196e+17

plt.figure(figsize=(6.4, 3.8))
plt.fill_between(xfull, p5[0], p95[0], color="0.85")
plt.plot(xfull, p5[0], color="0.4", lw=1, label="5th percentile")
plt.plot(xfull, p95[0], color="0.4", lw=1, label="95th percentile")
plt.plot(xfull, yfull[0], "ko", ms=5, label="surveillance data")
plt.xlabel("time (days)"); plt.ylabel("infected individuals I(t)")
plt.title(f"Reachable band before fitting ({contained:.0%} of data contained)")
plt.legend(); plt.grid(alpha=.3); plt.tight_layout(); plt.show()

Nearly all the observations fall inside the reachable band, so an epidemic of this shape is within the model's declared range and estimation is worth attempting. The band itself spans much of the population — far too wide to be an answer — which is the uncertainty the next sections set out to reduce.

4. Estimating from the full epidemic

With the whole curve in hand — growth, peak and decline — both factors are estimated by multistart least squares, twenty restarts.

# MATLAB: [Tfull,resFull] = gsua_pe(T, xfull, yfull, 'solver','lsqc','N',20,'margin',0.1)
pe_full = parameter_estimation(model_full, xfull, yfull, n=20,
                               solver="least_squares", margin=0.1, seed=0)
best_full = pe_full.x[np.argmin(pe_full.cost)]
ia_full = identifiability_analysis(model_full, pe_full.x, cost=pe_full.cost,
                                   cost_rtol=0.1, seed=0)
# margin = assumed relative standard deviation. MATLAB offsets this by one: margin=1.1 there.
pl_full = profile_likelihood(model_full, xfull, yfull, alpha=0.95, margin=0.1, params=FREE)

pd.DataFrame({
    "true":      truth[FREE],
    "estimated": best_full[FREE].round(4),
    "CI_low":    pl_full.range[FREE, 0].round(4),
    "CI_high":   pl_full.range[FREE, 1].round(4),
    "width":     (pl_full.range[FREE, 1] - pl_full.range[FREE, 0]).round(4),
}, index=["beta", "gamma"])

	true	estimated	CI_low	CI_high	width
beta	0.35	0.3442	0.3189	0.3556	0.0367
gamma	0.10	0.0996	0.0966	0.1429	0.0463

print(f"cost          : {pe_full.cost.min():.5g}")
print(f"R0 estimated  : {best_full[3] / best_full[4]:.3f}    (true 3.500)")

cost          : 3979.1
R0 estimated  : 3.454    (true 3.500)

5. Estimating from the early phase only

Now suppose the analysis had to be done during the outbreak, with only the first 25 days available — the situation every real-time epidemic assessment faces. Nothing about the model changes; only the observation window shrinks.

xearly = np.linspace(1, 25, 12)             # 25 days, 12 reports
model_early = build_model(xearly)           # same system, shorter domain

rng = np.random.default_rng(3)
clean_early = model_early.evaluate(truth, xearly)
yearly = clean_early + np.sqrt(np.maximum(clean_early, 1)) * rng.standard_normal(clean_early.shape)

pe_early = parameter_estimation(model_early, xearly, yearly, n=20,
                                solver="least_squares", margin=0.1, seed=0)
best_early = pe_early.x[np.argmin(pe_early.cost)]
ia_early = identifiability_analysis(model_early, pe_early.x, cost=pe_early.cost,
                                    cost_rtol=0.1, seed=0)
pl_early = profile_likelihood(model_early, xearly, yearly, alpha=0.95, margin=0.1, params=FREE)

pd.DataFrame({
    "true":      truth[FREE],
    "estimated": best_early[FREE].round(4),
    "CI_low":    pl_early.range[FREE, 0].round(4),
    "CI_high":   pl_early.range[FREE, 1].round(4),
    "width":     (pl_early.range[FREE, 1] - pl_early.range[FREE, 0]).round(4),
}, index=["beta", "gamma"])

	true	estimated	CI_low	CI_high	width
beta	0.35	0.4062	0.2336	0.4845	0.2509
gamma	0.10	0.1461	0.0200	0.2437	0.2237

The estimates still look plausible. Note the cost, though: the early-phase fit is better than the full-epidemic fit was — fewer points, all of them on a smooth exponential rise — while the intervals have widened.

6. Diagnosing identifiability

The correlation between the repeated estimates explains what happened.

pd.DataFrame({
    "full_epidemic": [ia_full.correlation[3, 4]],
    "early_phase":   [ia_early.correlation[3, 4]],
}, index=["corr(beta, gamma)"]).round(4)

	full_epidemic	early_phase
corr(beta, gamma)	0.913	1.0

In the early phase it is essentially $1$. During exponential growth the data constrain only the growth rate, roughly $\beta-\gamma$, so any pair with the right difference reproduces the observations equally well. It takes the peak — where susceptibles are depleted and the curve turns over — to separate the two.

7. The two windows side by side

pd.DataFrame({
    "full epidemic": [pe_full.cost.min(), ia_full.correlation[3, 4],
                      pl_full.range[3, 1] - pl_full.range[3, 0],
                      pl_full.range[4, 1] - pl_full.range[4, 0],
                      best_full[3] / best_full[4]],
    "early phase":   [pe_early.cost.min(), ia_early.correlation[3, 4],
                      pl_early.range[3, 1] - pl_early.range[3, 0],
                      pl_early.range[4, 1] - pl_early.range[4, 0],
                      best_early[3] / best_early[4]],
}, index=["best cost", "corr(beta, gamma)", "CI width beta",
          "CI width gamma", "R0 (true 3.5)"]).round(4)

	full epidemic	early phase
best cost	3979.0835	1089.0670
corr(beta, gamma)	0.9130	1.0000
CI width beta	0.0367	0.2509
CI width gamma	0.0463	0.2237
R0 (true 3.5)	3.4545	2.7809

fig, ax = plt.subplots(1, 2, figsize=(11, 3.8))

x = np.arange(2); w = 0.35
ax[0].bar(x - w/2, [pl_full.range[i, 1] - pl_full.range[i, 0] for i in FREE], w,
          label="full epidemic")
ax[0].bar(x + w/2, [pl_early.range[i, 1] - pl_early.range[i, 0] for i in FREE], w,
          label="early phase")
ax[0].set_xticks(x); ax[0].set_xticklabels([r"$\beta$", r"$\gamma$"])
ax[0].set_ylabel("95% CI width"); ax[0].set_title("Confidence interval width")
ax[0].legend(); ax[0].grid(alpha=.3, axis="y")

ax[1].plot(tdense, model_full.evaluate(best_full, tdense)[0], lw=1.6, label="fit to full epidemic")
ax[1].plot(tdense, model_full.evaluate(best_early, tdense)[0], "--", lw=1.6,
           label="fit to early phase, extrapolated")
ax[1].plot(xfull, yfull[0], "ko", ms=4, label="data")
ax[1].set_xlabel("time (days)"); ax[1].set_ylabel("infected individuals I(t)")
ax[1].set_title("Where the early-phase fit leads")
ax[1].legend(fontsize=8); ax[1].grid(alpha=.3)
plt.tight_layout(); plt.show()

8. What this example shows

The same model and the same estimator produced two very different states of knowledge, and the difference was the observation window rather than anything about the algorithm. The better-fitting dataset was the less informative one.

Fit quality measures agreement with the points you have; identifiability measures whether those points could have distinguished your parameters from the alternatives. They are different questions, and only the second tells you whether an estimate is worth reporting.

For a real outbreak the consequence follows directly: an $R_0$ estimated before the peak carries a confidence interval wide enough to change policy conclusions, and quoting the point estimate alone would hide that.

The companion user-defined example reaches the same conclusion from the opposite direction: there, adding information by fixing a known factor makes the fit slightly worse and the parameters recoverable.

Where a multistart run does spread across the factor space, identifiability_analysis adds detection of multiple global minima, noise_floor calibrates which fits to accept against the observation noise, and design_matrix(..., method="joint") propagates the accepted set without destroying its correlation structure.