Lab 2: Simulate an Epidemic

Abstract a person down to three states

The first act of modeling is deciding what to ignore. Real disease involves biology, geography, age, behavior — infinite detail. The classic SIR model throws almost all of it away and keeps three states:

• S — Susceptible: can catch it.
• I — Infectious: has it, can spread it.
• R — Recovered: had it, now immune (and no longer spreading).

Every person is exactly one letter, and the only allowed transitions are S → I → R. That is a radical abstraction, and its power is exactly that radicalism: it is simple enough to reason about, yet — as you are about to see — rich enough to reproduce real epidemic behavior. Start the file:

import random

POP_SIZE = 500
CONTACTS_PER_DAY = 6       # people each infectious person meets per day
P_TRANSMIT = 0.04          # chance a contact with S actually infects
DAYS_INFECTIOUS = 8        # days a person stays in state I
MAX_DAYS = 200

def new_population(pop_size, initial_infected=3, vaccinated_frac=0.0):
    states = ["S"] * pop_size
    timers = [0] * pop_size
    n_vac = int(pop_size * vaccinated_frac)
    for i in range(n_vac):               # pre-immune (vaccinated) people
        states[i] = "R"
    for i in range(n_vac, n_vac + initial_infected):
        states[i] = "I"
        timers[i] = DAYS_INFECTIOUS
    return states, timers

The four constants at the top are your model. Everything interesting in this lab comes from changing them.

Model one day of contagion

Now the simulation rule — what happens in a single day. Each infectious person meets a few random others; each susceptible contact might get infected; infectious people count down and eventually recover:

def step(states, timers, contacts, p_transmit):
    pop = len(states)
    infected_today = set()
    # 1. Infectious people expose random contacts
    for i, state in enumerate(states):
        if state == "I":
            for _ in range(contacts):
                j = random.randrange(pop)
                if states[j] == "S" and random.random() < p_transmit:
                    infected_today.add(j)
    # 2. Existing infections age; some recover
    for i, state in enumerate(states):
        if state == "I":
            timers[i] -= 1
            if timers[i] <= 0:
                states[i] = "R"
    # 3. Apply today's new infections (after, so they don't spread same day)
    for j in infected_today:
        states[j] = "I"
        timers[j] = DAYS_INFECTIOUS

Run it forward and watch the curve

Add a driver and a text plotter, then run:

def run(contacts=CONTACTS_PER_DAY, p_transmit=P_TRANSMIT,
        vaccinated_frac=0.0, seed=None):
    if seed is not None:
        random.seed(seed)
    states, timers = new_population(POP_SIZE, vaccinated_frac=vaccinated_frac)
    history = []
    for _ in range(MAX_DAYS):
        counts = (states.count("S"), states.count("I"), states.count("R"))
        history.append(counts)
        if counts[1] == 0:      # no one infectious: epidemic is over
            break
        step(states, timers, contacts, p_transmit)
    return history

def plot(history, every=4, width=60):
    peak = max(i for _, i, _ in history) or 1
    for day in range(0, len(history), every):
        s, i, r = history[day]
        bar = "#" * round(i / peak * width)
        print(f"day {day:3d} | S={s:3d} I={i:3d} R={r:3d} | {bar}")

if __name__ == "__main__":
    history = run(seed=42)
    plot(history)
    s, i, r = history[-1]
    print(f"\nEpidemic ended after {len(history)-1} days.")
    print(f"Total ever infected: {POP_SIZE - s}/{POP_SIZE} "
          f"({(POP_SIZE - s)/POP_SIZE:.0%})")

Predict the shape first, then run python3 epidemic.py. You will see the famous epidemic curve: slow start, explosive rise, a peak (around day 32 with seed=42), then decline as the disease runs out of susceptibles. About 78% of the population is infected — and, crucially, roughly 22% never get infected even though nobody protected them. Hold that thought; it is herd immunity, and you did not program it.

Find the tipping point

Now use the model as an instrument. The key quantity in epidemiology is R0 — the average number of people one infectious person infects in a fully susceptible population. In this model:

R0 = CONTACTS_PER_DAY * P_TRANSMIT * DAYS_INFECTIOUS

