Chance and Fizzles

How Epidemics Spread

Chapter 6 · Chance and Fizzles

A traveller steps off a train in a small town, carrying an infection that is, on paper, perfectly capable of causing an epidemic. Most of the time, nothing happens. The traveller recovers, perhaps passes it to one or two people, and the chain of cases quietly ends. Every so often the same arrival lights a fire. The models in How Epidemics Spread 3, How Epidemics Spread 4 and How Epidemics Spread 5 could not tell these two stories apart, because they tracked averages and gave the same answer every time. How Epidemics Spread 5 ended by suggesting that the deep troughs in its endemic runs would, in a real population, end in extinction. This chapter puts chance into the simulations and asks three questions: how often does a single case fizzle out, why does the answer depend on more than R0, and what do super-spreaders change?

What was run, and what to trust

Every Python snippet in this chapter and in the three exercise solutions was executed, and the output in the comments is what Python printed (random seeds are fixed, so you should get the same numbers). The criterion for branching processes (extinction certain if the mean number of offspring is 1 or less, uncertain if above 1, and the extinction probability found by starting at 0 and repeatedly applying the offspring generating function) comes from the Wikipedia article on branching processes, fetched with a tool that often returns truncated pages or model-written summaries, so treat the wording as a paraphrase. The super-spreading set-up (a gamma-distributed individual reproductive number, giving negative binomial offspring with mean R0 and dispersion k, with k = 1 the geometric and k going to infinity the Poisson case, and a reported SARS estimate of k = 0.16 with a 90% interval of 0.11 to 0.64) is from Lloyd-Smith and colleagues’ 2005 paper in Nature; I saw it only through search-result summaries, not the full paper, so check the paper before quoting the numbers. Those summaries conflict in emphasis on one point: some retellings speak of a “20/80 rule”, while the paper itself is reported as finding that no general 20/80 rule holds for directly transmitted diseases, so I do not claim one. The statement that a one-case outbreak fizzles with probability 1/R0 in the stochastic SIR model, and the formulas for the geometric and negative binomial cases, are my own algebra, checked here only against the simulations; they are inferences, not quotations. Every population size, R0, k, run count and cap in the code is assumed or illustrative, not an estimate for any real disease.

Putting Chance into the SIR Model

In a stochastic simulation, people are counted as whole individuals and each step is a random event. With S susceptibles and I infectious people, the next event is either an infection (at rate βSI/N) or a recovery (at rate γI). The chance that it is an infection is the first rate divided by the sum of the two. The code below ignores the exact time of each event, because we only want to know how many people end up infected. It starts with a single case in a population of 1,000, repeats the whole outbreak 2,000 times for each R0, and calls a run “fizzled” if 20 or fewer people were ever infected. That threshold of 20 is an arbitrary choice for this chapter. There is no calculus here, only a loop and a random number.

# Stochastic SIR: one random event at a time. Population, rates and run counts are assumed.
import random

gamma = 0.2                      # assumed: average infectious period 5 days
population = 1000

def outbreak(r0):
    """Start with one case; return how many people were ever infected."""
    beta = r0 * gamma
    S, I, ever = population - 1, 1, 1
    while I > 0:
        infect_rate = beta * S * I / population
        recover_rate = gamma * I
        if random.random() < infect_rate / (infect_rate + recover_rate):
            S -= 1
            I += 1
            ever += 1
        else:
            I -= 1
    return ever

random.seed(1)
runs = 2000
print("R0    fizzled (<=20 cases)   1/R0    mean size of the rest")
for r0 in [0.8, 1.5, 2.5, 4.0]:
    sizes = [outbreak(r0) for _ in range(runs)]
    small = [s for s in sizes if s <= 20]
    big = [s for s in sizes if s > 20]
    mean_big = sum(big) / len(big) / population if big else float("nan")
    print(f"{r0:3.1f} {len(small) / runs:20.1%} {min(1, 1 / r0):7.1%} {mean_big:20.1%}")
# R0    fizzled (<=20 cases)   1/R0    mean size of the rest
# 0.8                95.1%  100.0%                 4.3%
# 1.5                66.5%   66.7%                54.7%
# 2.5                41.1%   40.0%                89.2%
# 4.0                25.1%   25.0%                98.0%
    

