SIR Outbreak Size (Kermack–McKendrick)¶

Theory¶

Kermack and McKendrick (1927) showed that in a closed SIR model the final fraction of susceptibles $s_\infty = S(\infty)/N$ satisfies the final-size equation:

$$s_\infty = e^{-R_0\,(1 - s_\infty)}$$

where $R_0 = \beta / \gamma$. The attack rate (fraction ever infected) is:

$$\text{AR} = 1 - s_\infty$$

Key properties:

$R_0 \leq 1$: no epidemic, AR ≈ 0.
$R_0 > 1$: a major outbreak occurs; AR rises steeply with R₀.

An epidemic occurs only when $R_0 \cdot S(0)/N > 1$ (the epidemic threshold).

Goal: Compare the theoretical attack rate with the simulated final attack rate across a range of R₀ values, and show the full S/I/R time series for R₀ = 3.

In [1]:

Copied!





import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import brentq

from laser.core import PropertySet
from laser.core.utils import grid
from laser.generic.utils import ValuesMap
from laser.cohorts import Model
import laser.cohorts.SIR as SIR
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import brentq

from laser.core import PropertySet
from laser.core.utils import grid
from laser.generic.utils import ValuesMap
from laser.cohorts import Model
import laser.cohorts.SIR as SIR

In [2]:

Copied!





def solve_final_size(R0: float, s0: float = 1.0) -> float:
    """Return the expected attack rate for a given R0 and initial susceptible fraction.

    Solves the Kermack-McKendrick final-size equation
        s_inf = exp(-R0 * (1 - s_inf))
    for s_inf using Brent's method, then returns 1 - s_inf.

    Args:
        R0:  Basic reproduction number.
        s0:  Initial susceptible fraction (default 1.0 ≈ fully naive population).

    Returns:
        Attack rate AR = 1 - s_inf in [0, 1].
    """
    if R0 <= 1.0:
        return 0.0   # sub-threshold: no epidemic

    # f(s) = 0 at the non-trivial root  (s < s0)
    f = lambda s: s - np.exp(-R0 * (s0 - s))

    # Root lies in (0, s0); use a tiny lower bound to avoid the trivial root
    s_inf = brentq(f, 1e-9, s0 - 1e-9)
    return 1.0 - s_inf
def solve_final_size(R0: float, s0: float = 1.0) -> float:
    """Return the expected attack rate for a given R0 and initial susceptible fraction.

    Solves the Kermack-McKendrick final-size equation
        s_inf = exp(-R0 * (1 - s_inf))
    for s_inf using Brent's method, then returns 1 - s_inf.

    Args:
        R0:  Basic reproduction number.
        s0:  Initial susceptible fraction (default 1.0 ≈ fully naive population).

    Returns:
        Attack rate AR = 1 - s_inf in [0, 1].
    """
    if R0 <= 1.0:
        return 0.0   # sub-threshold: no epidemic

    # f(s) = 0 at the non-trivial root  (s < s0)
    f = lambda s: s - np.exp(-R0 * (s0 - s))

    # Root lies in (0, s0); use a tiny lower bound to avoid the trivial root
    s_inf = brentq(f, 1e-9, s0 - 1e-9)
    return 1.0 - s_inf

In [3]:

Copied!





# --- Sweep over R0 values ---
N      = 50_000
I0     = 10
nticks = 300
gamma  = 0.1      # recovery rate (per day); R0 = beta / gamma

R0_values = [0.5, 1.0, 1.5, 2.0, 3.0, 5.0]

theoretical_AR = []
simulated_AR   = []

