This skill covers the mathematical and computational workflow for building, fitting, and
analysing compartmental epidemic models in Python. Starting from classic SIR and SEIR
formulations, it extends to age-structured models, hospitalisation compartments (SEIHR),
real-time reproduction number estimation, and global sensitivity analysis via Partial Rank
Correlation Coefficients (PRCC).
All numerical integration is done with scipy.integrate.solve_ivp using the Radau solver,
which handles stiff ODE systems arising from wide parameter ranges common in outbreak modeling.
Setup
pip install numpy scipy pandas matplotlib numba
Core Functions
"""
epi_modeling.py
---------------
Compartmental epidemic model utilities: SIR, SEIR, SEIHR, Rt estimation,
parameter fitting, age-structured models, and sensitivity analysis.
"""from __future__ import annotations
import warnings
from typing importCallable, Dict, List, Optional, Tuple, Unionimport numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import matplotlib.dates as mdates
from scipy.integrate import solve_ivp
from scipy.optimize import minimize, differential_evolution
from scipy.stats import pearsonr, spearmanr
warnings.filterwarnings("ignore", category=RuntimeWarning)
# ---------------------------------------------------------------------------# 1. SEIR Model# ---------------------------------------------------------------------------def () -> :
() -> []:
S, E, I, R = y
force_of_infection = beta * S * I / N
dS = mu * (N - S) - force_of_infection
dE = force_of_infection - (sigma + mu) * E
dI = sigma * E - (gamma + mu) * I
dR = gamma * I - mu * R
[dS, dE, dI, dR]
seir_ode.__doc__ = (
)
seir_ode
() -> :
():
S, I, R = y
foi = beta * S * I / N
[-foi, foi - gamma * I, gamma * I]
sir_ode
() -> :
():
S, E, I, H, R = y
foi = beta * S * I / N
dS = -foi
dE = foi - sigma * E
dI = sigma * E - gamma * I
dH = hosp_rate * gamma * I - gamma_h * H
dR = ( - hosp_rate) * gamma * I + gamma_h * H
[dS, dE, dI, dH, dR]
seihr_ode
() -> :
T = (incidence)
t_eval = np.arange(T)
I0 =
S0 = N - I0 - initial_exposed
param_names = (param_bounds.keys())
bounds_list = [param_bounds[k] k param_names]
():
kwargs = ((param_names, params_vec))
kwargs[] = N
:
ode = model_func(**kwargs)
TypeError:
y0_len = ode.__code__.co_varcount -
n_states = (ode, , ).lower()
(model_func, , ).lower():
n_states =
y0_map = {
: [S0, I0, ],
: [S0, initial_exposed, I0, ],
: [S0, initial_exposed, I0, , ],
}
y0 = y0_map.get(n_states, [S0, initial_exposed, I0, ])
:
sol = solve_ivp(
ode, [, T - ], y0,
t_eval=t_eval, method=,
max_step=, dense_output=,
)
sol.success:
cumulative = sol.y[compartment_idx - ] + sol.y[compartment_idx]
modelled = np.maximum(np.diff(cumulative, prepend=cumulative[]), )
(np.sqrt(np.mean((modelled - incidence) ** )))
Exception:
method == :
result = differential_evolution(
_residuals, bounds_list,
maxiter=, tol=, seed=,
workers=, polish=,
)
:
x0 = [(lo + hi) / lo, hi bounds_list]
result = minimize(_residuals, x0, method=,
options={: , : })
fitted_params = ((param_names, result.x))
fitted_params[] = N
ode_fitted = model_func(**fitted_params)
n_states_f =
y0_final = [S0, initial_exposed, I0, ]
sol_final = solve_ivp(
ode_fitted, [, T - ], y0_final,
t_eval=t_eval, method=, max_step=,
)
cumulative_f = sol_final.y[compartment_idx - ] + sol_final.y[compartment_idx]
fitted_inc = np.maximum(np.diff(cumulative_f, prepend=cumulative_f[]), )
R0 =
fitted_params fitted_params:
R0 = fitted_params[] / fitted_params[]
{
: {k: v k, v fitted_params.items() k != },
: sol_final,
: fitted_inc,
: result.fun,
: R0,
}
() -> pd.DataFrame:
scipy.stats gamma gamma_dist
T = (incidence)
cv2 = (si_sd / si_mean) **
k = / cv2
theta = si_mean * cv2
rt_vals = np.full(T, np.nan)
rt_lower = np.full(T, np.nan)
rt_upper = np.full(T, np.nan)
r_vals = np.full(T, np.nan)
inc_smooth = pd.Series(incidence).rolling(window, center=,
min_periods=).mean().values
t (window, T - window):
segment = inc_smooth[t - window:t + window + ]
np.(segment <= ) np.(np.isnan(segment)):
log_inc = np.log(segment + )
time_idx = np.arange((segment)) - window
slope, intercept = np.polyfit(time_idx, log_inc, )
r = slope
rt = ( + r * theta) ** k
boot_rt = []
rng = np.random.default_rng(t)
_ ():
noise = rng.normal(, * np.(r) + , (segment))
log_b = log_inc + noise
b_slope, _ = np.polyfit(time_idx, log_b, )
boot_rt.append(( + b_slope * theta) ** k)
rt_lower[t] = (, np.quantile(boot_rt, quantiles[]))
rt_upper[t] = np.quantile(boot_rt, quantiles[])
rt_vals[t] = rt
r_vals[t] = r
df = pd.DataFrame({
: rt_vals,
: rt_lower,
: rt_upper,
: r_vals,
: ~np.isnan(rt_vals),
})
df
() -> plt.Figure:
n_panels = rt_df
fig, axes = plt.subplots(n_panels, , figsize=(, * n_panels),
sharex=)
n_panels == :
axes = [axes]
ax = axes[]
ax.bar(dates, cases, color=, alpha=, label=, width=)
model_output :
ax.plot(dates, model_output, color=, linewidth=,
label=)
ax.set_ylabel()
ax.set_title(title, fontsize=)
ax.legend()
ax.xaxis.set_major_formatter(mdates.DateFormatter())
ax.xaxis.set_major_locator(mdates.MonthLocator())
rt_df :
ax2 = axes[]
valid = rt_df[]
ax2.plot(np.array(dates)[valid], rt_df.loc[valid, ],
color=, linewidth=, label=)
ax2.fill_between(
np.array(dates)[valid],
rt_df.loc[valid, ],
rt_df.loc[valid, ],
color=, alpha=, label=,
)
ax2.axhline(, color=, linewidth=, linestyle=)
ax2.set_ylabel()
ax2.set_ylim(, (, rt_df[].quantile() * ))
ax2.legend()
fig.autofmt_xdate()
fig.tight_layout()
output_path:
fig.savefig(output_path, dpi=, bbox_inches=)
()
fig
() -> :
n = (N_by_age)
N_total = N_by_age.()
():
y = np.reshape(y, (, n))
S, E, I, R = y
foi = beta * (contact_matrix @ (I / N_by_age))
dS = -foi * S
dE = foi * S - sigma * E
dI = sigma * E - gamma * I
dR = gamma * I
np.concatenate([dS, dE, dI, dR])
ode
"""
Fit a compartmental model to observed incidence data using least-squares
optimisation.
Parameters
----------
incidence : np.ndarray
Daily or weekly new case counts (length T).
dates : array-like
Date array corresponding to incidence (length T).
model_func : Callable
A function that accepts (beta, ...) and returns a scipy solve_ivp ODE
callable. Must be consistent with ``param_bounds`` keys.
N : int
Population size.
param_bounds : dict
Dictionary of {parameter_name: (lower_bound, upper_bound)}.
Keys determine order of optimisation variables.
compartment_idx : int
Index of the ODE state vector corresponding to the observed count
(e.g. 2 for I in SEIR with states [S,E,I,R]).
initial_exposed : int
Assumed initial exposed individuals (E0).
method : str
Optimisation method: 'differential_evolution' (global, robust) or
'Nelder-Mead' (fast, local).
Returns
-------
dict with keys:
- 'params' : fitted parameter dict
- 'solution' : solve_ivp solution object
- 'fitted_incidence': np.ndarray of modelled new cases
- 'rmse' : root mean squared error
- 'R0' : basic reproduction number (beta/gamma if present)
"""
len
# Initial conditions: one infectious seed
1
list
for
in
def
_residuals
params_vec
dict
zip
"N"
try
except
return
1e12
2
# rough estimate; override below
# Determine y0 from compartment count via a probe call
4
if
"seir"
in
getattr
"__doc__"
""
else
3
if
"seihr"
in
getattr
"__name__"
""
5
3
0
4
0
5
0
0
0
try
0
1
"Radau"
1.0
False
if
not
return
1e12
# Convert cumulative to incidence via daily difference
"""
Estimate the time-varying reproduction number Rt using the Wallinga-Lipsitch
exponential growth method with a discretised serial interval distribution.
Parameters
----------
incidence : np.ndarray
Daily new case counts (length T). Should be smoothed (7-day rolling average).
si_mean : float
Mean serial interval in days.
si_sd : float
Standard deviation of serial interval in days.
window : int
Rolling window size in days for growth rate estimation.
quantiles : tuple
Lower and upper quantile for the uncertainty interval (default 95% CI).
Returns
-------
pd.DataFrame with columns: ['Rt', 'Rt_lower', 'Rt_upper', 'r', 'valid']
Notes
-----
The relationship Rt = 1 / M(-r) where M is the moment generating function of
the serial interval distribution is used (Wallinga & Lipsitch 2007).
For a gamma-distributed SI with mean mu and sd sigma:
Rt = (1 + r * sigma^2 / mu)^(mu^2 / sigma^2)
"""
from
import
as
len
# Fit gamma distribution to serial interval
2
1.0
# shape
# scale
True
3
for
in
range
1
if
any
0
or
any
continue
# Log-linear regression to estimate instantaneous growth rate r
"""
Build an age-structured SEIR ODE system.
Parameters
----------
contact_matrix : np.ndarray
(n_age x n_age) POLYMOD-style contact matrix. Entry C[i,j] is the
average number of contacts individuals in age group i make with
age group j per day.
beta : float
Per-contact transmission probability.
sigma : float
Incubation rate (scalar, shared across age groups).
gamma : float
Recovery rate (scalar, shared across age groups).
N_by_age : np.ndarray
Population size in each age group (length n_age).
Returns
-------
Callable
ODE function f(t, y) with y flattened as [S0,S1,...,E0,E1,...,I0,I1,...,R0,R1,...].
"""
len
sum
def
ode
t, y
4
# Force of infection for each age group
return
return
Example 1: Fit SEIR to COVID-19 Data and Estimate Rt
"""
example_covid_seir_rt.py
-------------------------
Fit a SEIR model to a COVID-19-like incidence curve and estimate Rt
using the sliding-window exponential growth method.
"""import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from epi_modeling import (
build_seir_model,
fit_model_to_data,
estimate_rt,
plot_epidemic_curve,
)
# ---- 1. Generate synthetic "observed" COVID-like incidence ----# In production, replace this with real data from:# WHO: https://covid19.who.int/data# CDC: https://data.cdc.gov
rng = np.random.default_rng(7)
N_pop = 5_000_000# True parameters (to be recovered by fitting)
BETA_TRUE = 0.28
SIGMA_TRUE = 1 / 5.1# mean incubation 5.1 days
GAMMA_TRUE = 1 / 10.0# mean infectious period 10 daysfrom scipy.integrate import solve_ivp
true_ode = build_seir_model(BETA_TRUE, SIGMA_TRUE, GAMMA_TRUE, N_pop)
T_days = 120
sol_true = solve_ivp(
true_ode, [0, T_days - 1],
[N_pop - 20, 10, 10, 0],
t_eval=np.arange(T_days), method="Radau",
)
S_true, E_true, I_true, R_true = sol_true.y
new_infections_true = np.maximum(np.diff(SIGMA_TRUE * E_true, prepend=0), 0)
observed = rng.poisson(np.clip(new_infections_true, 0, None)).astype(float)
start_date = pd.Timestamp("2023-01-01")
dates = pd.date_range(start_date, periods=T_days, freq="D")
# ---- 2. Fit SEIR model ----print("Fitting SEIR model (this may take ~30-60 seconds) ...")
fit_result = fit_model_to_data(
incidence=observed,
dates=dates,
model_func=build_seir_model,
N=N_pop,
param_bounds={
"beta": (0.05, 1.0),
"sigma": (1/14, 1/2),
"gamma": (1/21, 1/4),
},
compartment_idx=2,
initial_exposed=10,
method="differential_evolution",
)
fp = fit_result["params"]
print(f"\n--- Fitted Parameters ---")
print(f" beta : {fp['beta']:.4f} (true: {BETA_TRUE:.4f})")
print(f" sigma : {fp['sigma']:.4f} (true: {SIGMA_TRUE:.4f})")
print(f" gamma : {fp['gamma']:.4f} (true: {GAMMA_TRUE:.4f})")
print(f" R0 : {fit_result['R0']:.2f} (true: {BETA_TRUE/GAMMA_TRUE:.2f})")
print(f" RMSE : {fit_result['rmse']:.1f} cases/day")
# ---- 3. Estimate Rt ----
rt_df = estimate_rt(
incidence=observed,
si_mean=7.5, # mean serial interval (days)
si_sd=3.4,
window=7,
)
rt_df.index = range(len(rt_df))
# ---- 4. Plot epidemic curve + Rt ----
fig = plot_epidemic_curve(
dates=dates,
cases=observed,
model_output=fit_result["fitted_incidence"],
rt_df=rt_df,
title="COVID-19-like Outbreak: SEIR Fit and Rt Estimation",
output_path="seir_covid_fit.png",
)
plt.show()
# ---- 5. Peak incidence statistics ----
peak_day = np.argmax(fit_result["fitted_incidence"])
print(f"\nPeak incidence (model): {fit_result['fitted_incidence'][peak_day]:.0f} cases "f"on {dates[peak_day].strftime('%Y-%m-%d')} (day {peak_day})")
valid_rt = rt_df[rt_df["valid"]]
rt_above_1 = (valid_rt["Rt"] > 1).sum()
print(f"Days with Rt > 1: {rt_above_1} / {len(valid_rt)}")
"""
sensitivity_analysis.py
-----------------------
Partial Rank Correlation Coefficient (PRCC) global sensitivity analysis
for epidemic model outputs (e.g., peak prevalence, final attack rate).
"""import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy.stats import spearmanr, rankdata
from scipy.integrate import solve_ivp
from epi_modeling import build_seir_model
defprcc_sensitivity(
param_distributions: dict,
model_func: callable,
output_func: callable,
N: int,
n_samples: int = 1000,
seed: int = 42,
) -> pd.DataFrame:
"""
Compute PRCC indices for each parameter against a scalar model output.
Parameters
----------
param_distributions : dict
{param_name: (low, high)} uniform distribution bounds.
model_func : callable
Function(params_dict, N) -> solve_ivp solution.
output_func : callable
Function(solution) -> scalar output of interest.
N : int
Population size.
n_samples : int
LHS sample size.
seed : int
Returns
-------
pd.DataFrame with columns: param, prcc, p_value.
"""
rng = np.random.default_rng(seed)
param_names = list(param_distributions.keys())
k = len(param_names)
# Latin Hypercube Sampling
lhs = np.zeros((n_samples, k))
for j, name inenumerate(param_names):
lo, hi = param_distributions[name]
perm = rng.permutation(n_samples)
u = (perm + rng.uniform(size=n_samples)) / n_samples
lhs[:, j] = lo + u * (hi - lo)
# Evaluate model
outputs = np.full(n_samples, np.nan)
for i inrange(n_samples):
params = dict(zip(param_names, lhs[i]))
try:
sol = model_func(params, N)
outputs[i] = output_func(sol)
except Exception:
pass
valid = ~np.isnan(outputs)
lhs_v = lhs[valid]
out_v = outputs[valid]
# Rank-transform
ranked = np.column_stack([rankdata(lhs_v[:, j]) for j inrange(k)] +
[rankdata(out_v)])
# Partial correlations via residuals
prcc_vals, pvals = [], []
for j inrange(k):
other = [jj for jj inrange(k) if jj != j]
X = ranked[:, j]
Y = ranked[:, k]
Z = ranked[:, other]
# Regress X on Z
coef_x = np.linalg.lstsq(np.column_stack([np.ones(len(X)), Z]), X, rcond=None)[0]
res_x = X - np.column_stack([np.ones(len(X)), Z]) @ coef_x
# Regress Y on Z
coef_y = np.linalg.lstsq(np.column_stack([np.ones(len(Y)), Z]), Y, rcond=None)[0]
res_y = Y - np.column_stack([np.ones(len(Y)), Z]) @ coef_y
r, p = spearmanr(res_x, res_y)
prcc_vals.append(r)
pvals.append(p)
df = pd.DataFrame({
"param": param_names,
"prcc": prcc_vals,
"p_value": pvals,
}).sort_values("prcc", key=abs, ascending=False)
return df
if __name__ == "__main__":
N = 1_000_000defrun_seir(params, N):
ode = build_seir_model(params["beta"], params["sigma"], params["gamma"], N)
return solve_ivp(ode, [0, 200], [N - 11, 10, 1, 0],
t_eval=np.arange(201), method="Radau")
defpeak_prevalence(sol):
return sol.y[2].max() / N * 100
results = prcc_sensitivity(
param_distributions={
"beta": (0.1, 0.6),
"sigma": (1/14, 1/3),
"gamma": (1/21, 1/3),
},
model_func=run_seir,
output_func=peak_prevalence,
N=N,
n_samples=500,
)
print("\nPRCC Sensitivity Analysis — Peak Infectious Prevalence")
print(results.to_string(index=False))
fig, ax = plt.subplots(figsize=(7, 4))
colors = ["firebrick"if v > 0else"steelblue"for v in results["prcc"]]
ax.barh(results["param"], results["prcc"], color=colors)
ax.axvline(0, color="black", linewidth=0.8)
ax.set_xlabel("PRCC")
ax.set_title("Sensitivity Analysis: Peak Infectious Prevalence")
fig.tight_layout()
fig.savefig("prcc_sensitivity.png", dpi=150)
print("Saved: prcc_sensitivity.png")
plt.show()
Tips and Best Practices
Solver choice: Use method="Radau" for stiff systems (large populations, wide R0 ranges).
LSODA is also stiff-capable. RK45 is fast but may fail for stiff problems.
Initial conditions: Seed with a small number of exposed individuals (E0 > 0) to allow
the latent period to seed infectious cases naturally; avoids discontinuities.
Identifiability: SEIR with only I(t) observed is structurally identifiable only when
sigma is fixed from literature. Fitting all three parameters simultaneously risks
converging to incorrect minima.
Data smoothing: Apply a 7-day rolling average to raw case counts before Rt estimation
to suppress weekly reporting artefacts.
Serial interval vs generation time: The serial interval (time between symptom onsets)
is observable; the generation time (time between infections) is not. For Rt estimation,
use the serial interval.
Age-structured models: Use POLYMOD contact matrices (Mossong et al., 2008) or
country-specific matrices from the socialmixr R package (exportable to CSV).
References
Kermack & McKendrick (1927). A contribution to the mathematical theory of epidemics.
Proc. R. Soc. A, 115, 700–721.
Wallinga & Lipsitch (2007). How generation intervals shape the relationship between
growth rates and reproductive numbers. Proc. R. Soc. B, 274, 599–604.
Cori et al. (2013). A new framework and software to estimate time-varying reproduction
numbers during epidemics. Am. J. Epidemiol., 178, 1505–1512.
Marino et al. (2008). A methodology for performing global uncertainty and sensitivity
analysis in systems biology. J. Theor. Biol., 254, 178–196.
