OpenSim.zip

OS-029 · Population & epidemics

SIR Epidemic Model Simulation

The SIR model splits a population into susceptible, infected and recovered people and moves them between the groups at two rates. Three differential equations, two parameters — and already the familiar rise, peak and fall of an outbreak.

Interactive model

OS-029 · sir-model · deterministic
S — susceptibleI — infectedR — recoveredS = N/R₀peak day
Parameters
Closed population, no births or deaths
Effective contacts per infected person per day
1/γ is the mean infectious period
Basic reproduction number R₀
3.00
Peak infected
30,050
Peak day
49.9
Ever infected by end
94.0%
Final size (theory, t → ∞)
94.0%
max |S + I + R − N|
2.0e-10

Example run: default parameters.

What the simulation shows

Everyone starts susceptible except a handful of infected people. Infection moves people from S to I at a rate that depends on how many susceptible and infected people there are; recovery moves them from I to R at a constant per-person rate. Nobody enters or leaves the population, and recovered people stay immune.

  • R₀ decides whether there is an epidemic at all. With β/γ below 1 each infected person causes fewer than one new infection and the curve of I only falls. Push β above γ and a wave appears.
  • The peak happens when S reaches N/R₀. The dashed blue line marks that level. Once susceptibles drop below it, each case replaces itself less than once and infections decline — even though many people are still susceptible.
  • Not everyone gets infected. The epidemic overshoots the threshold and then burns out with some people never infected. The lab compares the simulated total with the final-size equation below.

Try doubling the infectious period (halve γ) with β fixed, then halve β instead: both change R₀, but they change the timing differently.

The equations

With population N = S + I + R, transmission rate β and recovery rate γ (both per day):

dS/dt = −β · S · I / N
dI/dt = β · S · I / N − γ · I
dR/dt = γ · I

The three right-hand sides sum to zero, so N stays constant; the lab reports the largest numerical drift from that as a check. The mean infectious period is 1/γ days.

Two quantities follow directly from the equations:

  • Basic reproduction number R₀ = β/γ — new infections caused by one case in a fully susceptible population.
  • Peak condition — dI/dt = 0 when βS/N = γ, that is when S = N/R₀.

Dividing dS/dt by dR/dt and integrating gives the final-size relation for an epidemic that starts with R(0) = 0 and ends with no infected people:

S∞ = S0 · exp(−R₀ · (N − S∞) / N)

For R₀ = 3 the solution is about 94% ever infected, well above the 1 − 1/R₀ ≈ 67% herd-immunity threshold.

How the simulation solves it

The system has no closed-form solution for I(t), so the model steps it forward with the classical fourth-order Runge–Kutta method (RK4) at a fixed step of 0.05 days, recording the state once per day. RK4 evaluates the derivatives four times per step and combines them; its error per unit time shrinks with the fourth power of the step size. For these smooth curves a 0.05-day step is far more accurate than needed — halving it changes the reported peak by less than one person. Euler vs Runge–Kutta shows why the simpler Euler method is a poor choice here.

The peak is taken from the RK4 steps themselves, so it is resolved to 0.05 days.

Example result

With β = 0.3, γ = 0.1 (R₀ = 3.00) in a population of 100,000 starting from 10 cases, infections peak around day 49.9 at about 30,050 people. By day 180, 94.0% of the population has been infected, matching the final-size equation (94.0%). The sum S + I + R never drifts more than 2.0e-10 people from N.

Output of sir-model with population=100000 initial_infected=10 beta=0.3 gamma=0.1 days=180
Basic reproduction number R₀3.00
Peak infected30,050
Peak day49.9
Ever infected by end94.0%
Final size (theory, t → ∞)94.0%
max |S + I + R − N|2.0e-10

Implementation

The page runs simulation.js from the package below; the Python file is a line-by-line port with the same random-number generator. The excerpts show the core of each.

JavaScript — simulation.js
// dS/dt = −βSI/N,  dI/dt = βSI/N − γI,  dR/dt = γI
export function derivatives(s, i, beta, gamma, n) {
  const infection = (beta * s * i) / n;
  const recovery = gamma * i;
  return [-infection, infection - recovery, recovery];
}

