Solving the SIR model from scratch

The model

Split a population into the susceptible $S$ , the infectious $I$ and the recovered $R$ . Infections happen when susceptible and infectious people meet, at a rate $\beta$ per day, and the infectious recover at a rate $\gamma$ . Working in shares of the population (so $S + I + R = 1$ ) gives

\begin{aligned} \frac{dS}{dt} &= -\beta S I \\[1.2ex] \frac{dI}{dt} &= \beta S I - \gamma I \\[1.2ex] \frac{dR}{dt} &= \gamma I \end{aligned}

The second equation already holds the main result. $I$ grows only while $\beta S > \gamma$ , that is while $R_0 S > 1$ with the basic reproduction number $R_0 = \beta / \gamma$ . Every infection uses up susceptibles, so $S$ falls until the condition fails and the curve turns over.

Try it: an outbreak in 160 daysSIR

Transmission rate β

Recovery rate γ

R₀ = β / γ–

Peak infected–

Peak day–

Ever infected–

Computed in your browser with RK4, step 0.25 days, starting from 0.1% infected.

Integrating it

There is no closed-form solution for $I(t)$ , so we integrate numerically. SciPy’s solve_ivp does it in a few lines. The figure above uses a hand-written fourth-order Runge–Kutta step instead, which is the same idea in about fifteen lines of JavaScript.

1
import numpy as np
2
from scipy.integrate import solve_ivp
3

4
def sir(t, y, beta, gamma):
5
    s, i, r = y                     # shares of the population
6
    return [-beta * s * i, beta * s * i - gamma * i, gamma * i]
7

8
beta, gamma = 0.32, 0.10            # per day
9
y0 = [0.999, 0.001, 0.0]
10
t = np.linspace(0, 160, 641)
11

12
sol = solve_ivp(sir, (0, 160), y0, args=(beta, gamma), t_eval=t)
13
s, i, r = sol.y
14
print(f"R0 = {beta / gamma:.2f}, peak on day {t[i.argmax()]:.0f}")

How big does it get?

Divide the first equation by the third and integrate. With almost everyone susceptible at the start, the share that is still susceptible at the end, $s_\infty$ , solves

s_\infty = \exp\!\big(-R_0\,(1 - s_\infty)\big)

For $R_0 = 3.2$ that leaves about 5% untouched, so roughly 95% are infected at some point. Drag $\gamma$ up until $R_0$ drops below 1 and the outbreak never takes off.