The SIR Model

How Epidemics Spread

Chapter 3 · The SIR Model

How Epidemics Spread 2 ended with a crude simulation: new cases were R0 times the fraction still susceptible times the cases one generation earlier. It showed the right shape, a climb, a peak and a fall, but it had no real notion of who was ill now, and it ignored how long people stay infectious. This chapter replaces it with the model that nearly every later model grows from, the SIR model. Everyone is in exactly one of three groups: Susceptible (can catch it), Infected, meaning currently infectious (can pass it on), or Recovered (no longer infectious and, in this model, immune, or removed). The model is a pair of rules for how many people move from S to I and from I to R each day. The textbook calls such rules differential equations, but you will not need calculus: we simply apply the rules over small steps of time, one step after another, in Python. That is called stepping the equations (the method used here is Euler’s method).

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’s equations, the meaning of the two rate constants, the relation R0 = beta/gamma, the 1927 origin (Kermack and McKendrick), the average infectious period 1/gamma and the final-size equation come from the Wikipedia article on compartmental models in epidemiology, fetched with a tool that often returns truncated pages or model-written summaries; the page was far longer than the part the tool read, and the summary noted that the article writes the equations with and without dividing by the population in different sections, so treat the wording as a paraphrase. The statement that infections peak when S falls to N/R0 is not quoted in the article: it is worked out below from the equations. All rates, populations and R0 values in the code are assumed or illustrative; the boarding-school fit in the third exercise is a classroom fit, not a published estimate.

Three Boxes and Two Arrows

Picture three boxes in a row: S, then I, then R, with an arrow from S to I and another from I to R. People only move forwards. The population is closed: nobody is born, dies or arrives (assumed, and relaxed in How Epidemics Spread 5), so S + I + R stays equal to the total N all the time. In How Epidemics Spread 1 the boarding school counted boys “in bed”, and a rough match is that the in-bed count is the I box and the convalescent boys are heading for R (an assumption, since being in bed is not quite the same as being infectious).

BoxWho is in itHow people enterHow people leave
S susceptiblecan catch the infectionstart here (nobody enters later)get infected, move to I
I infectedinfectious nowinfected by someone in Irecover (or are removed), move to R
R recoveredno longer infectious, assumed immunerecovery from Inobody leaves

The Two Rules

Two numbers control the flows. The infection rate constant beta (written β) says how fast infectious people infect susceptible ones: the Wikipedia article gives it in units of people infected per day per infectious person. The recovery rate constant gamma (γ) is the fraction of the infectious who recover each day. With the population N, the rules for a small step of time are:

  • New infections in the step = β × S × I / N × (step length). The product S × I / N is large when there are many people who can pass it on and many who can catch it; it shrinks as S falls or as I falls.
  • New recoveries in the step = γ × I × (step length).
  • S falls by the new infections; I gains the new infections and loses the new recoveries; R gains the new recoveries.

If the infectious recover at the rate γ = 0.2 per day, then on average a person stays infectious for 1/γ = 5 days. (The article states that the mean infectious period is 1/γ; this holds when the time to recover varies in the particular way built into the model, which is a simplification of real diseases.) Compare this with How Epidemics Spread 2:

Generation model (How Epidemics Spread 2)SIR model
Timejumps of one generation (a fixed number of days)small steps, such as a day or a tenth of a day
Who is infectiousonly the latest generationthe whole box I, with people leaving it at rate γ
What it needsR0 and a generation timeβ and γ (so R0 and the infectious period)
Where it stopsa rough limit on the susceptible fractioninfections fall as S and I both shrink
Originthis course’s simplificationKermack and McKendrick, 1927 (per the article)

Stepping the Model in Python

The code below follows the rules exactly, one day at a time, for a closed community of 10,000 with a single infectious person at the start. The values β = 0.5 and γ = 0.2 are assumed for illustration; they give R0 = 2.5. Each loop turns the rules into three lines of arithmetic: the new infections, the new recoveries, and the three updates. It also records the largest I along the way.

