Adding Realism

How Epidemics Spread

Chapter 5 · Adding Realism

The model in How Epidemics Spread 3 had three boxes, and How Epidemics Spread 4 used it to find the herd immunity threshold. Real outbreaks break its rules in familiar ways. People are infected a few days before they can infect anyone else. Some die. Immunity can fade. Babies are born, and each new baby is a fresh susceptible. This chapter adds those four things one at a time and watches what each does. Two lessons emerge: some changes alter only the timing of an epidemic, and some change the long-run outcome, including the possibility that a disease never leaves and becomes endemic.

What was run, and what to trust

Every Python snippet in this chapter and in the three exercise solutions was executed, and the output in the comments is what Python printed. The model equations with births and deaths (R0 = β/(μ + γ), a disease-free state, and an endemic state with the susceptible share at (γ + μ)/β) come from the Wikipedia article on compartmental models in epidemiology, fetched with a tool that often returns truncated pages or model-written summaries; this fetch covered only the first 100,000 of about 229,000 characters and stopped before the SEIR section, so treat the wording as a paraphrase and the article as unchecked for SEIR and SIRS equations. I found no conflict between sources, but I only had the one. The SEIR and SIRS rules, the death rule, and the formulas for the endemic level with waning immunity or newborn vaccination are my own algebra on the standard set-up, checked here only against the simulations; they are inferences, not quotations. Every population size, R0, latent period, infectious period, fatality share, lifetime and waning time in the code is assumed or illustrative, not an estimate for any real disease or place.

A Latent Period: the E Box

After catching a disease, a person is usually not infectious straight away. The Wikipedia article notes that SEIS and SEIR models include a latent period during which the person is not infectious. The extra compartment is called E, for exposed, and the model is called SEIR: people move S to E at the infection rate (exactly as before), then E to I at a rate σ = 1 / (latent period), then I to R at γ. The code below runs R0 = 2.5 with a 5-day infectious period and latent periods of 0 (the plain SIR of earlier chapters), 2, 4 and 8 days, and reports the day of the peak, the peak number infectious, and the share ever infected. The step is 0.05 days, and this is just numerical stepping, no calculus.

# SEIR: adding a latent period. R0 = 2.5 and a 5-day infectious period are assumed.
population = 10_000
gamma = 0.2            # assumed: average infectious period 5 days
r0 = 2.5
beta = r0 * gamma
dt = 0.05

def run(latent_days):
    """latent_days = 0 means no latent compartment (plain SIR)."""
    S, E, I, R = population - 1.0, 0.0, 1.0, 0.0
    if latent_days > 0:
        sigma = 1 / latent_days
    peak_I, peak_t, t = I, 0.0, 0.0
    while (E + I) >= 0.01 and t < 2000:
        new_exposed = beta * S * I / population * dt
        if latent_days > 0:
            new_infectious = sigma * E * dt
        else:
            new_infectious = new_exposed          # skip the E stage
            E = 0.0
        new_rec = gamma * I * dt
        S -= new_exposed
        if latent_days > 0:
            E += new_exposed - new_infectious
        I += new_infectious - new_rec
        R += new_rec
        t += dt
        if I > peak_I:
            peak_I, peak_t = I, t
    return peak_t, peak_I, R / population

print("latent days   peak day   peak infectious   share ever infected")
for latent in [0, 2, 4, 8]:
    day, peak, share = run(latent)
    print(f"{latent:11d} {day:10.0f} {peak:17.0f} {share:21.1%}")
# latent days   peak day   peak infectious   share ever infected
#           0         32              2343                 89.3%
#           2         55              1653                 89.3%
#           4         72              1279                 89.3%
#           8        102               886                 89.3%
    