At R0 = 2.5, the plain model of How Epidemics Spread 3 says the epidemic always happens and reaches 89.3 per cent of the population. The simulation says that about 41 runs in 100 fizzle, and the others reach 89.2 per cent on average, almost exactly the deterministic figure. At R0 = 1.5 and 4.0 the picture is the same: a mixture of fizzles and full epidemics, and the share that fizzle is almost exactly 1/R0 (66.5 against 66.7 per cent, and 25.1 against 25.0 per cent). At R0 = 0.8 nearly every run fizzles; the 4.9 per cent that passed 20 cases are unlucky long chains, and they average only 4.3 per cent of the population. Run the deterministic rule (repeating z = 1 − exp(−R0 × z)) and you get about 58.3 per cent for R0 = 1.5 and 98.0 per cent for R0 = 4.0; the simulated means are 54.7 and 98.0, so the model gets the size of a real outbreak roughly right, and the fizzle rate is the part the earlier chapters could not see. That 1/R0 rule is my inference, and the next section explains why it should hold.

	Deterministic (How Epidemics Spread 3)	Stochastic (this chapter)
Unit	a fraction of the population	whole people
Same inputs, two runs	identical output	different output every run (unless the seed is fixed)
One case, R0 = 2.5	always an epidemic	about 4 runs in 10 fizzle
Result reported	one curve	a spread of outcomes, a fizzle rate, a mean size
Cost	one run	hundreds or thousands of runs

Branching Processes: the Early Days

Why 1/R0? At the start, almost everyone is susceptible, so each case behaves independently of the others. A case infects some number of people, each of whom infects some number, and so on. That is a branching process. If the average number of offspring is 1 or less, the chain dies out for certain; if it is above 1, it may still die out, with a probability q below 1 (this is the criterion from the article). To find q, write the probability of having k offspring, build the function f(q) = sum of P(k) × q^k, and repeat q ← f(q) starting from q = 0 until the number stops changing. The method is just repeated substitution.

The offspring distribution matters. Two natural choices both have mean R0. With Poisson offspring (many independent contacts, each with a small chance of transmission) f(q) = exp(R0(q − 1)). In the stochastic SIR model, each step before the susceptibles run low is a race between an infection and the case’s own recovery, so the number of people one case infects is geometric: f(q) = 1 / (1 + R0(1 − q)). Solving q = f(q) for that case gives q = 1/R0, which explains the fizzle rates above (my algebra). The code checks both against direct simulation, calling a chain “taken off” once it has reached 500 people (an assumed cap).

# Branching process: each case independently infects a random number of others.
# Poisson and geometric offspring, both with mean R0. R0 values and run counts are assumed.
import math
import random

def poisson(lam):
    """Knuth's method, in chunks so exp(-lam) never underflows."""
    total = 0
    while lam > 30:
        total += poisson(30)
        lam -= 30
    limit, k, p = math.exp(-lam), 0, 1.0
    while True:
        p *= random.random()
        if p <= limit:
            return total + k
        k += 1

def geometric(mean):
    """P(k) = (1 - p) * p**k, with p chosen so the mean is `mean`."""
    p = mean / (1 + mean)
    return int(math.log(1 - random.random()) / math.log(p))

def extinct_probability(generating, steps=2000):
    """Smallest solution of q = generating(q), found by repeating from q = 0."""
    q = 0.0
    for _ in range(steps):
        q = generating(q)
    return q

def escapes(draw, cap=500):
    """True if the chain of cases reaches `cap` people before dying out."""
    active, total = 1, 1
    while active > 0 and total < cap:
        new = sum(draw() for _ in range(active))
        active = new
        total += new
    return active > 0

