SIR Model with Vaccination

An SIR epidemic at the threshold: vaccinating past the herd-immunity line pushes effective R0 below 1 and the outbreak fizzles.

Level: Intermediate

populationreinforcing-loopbalancing-loopcontrol

  • Stocks: susceptible, infected, recovered
  • Flows: infections, recoveries
  • Feedback Loops: disease spread (reinforcing), susceptible depletion / herd immunity (balancing)
  • Probes: susceptible, infected, recovered, r_eff

Leverage Points

Discover how small, well-focused actions can produce significant, lasting improvements in complex systems.

Take the Quiz
simulation.py

SIR epidemics and the herd-immunity threshold

A disease moves through a population: susceptible people catch it, infected people recover and stay immune. Whether it erupts into an outbreak or fizzles out turns on a single number — the effective reproduction number R_eff = R₀ · (fraction still susceptible). Above 1, each case sparks more than one new case and infections climb; below 1 the chain of transmission dies out.

Here R₀ = β/γ is pinned just above 1, so with no protection it is a real (if slow) outbreak. The one knob is how much of the population starts immune from vaccination. Cross the herd-immunity threshold 1 − 1/R₀ and R_eff begins below 1, so infections fall from day one. Default: no vaccine (the curve peaks). Contrast: enough coverage to fizzle.


from tys import probe, progress

Simulate an SIR outbreak with a fraction of the population pre-vaccinated.

def simulate(cfg: dict):
    import simpy
    env = simpy.Environment()

Parameters

    population = cfg["population"]
    infected = cfg["initial_infected"]
    beta = cfg["beta"]                                   # infections per contact
    gamma = cfg["gamma"]                                 # recovery rate / day
    vacc_fraction = cfg.get("vaccinated_fraction", 0.0)  # the knob: share immune at t=0
    days = cfg["sim_days"]

A fraction of the population starts immune (pre-epidemic vaccination), which shrinks the susceptible pool the disease can spread through.

    vaccinated = vacc_fraction * population
    susceptible = population - infected - vaccinated
    recovered = vaccinated

    basic_r0 = beta / gamma
    r_eff_start = basic_r0 * susceptible / population
    peak_infected = infected
    crossed = False

    done = env.event()

Record S, I, R and the effective reproduction number each day. R_eff is the leading indicator: it crosses 1.0 exactly at the infection peak.

    def recorder():
        while True:
            probe("susceptible", env.now, susceptible)
            probe("infected", env.now, infected)
            probe("recovered", env.now, recovered)
            probe("r_eff", env.now, basic_r0 * susceptible / population)
            yield env.timeout(1)
    env.process(recorder())

Standard SIR update each day.

    def dynamics():
        nonlocal susceptible, infected, recovered, peak_infected, crossed
        for d in range(days):
            new_infections = beta * susceptible * infected / population
            new_recoveries = gamma * infected

            susceptible = max(susceptible - new_infections, 0.0)
            infected = max(infected + new_infections - new_recoveries, 0.0)
            recovered = min(recovered + new_recoveries, population)
            peak_infected = max(peak_infected, infected)

            if not crossed and basic_r0 * susceptible / population < 1.0:
                crossed = True
                progress(int(100 * d / days), "R_eff fell below 1 - epidemic turning over")

            yield env.timeout(1)

        progress(100)
        done.succeed({
            "peak_infected": peak_infected,
            "final_infected": infected,
            "final_recovered": recovered,
            "basic_r0": basic_r0,
            "r_eff_start": r_eff_start,

An outbreak means infections grew before they shrank.

            "outbreak": peak_infected > cfg["initial_infected"] * 1.05,
        })

    env.process(dynamics())
    env.run(until=done)
    return done.value


def requirements():
    return {
        "external": ["simpy==4.1.1"],
    }