// One classical Runge–Kutta (RK4) step of size h.
export function rk4Step(state, h, beta, gamma, n) {
  const [s, i, r] = state;
  const k1 = derivatives(s, i, beta, gamma, n);
  const k2 = derivatives(s + (h / 2) * k1[0], i + (h / 2) * k1[1], beta, gamma, n);
  const k3 = derivatives(s + (h / 2) * k2[0], i + (h / 2) * k2[1], beta, gamma, n);
  const k4 = derivatives(s + h * k3[0], i + h * k3[1], beta, gamma, n);
  return [
    s + (h / 6) * (k1[0] + 2 * k2[0] + 2 * k3[0] + k4[0]),
    i + (h / 6) * (k1[1] + 2 * k2[1] + 2 * k3[1] + k4[1]),
    r + (h / 6) * (k1[2] + 2 * k2[2] + 2 * k3[2] + k4[2]),
  ];
}
Python — simulation.py
def derivatives(s, i, beta, gamma, n):
    """dS/dt = -beta*S*I/N, dI/dt = beta*S*I/N - gamma*I, dR/dt = gamma*I"""
    infection = beta * s * i / n
    recovery = gamma * i
    return (-infection, infection - recovery, recovery)


def rk4_step(state, h, beta, gamma, n):
    s, i, r = state
    k1 = derivatives(s, i, beta, gamma, n)
    k2 = derivatives(s + h / 2 * k1[0], i + h / 2 * k1[1], beta, gamma, n)
    k3 = derivatives(s + h / 2 * k2[0], i + h / 2 * k2[1], beta, gamma, n)
    k4 = derivatives(s + h * k3[0], i + h * k3[1], beta, gamma, n)
    return (
        s + h / 6 * (k1[0] + 2 * k2[0] + 2 * k3[0] + k4[0]),
        i + h / 6 * (k1[1] + 2 * k2[1] + 2 * k3[1] + k4[1]),
        r + h / 6 * (k1[2] + 2 * k2[2] + 2 * k3[2] + k4[2]),
    )

Run the full model from a terminal:

Shell
node simulation.js --population 100000 --initial-infected 10
python simulation.py --population 100000 --initial-infected 10

Download the model

sir-model.zip6 files · 14.0 KBDownload .zip
Files in sir-model.zip
FileContentsSize
README.mdHow to run it, parameters, example result2.2 KB
model.jsonParameters with ranges and defaults, outputs, the seed-42 example2.3 KB
simulation.jsThe model used on this page plus a command-line runner (Node 18+)7.9 KB
simulation.pyPython 3.9+ port, standard library only6.7 KB
example.csvResult table from the example run10.2 KB
LICENSE.txtMIT License1.0 KB

Assumptions & limitations

  • Homogeneous mixing. Every infected person is equally likely to meet every susceptible one. Real contact is structured by households, age, schools and geography; an epidemic on a network shows what changes when it is not.
  • Deterministic and continuous. The equations treat people as a continuous quantity and have no chance events, so an outbreak starting from 10 cases always takes off when R₀ > 1. In small populations, early chance extinction is common.
  • Exponential infectious period. A constant recovery rate γ means infectious periods are exponentially distributed; real diseases have more regular durations and often a latent (exposed) stage, handled by SEIR-type models.
  • Fixed parameters. β and γ do not change with behavior, seasons, interventions or vaccination. There are no births, deaths or waning immunity.
  • An educational model. It captures the mechanism of an outbreak, not the course of any real disease. It is not a forecasting tool.

References

  1. Kermack, W. O. & McKendrick, A. G. (1927). A contribution to the mathematical theory of epidemics. Proceedings of the Royal Society A, 115(772), 700–721. doi.org/10.1098/rspa.1927.0118
  2. Keeling, M. J. & Rohani, P. (2008). Modeling Infectious Diseases in Humans and Animals, chapter 2. Princeton University Press.
  3. Hethcote, H. W. (2000). The mathematics of infectious diseases. SIAM Review, 42(4), 599–653. doi.org/10.1137/S0036144500371907