Simulating the SIR Model

Module 1: Low-Dimensional Dynamics

Reading

Compulsory: (castro2020turning?).
Course material: Compartmental Models, The Reproduction Number, The Final Size of an Epidemic.

Core activity

Simulate SIR and SIS (this page and the next), fit the SIR model to incidence data, and test how well the peak can be predicted from early data.

The SIR model is a small system of three ODEs (Kermack and McKendrick 1927; Keeling and Rohani 2008). The solve_ivp workflow from Solving Initial Value Problems applies unchanged: the state is now a vector.

\[ \dot S = -\beta SI, \qquad \dot I = \beta SI - \gamma I, \qquad \dot R = \gamma I . \]

The state has three components, but $S + I + R = 1$ is conserved, so only two are independent.

Prediction

Exercise

For $\beta = 0.3$, $\gamma = 0.1$, $S(0) = 0.99$, $I(0) = 0.01$, use the formulas from The Reproduction Number and The Final Size of an Epidemic to predict the peak prevalence $I_{\max}$ and the final size before you simulate.

Solution

Here $R_0 = 3$. The peak prevalence is $I_{\max} = I_0 + S_0 - \frac{1}{R_0}\bigl(1 + \ln(R_0 S_0)\bigr) \approx 0.304$, reached when $S = 1/R_0 \approx 0.333$. The final size relation $\ln(S_0/S_\infty) = R_0 (1 - S_\infty)$ gives $S_\infty \approx 0.059$, so about 93% of the susceptibles are infected and $R_\infty \approx 0.94$.

from amlab.epidemics import final_size, peak_prevalence

print(f"I_max = {peak_prevalence(3.0, 0.99, 0.01):.3f}")
print(f"final size = {final_size(3.0, 0.99):.3f}")

I_max = 0.304
final size = 0.931

Implementation

Exercise

Complete the right-hand side of the SIR model. The solver only needs one function that returns the three derivatives.

def sir_model(t, y, beta, gamma):
    s, i, r = y
    dsdt = # TODO
    didt = # TODO
    drdt = # TODO
    return [dsdt, didt, drdt]

Solution

def sir_model(t, y, beta, gamma):
    s, i, r = y
    dsdt = -beta * s * i
    didt = beta * s * i - gamma * i
    drdt = gamma * i
    return [dsdt, didt, drdt]

The solver call follows the usual pattern:

solution = solve_ivp(
    sir_model,
    t_span=(0, 160),
    y0=[0.99, 0.01, 0.0],
    args=(0.3, 0.1),
    t_eval=np.linspace(0, 160, 1000),
)

The reference implementation amlab/odes_1d/sir_model.py provides sir_model and the plotting helper plot_sir_model.

Figure 1: SIR dynamics for beta = 0.3, gamma = 0.1 (R0 = 3), S(0) = 0.99, I(0) = 0.01.

Validation

Check three facts against the theory: $S + I + R$ stays equal to 1, the peak of $I$ matches the formula, and the peak happens when $S = 1/R_0$.

import numpy as np
from scipy.integrate import solve_ivp

from amlab.epidemics import peak_prevalence
from amlab.odes_1d.sir_model import sir_model

sol = solve_ivp(sir_model, (0, 160), [0.99, 0.01, 0.0], args=(0.3, 0.1),
                t_eval=np.linspace(0, 160, 4001), rtol=1e-9, atol=1e-12)
s, i, r = sol.y
print(f"max |S+I+R-1| = {np.max(np.abs(s + i + r - 1)):.1e}")
print(f"I_max simulated = {i.max():.4f}, formula = {peak_prevalence(3.0, 0.99, 0.01):.4f}")
print(f"S at the peak = {s[i.argmax()]:.4f}, 1/R0 = {1 / 3:.4f}")

max |S+I+R-1| = 6.7e-16
I_max simulated = 0.3038, formula = 0.3038
S at the peak = 0.3331, 1/R0 = 0.3333

What to look for

The susceptible fraction decreases monotonically.
The infected fraction rises, peaks and then falls.
The recovered fraction grows monotonically.
Changing $\beta$ mostly affects how fast the outbreak grows.
Changing $\gamma$ mostly affects how quickly the outbreak resolves.

Experiment

Exercise

Start from $\beta = 0.3$, $\gamma = 0.1$, $I(0) = 0.01$ and change one thing at a time.

Increase $\beta$ with $\gamma$ fixed. What happens to the height and the time of the peak?
Increase $\gamma$ with $\beta$ fixed. Does the epidemic end sooner? Is it smaller?
Change $I(0)$ from $10^{-2}$ to $10^{-5}$. Which quantities change: peak height, peak time, final size?

Solution

