Compartmental transmission model

Build an infectious-disease model with ODEModel, a susceptible-exposed-infectious-recovered model of a vaccination program.

Keywords

compartmental model, SEIR model, dynamic transmission model, infectious disease modeling, vaccination cost-effectiveness, ordinary differential equation model, health economics, Python

This tutorial shows how to build an infectious-disease cost-effectiveness model with ODEModel, the engine for systems of ordinary differential equations, using a susceptible-exposed-infectious-recovered (SEIR) model of a vaccination program. In this model one person’s outcome depends on the others. Vaccinating some susceptible people protects the rest by lowering the force of infection, an effect that neither a cohort model nor a microsimulation of independent individuals can represent. Full script: examples/seir_vaccination.py.

Specifying the compartments and their dynamics

The model follows a closed population of 100,000 through an epidemic of a directly transmitted, non-fatal infection. Five compartments track the susceptible (S), the exposed but not yet infectious (E), the infectious (I), the recovered and immune (R), and the vaccinated and immune (V). The dynamics_and_rewards function returns an ODESpec: the right-hand side of the system (the rate of change of each compartment), the initial compartment sizes, and the reward rates.

The coupling between people enters through the force of infection, the rate at which a susceptible person becomes infected. It is proportional to the current infectious prevalence, $\lambda(t) = \beta\, I / N$, so each new infection raises the risk for everyone still susceptible. The transmission rate follows from the basic reproduction number as $\beta = R_0\,\gamma$, where $R_0$ is the number of secondary infections one case produces in a fully susceptible population. A cohort model cannot represent this dependence.

Rewards use the two channels the engine offers, matched to continuous time. Quality-adjusted life-years accrue on compartment occupancy through state_effect. Everyone healthy contributes one per year and the infectious contribute less while ill, so averted illness raises the effect. Costs fall on flows rather than states. event_rates returns the per-year rate of two events, doses administered and new infections, and event_cost charges a one-time amount to each. The engine integrates the discounted reward flows alongside the compartments.

import numpy as np
import pandas as pd
from heormodel.models import ODEModel, ODESpec

STATES = ("S", "E", "I", "R", "V")
INTERVENTIONS = ("No vaccination", "Vaccination program")
POPULATION, INITIAL_INFECTIOUS = 100_000.0, 10.0

def seir(p, intervention):
    beta = p["R0"] * p["gamma"]  # transmission rate from the reproduction number
    nu = p["nu"] if intervention == "Vaccination program" else 0.0
    def derivatives(t, y):
        s, e, i, r, v = y
        foi = beta * i / (s + e + i + r + v)  # force of infection
        new_infections, doses = foi * s, nu * s
        return np.array([-new_infections - doses, new_infections - p["sigma"] * e,
                         p["sigma"] * e - p["gamma"] * i, p["gamma"] * i, doses])
    def event_rates(t, y):
        s, e, i, r, v = y
        foi = beta * i / (s + e + i + r + v)
        return np.array([nu * s, foi * s])  # doses per year, new infections per year
    return ODESpec(
        derivatives=derivatives,
        initial=np.array([POPULATION - INITIAL_INFECTIOUS, 0.0, INITIAL_INFECTIOUS, 0.0, 0.0]),
        state_cost=np.zeros(5), state_effect=np.array([1.0, 1.0, p["u_I"], 1.0, 1.0]),
        event_rates=event_rates,
        event_cost=np.array([p["c_vacc"], p["c_case"]]), event_effect=np.zeros(2))

engine = ODEModel(states=STATES, interventions=INTERVENTIONS,
                  dynamics_and_rewards=seir, horizon=10.0, discount_rate=0.03)

Reading the epidemic curves

trajectory integrates the compartments for one parameter set and returns them over time. It is useful for inspecting the model before relying on its cost-effectiveness output. Plotting the two arms side by side shows what the program does. Without vaccination, infection depletes the susceptible pool and produces a wave of recovered individuals. Under the program, susceptible people move to the vaccinated compartment before the epidemic grows.

import matplotlib.pyplot as plt

base = pd.Series(dict(R0=1.8, sigma=4.0, gamma=6.0, nu=1.0, c_vacc=200.0, c_case=150.0, u_I=0.7))
no_vacc, vacc = engine.trajectory(base, "No vaccination"), engine.trajectory(base, "Vaccination program")