for R0 in R0_values:
    beta_val = R0 * gamma

    scenario = grid(M=1, N=1)
    scenario["S"] = N - I0
    scenario["I"] = I0
    scenario["R"] = 0

    params = PropertySet({"nticks": nticks, "beta": beta_val, "r_recovery": gamma})
    model  = Model(scenario, params)

    betas        = ValuesMap.from_scalar(beta_val, nticks, len(scenario))
    r_recoveries = ValuesMap.from_scalar(gamma,    nticks, len(scenario))

    model.components = [
        SIR.Susceptible(model),
        SIR.Infectious(model,  r_recovery=r_recoveries),
        SIR.Recovered(model),
        SIR.Transmission(model, beta=betas),
    ]
    model.run()

    # Attack rate = everyone who was ever infected = I_final + R_final
    I_final = model.states.I[-1, 0]
    R_final = model.states.R[-1, 0]
    sim_AR  = (I_final + R_final) / N

    s0 = (N - I0) / N
    th_AR = solve_final_size(R0, s0=s0)

    theoretical_AR.append(th_AR)
    simulated_AR.append(sim_AR)
    print(f"R₀ = {R0:.1f}  |  theory AR = {th_AR:.3f}  |  sim AR = {sim_AR:.3f}")
# --- Sweep over R0 values ---
N      = 50_000
I0     = 10
nticks = 300
gamma  = 0.1      # recovery rate (per day); R0 = beta / gamma

R0_values = [0.5, 1.0, 1.5, 2.0, 3.0, 5.0]

theoretical_AR = []
simulated_AR   = []

for R0 in R0_values:
    beta_val = R0 * gamma

    scenario = grid(M=1, N=1)
    scenario["S"] = N - I0
    scenario["I"] = I0
    scenario["R"] = 0

    params = PropertySet({"nticks": nticks, "beta": beta_val, "r_recovery": gamma})
    model  = Model(scenario, params)

    betas        = ValuesMap.from_scalar(beta_val, nticks, len(scenario))
    r_recoveries = ValuesMap.from_scalar(gamma,    nticks, len(scenario))

    model.components = [
        SIR.Susceptible(model),
        SIR.Infectious(model,  r_recovery=r_recoveries),
        SIR.Recovered(model),
        SIR.Transmission(model, beta=betas),
    ]
    model.run()

    # Attack rate = everyone who was ever infected = I_final + R_final
    I_final = model.states.I[-1, 0]
    R_final = model.states.R[-1, 0]
    sim_AR  = (I_final + R_final) / N

    s0 = (N - I0) / N
    th_AR = solve_final_size(R0, s0=s0)

    theoretical_AR.append(th_AR)
    simulated_AR.append(sim_AR)
    print(f"R₀ = {R0:.1f}  |  theory AR = {th_AR:.3f}  |  sim AR = {sim_AR:.3f}")

R₀ = 0.5  |  theory AR = 0.000  |  sim AR = 0.000
R₀ = 1.0  |  theory AR = 0.000  |  sim AR = 0.000

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
Cell In[3], line 40
     36     R_final = model.states.R[-1, 0]
     37     sim_AR  = (I_final + R_final) / N
     38 
     39     s0 = (N - I0) / N
---> 40     th_AR = solve_final_size(R0, s0=s0)
     41 
     42     theoretical_AR.append(th_AR)
     43     simulated_AR.append(sim_AR)

Cell In[2], line 22, in solve_final_size(R0, s0)
     18     # f(s) = 0 at the non-trivial root  (s < s0)
     19     f = lambda s: s - np.exp(-R0 * (s0 - s))
     20 
     21     # Root lies in (0, s0); use a tiny lower bound to avoid the trivial root
---> 22     s_inf = brentq(f, 1e-9, s0 - 1e-9)
     23     return 1.0 - s_inf

File ~/projects/laser-fresh/laser.cohorts/.docs/lib/python3.12/site-packages/scipy/optimize/_zeros_py.py:846, in brentq(f, a, b, args, xtol, rtol, maxiter, full_output, disp)
    844     raise ValueError(f"rtol too small ({rtol:g} < {_rtol:g})")
    845 f = _wrap_nan_raise(f)
--> 846 r = _zeros._brentq(f, a, b, xtol, rtol, maxiter, args, full_output, disp)
    847 return results_c(full_output, r, "brentq")

ValueError: f(a) and f(b) must have different signs