# SIR stepped one day at a time. All numbers are assumed, for illustration. population = 10_000 beta = 0.5 # assumed: contacts that could transmit, per infectious person per day gamma = 0.2 # assumed: 1/5 of the infectious recover each day (average 5 days) dt = 1.0 # step length in days S, I, R = population - 1.0, 1.0, 0.0 peak_I, peak_day, S_at_peak = I, 0, S print("day S I R") for day in range(0, 121): if day % 10 == 0: print(f"{day:3d} {S:8.1f} {I:8.1f} {R:8.1f}") if I > peak_I: peak_I, peak_day, S_at_peak = I, day, S new_infections = beta * S * I / population * dt new_recoveries = gamma * I * dt S -= new_infections I += new_infections - new_recoveries R += new_recoveries print(f"R0 = beta/gamma = {beta / gamma}") print(f"peak: {peak_I:.0f} infectious on day {peak_day}, when S was {S_at_peak:.0f}") print(f"S = N/R0 predicts the peak at S = {population / (beta / gamma):.0f}") print(f"total never infected at the end: {S:.0f}; S+I+R = {S + I + R:.1f}") # day S I R # 0 9999.0 1.0 0.0 # 10 9977.7 13.8 8.5 # 20 9693.4 182.7 123.9 # 30 6983.1 1610.1 1406.8 # 40 2161.6 2009.9 5828.5 # 50 1129.7 539.4 8330.9 # 60 974.5 108.6 8916.9 # 70 946.7 20.9 9032.4 # 80 941.5 4.0 9054.6 # 90 940.5 0.8 9058.8 # 100 940.3 0.1 9059.6 # 110 940.3 0.0 9059.7 # 120 940.3 0.0 9059.7 # R0 = beta/gamma = 2.5 # peak: 2487 infectious on day 36, when S was 3570 # S = N/R0 predicts the peak at S = 4000 # total never infected at the end: 940; S+I+R = 10000.0

Read the table like a story. For ten days almost nothing seems to happen: 13.8 people are infectious on day 10, then 183 on day 20 and 1,610 on day 30. That quiet start is the exponential growth of How Epidemics Spread 2, and it is why outbreaks often look harmless until they are not. The infectious number peaks at 2,487 on day 36; by day 60 only about 109 remain, and the epidemic ends with 940 people never infected. Two checks are built in. First, S + I + R is 10,000 to the last decimal (a sign the arithmetic moves people between boxes without losing any). Second, the peak line and the prediction line compare the S at the peak with N/R0 = 4,000. They do not agree exactly (3,570 against 4,000), and the next sections explain both the prediction and the gap.

Why R0 = β / γ

Take one infectious person in a population where everyone else is susceptible, so S/N is about 1. Each day they pass the infection to β people. They stay infectious for 1/γ days on average. So over their whole infectious period they infect β × (1/γ) = β/γ people. That is exactly R0, the expected number of cases one case generates in a fully susceptible population: the same definition as in How Epidemics Spread 2, now built from two things you can measure separately, how infectious the disease is per day and how long it lasts. The same reasoning shows how the effective reproduction number falls: when only a fraction S/N of the population can still be infected, each case infects R0 × S/N people.

That gives the peak rule. Rewrite the change in I per day as (R0 × S/N − 1) × γ × I. The number of infectious people keeps growing while R0 × S/N is above 1 and starts to shrink when it drops below 1, so the peak is at S = N/R0. For R0 = 2.5 and N = 10,000 that is S = 4,000. The article states this threshold only as a herd-immunity threshold; reading it as the peak condition is this course’s own inference from the equation. It is the seed of How Epidemics Spread 4. Notice what the peak is not: it is not the end. After I peaks, I still falls only gradually while it keeps infecting people, so S continues to drop below 4,000. In the run above, S was 2,162 by day 40 and ended at 940.

How Small Must the Step Be?

The run above used a one-day step, and the peak S came out at 3,570 rather than 4,000. Stepping the rules is an approximation: during a whole day the rule assumes the rate of infection stays fixed even though S and I change as the day goes on. Shorter steps follow the real curve more closely. The code below runs the same model with steps from 5 days down to 0.01 days (1/100 of a day). It also stops the day’s new infections from exceeding S, a guard that matters with very large steps.

