SIR and SEIR Epidemic Models — Population biology, R₀ and herd immunity

1. The SIR compartments and assumptions

The population of size N is divided into three compartments. Every individual belongs to exactly one at any given time:

The classic model rests on three idealisations:

• Homogeneous mixing: every susceptible individual is equally likely to contact any infectious individual — the mass-action assumption.
• Closed population: N = S + I + R is constant (no births, deaths, or migration).
• Recovery grants permanent immunity: once removed, individuals never return to S.

2. System of ODEs — derivation

Denote by β the transmission rate (contacts per unit time × probability of transmission per contact) and by γ the recovery rate (1/γ = mean infection duration). The flow is:

The incidence term β · S · I / N is frequency-dependent: the N in the denominator means the per-capita contact rate is fixed regardless of population size — the appropriate form for directly-transmitted diseases where contacts are limited by time, not population density.

For density-dependent transmission (e.g., vector-borne disease, where more hosts means more contacts) the denominator is omitted: β · S · I.

3. Basic reproduction number R₀

R₀ (pronounced "R-naught") is the average number of secondary infections produced by one infectious individual introduced into a fully susceptible population. It is the single most useful epidemiological quantity:

R₀ should not be confused with the effective reproduction number Rₜ, which accounts for partial immunity and interventions: Rₜ = R₀ · S/N. The epidemic declines when Rₜ drops below 1.

4. Herd immunity threshold

At the epidemic peak, dI/dt = 0, which gives the critical susceptible fraction:

If the immune fraction exceeds h — whether by vaccination or prior infection — a single introduction cannot spark an epidemic. This is population-level (herd) immunity.

5. SEIR model — the exposed compartment

Many diseases have a latent period: individuals are infected but not yet contagious (COVID-19: ~5 days; measles ~10 days). SEIR adds an Exposed compartment to capture this:

The latent period delays the epidemic peak — this is why contact tracing works: exposed individuals can be identified and isolated before they become infectious.

6. JavaScript implementation with RK4

We use a classic fourth-order Runge-Kutta (RK4) ODE solver. The state vector is [S, E, I, R].

// ── SEIR model parameters ─────────────────────────────────────────
const params = {
  N    : 1_000_000,  // total population
  beta : 0.3,        // transmission rate (day⁻¹)
  sigma: 1 / 5,     // incubation rate (1/latent period)
  gamma: 1 / 10,    // recovery rate (1/infectious period)
};

// ── SEIR right-hand side ──────────────────────────────────────────
function seir([S, E, I, R], { N, beta, sigma, gamma }) {
  const force = beta * I / N;  // force of infection (per susceptible)
  return [
    -force * S,               // dS/dt
     force * S - sigma * E,   // dE/dt
     sigma * E - gamma * I,   // dI/dt
     gamma * I                // dR/dt
  ];
}

// ── RK4 single step ──────────────────────────────────────────────
function rk4Step(f, y, dt, p) {
  const add = (a, b, s) => a.map((v, i) => v + b[i] * s);
  const k1 = f(y,                dt, p);
  const k2 = f(add(y, k1, dt/2), dt, p);
  const k3 = f(add(y, k2, dt/2), dt, p);
  const k4 = f(add(y, k3, dt),    dt, p);
  return y.map((v, i) =>
    v + (k1[i] + 2*k2[i] + 2*k3[i] + k4[i]) * dt / 6
  );
}

// Fix rk4Step to pass correct signature
function rk4StepFixed(y, dt, p) {
  const add = (a, b, s) => a.map((v, i) => v + b[i] * s);
  const k1 = seir(y,                p);
  const k2 = seir(add(y, k1, dt/2), p);
  const k3 = seir(add(y, k2, dt/2), p);
  const k4 = seir(add(y, k3, dt),    p);
  return y.map((v, i) =>
    v + (k1[i] + 2*k2[i] + 2*k3[i] + k4[i]) * dt / 6
  );
}

// ── Run simulation for 365 days ───────────────────────────────────
function runSEIR(p, days = 365, dt = 0.25) {
  const results = [];
  let y = [p.N - 1, 0, 1, 0];  // S, E, I=1 seed, R
  for (let t = 0; t <= days; t += dt) {
    results.push({ t, S: y[0], E: y[1], I: y[2], R: y[3] });
    y = rk4StepFixed(y, dt, p);
  }
  return results;
}

// Compute R₀ and herd immunity threshold
const R0     = params.beta / params.gamma;        // 0.3 / 0.1 = 3
const herd   = 1 - 1 / R0;                       // ≈ 0.667 = 66.7%
const result = runSEIR(params);
const peak   = result.reduce((mx, r) => r.I > mx.I ? r : mx, result[0]);
console.log(`R₀ = ${R0.toFixed(2)} | herd immunity = ${(herd*100).toFixed(1)}%`);
console.log(`Peak I = ${(peak.I/params.N*100).toFixed(2)}% of N at day ${peak.t.toFixed(0)}`);

Plotting with Canvas 2D

function plotSEIR(canvas, results, N) {
  const ctx = canvas.getContext('2d');
  const { width: W, height: H } = canvas;
  ctx.clearRect(0, 0, W, H);

  const pad = { l: 50, r: 20, t: 20, b: 40 };
  const cw = W - pad.l - pad.r;
  const ch = H - pad.t - pad.b;

  const px = t  => pad.l + (t / results.at(-1).t) * cw;
  const py = v  => pad.t + ch - (v / N) * ch;

  const lines = [
    { key: 'S', color: '#60a5fa', label: 'Susceptible'  },
    { key: 'E', color: '#f59e0b', label: 'Exposed'      },
    { key: 'I', color: '#f97316', label: 'Infectious'   },
    { key: 'R', color: '#22c55e', label: 'Removed'      },
  ];

  for (const { key, color } of lines) {
    ctx.beginPath();
    ctx.strokeStyle = color;
    ctx.lineWidth   = 2;
    results.forEach((r, idx) => {
      idx === 0 ? ctx.moveTo(px(r.t), py(r[key]))
                 : ctx.lineTo(px(r.t), py(r[key]));
    });
    ctx.stroke();
  }
}

7. Extensions: vital dynamics, stochasticity, spatial spread

SIR with vital dynamics

Adding birth rate μ (equal to death rate for constant N) produces an endemic equilibrium — the pathogen persists indefinitely and oscillates. The endemic equilibrium is S* = N/R₀, I* = μN(1 − 1/R₀)/γ.

SIRS model — waning immunity

If immunity wanes at rate ξ, recovered individuals return to S. This is relevant for coronaviruses and influenza: R → S at rate ξ·R. The system can cycle, producing the seasonal waves observed in influenza.

Stochastic models

Deterministic models give average trajectories; real outbreaks involve randomness. Two approaches:

• Gillespie algorithm: exact stochastic simulation. Draw the time to next event from an exponential distribution; sample which event (infection, recovery) proportionally to its rate. Exact but slow for large N.
• Agent-based (IBM): simulate each individual, assign explicit contact network. Captures household structure, spatial heterogeneity and super-spreader events.

Spatial spread and metapopulation models

Real epidemics spread through a network of interconnected subpopulations (cities, regions). A metapopulation model applies SIR dynamics within each patch and adds inter-patch mixing via a mobility matrix — used for airline-borne pandemic modelling (GLEaM model, MIT/Harvard).