random.seed(2)
runs = 4000
print("R0    Poisson: formula  simulated    geometric: formula  simulated")
for r0 in [0.9, 1.2, 1.5, 2.5, 4.0]:
    q_pois = extinct_probability(lambda q: math.exp(r0 * (q - 1)))
    q_geom = extinct_probability(lambda q: 1 / (1 + r0 * (1 - q)))
    sim_p = 1 - sum(escapes(lambda: poisson(r0)) for _ in range(runs)) / runs
    sim_g = 1 - sum(escapes(lambda: geometric(r0)) for _ in range(runs)) / runs
    print(f"{r0:3.1f} {q_pois:17.3f} {sim_p:10.3f} {q_geom:20.3f} {sim_g:10.3f}")
# R0    Poisson: formula  simulated    geometric: formula  simulated
# 0.9             1.000      1.000                1.000      0.999
# 1.2             0.686      0.681                0.833      0.837
# 1.5             0.417      0.430                0.667      0.663
# 2.5             0.107      0.106                0.400      0.404
# 4.0             0.020      0.017                0.250      0.247
    

The formula and the simulation agree (to within noise from 4,000 runs each; the 0.999 at R0 = 0.9 is a rare chain that reached the cap of 500 before dying out). The surprise is how far apart the two columns are. With R0 = 2.5 the Poisson case fizzles only 10.7 per cent of the time, while the geometric case fizzles 40.0 per cent of the time, nearly four times as often. Both have exactly the same average number of secondary cases. The difference is the spread: the geometric distribution has many more zeros and a few large counts. R0 alone does not say how likely an outbreak is to fizzle. The spread of the offspring counts matters too.

Super-Spreaders: When Transmission Is Uneven

Real transmission is often more uneven than either of those. The 2005 paper models each person’s own reproductive number as a random value with mean R0, gamma-distributed with a dispersion parameter k; the resulting offspring counts are negative binomial. Small k means very uneven: most cases infect nobody and a few infect many. k = 1 gives the geometric case and k going to infinity gives Poisson. The paper reports k of about 0.16 for SARS. The code samples from this model with R0 = 2.5 and several assumed values of k, and reports four things: the share of cases who infect nobody, the fizzle probability from the formula q = (1 + R0(1 − q)/k)^−k (same repeated substitution; my derivation from the standard generating function, checked below), the fizzle rate from direct simulation, and the share of all onward infections caused by the 20 per cent most infectious cases.

# Super-spreading: negative binomial offspring with mean R0 and dispersion k (small k = more uneven).
# R0 = 2.5 and every k below are assumed, except that k = 0.16 echoes a published SARS estimate.
import math
import random

r0 = 2.5

def poisson(lam):
    total = 0
    while lam > 30:
        total += poisson(30)
        lam -= 30
    limit, k, p = math.exp(-lam), 0, 1.0
    while True:
        p *= random.random()
        if p <= limit:
            return total + k
        k += 1

def offspring(k):
    """Gamma-Poisson mixture: each case gets its own mean, then a Poisson draw."""
    if k == float("inf"):
        return poisson(r0)
    return poisson(random.gammavariate(k, r0 / k))

def generating(q, k):
    if k == float("inf"):
        return math.exp(r0 * (q - 1))
    return (1 + r0 * (1 - q) / k) ** (-k)

def extinct_probability(k):
    q = 0.0
    for _ in range(3000):
        q = generating(q, k)
    return q

def escapes(k, cap=500):
    active, total = 1, 1
    while active > 0 and total < cap:
        active = sum(offspring(k) for _ in range(active))
        total += active
    return active > 0