The final size is the same in all four runs: 89.3 per cent, the figure for R0 = 2.5 in How Epidemics Spread 3 and How Epidemics Spread 4. What changes is the schedule. A latent period of 8 days moves the peak from day 32 to day 102 and cuts the peak number infectious from 2,343 to 886. In these runs R0 is untouched (it is still β/γ, the number of people one case infects), but each step along the chain now takes longer, so the epidemic grows more slowly, spreads the same total over a longer time, and at any instant fewer people are infectious (many are waiting in E). That is my reading of the output, not a result from the article. Reading for practice: the same number of eventual cases can mean a very different load at the peak, which is why timing matters for hospitals.

	SIR (How Epidemics Spread 3)	SEIR (this chapter)
Compartments	S, I, R	S, E, I, R
Path of a case	infected and infectious at once	infected, wait (latent), then infectious
Peak day (R0 = 2.5)	day 32	day 55 / 72 / 102 for 2 / 4 / 8 days latent
Final share infected	89.3%	89.3% in each of those runs
Best for	diseases with little delay	diseases with an incubation gap before infectiousness

Adding Deaths

So far people left the infectious stage into R and the population stayed the same size. A simple way to add deaths is to say that a share f of the people leaving the infectious stage die rather than recover. The code does that, with the SEIR rules above (3-day latent, 5-day infectious). It also divides by the number still alive when working out how often a susceptible meets an infectious person. That is a modelling choice, assumed. The fatality shares 0.1, 1 and 5 per cent are illustrative.

# Deaths: a share f of the people who leave the infectious stage die instead of recovering.
# R0 = 2.5, a 3-day latent and 5-day infectious period, and every f below are assumed.
population = 10_000
gamma, sigma, r0, dt = 0.2, 1 / 3, 2.5, 0.05
beta = r0 * gamma

def run(f):
    S, E, I, R, D = population - 1.0, 0.0, 1.0, 0.0, 0.0
    t = 0.0
    while (E + I) >= 0.01 and t < 2000:
        alive = population - D              # the dead no longer take part in mixing
        new_exposed = beta * S * I / alive * dt
        new_infectious = sigma * E * dt
        leaving = gamma * I * dt
        S -= new_exposed
        E += new_exposed - new_infectious
        I += new_infectious - leaving
        R += (1 - f) * leaving
        D += f * leaving
        t += dt
    return D, (R + D) / population

print("fatality share f   deaths per 10,000   share ever infected")
for f in [0.001, 0.01, 0.05]:
    deaths, ever = run(f)
    print(f"{f:15.1%} {deaths:19.1f} {ever:21.1%}")
# fatality share f   deaths per 10,000   share ever infected
#            0.1%                 8.9                 89.3%
#            1.0%                89.5                 89.5%
#            5.0%               450.2                 90.0%
    

Deaths scale almost directly with f: 8.9, 89.5 and 450.2 per 10,000. The share ever infected barely moves (89.3, 89.5 and 90.0 per cent). The 90.0 comes from my choice to divide by the survivors: as people die, each remaining susceptible meets a slightly larger share of infectious people, so a few more get infected. The point is that this kind of death changes the toll but hardly the course of the epidemic. It is also a different thing from a deadly disease that kills people before they can infect others, which would shorten the infectious period and so lower R0, an effect not modelled here.

Births, Waning Immunity and the Endemic State

The big change comes from a longer horizon. Over decades, people are born and die for other reasons too. The article's version of the model has a birth rate and a death rate; with them, R0 becomes β/(μ + γ), because people leave the infectious state by dying as well as recovering. It gives two long-run outcomes: if R0 is 1 or less, the disease dies out and the population returns to the all-susceptible, disease-free state; if R0 is above 1 and some infection is present, the system tends to an endemic state in which the disease stays permanently, with the susceptible share settling at (γ + μ)/β, which is 1/R0. The disease stays because births keep refilling S.

The code compares two worlds with R0 = 2.5, a population of 1,000,000 and people living on average 70 years (birth rate equal to death rate, so the population stays constant). In the first, immunity lasts for life, so only births refill S. In the second, immunity wanes: recovered people return to S after 2 years on average, giving the SIRS model (mentioned briefly in the article, with no equations in the portion I saw). It prints, for each 25-year window, the lowest, highest and mean number infectious per 100,000, plus the endemic level my algebra predicts. Setting the rates of change to zero and solving is just algebra: for births only it gives μ(R0 − 1)/β, the article’s endemic infectious share; with waning ω, it gives (1 − 1/R0) / (1 + γ/(μ + ω)).