In [4]:

Copied!





fig, ax = plt.subplots(figsize=(6, 5))

ax.scatter(theoretical_AR, simulated_AR, s=80, zorder=3,
           c=R0_values, cmap="plasma", label="(theory, sim) pairs")

# Annotate each point with its R0
for th, sim, R0 in zip(theoretical_AR, simulated_AR, R0_values):
    ax.annotate(f"R₀={R0}", (th, sim), textcoords="offset points",
                xytext=(6, 3), fontsize=8)

# y = x reference line
lim = max(max(theoretical_AR), max(simulated_AR)) * 1.05
ax.plot([0, lim], [0, lim], "k--", linewidth=0.8, label="y = x (perfect agreement)")

ax.set_xlim(0, lim)
ax.set_ylim(0, lim)
ax.set_xlabel("Theoretical attack rate")
ax.set_ylabel("Simulated attack rate")
ax.set_title("Kermack–McKendrick final-size: theory vs. simulation")
ax.legend(fontsize=8)
plt.tight_layout()
plt.show()
fig, ax = plt.subplots(figsize=(6, 5))

ax.scatter(theoretical_AR, simulated_AR, s=80, zorder=3,
           c=R0_values, cmap="plasma", label="(theory, sim) pairs")

# Annotate each point with its R0
for th, sim, R0 in zip(theoretical_AR, simulated_AR, R0_values):
    ax.annotate(f"R₀={R0}", (th, sim), textcoords="offset points",
                xytext=(6, 3), fontsize=8)

# y = x reference line
lim = max(max(theoretical_AR), max(simulated_AR)) * 1.05
ax.plot([0, lim], [0, lim], "k--", linewidth=0.8, label="y = x (perfect agreement)")

ax.set_xlim(0, lim)
ax.set_ylim(0, lim)
ax.set_xlabel("Theoretical attack rate")
ax.set_ylabel("Simulated attack rate")
ax.set_title("Kermack–McKendrick final-size: theory vs. simulation")
ax.legend(fontsize=8)
plt.tight_layout()
plt.show()

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
File ~/projects/laser-fresh/laser.cohorts/.docs/lib/python3.12/site-packages/matplotlib/axes/_axes.py:4761, in Axes._parse_scatter_color_args(c, edgecolors, kwargs, xsize, get_next_color_func)
   4760 try:  # Is 'c' acceptable as PathCollection facecolors?
-> 4761     colors = mcolors.to_rgba_array(c)
   4762 except (TypeError, ValueError) as err:

File ~/projects/laser-fresh/laser.cohorts/.docs/lib/python3.12/site-packages/matplotlib/colors.py:515, in to_rgba_array(c, alpha)
    514 else:
--> 515     rgba = np.array([to_rgba(cc) for cc in c])
    517 if alpha is not None:

File ~/projects/laser-fresh/laser.cohorts/.docs/lib/python3.12/site-packages/matplotlib/colors.py:317, in to_rgba(c, alpha)
    316 if rgba is None:  # Suppress exception chaining of cache lookup failure.
--> 317     rgba = _to_rgba_no_colorcycle(c, alpha)
    318     try:

File ~/projects/laser-fresh/laser.cohorts/.docs/lib/python3.12/site-packages/matplotlib/colors.py:401, in _to_rgba_no_colorcycle(c, alpha)
    400 if not np.iterable(c):
--> 401     raise ValueError(f"Invalid RGBA argument: {orig_c!r}")
    402 if len(c) not in [3, 4]:

ValueError: Invalid RGBA argument: np.float64(0.5)

The above exception was the direct cause of the following exception:

ValueError                                Traceback (most recent call last)
Cell In[4], line 3
      1 fig, ax = plt.subplots(figsize=(6, 5))
      2 
----> 3 ax.scatter(theoretical_AR, simulated_AR, s=80, zorder=3,
      4            c=R0_values, cmap="plasma", label="(theory, sim) pairs")
      5 
      6 # Annotate each point with its R0