With $\beta = 0.45$ ($R_0 = 4.5$) the peak rises from 0.30 to 0.45 and comes earlier, near day 17 instead of day 27. With $\gamma = 0.15$ ($R_0 = 2$) the peak drops to 0.16 and the final size to 0.80; the epidemic is smaller, and it dies out a little sooner. With $I(0) = 10^{-5}$ the peak height (0.30) and the final size (0.94) barely change, but the peak moves from day 27 to about day 61.

for beta, gamma, i0 in [(0.3, 0.1, 1e-2), (0.45, 0.1, 1e-2), (0.3, 0.15, 1e-2), (0.3, 0.1, 1e-5)]:
    run = solve_ivp(sir_model, (0, 300), [1 - i0, i0, 0.0], args=(beta, gamma),
                    t_eval=np.linspace(0, 300, 30001), rtol=1e-9, atol=1e-12)
    i_run = run.y[1]
    print(f"beta={beta}, gamma={gamma}, I0={i0:.0e}: peak {i_run.max():.2f} "
          f"on day {run.t[i_run.argmax()]:.0f}, final R {run.y[2, -1]:.2f}")

beta=0.3, gamma=0.1, I0=1e-02: peak 0.30 on day 27, final R 0.94
beta=0.45, gamma=0.1, I0=1e-02: peak 0.45 on day 17, final R 0.99
beta=0.3, gamma=0.15, I0=1e-02: peak 0.16 on day 29, final R 0.80
beta=0.3, gamma=0.1, I0=1e-05: peak 0.30 on day 61, final R 0.94

Interpretation

$\beta$ and $\gamma$ enter the shape of the outbreak mainly through $R_0 = \beta/\gamma$ (peak height, final size) and $\beta - \gamma$ (early growth rate). The initial seed $I(0)$ shifts the curve in time but hardly changes its shape.

Optional extension

Add Streamlit sliders for $\beta$, $\gamma$ and $I(0)$, following Budworm: a Streamlit App. Show the peak prevalence formula next to the plot.

Run locally

To run the standalone script and show the plot:

python amlab/odes_1d/sir_model.py

Keeling, Matt J., and Pejman Rohani. 2008. Modeling Infectious Diseases in Humans and Animals. Princeton University Press.

Kermack, William O., and Anderson G. McKendrick. 1927. “A Contribution to the Mathematical Theory of Epidemics.” Proceedings of the Royal Society A 115 (772): 700–721. https://doi.org/10.1098/rspa.1927.0118.