# Births and waning immunity: SIR with births vs SIRS with births (all values assumed).
population = 1_000_000
gamma, r0, dt = 0.2, 2.5, 0.25
mu = 1 / (70 * 365)                  # birth rate = death rate; mean lifetime 70 years
beta = r0 * (gamma + mu)             # so that beta / (gamma + mu) = 2.5

def run(omega, years=150):
    S, I, R = 0.9 * population, 100.0, 0.1 * population
    steps_per_year = int(365 / dt)
    windows = {}                      # 25-year windows: lowest, highest, mean infectious
    for year in range(years):
        for _ in range(steps_per_year):
            infect = beta * S * I / population
            S += (mu * population - mu * S - infect + omega * R) * dt
            I += (infect - (gamma + mu) * I) * dt
            R += (gamma * I - (mu + omega) * R) * dt
            window = year // 25
            lo, hi, total, n = windows.get(window, (1e18, 0.0, 0.0, 0))
            windows[window] = (min(lo, I), max(hi, I), total + I, n + 1)
    return windows, S / population

for label, omega in [("births only (lifelong immunity)", 0.0), ("waning, mean 2 years", 1 / (2 * 365))]:
    print(label)
    windows, s_end = run(omega)
    for w, (lo, hi, total, n) in windows.items():
        print(f"  years {25 * w:3d}-{25 * w + 24:3d}: infectious per 100,000: lowest {lo / 10:9.3f}, highest {hi / 10:9.3f}, mean {total / n / 10:7.1f}")
    print(f"  susceptible share at the end: {s_end:.3f}   (1/R0 = {1 / r0:.3f})")
    if omega == 0:
        formula = mu * (r0 - 1) / beta           # (mu / beta)(R0 - 1)
    else:
        formula = (1 - 1 / r0) / (1 + gamma / (mu + omega))
    print(f"  formula for the endemic level: {formula * 100_000:.1f} per 100,000")
# births only (lifelong immunity)
#   years   0- 24: infectious per 100,000: lowest     0.000, highest 17832.605, mean    42.3
#   years  25- 49: infectious per 100,000: lowest     0.000, highest     0.000, mean     0.0
#   years  50- 74: infectious per 100,000: lowest     0.000, highest  4139.665, mean    20.1
#   years  75- 99: infectious per 100,000: lowest     0.000, highest  1842.713, mean    13.4
#   years 100-124: infectious per 100,000: lowest     0.000, highest  1041.567, mean    10.1
#   years 125-149: infectious per 100,000: lowest     0.000, highest   669.340, mean    14.9
#   susceptible share at the end: 0.346   (1/R0 = 0.400)
#   formula for the endemic level: 11.7 per 100,000
# waning, mean 2 years
#   years   0- 24: infectious per 100,000: lowest     0.000, highest 18157.649, mean   446.2
#   years  25- 49: infectious per 100,000: lowest   419.843, highest   419.886, mean   419.9
#   years  50- 74: infectious per 100,000: lowest   419.813, highest   419.843, mean   419.8
#   years  75- 99: infectious per 100,000: lowest   419.792, highest   419.813, mean   419.8
#   years 100-124: infectious per 100,000: lowest   419.778, highest   419.792, mean   419.8
#   years 125-149: infectious per 100,000: lowest   419.767, highest   419.778, mean   419.8
#   susceptible share at the end: 0.400   (1/R0 = 0.400)
#   formula for the endemic level: 419.7 per 100,000
    

The two worlds behave very differently. With waning immunity the system settles within 25 years at 419.8 infectious per 100,000, a steady background level that matches the formula (419.7) and keeps the susceptible share at exactly 1/R0 = 0.400. With lifelong immunity the disease does not settle quickly. The first wave peaks at 17,832.6 per 100,000, and then infection fades to a lowest and highest of under 0.0005 per 100,000 in the years 25 to 49, which in a population of one million is far less than a single person. Only as births rebuild S does it rebound (to 4,139.7, then 1,842.7, 1,041.6 and 669.3 per 100,000 in later windows). The waves shrink and slowly home in on the formula’s 11.7 per 100,000 (the final-window mean is 14.9 and the susceptible share has reached 0.346, not yet 0.400).

