SIR models

How outbreaks grow, peak, and end

Last time

Lineages split (birth) and die (death) at random

Today: The same model, but the “lineages” are infections

  • An infected person passes the virus on: “birth” of an infection
  • An infected person recovers: “death” of an infection

The SIR model

Three groups of people1

  • \(N\) people: \(S + I + R = N\)
  • Each infected person recovers at rate \(\gamma\) (on average, after \(1 / \gamma\) days)
  • Each infected person passes on the virus at rate \(\beta\), but only a fraction, \(S / N\), of the people they meet are susceptible
    • So, transmission rate is \(\beta \frac{S}{N}\) per infected person

Compared to a birth-death model

Birth-death SIR
What multiplies? Lineages, \(n\) Infections, \(I\)
Total rate of “births” \(\lambda n\) \(\beta \frac{S}{N} \cdot I\)
Total rate of “deaths” \(\mu n\) \(\gamma I\)

At the start of an outbreak, \(S / N \approx 1\), so it is a birth-death process with \(\lambda = \beta\) and \(\mu = \gamma\)

The difference: Each infection “uses up” a susceptible person, so the outbreak slows down as \(S\) shrinks

The basic reproduction number

\(R_0\): The average number of people one infected person infects, when everyone else is susceptible

Infects at rate \(\beta\), for an average of \(1 / \gamma\) days:

\[ R_0 = \frac{\beta}{\gamma} \]

Question

What do you think happens when \(R_0 < 1\)? When \(R_0 > 1\)?

Simulating an outbreak one event at a time

The Gillespie algorithm, again1

  1. Calculate the total (population) rates: infection \(\beta I S / N\) and recovery \(\gamma I\)
  2. Draw an exponential waiting time using the total rate
  3. Choose which event happened, with probability proportional to its rate
    • Infection: \(S - 1\), \(I + 1\)
    • Recovery: \(I - 1\)
  4. Repeat until \(I = 0\)

Same recipe as simulate_birth_death, with different rates

Minor and major outbreaks

\(N = 763\), one infected person, \(\beta = 1\), \(\gamma = 0.5\) per day (\(R_0 = 2\))

How often does an outbreak fizzle out?

\(R_0\) Fraction of minor outbreaks (simulated)
1.5 0.65
2.0 0.52
3.0 0.34

Question

Last time, a clade starting with one lineage eventually died out with probability \(\mu / \lambda\).1 What should the probability be here?

The average outbreak

Euler’s method

In a short time \(\Delta t\), the expected numbers of events are:

\[ \text{new infections} = \frac{\beta I S}{N} \Delta t \qquad \text{recoveries} = \gamma I \Delta t \]

Update, and repeat, one time step at a time:

\[ \begin{aligned} S &\leftarrow S - \text{new infections} \\ I &\leftarrow I + \text{new infections} - \text{recoveries} \\ R &\leftarrow R + \text{recoveries} \end{aligned} \]

  • Solves the SIR differential equations, e.g., \(\frac{dI}{dt} = \frac{\beta I S}{N} - \gamma I\)
  • Smaller \(\Delta t\) → closer to the exact solution (but more steps)

Random outbreaks vs. the average

  • Major outbreaks have the same shape
    • Randomness early on mostly shifts when they take off

Why does an outbreak peak?

Each infected person infects at rate \(\beta S / N\) and recovers at rate \(\gamma\)

Question

What is \(S / N\) when \(I\) stops growing? (You’ll find out in the capstone!)

The boarding school

Influenza at a boarding school, 19781

  • A strain of H1N1 that hadn’t circulated since the 1950s: no one was immune
  • 512 of the school’s 763 students were confined to bed within two weeks

Fitting the model to the data

For each \(\beta\) and \(\gamma\), how far are the predictions from the counts?

\[ \text{sum of squares} = \sum_{\text{days}} \left(\text{observed} - \text{predicted}\right)^2 \]

  • Least squares: The best fit has the smallest sum of squares (like linregress in Capstone 2)
  • Grid search: Try every combination of \(\beta\) and \(\gamma\) from two lists, and keep the best

Watch for

How contagious was this virus? And why did only 512 students end up in bed?

Summary

  • The SIR model: infections at rate \(\beta I S / N\), recoveries at rate \(\gamma I\)
  • \(R_0 = \beta / \gamma\): the average number of people one case infects
  • Near its start, an outbreak is a birth-death process, so many outbreaks fizzle out by chance
  • Euler’s method calculates the average course of an outbreak, one small time step at a time
  • Fitting the model to data tells us how contagious a virus is

Capstone 3

  1. simulate_sir: The Gillespie algorithm for outbreaks
  2. fraction_minor_outbreaks: How often does an outbreak fizzle out?
  3. simulate_sir_euler and plot_sir: Explore the average outbreak, with sliders
  4. sum_of_squares and fit_sir: Fit the model to the boarding school
  5. Extra credit: How many students would the school have needed to vaccinate?

Capstone 3

  1. Log into molab
  2. Open Capstone 3 in molab
  3. Click “Fork” button

References

Anonymous. 1978. “Influenza in a Boarding School.” British Medical Journal 1 (6112): 587. https://www.ncbi.nlm.nih.gov/pmc/articles/PMC1603269/.
Gillespie, Daniel T. 1977. “Exact Stochastic Simulation of Coupled Chemical Reactions.” The Journal of Physical Chemistry 81 (25): 2340–61. https://doi.org/10.1021/j100540a008.
Kendall, David G. 1948. “On the Generalized ‘Birth-and-Death’ Process.” Annals of Mathematical Statistics 19 (1): 1–15. https://doi.org/10.1214/aoms/1177730285.
Kermack, William Ogilvy, and Anderson G. McKendrick. 1927. “A Contribution to the Mathematical Theory of Epidemics.” Proceedings of the Royal Society of London. Series A 115 (772): 700–721. https://doi.org/10.1098/rspa.1927.0118.

Acknowledgements

I drafted and iterated on these slides in collaboration with Anthropic’s Claude (Opus 5.5 model).