File ~/projects/laser-fresh/laser.cohorts/.docs/lib/python3.12/site-packages/matplotlib/_api/deprecation.py:453, in make_keyword_only.<locals>.wrapper(*args, **kwargs)
    447 if len(args) > name_idx:
    448     warn_deprecated(
    449         since, message="Passing the %(name)s %(obj_type)s "
    450         "positionally is deprecated since Matplotlib %(since)s; the "
    451         "parameter will become keyword-only in %(removal)s.",
    452         name=name, obj_type=f"parameter of {func.__name__}()")
--> 453 return func(*args, **kwargs)

File ~/projects/laser-fresh/laser.cohorts/.docs/lib/python3.12/site-packages/matplotlib/__init__.py:1524, in _preprocess_data.<locals>.inner(ax, data, *args, **kwargs)
   1521 @functools.wraps(func)
   1522 def inner(ax, *args, data=None, **kwargs):
   1523     if data is None:
-> 1524         return func(
   1525             ax,
   1526             *map(cbook.sanitize_sequence, args),
   1527             **{k: cbook.sanitize_sequence(v) for k, v in kwargs.items()})
   1529     bound = new_sig.bind(ax, *args, **kwargs)
   1530     auto_label = (bound.arguments.get(label_namer)
   1531                   or bound.kwargs.get(label_namer))

File ~/projects/laser-fresh/laser.cohorts/.docs/lib/python3.12/site-packages/matplotlib/axes/_axes.py:4954, in Axes.scatter(self, x, y, s, c, marker, cmap, norm, vmin, vmax, alpha, linewidths, edgecolors, colorizer, plotnonfinite, **kwargs)
   4951 if edgecolors is None:
   4952     orig_edgecolor = kwargs.get('edgecolor', None)
   4953 c, colors, edgecolors = \
-> 4954     self._parse_scatter_color_args(
   4955         c, edgecolors, kwargs, x.size,
   4956         get_next_color_func=self._get_patches_for_fill.get_next_color)
   4958 if plotnonfinite and colors is None:
   4959     c = np.ma.masked_invalid(c)

File ~/projects/laser-fresh/laser.cohorts/.docs/lib/python3.12/site-packages/matplotlib/axes/_axes.py:4767, in Axes._parse_scatter_color_args(c, edgecolors, kwargs, xsize, get_next_color_func)
   4765 else:
   4766     if not valid_shape:
-> 4767         raise invalid_shape_exception(c.size, xsize) from err
   4768     # Both the mapping *and* the RGBA conversion failed: pretty
   4769     # severe failure => one may appreciate a verbose feedback.
   4770     raise ValueError(
   4771         f"'c' argument must be a color, a sequence of colors, "
   4772         f"or a sequence of numbers, not {c!r}") from err

ValueError: 'c' argument has 6 elements, which is inconsistent with 'x' and 'y' with size 2.

No description has been provided for this image

In [5]:

Copied!





# --- Full SIR time series for R0 = 3 ---
R0_demo  = 3.0
beta_demo = R0_demo * gamma

scenario_demo = grid(M=1, N=1)
scenario_demo["S"] = N - I0
scenario_demo["I"] = I0
scenario_demo["R"] = 0

params_demo = PropertySet({"nticks": nticks, "beta": beta_demo, "r_recovery": gamma})
model_demo  = Model(scenario_demo, params_demo)

betas_demo       = ValuesMap.from_scalar(beta_demo, nticks, len(scenario_demo))
r_recoveries_demo = ValuesMap.from_scalar(gamma,    nticks, len(scenario_demo))

model_demo.components = [
    SIR.Susceptible(model_demo),
    SIR.Infectious(model_demo,  r_recovery=r_recoveries_demo),
    SIR.Recovered(model_demo),
    SIR.Transmission(model_demo, beta=betas_demo),
]
model_demo.run()

t = np.arange(nticks + 1)
S_d = model_demo.states.S[:, 0]
I_d = model_demo.states.I[:, 0]
R_d = model_demo.states.R[:, 0]