fig, axes = plt.subplots(1, 2, figsize=(9.5, 3.6), sharey=True)

for ax, traj, title in ((axes[0], no_vacc, "No vaccination"), (axes[1], vacc, "Vaccination program")):
    for comp in ("S", "I", "R", "V"):
        ax.plot(traj["time"], traj[comp], label=comp)
    ax.set_title(title); ax.set_xlabel("Year")

axes[0].set_ylabel("People"); axes[1].legend(loc="center right")

plt.show()

Vaccination reduces the total number ever infected over the ten-year horizon from about 73,000 to under 100, nearly eliminating the epidemic.

Computing the base-case cost-effectiveness results

Integrating each arm once at the point estimates gives the deterministic cost-effectiveness result. Vaccination is not cost-saving here. The doses reach nearly everyone while the treatment cost per infection is modest, so the program costs more than the treatment it averts. It is cost-effective, with an incremental cost-effectiveness ratio around 3,200 per quality-adjusted life-year, far below a 50,000 threshold.

from heormodel.cea import icer_table

draws0 = pd.DataFrame([base.to_dict()], index=pd.RangeIndex(1, name="iteration"))

icer_table(engine.evaluate(draws0)).round(2)

	cost	effect	inc_cost	inc_effect	icer	status
intervention
No vaccination	9512537.48	860808.38	NaN	NaN	NaN	ND
Vaccination program	19412280.17	863935.86	9899742.69	3127.48	3165.41	ND

Running the probabilistic sensitivity analysis

The epidemic and cost parameters go into a ParameterSet. The reproduction number, the transition rates, and the costs take log-normal distributions, which keep every draw positive. The infectious-state utility takes a beta distribution, which keeps it between zero and one. The reproduction number carries the widest spread because it is the least certain and, through the force of infection, the parameter the epidemic is most sensitive to.

As in the other engines, we define the parameter set, sample it, run the analysis with run_psa, and read the incremental cost-effectiveness ratio table.

from heormodel.params import Beta, LogNormal, ParameterSet
from heormodel.run import SeedManager, run_psa

params = ParameterSet({
    "R0": LogNormal(np.log(1.8), 0.30), "sigma": LogNormal(np.log(4.0), 0.1),
    "gamma": LogNormal(np.log(6.0), 0.1), "nu": LogNormal(np.log(1.0), 0.15),
    "c_vacc": LogNormal(np.log(200.0), 0.25), "c_case": LogNormal(np.log(150.0), 0.25),
    "u_I": Beta(8.4, 3.6)})

draws = params.sample(1000, seed=SeedManager(20260711).generator())

outcomes = run_psa(engine, draws).outcomes

icer_table(outcomes).round(2)

	cost	cost_lo	cost_hi	effect	effect_lo	effect_hi	inc_cost	inc_cost_lo	inc_cost_hi	inc_effect	inc_effect_lo	inc_effect_hi	icer	icer_lo	icer_hi	status
intervention
No vaccination	8547209.69	35769.97	18139751.27	861130.58	857117.30	863928.37	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	ND
Vaccination program	20001324.57	12252301.29	30963512.04	863933.24	863913.01	863938.62	11454114.88	-1213791.31	25880147.58	2802.66	9.37	6811.87	4086.87	-274.22	2202925.09	ND

Cost-effectiveness acceptability curve (CEAC)

The cost-effectiveness acceptability curve reports, at each willingness-to-pay threshold, the probability that each arm is the best alternative across the draws. Vaccination becomes the preferred arm once the threshold passes its incremental cost-effectiveness ratio, and the crossover is where the two curves meet.

from heormodel.cea import ceac
from heormodel.report import plot_ceac
plot_ceac(ceac(outcomes, np.linspace(0, 30_000, 61)))
plt.show()

The value-of-information functions read this engine’s output like any other. The expected value of perfect information (EVPI) puts a monetary value on resolving the remaining uncertainty at the willingness-to-pay threshold.

from heormodel.voi import evpi, evppi_ranking
print(f"EVPI at WTP 50,000: {evpi(outcomes, 50_000.0):,.0f}")

EVPI at WTP 50,000: 1,938,103

evppi_ranking attributes that value to individual parameters, so a study can target whichever one the decision is most sensitive to.

evppi_ranking(outcomes, draws, 50_000.0).round(0).head()