# How much does the step length matter? Same model, smaller and smaller steps. population, beta, gamma = 10_000, 0.5, 0.2 # assumed values def run(dt, days=150): S, I, R = population - 1.0, 1.0, 0.0 peak_I, peak_t = I, 0.0 for k in range(int(days / dt)): new_inf = min(S, beta * S * I / population * dt) new_rec = gamma * I * dt S -= new_inf I += new_inf - new_rec R += new_rec if I > peak_I: peak_I, peak_t = I, (k + 1) * dt return peak_I, peak_t, S print("step (days) peak infectious peak time (days) never infected") for dt in [5.0, 2.0, 1.0, 0.5, 0.1, 0.01]: peak_I, peak_t, S = run(dt) print(f"{dt:10.2f} {peak_I:16.0f} {peak_t:17.1f} {S:15.0f}") # step (days) peak infectious peak time (days) never infected # 5.00 3272 50.0 206 # 2.00 2636 40.0 791 # 1.00 2487 36.0 940 # 0.50 2411 34.0 1009 # 0.10 2350 32.4 1061 # 0.01 2337 32.1 1072

The 5-day step is clearly wrong: it reports a peak of 3,272 on day 50 when the answer is closer to 2,337 on day 32. The 1-day step used above overshoots the peak by about 6% (2,487 against 2,337) and by 4 days. Halving to 0.5 days and then 0.1 days brings the answers close to the 0.01 results (peak around 2,340 on day 32, about 1,070 never infected), and the peaks in the last two rows differ by less than 1%. The lesson is practical: when a result barely changes as the step is halved, it is converged; when it keeps changing, the step is too coarse. Always repeat a simulation with a smaller step before you trust a number. (The simple method used here has a known weakness: it needs a small step. Better methods exist, and calculus-based treatments use them, but this one is enough for everything in this course.)

Where Does the Epidemic End?

Look again at the first table: 940 people were never infected. An SIR epidemic does not stop because the susceptible run out. It stops because, as S falls, each case infects fewer than one other (R0 × S/N drops below 1), so I declines to nothing while some people remain in S. The article gives a final-size equation for the share who escape, and the code below uses a rearranged version: if z is the share of the population ever infected (starting with almost nobody infected), then z = 1 − e−R0×z, which you can solve with no calculus by guessing z, plugging it into the right-hand side, and repeating until z stops changing. The code compares that answer with the full step-by-step simulation, using a step of 0.05 days. (The equation as printed in the article uses its own notation; the rearrangement here is mine, for the case of almost everyone starting susceptible.)

# Final size: simulate, then check against a number found by repeated guessing. import math population, gamma, dt = 10_000, 0.2, 0.05 # assumed values def simulate(r0): beta = r0 * gamma S, I, R = population - 1.0, 1.0, 0.0 while I >= 0.01: new_inf = beta * S * I / population * dt new_rec = gamma * I * dt S -= new_inf I += new_inf - new_rec R += new_rec return (population - S) / population def final_size_guess(r0): z = 0.5 # start with any guess for the share ever infected for _ in range(200): z = 1 - math.exp(-r0 * z) # a self-consistency equation; repeat until steady return z print(" R0 simulated share infected self-consistency answer") for r0 in [0.8, 1.0, 1.2, 1.5, 2.0, 2.5, 3.0, 5.0]: print(f"{r0:5.1f} {simulate(r0):14.1%} {final_size_guess(r0):24.1%}") # R0 simulated share infected self-consistency answer # 0.8 0.0% 0.0% # 1.0 1.4% 1.0% # 1.2 31.4% 31.4% # 1.5 58.3% 58.3% # 2.0 79.7% 79.7% # 2.5 89.3% 89.3% # 3.0 94.1% 94.0% # 5.0 99.3% 99.3%

The two columns agree to within about a tenth of a percentage point for every R0 of 1.2 and above (for example 58.3% at R0 = 1.5 by both methods), which is a good check that both are right. Three points stand out. First, R0 = 1.2 infects 31% in the end, and R0 = 2 infects 80%: the share rises steeply near 1 and then flattens, so the extra infected from doubling R0 from 2.5 to 5 is only about 10 percentage points. Second, at R0 = 0.8 the final number is effectively zero: with R0 below 1 there is no epidemic at all, just a few cases. Third, at R0 exactly 1 the simulation reaches 1.4% while the repeated-guess method shows 1.0% after only 200 repeats: convergence is very slow at R0 = 1, so that row mostly shows the limits of the quick method, not a real difference. At R0 = 1 the true answer for a large population is a tiny share.

Early Growth, Revisited

How Epidemics Spread 2 turned R0 into a doubling time using a fixed generation time. The SIR model gives a different answer with no extra assumption. At the start almost everyone is susceptible (S/N is about 1), so the rule for I becomes: change per day = (β − γ) × I = γ × (R0 − 1) × I. That is exponential growth at the rate γ × (R0 − 1) per day, and so a doubling time of ln(2) / (γ × (R0 − 1)). The code checks this on a large population (so S really is near N) and compares it with the fixed-generation idea.

# Early growth in the SIR model compared with the generation-time idea. import math population, gamma, dt = 1_000_000, 0.2, 0.01 # assumed values r0 = 2.5 beta = r0 * gamma S, I = population - 10.0, 10.0 t, history = 0.0, [] while t < 10: history.append((t, I)) new_inf = beta * S * I / population * dt S -= new_inf I += new_inf - gamma * I * dt t += dt t1, i1 = history[int(2 / dt)] t2, i2 = history[int(8 / dt)] measured = (t2 - t1) * math.log(2) / math.log(i2 / i1) rate = gamma * (r0 - 1) print(f"measured doubling time, days 2 to 8: {measured:.2f} days") print(f"growth rate gamma*(R0-1) = {rate:.2f} per day -> doubling {math.log(2) / rate:.2f} days") print(f"fixed 5-day generation idea (Chapter 2 style): {5 * math.log(2) / math.log(r0):.2f} days") # measured doubling time, days 2 to 8: 2.31 days # growth rate gamma*(R0-1) = 0.30 per day -> doubling 2.31 days # fixed 5-day generation idea (Chapter 2 style): 3.78 days

The measured doubling time, 2.31 days, matches the formula exactly. The fixed 5-day generation of How Epidemics Spread 2 would give 3.78 days for the same R0 and the same average infectious period. The two models make different assumptions about timing: in the SIR model, people stop being infectious at random moments, and many pass the infection on early in their infectious period. So the same R0 and the same mean infectious time can give a faster or slower epidemic depending on how the timing is spread out, which is the point How Epidemics Spread 8 returns to when it covers serial intervals. Do not treat either number as the truth for a real disease.

Common mistakes
Using β as if it were R0 (R0 is β divided by γ); forgetting to divide by N in the infection rule, which makes the epidemic depend on the population size; updating S, I and R one after another using already-changed values instead of working out both flows first; trusting a result from one step size without halving the step; reading “the epidemic peaks when S = N/R0” as “the epidemic ends when S = N/R0”; expecting everyone to be infected in the end (a share always escapes unless R0 is huge); and forgetting that the model assumes a closed population, equal mixing, and permanent immunity, none of which is exactly true of a real disease.
Exercise 1

Keep R0 = 2.5 in a population of 10,000 with one starting case, but give the disease infectious periods of 2, 5 and 10 days (assumed). Work out β and γ for each, run the model, and report the peak day, the peak number infectious and the number never infected. What changes and what stays the same?

📄 View solution
Exercise 2

With a strongly infectious disease (assumed β = 2.0, γ = 0.2, so R0 = 10), run the model for 60 days with steps of 1, 0.5, 0.1 and 0.01 days, once with and once without the guard that stops new infections exceeding S. Report the lowest S reached in each run, and explain what goes wrong with a step that is too large.

📄 View solution
Exercise 3

Use the 14 daily in-bed counts from the boarding school in How Epidemics Spread 1 (763 boys, 3 infectious on day 1, assumed to match the I box). Search a grid of β and γ values for the pair that gives the smallest sum of squared differences from the data, then print the observed and model numbers side by side, the fitted R0, and the share of the school the model infects against the 512 recorded cases.

📄 View solution

Quick Reference

  • Boxes: S (susceptible), I (infectious), R (recovered or removed); S + I + R = N in a closed population
  • New infections per step: β × S × I / N × step
  • New recoveries per step: γ × I × step
  • R0 = β / γ; the average infectious period is 1/γ
  • Peak of I: when R0 × S/N = 1, that is S = N/R0 (worked out from the equation)
  • Early growth rate: γ × (R0 − 1) per day
  • Final size: z = 1 − e−R0×z for the share ever infected (almost everyone starting susceptible); never 100%
  • Check: halve the step; if the answer moves, the step was too big
Coming next
How Epidemics Spread 4 is about herd immunity: the threshold 1 - 1/R0, vaccination, imperfect vaccines, and what coverage is needed.