fig, ax = plt.subplots(figsize=(9, 4))
ax.plot(t, S_d, color="steelblue", label="S")
ax.plot(t, I_d, color="firebrick", label="I")
ax.plot(t, R_d, color="seagreen",  label="R")
ax.set_xlabel("Day")
ax.set_ylabel("Count")
ax.set_title(f"SIR model — full time series  (R₀ = {R0_demo})")
ax.legend()
plt.tight_layout()
plt.show()

final_AR = (I_d[-1] + R_d[-1]) / N
print(f"Final attack rate (simulated): {final_AR:.3f}")
print(f"Theoretical attack rate      : {solve_final_size(R0_demo, s0=(N-I0)/N):.3f}")
# --- Full SIR time series for R0 = 3 ---
R0_demo  = 3.0
beta_demo = R0_demo * gamma

scenario_demo = grid(M=1, N=1)
scenario_demo["S"] = N - I0
scenario_demo["I"] = I0
scenario_demo["R"] = 0

params_demo = PropertySet({"nticks": nticks, "beta": beta_demo, "r_recovery": gamma})
model_demo  = Model(scenario_demo, params_demo)

betas_demo       = ValuesMap.from_scalar(beta_demo, nticks, len(scenario_demo))
r_recoveries_demo = ValuesMap.from_scalar(gamma,    nticks, len(scenario_demo))

model_demo.components = [
    SIR.Susceptible(model_demo),
    SIR.Infectious(model_demo,  r_recovery=r_recoveries_demo),
    SIR.Recovered(model_demo),
    SIR.Transmission(model_demo, beta=betas_demo),
]
model_demo.run()

t = np.arange(nticks + 1)
S_d = model_demo.states.S[:, 0]
I_d = model_demo.states.I[:, 0]
R_d = model_demo.states.R[:, 0]

fig, ax = plt.subplots(figsize=(9, 4))
ax.plot(t, S_d, color="steelblue", label="S")
ax.plot(t, I_d, color="firebrick", label="I")
ax.plot(t, R_d, color="seagreen",  label="R")
ax.set_xlabel("Day")
ax.set_ylabel("Count")
ax.set_title(f"SIR model — full time series  (R₀ = {R0_demo})")
ax.legend()
plt.tight_layout()
plt.show()

final_AR = (I_d[-1] + R_d[-1]) / N
print(f"Final attack rate (simulated): {final_AR:.3f}")
print(f"Theoretical attack rate      : {solve_final_size(R0_demo, s0=(N-I0)/N):.3f}")

Final attack rate (simulated): 0.928

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
Cell In[5], line 42
     38 plt.show()
     39 
     40 final_AR = (I_d[-1] + R_d[-1]) / N
     41 print(f"Final attack rate (simulated): {final_AR:.3f}")
---> 42 print(f"Theoretical attack rate      : {solve_final_size(R0_demo, s0=(N-I0)/N):.3f}")

Cell In[2], line 22, in solve_final_size(R0, s0)
     18     # f(s) = 0 at the non-trivial root  (s < s0)
     19     f = lambda s: s - np.exp(-R0 * (s0 - s))
     20 
     21     # Root lies in (0, s0); use a tiny lower bound to avoid the trivial root
---> 22     s_inf = brentq(f, 1e-9, s0 - 1e-9)
     23     return 1.0 - s_inf

File ~/projects/laser-fresh/laser.cohorts/.docs/lib/python3.12/site-packages/scipy/optimize/_zeros_py.py:846, in brentq(f, a, b, args, xtol, rtol, maxiter, full_output, disp)
    844     raise ValueError(f"rtol too small ({rtol:g} < {_rtol:g})")
    845 f = _wrap_nan_raise(f)
--> 846 r = _zeros._brentq(f, a, b, xtol, rtol, maxiter, args, full_output, disp)
    847 return results_c(full_output, r, "brentq")

ValueError: f(a) and f(b) must have different signs

In [ ]: