Epidemic modeling: the SIR model

How does a disease spread through a population — and what makes it explode into an epidemic versus fizzle out? The SIR model is the classic first answer. It splits everyone into three groups:

S — Susceptible (can still catch it)
I — Infected (have it now, and can pass it on)
R — Recovered (had it, now immune)

Two numbers drive everything: how fast infected people pass it on, and how fast they recover. Their ratio is the famous basic reproduction number R0 — the average number of people one infected person infects in a fully susceptible population.

Turn the sliders below and watch the epidemic curve respond live, computed in your browser as WebAssembly — no Julia server.

begin
    using PlutoUI, WasmMakie
end

The dynamics

Each day, a fraction of susceptible-meets-infected pairings produces new infections, and a fraction of the infected recover. Writing S, I, R as fractions of the population:

$$\Delta S = -\beta\, S\, I, \qquad \Delta I = \beta\, S\, I - \gamma\, I, \qquad \Delta R = \gamma\, I.$$

Here $\beta$ is the infection rate and $\gamma$ the recovery rate, with $R_0 = \beta / \gamma$. We step this forward one day at a time (Euler's method) — the same loop the original lecture runs, here compiled to WebAssembly.

R0 (basic reproduction number) = 25 ÷10

infectious period (days) = 7

initial infected (%) = 1

days to simulate = 160

let
    # derive the rates from the (integer) sliders inside this cell, then step the
    # SIR equations forward one day at a time. Self-contained so every value the
    # wasm island needs is computed right here.
    gamma = 1.0 / Float64(infdays)
    r0 = Float64(r0x10) / 10.0
    beta = r0 * gamma
    i0 = Float64(i0pct) / 100.0

    S = Vector{Float64}(undef, ndays + 1)
    I = Vector{Float64}(undef, ndays + 1)
    R = Vector{Float64}(undef, ndays + 1)
    days_axis = Vector{Float64}(undef, ndays + 1)
    s = 1.0 - i0
    inf = i0
    rec = 0.0
    S[1] = s
    I[1] = inf
    R[1] = rec
    days_axis[1] = 0.0
    for d in 1:ndays
        newinf = beta * s * inf
        newrec = gamma * inf
        s = s - newinf
        inf = inf + newinf - newrec
        rec = rec + newrec
        S[d + 1] = s
        I[d + 1] = inf
        R[d + 1] = rec
        days_axis[d + 1] = Float64(d)
    end

    fig = Figure(size = (600, 360))
    ax = Axis(fig[1, 1])
    lines!(ax, days_axis, S)   # susceptible (falls)
    lines!(ax, days_axis, I)   # infected (the epidemic curve)
    lines!(ax, days_axis, R)   # recovered (rises)
    fig
end

sir_stats = let
    # rerun the simulation in a SEPARATE bond-dependent cell, tracking the peak, and
    # return the readout numbers as a tuple (the markdown cell below interpolates them;
    # values computed inside a markdown cell's own `let` bake to the slider defaults).
    gamma = 1.0 / Float64(infdays)
    r0 = Float64(r0x10) / 10.0
    beta = r0 * gamma
    i0 = Float64(i0pct) / 100.0
    s = 1.0 - i0
    inf = i0
    rec = 0.0
    peak = inf
    peakday = 0
    for d in 1:ndays
        newinf = beta * s * inf
        newrec = gamma * inf
        s = s - newinf
        inf = inf + newinf - newrec
        rec = rec + newrec
        if inf > peak
            peak = inf
            peakday = d
        end
    end
    (floor(r0 * 100.0) / 100.0, floor(peak * 1000.0) / 10.0, peakday, floor(rec * 1000.0) / 10.0)
end;

R0 = 2.5. Peak infection reaches about 24.8% of the population at once, around day 24, and by the end 90.3% have been infected at some point. Push R0 below 1 (low infectiousness, fast recovery) and the outbreak never gets going.

R0 and the epidemic threshold

The single most important quantity is R0:

R0 > 1 — each case causes more than one new case on average → the infected curve climbs, peaks, then falls as the supply of susceptibles runs out. An epidemic.
R0 < 1 — each case causes fewer than one → the disease dies out almost immediately.

Notice the infected curve peaks and then declines even though no one intervened: the epidemic runs out of susceptible people to infect. The final recovered fraction — how many got infected in total — is the model's grim headline number, and lowering R0 (distancing, vaccination, faster recovery) is what bends that curve down.

Appendix

The MIT lecture integrates the continuous SIR differential equations with DifferentialEquations.jl and plots with Plots.jl. WebAssembly can't run that stack in the browser, so here we take simple one-day Euler steps by hand and draw with WasmMakie. For the gentle curves of an epidemic this discrete version tracks the continuous one closely, and every concept — beta, gamma, R0, the peak, herd immunity — is identical.