R0        2358086.0
u_I             0.0
c_case          0.0
gamma           0.0
sigma           0.0
Name: evppi, dtype: float64

The ordinary differential equation engine returns the same Outcomes structure as the other engines, so the cost-effectiveness and value-of-information code is unchanged.

Keywords: compartmental model, SEIR model, dynamic transmission model, infectious disease modeling, vaccination cost-effectiveness, ordinary differential equation model, health economics, Python

--- title: Compartmental transmission model description: "Build an infectious-disease model with `ODEModel`, a susceptible-exposed-infectious-recovered model of a vaccination program." keywords: [compartmental model, SEIR model, dynamic transmission model, infectious disease modeling, vaccination cost-effectiveness, ordinary differential equation model, health economics, Python] ---  <a href="https://colab.research.google.com/github/pedroliman/heormodel/blob/main/docs/_notebooks/seir-vaccination.ipynb"><img src="https://colab.research.google.com/assets/colab-badge.svg" alt="Open In Colab"></a>  This tutorial shows how to build an infectious-disease cost-effectiveness model with `ODEModel`, the engine for systems of ordinary differential equations, using a susceptible-exposed-infectious-recovered (SEIR) model of a vaccination program. In this model one person's outcome depends on the others. Vaccinating some susceptible people protects the rest by lowering the force of infection, an effect that neither a cohort model nor a microsimulation of independent individuals can represent. Full script: [`examples/seir_vaccination.py`](https://github.com/pedroliman/heormodel/blob/main/examples/seir_vaccination.py). ## Specifying the compartments and their dynamics The model follows a closed population of 100,000 through an epidemic of a directly transmitted, non-fatal infection. Five compartments track the susceptible (S), the exposed but not yet infectious (E), the infectious (I), the recovered and immune (R), and the vaccinated and immune (V). The `dynamics_and_rewards` function returns an `ODESpec`: the right-hand side of the system (the rate of change of each compartment), the initial compartment sizes, and the reward rates. The coupling between people enters through the force of infection, the rate at which a susceptible person becomes infected. It is proportional to the current infectious prevalence, $\lambda(t) = \beta\, I / N$, so each new infection raises the risk for everyone still susceptible. The transmission rate follows from the basic reproduction number as $\beta = R_0\,\gamma$, where $R_0$ is the number of secondary infections one case produces in a fully susceptible population. A cohort model cannot represent this dependence. Rewards use the two channels the engine offers, matched to continuous time. Quality-adjusted life-years accrue on compartment occupancy through `state_effect`. Everyone healthy contributes one per year and the infectious contribute less while ill, so averted illness raises the effect. Costs fall on flows rather than states. `event_rates` returns the per-year rate of two events, doses administered and new infections, and `event_cost` charges a one-time amount to each. The engine integrates the discounted reward flows alongside the compartments. ```{python} import numpy as np import pandas as pd from heormodel.models import ODEModel, ODESpec STATES = ("S", "E", "I", "R", "V") INTERVENTIONS = ("No vaccination", "Vaccination program") POPULATION, INITIAL_INFECTIOUS = 100_000.0, 10.0 def seir(p, intervention): beta = p["R0"] * p["gamma"] # transmission rate from the reproduction number nu = p["nu"] if intervention == "Vaccination program" else 0.0 def derivatives(t, y): s, e, i, r, v = y foi = beta * i / (s + e + i + r + v) # force of infection new_infections, doses = foi * s, nu * s return np.array([-new_infections - doses, new_infections - p["sigma"] * e, p["sigma"] * e - p["gamma"] * i, p["gamma"] * i, doses]) def event_rates(t, y): s, e, i, r, v = y foi = beta * i / (s + e + i + r + v) return np.array([nu * s, foi * s]) # doses per year, new infections per year return ODESpec( derivatives=derivatives, initial=np.array([POPULATION - INITIAL_INFECTIOUS, 0.0, INITIAL_INFECTIOUS, 0.0, 0.0]), state_cost=np.zeros(5), state_effect=np.array([1.0, 1.0, p["u_I"], 1.0, 1.0]), event_rates=event_rates, event_cost=np.array([p["c_vacc"], p["c_case"]]), event_effect=np.zeros(2)) engine = ODEModel(states=STATES, interventions=INTERVENTIONS, dynamics_and_rewards=seir, horizon=10.0, discount_rate=0.03) ``` ## Reading the epidemic curves `trajectory` integrates the compartments for one parameter set and returns them over time. It is useful for inspecting the model before relying on its cost-effectiveness output. Plotting the two arms side by side shows what the program does. Without vaccination, infection depletes the susceptible pool and produces a wave of recovered individuals. Under the program, susceptible people move to the vaccinated compartment before the epidemic grows. ```{python} import matplotlib.pyplot as plt base = pd.Series(dict(R0=1.8, sigma=4.0, gamma=6.0, nu=1.0, c_vacc=200.0, c_case=150.0, u_I=0.7)) no_vacc, vacc = engine.trajectory(base, "No vaccination"), engine.trajectory(base, "Vaccination program") fig, axes = plt.subplots(1, 2, figsize=(9.5, 3.6), sharey=True) for ax, traj, title in ((axes[0], no_vacc, "No vaccination"), (axes[1], vacc, "Vaccination program")): for comp in ("S", "I", "R", "V"): ax.plot(traj["time"], traj[comp], label=comp) ax.set_title(title); ax.set_xlabel("Year") axes[0].set_ylabel("People"); axes[1].legend(loc="center right") plt.show() ``` Vaccination reduces the total number ever infected over the ten-year horizon from about 73,000 to under 100, nearly eliminating the epidemic. ## Computing the base-case cost-effectiveness results Integrating each arm once at the point estimates gives the deterministic cost-effectiveness result. Vaccination is not cost-saving here. The doses reach nearly everyone while the treatment cost per infection is modest, so the program costs more than the treatment it averts. It is cost-effective, with an incremental cost-effectiveness ratio around 3,200 per quality-adjusted life-year, far below a 50,000 threshold. ```{python} from heormodel.cea import icer_table draws0 = pd.DataFrame([base.to_dict()], index=pd.RangeIndex(1, name="iteration")) icer_table(engine.evaluate(draws0)).round(2) ``` ## Running the probabilistic sensitivity analysis The epidemic and cost parameters go into a `ParameterSet`. The reproduction number, the transition rates, and the costs take log-normal distributions, which keep every draw positive. The infectious-state utility takes a beta distribution, which keeps it between zero and one. The reproduction number carries the widest spread because it is the least certain and, through the force of infection, the parameter the epidemic is most sensitive to. As in the other engines, we define the parameter set, sample it, run the analysis with `run_psa`, and read the incremental cost-effectiveness ratio table. ```{python} from heormodel.params import Beta, LogNormal, ParameterSet from heormodel.run import SeedManager, run_psa params = ParameterSet({ "R0": LogNormal(np.log(1.8), 0.30), "sigma": LogNormal(np.log(4.0), 0.1), "gamma": LogNormal(np.log(6.0), 0.1), "nu": LogNormal(np.log(1.0), 0.15), "c_vacc": LogNormal(np.log(200.0), 0.25), "c_case": LogNormal(np.log(150.0), 0.25), "u_I": Beta(8.4, 3.6)}) draws = params.sample(1000, seed=SeedManager(20260711).generator()) outcomes = run_psa(engine, draws).outcomes icer_table(outcomes).round(2) ``` ## Cost-effectiveness acceptability curve (CEAC) The cost-effectiveness acceptability curve reports, at each willingness-to-pay threshold, the probability that each arm is the best alternative across the draws. Vaccination becomes the preferred arm once the threshold passes its incremental cost-effectiveness ratio, and the crossover is where the two curves meet. ```{python} from heormodel.cea import ceac from heormodel.report import plot_ceac plot_ceac(ceac(outcomes, np.linspace(0, 30_000, 61))) plt.show() ``` The value-of-information functions read this engine's output like any other. The expected value of perfect information (EVPI) puts a monetary value on resolving the remaining uncertainty at the willingness-to-pay threshold. ```{python} from heormodel.voi import evpi, evppi_ranking print(f"EVPI at WTP 50,000: {evpi(outcomes, 50_000.0):,.0f}") ``` `evppi_ranking` attributes that value to individual parameters, so a study can target whichever one the decision is most sensitive to. ```{python} evppi_ranking(outcomes, draws, 50_000.0).round(0).head() ``` The ordinary differential equation engine returns the same `Outcomes` structure as the other engines, so the cost-effectiveness and value-of-information code is unchanged.