--- title: "Simulating the SIR Model" subtitle: "Module 1: Low-Dimensional Dynamics" format: html --- ::: {.callout-note} ## Reading - **Compulsory:** @castro2020turning. - **Course material:** [Compartmental Models](compartmental-models.qmd), [The Reproduction Number](reproduction-number.qmd), [The Final Size of an Epidemic](final-size.qmd). ::: ::: {.callout-important} ## Core activity Simulate SIR and SIS (this page and the next), fit the SIR model to incidence data, and test how well the peak can be predicted from early data. ::: The SIR model is a small system of three ODEs [@kermack1927contribution; @keeling2008modeling]. The `solve_ivp` workflow from [Solving Initial Value Problems](solve-ivp.qmd) applies unchanged: the state is now a vector. $$ \dot S = -\beta SI, \qquad \dot I = \beta SI - \gamma I, \qquad \dot R = \gamma I . $$ The state has three components, but $S + I + R = 1$ is conserved, so only two are independent. ## Prediction ::: {.callout-note} ## Exercise For $\beta = 0.3$, $\gamma = 0.1$, $S(0) = 0.99$, $I(0) = 0.01$, use the formulas from [The Reproduction Number](reproduction-number.qmd) and [The Final Size of an Epidemic](final-size.qmd) to predict the peak prevalence $I_{\max}$ and the final size before you simulate. ::: ::: {.callout-tip collapse="true"} ## Solution Here $R_0 = 3$. The peak prevalence is $I_{\max} = I_0 + S_0 - \frac{1}{R_0}\bigl(1 + \ln(R_0 S_0)\bigr) \approx 0.304$, reached when $S = 1/R_0 \approx 0.333$. The final size relation $\ln(S_0/S_\infty) = R_0 (1 - S_\infty)$ gives $S_\infty \approx 0.059$, so about 93% of the susceptibles are infected and $R_\infty \approx 0.94$. ```{python} from amlab.epidemics import final_size, peak_prevalence print(f"I_max = {peak_prevalence(3.0, 0.99, 0.01):.3f}") print(f"final size = {final_size(3.0, 0.99):.3f}") ``` ::: ## Implementation ::: {.callout-note} ## Exercise Complete the right-hand side of the SIR model. The solver only needs one function that returns the three derivatives. ```python def sir_model(t, y, beta, gamma): s, i, r = y dsdt = # TODO didt = # TODO drdt = # TODO return [dsdt, didt, drdt] ``` ::: ::: {.callout-tip collapse="true"} ## Solution ```python def sir_model(t, y, beta, gamma): s, i, r = y dsdt = -beta * s * i didt = beta * s * i - gamma * i drdt = gamma * i return [dsdt, didt, drdt] ``` ::: The solver call follows the usual pattern: ```python solution = solve_ivp( sir_model, t_span=(0, 160), y0=[0.99, 0.01, 0.0], args=(0.3, 0.1), t_eval=np.linspace(0, 160, 1000), ) ``` The reference implementation `amlab/odes_1d/sir_model.py` provides `sir_model` and the plotting helper `plot_sir_model`. ```{python} #| label: fig-sir #| fig-cap: "SIR dynamics for beta = 0.3, gamma = 0.1 (R0 = 3), S(0) = 0.99, I(0) = 0.01." #| fig-width: 8 #| fig-height: 4 #| echo: false import matplotlib.pyplot as plt from amlab.odes_1d.sir_model import plot_sir_model plot_sir_model(beta=0.3, gamma=0.1, s0=0.99, i0=0.01, r0=0.0, t_show=160) plt.show() plt.close() ``` ## Validation Check three facts against the theory: $S + I + R$ stays equal to 1, the peak of $I$ matches the formula, and the peak happens when $S = 1/R_0$. ```{python} import numpy as np from scipy.integrate import solve_ivp from amlab.epidemics import peak_prevalence from amlab.odes_1d.sir_model import sir_model sol = solve_ivp(sir_model, (0, 160), [0.99, 0.01, 0.0], args=(0.3, 0.1), t_eval=np.linspace(0, 160, 4001), rtol=1e-9, atol=1e-12) s, i, r = sol.y print(f"max |S+I+R-1| = {np.max(np.abs(s + i + r - 1)):.1e}") print(f"I_max simulated = {i.max():.4f}, formula = {peak_prevalence(3.0, 0.99, 0.01):.4f}") print(f"S at the peak = {s[i.argmax()]:.4f}, 1/R0 = {1 / 3:.4f}") ``` ## What to look for - The susceptible fraction decreases monotonically. - The infected fraction rises, peaks and then falls. - The recovered fraction grows monotonically. - Changing $\beta$ mostly affects how fast the outbreak grows. - Changing $\gamma$ mostly affects how quickly the outbreak resolves. ## Experiment ::: {.callout-note} ## Exercise Start from $\beta = 0.3$, $\gamma = 0.1$, $I(0) = 0.01$ and change one thing at a time. - Increase $\beta$ with $\gamma$ fixed. What happens to the height and the time of the peak? - Increase $\gamma$ with $\beta$ fixed. Does the epidemic end sooner? Is it smaller? - Change $I(0)$ from $10^{-2}$ to $10^{-5}$. Which quantities change: peak height, peak time, final size? ::: ::: {.callout-tip collapse="true"} ## Solution With $\beta = 0.45$ ($R_0 = 4.5$) the peak rises from 0.30 to 0.45 and comes earlier, near day 17 instead of day 27. With $\gamma = 0.15$ ($R_0 = 2$) the peak drops to 0.16 and the final size to 0.80; the epidemic is smaller, and it dies out a little sooner. With $I(0) = 10^{-5}$ the peak height (0.30) and the final size (0.94) barely change, but the peak moves from day 27 to about day 61. ```{python} for beta, gamma, i0 in [(0.3, 0.1, 1e-2), (0.45, 0.1, 1e-2), (0.3, 0.15, 1e-2), (0.3, 0.1, 1e-5)]: run = solve_ivp(sir_model, (0, 300), [1 - i0, i0, 0.0], args=(beta, gamma), t_eval=np.linspace(0, 300, 30001), rtol=1e-9, atol=1e-12) i_run = run.y[1] print(f"beta={beta}, gamma={gamma}, I0={i0:.0e}: peak {i_run.max():.2f} " f"on day {run.t[i_run.argmax()]:.0f}, final R {run.y[2, -1]:.2f}") ``` ::: ## Interpretation $\beta$ and $\gamma$ enter the shape of the outbreak mainly through $R_0 = \beta/\gamma$ (peak height, final size) and $\beta - \gamma$ (early growth rate). The initial seed $I(0)$ shifts the curve in time but hardly changes its shape. ::: {.callout-tip} ## Optional extension Add Streamlit sliders for $\beta$, $\gamma$ and $I(0)$, following [Budworm: a Streamlit App](budworm-streamlit.qmd). Show the peak prevalence formula next to the plot. ::: ## Run locally To run the standalone script and show the plot: ```bash python amlab/odes_1d/sir_model.py ``` ## Next