Outbreak (R_eff > 1).yaml
# outbreak.yaml — the default. R0 = beta/gamma = 1.2, just above 1, and no one
# is vaccinated, so effective R0 starts at 1.2: a real (if slow) outbreak.
population: 1000
initial_infected: 5
beta: 0.24
gamma: 0.2
vaccinated_fraction: 0.0
sim_days: 180
Charts (Outbreak (R_eff > 1))

susceptible

Samples181 @ 0.00–180.00
Valuesmin 668.00, mean 781.12, median 732.83, max 995.00, σ 113.37

infected

Samples181 @ 0.00–180.00
Valuesmin 0.50, mean 9.16, median 8.42, max 19.17, σ 6.37

recovered

Samples181 @ 0.00–180.00
Valuesmin 0.00, mean 209.72, median 254.38, max 331.51, σ 116.80

r_eff

Samples181 @ 0.00–180.00
Valuesmin 0.80, mean 0.94, median 0.88, max 1.19, σ 0.14
Final Results (Outbreak (R_eff > 1))
MetricValue
peak_infected19.17
final_infected0.50
final_recovered331.51
basic_r01.20
r_eff_start1.19
outbreaktrue
Herd Immunity (R_eff < 1).yaml
# herd-immunity.yaml — the contrast: only vaccinated_fraction changes (0 -> 0.30).
# That is past the herd-immunity threshold 1 - 1/R0 = 0.167, so effective R0
# starts at 0.83 < 1 and infections fall from day one.
population: 1000
initial_infected: 5
beta: 0.24
gamma: 0.2
vaccinated_fraction: 0.30
sim_days: 180
Charts (Herd Immunity (R_eff < 1))

susceptible

Samples181 @ 0.00–180.00
Valuesmin 672.17, mean 675.51, median 672.84, max 695.00, σ 5.31

infected

Samples181 @ 0.00–180.00
Valuesmin 0.00, mean 0.77, median 0.17, max 5.00, σ 1.19

recovered

Samples181 @ 0.00–180.00
Valuesmin 300.00, mean 323.72, median 326.99, max 327.83, σ 6.50

r_eff

Samples181 @ 0.00–180.00
Valuesmin 0.81, mean 0.81, median 0.81, max 0.83, σ 0.01
Final Results (Herd Immunity (R_eff < 1))
MetricValue
peak_infected5.00
final_infected0.00
final_recovered327.83
basic_r01.20
r_eff_start0.83
outbreakfalse
FAQ
What sets the tipping point between an outbreak and a fizzle?
The effective reproduction number R_eff = R0 · (susceptible fraction), where R0 = beta/gamma is the average new cases one infection makes in a fully susceptible population. When R_eff > 1 each case more than replaces itself and infections climb; when R_eff < 1 the chain shrinks every generation. The whole epidemic turns over at the moment R_eff crosses 1.0 — which is exactly when the infected curve peaks.
How does vaccination stop the spread without vaccinating everyone?
Vaccinated people are removed from the susceptible pool, so they shrink the susceptible fraction R_eff depends on. You only need enough of them to drag R_eff below 1 from the start. That critical share is the herd-immunity threshold, 1 − 1/R0. Below it the disease still spreads (just slower); above it, it cannot get going — the unvaccinated are protected by the immune wall around them.
Why is the threshold here only about 17%?
Because R0 is just 1.2: 1 − 1/1.2 ≈ 0.167. The threshold scales with how contagious the disease is. Measles, with R0 around 15, needs roughly 1 − 1/15 ≈ 95% coverage — which is why measles outbreaks return as soon as local coverage slips below that line, while a mildly contagious disease needs far less.
Which probe shows the story best?
r_eff. In the outbreak run it starts above 1 and slides down through 1.0 right as infections peak; in the herd-immunity run it starts below 1 and never crosses, so infections only ever decline. The 1.0 line is the threshold made visible — the single chart where the two configs split.
What misconception does this debunk?
That you must vaccinate everyone to stop a disease. You only need to cross the threshold. It also shows R0 alone does not decide an outbreak — the susceptible fraction does, which is why the same disease can rage in one community and fizzle in a better-vaccinated one.