The deterministic model keeps a vanishing fraction of a person infected and lets the disease recover from it. A real population has whole people. My inference (not a result from the article) is that those deep troughs would end in real extinction, with the disease coming back only if someone brought it in from outside. This is where chance matters, and How Epidemics Spread 6 is about it.

	Epidemic (earlier chapters)	Endemic (this chapter)
Time scale	weeks to months	years to decades
Susceptible supply	only the starting population	refilled by births and by waning immunity
Long-run susceptible share	falls and stays low	settles near 1/R0 (0.400 for R0 = 2.5)
Shape	one wave, then it ends	waves that shrink, or a steady background level
Effective reproduction number	falls below 1 as S falls	hovers around 1 (R0 × S/N = 1 at the endemic state)

Why Vaccination Looks Different in an Endemic World

How Epidemics Spread 4 vaccinated people before one outbreak. When births keep arriving, vaccination has to be an ongoing programme, and the same threshold appears in a new form: if a share p of newborns is protected, the formula for the endemic level becomes μ(R0(1 − p) − 1)/β, which falls to zero when p reaches 1 − 1/R0. This is my own derivation, checked by Exercise 3. Notice how a plain 1 − 1/R0 still governs elimination, but here it applies to each new generation of susceptibles.

Common mistakes

Believing a latent period changes the final size (in these runs it only changed the timing); counting people in E as infectious (they are not, and the peak reported is of I only); forgetting that R0 changes to β/(μ + γ) once deaths or births are in the model; reading an endemic level of 11.7 per 100,000 as a promise that a real town always has cases (very small populations die out by chance); mixing up a low endemic level with a harmless disease; and using one set of assumed numbers as if it were data (every rate here is illustrative).

Exercise 1

Run the SEIR model with R0 = 3 (assumed) and latent periods of 2, 4 and 8 days with a 5-day infectious period, population 10,000. Print the day of the peak, the peak number infectious and the share ever infected for each, and say which of the three numbers changes with the latent period and which does not.

📄 View solution

Exercise 2

Using the SIRS model with births in a population of 1,000,000 (R0 = 2.5, 70-year lifetime, assumed), let immunity last 0.5, 1 and 3 years. Simulate 60 years and compare the infectious level per 100,000 at the end with the formula (1 − 1/R0) / (1 + γ/(μ + ω)).

📄 View solution

Exercise 3

Start from the unvaccinated endemic state (R0 = 2.5, lifelong immunity, assumed). Switch on a perfect vaccine for 0, 30, 50, 60 and 70 per cent of newborns, run 100 years, and compare the final infectious level per 100,000 with the prediction μ(R0(1 − p) − 1)/β. Which coverages eliminate the disease?

📄 View solution

Quick Reference

SEIR: S to E at the infection rate, E to I at σ = 1 / latent period, I to R at γ
Latent period effect: later and lower peak; final size unchanged in these runs (R0 = 2.5: 89.3%)
Deaths: a share f leaves I by dying; the toll scales with f, the epidemic barely changes
With births and deaths: R0 = β/(μ + γ); disease dies out if R0 ≤ 1, otherwise tends to an endemic state
Endemic state: susceptible share (γ + μ)/β = 1/R0; infectious share μ(R0 − 1)/β with lifelong immunity
SIRS (waning ω): infectious share (1 − 1/R0) / (1 + γ/(μ + ω))
Newborn vaccination p: level μ(R0(1 − p) − 1)/β, zero once p ≥ 1 − 1/R0
Assumes: well-mixed population, constant rates, illustrative numbers

Coming next

How Epidemics Spread 6 is about chance and small outbreaks: stochastic simulation, branching processes, super-spreaders, why outbreaks sometimes fizzle.