Theory says R0 = 1 is a knife-edge: below it epidemics die out, above it they explode. Test whether your simulation agrees. Add a sweep (a separate script, or a new block) — averaging several runs because the model is stochastic:

from epidemic import run, POP_SIZE, CONTACTS_PER_DAY, DAYS_INFECTIOUS

def avg_total_infected(p, trials=10):
    return sum(POP_SIZE - run(p_transmit=p)[-1][0]
               for _ in range(trials)) / trials

print(f"{'p':>5} {'R0':>5} {'avg total infected':>20}")
for p in [0.01, 0.02, 0.03, 0.04, 0.06]:
    r0 = CONTACTS_PER_DAY * p * DAYS_INFECTIOUS
    tot = avg_total_infected(p)
    print(f"{p:5.2f} {r0:5.2f} {tot:12.0f} ({tot/POP_SIZE:4.0%})")

Typical output:

    p    R0   avg total infected
 0.01  0.48            4 (  1%)
 0.02  0.96           19 (  4%)
 0.03  1.44          165 ( 33%)
 0.04  1.92          348 ( 70%)
 0.06  2.88          460 ( 92%)

Look at what happens between R0 = 0.96 and R0 = 1.44: the outcome leaps from 4% to 33% of the population. That is a tipping point — a small, smooth change in a parameter producing a large, abrupt change in the system’s behavior. Nothing in your code says “explode above 1”; the threshold emerges from the interaction of contacts, transmission, and recovery. Recognizing that a system has such thresholds — and roughly where they sit — is the essence of systems thinking.

Discover herd immunity

Return to the 22%-never-infected observation and probe it deliberately. Pre-immunize part of the population (the vaccinated_frac parameter already supports this) and watch what happens to the people who are not vaccinated:

from epidemic import run, POP_SIZE

for v in [0.0, 0.2, 0.4, 0.6, 0.8]:
    vaccinated = int(POP_SIZE * v)
    trials = [run(p_transmit=0.06, vaccinated_frac=v) for _ in range(15)]
    peak = sum(max(i for _, i, _ in h) for h in trials) / len(trials)
    # subtract the pre-immune so we count only NEW infections
    newly = sum(POP_SIZE - h[-1][0] - vaccinated for h in trials) / len(trials)
    print(f"vaccinated={v:.0%}  avg peak infectious={peak:4.0f}  "
          f"avg newly infected={newly:4.0f}")

Typical output (with a nasty R0 ≈ 2.9 disease):

vaccinated= 0%  avg peak infectious= 257  avg newly infected= 468
vaccinated=20%  avg peak infectious= 154  avg newly infected= 344
vaccinated=40%  avg peak infectious=  59  avg newly infected= 208
vaccinated=60%  avg peak infectious=  11  avg newly infected=  34
vaccinated=80%  avg peak infectious=   5  avg newly infected=   8

Notice the nonlinearity. Going from 40% to 60% vaccinated collapses the peak from 59 to 11 and newly infected from 208 to 34 — far more than a proportional dent. Somewhere near 60% coverage the effective reproduction number drops below 1 and the epidemic can no longer sustain itself. Beyond that threshold, even unvaccinated people are largely protected because the chains of transmission break. That is herd immunity, and it is a pure systems-thinking phenomenon: a property of the network, not of any individual.

See the feedback loop

Why does every curve eventually fall even with nobody recovering-then-reinfecting? Add a one-line diagnostic to run’s loop to print, each day, the current number of susceptibles alongside new infections, and you will see the mechanism directly:

    flowchart LR
    I[More infectious people] --> N[More new infections]
    N --> R[Susceptible pool shrinks]
    R --> F[Fewer S left to infect]
    F --> D[Infection rate falls]
    D --> I
  

This is a balancing feedback loop. Early on, more infections beget more infections — a reinforcing spiral that produces the explosive rise. But every infection permanently removes a susceptible, and once susceptibles get scarce the same loop runs in reverse and starves the epidemic. The peak is exactly the moment the loop flips sign. You built three simple rules; the loop, the peak, and the tipping point are all things the rules do without being told to.