random.seed(3)
runs, sample = 3000, 100_000
print("   k   infect nobody   fizzle: formula  simulated   share of spread from top 20%")
for k in [0.1, 0.3, 1.0, 10.0, float("inf")]:
    draws = sorted((offspring(k) for _ in range(sample)), reverse=True)
    top_share = sum(draws[: sample // 5]) / sum(draws)
    nobody = sum(1 for d in draws if d == 0) / sample
    sim = 1 - sum(escapes(k) for _ in range(runs)) / runs
    label = "inf" if k == float("inf") else f"{k:.1f}"
    print(f"{label:>4} {nobody:14.1%} {extinct_probability(k):16.3f} {sim:10.3f} {top_share:24.1%}")
#    k   infect nobody   fizzle: formula  simulated   share of spread from top 20%
#  0.1          72.2%            0.861      0.859                    96.5%
#  0.3          51.1%            0.675      0.678                    79.0%
#  1.0          28.6%            0.400      0.398                    58.0%
# 10.0          10.7%            0.144      0.152                    41.3%
#  inf           8.3%            0.107      0.113                    38.9%
    

As k falls, everything shifts in the same direction. At k = 0.1, 72.2 per cent of cases infect nobody and the top fifth of cases cause 96.5 per cent of the onward spread, yet the average is still 2.5. A single introduction fizzles 86 times in 100, against about 11 times in 100 for Poisson. The formula and simulation agree throughout (0.861 and 0.859 at k = 0.1, 0.400 and 0.398 at k = 1). Notice also the last column: even with Poisson offspring the top 20 per cent cause 38.9 per cent of spread, simply because some cases happen to have more offspring than others. These figures come from my simulation with assumed parameters; they are not the paper’s results, and they are not a general 20/80 rule.

The practical lesson is modest but important. When transmission is uneven, a single introduction is much less likely to start an epidemic, but when the arrival is a super-spreading event the same R0 can be reached in one step. Control measures aimed at the unusually infectious events can then matter more than measures aimed at everyone equally. That last sentence is an inference from the model, not a result I tested; How Epidemics Spread 7 and How Epidemics Spread 9 will test targeted control properly.

	Poisson (k = infinity)	Geometric (k = 1)	Very uneven (k = 0.1)
Mean offspring (R0)	2.5	2.5	2.5
Cases who infect nobody	8.3%	28.6%	72.2%
One case fizzles	about 11%	40%	about 86%
Share of spread from top 20% of cases	38.9%	58.0%	96.5%
Typical story	almost every introduction grows	a coin flip with a bias	many dead ends, a few explosive clusters

Common mistakes

Believing that R0 > 1 guarantees an epidemic (it only makes one possible; the fizzle chance can be large); reading “the epidemic fizzled” as proof the disease is mild or that a control measure worked (many untreated introductions fizzle anyway); using a branching-process formula late in an outbreak, when the susceptibles are no longer plentiful and cases are no longer independent; treating the 20/80 idea as a law (the paper reports it does not hold in general); running one simulation and quoting its result (a single run is one random story, so report a rate from many runs); and using the illustrative numbers here (R0, k, the 500-person cap, the 20-case threshold) as if they were estimates for a real disease.

Exercise 1

Run the stochastic SIR model with R0 = 2 (assumed) for populations of 50, 200 and 1,000, 2,000 runs each. Print the share of runs that fizzle (10 or fewer cases) and the mean share infected in the rest, and compare them with 1/R0 and with the deterministic final share.

📄 View solution

Exercise 2

With Poisson offspring and R0 = 1.1, 1.5, 2 and 3 (assumed), find the extinction probability q by repeated substitution. Print the chance that 1, 3 and 5 independent introductions take off (1 − q to the power of the number of seeds), and check the 1-seed and 3-seed values by simulation.

📄 View solution

Exercise 3

For R0 = 1.2 and 2.5 and k = 0.16, 1 and infinity (all assumed), print the chance a case infects nobody, the fizzle probability from the formula, and the fizzle rate from 3,000 simulated chains. For which R0 does overdispersion change the fizzle chance more?

📄 View solution

Quick Reference

Stochastic model: whole people, one random event at a time; report rates from many runs, not one run
Branching process: early on each case acts independently; extinction certain if mean offspring is 1 or less
Extinction probability q: smallest solution of q = f(q); repeat q ← f(q) from q = 0
Poisson offspring: f(q) = exp(R0(q − 1))
Geometric offspring (stochastic SIR, one case): q = 1/R0 (my algebra, matched by the runs)
Negative binomial (R0, k): q = (1 + R0(1 − q)/k)^−k; small k means many dead ends and a few big spreaders; SARS reported near k = 0.16
With several introductions: chance at least one takes off = 1 − q^seeds
Assumes: independent cases, illustrative numbers, an arbitrary cap and threshold

Coming next

How Epidemics Spread 7 is about networks and contacts: random, clustered and hub networks, simulating spread on a graph, targeted control.