import matplotlib.pyplot as plt
import numpy as np
from scipy.integrate import odeint
plt.style.use('custom.mplstyle')
%config InlineBackend.figure_format = 'retina'
SIR models describe the spread of a pathogen in a population. The SIR model breaks down the overall population into three subsets:
The system of differential equations for this model are:
$$ \begin{align*} \frac{dS(t)}{dt} &= -\beta S(t) I(t) \\ \frac{dI(t)}{dt} &= \beta S(t) I(t) - \nu I(t) \\ \frac{dR(t)}{dt} &= \nu I(t) \end{align*} $$where $\beta > 0$ is the rate at which a susceptible person meeting an infected individual becomes infected. Note that such encounters occur at rate $S(t) I(t)$. $\nu > 0$ is the rate at which infectious people recover. Here, we are letting $S$, $I$, and $R$ be the fractions of the population which are susceptible, infectious, or recovered, respectively: $S(t) + I(t) + R(t) = 1$. Because of this relationship, we consider the equation for $dR(t)/dt$ to be superfluous, as at each time point: $R(t) = 1- S(t) -I(t)$. Additionally, by modeling the fractions of the overall population, we don't need to worry about what the size of the total population size, making our lives a bit easier.
The ratio of $\beta$ and $\nu$ constitutes what is called $R_0$, the basic reproduction number:
$$ R_0 = \frac{\beta}{\nu} $$$R_0$ describes how many people an infectious person infects on average; you probably have heard about $R_0$ in the context of SARS-CoV-2 infections.
your answer:
For a given pathogen, the exact value of $R_0$ can change by changing the encounter rate within a popualtion. I.e., if we limit encounters, the basic reproduction number goes down. That is why lockdowns have been used to reduce infections.
Different pathogens can also have different intrinsic pathogenicity and transmissibilities. For example, the Omicron variant of SARS-CoV-2 is highly transmissable (i.e., a relatively large $R_0$).
Estimating $R_0$ is not straightforward and can be sensitive to data sampling. Here are some of the $R_0$ estimates for various pathogens:
Note that $dS(t)/dt$ and $dI(t)/dt$ are coupled in a nonlinear way due to the presence of $\beta S(t) I(t)$. Systems of differential equations coupled nonlinearly rarely have analytical solutions. We can play the game of looking at how the populations evolve in limiting cases which make the solution tractable.
your answer:
your answer:
The solutions from Questions 2,3 are accurate only in the limiting regimes. Here, we'll solve the equations irrespective of any regime using numerical integration. ___
Because $R(t) = 1 - S(t) - I(t)$, we need only solve for the evolution of $S(t)$ and $I(t)$. As before (see tutorial 2),
$$ \begin{align*} S(t + \Delta t) &= S(t) + \frac{dS(t)}{dt} \Delta t + O(\Delta t^2) \\ I(t + \Delta t) &= I(t) + \frac{dI(t)}{dt} \Delta t + O(\Delta t^2) \end{align*} $$Applying the approximation and plugging in the values for the derivatives, we have
$$ \begin{align*} S(t + \Delta t) &\approx S(t) - \beta S(t) I(t) \Delta t \\ I(t + \Delta t) &\approx I(t) + \beta S(t) I(t) \Delta t - \nu I(t) \Delta t \end{align*} $$def euler_sir(s0, i0, r0, nu, T, dt):
# Your code here.
pass
pass
pass
We stated above the Euler's method employs a first-order Taylor series approximation to solve differential equations. The SIR model is seemingly simple, but it is much more complicated than solving the differential equation for exponential growth. Let's now use scipy
which has functions to integrate differential equations more accurately. Specifically let's use scipy.integrate.odeint. (We don't concern ourselves with the implementation here. Though for those who are curious, this function uses multistep methods whereas Euler's method uses only one step. So this scipy
function should be more accurate.)
Reading the documentation, we see an example that is illuminating and follow it. We need to define a function which outputs an array-like object which contains the derivatives.
def sir_ode(current_vals, T, beta, nu):
pass
scipy
ode integrator using your function, obtaining trajectories for $S$ and $I$ using the following parameters:¶np.linspace(0, 100, 101)
(specifies that S(t), I(t) is evaluated at $t \in [0, 100]$ (days) $ \subset \mathbb{Z}$)pass
scipy
as solid lines. When plotting the results from using the Euler method, plot every 100th point since we had $\Delta t = 0.01$. Set the size of the markers of the scatter to be 5: s=5
. Label all plotted objects. Because we now have six plotted objects, the legend is large. We move the legend outside the axes: plt.legend(loc='center right', bbox_to_anchor=(1.3, 0.5))
.¶pass
We see that, at least for this choice of parameters, Euler's method produces an approximation which matches the results given by the scipy
integrator (at least visually). Let's continue to understand the trajectories of $S(t)$ and $I(t)$ using the function from scipy
since it's much quicker and more accurate in principle.
Let's first find the relation between $R_0$, $I(t)$ and $S(t)$. We compute $dI/dS$ to see how $R_0$ appears in our model naturally.
$$ \begin{align*} \frac{dI(t)}{dS(t)} &= \frac{dI(t)}{dt} \left( \frac{dS(t)}{dt} \right)^{-1} \\ &= \frac{\beta S(t) I(t) - \nu I(t)}{-\beta S(t) I(t)} \\ &= -1 + \frac{1}{R_0 S(t)} \end{align*} $$This is a simple ordinary differential equation which we can solve using separation of variables.
$$ I(t) = -S(t) + \frac{\ln S(t)}{R_0} + C \quad \forall t $$where $C$ is some constant of integration found from initial conditions.
your answer:
your answer:
Essentially, this shows that we can determine the final state of an outbreak based on $R_0$ and the initial conditions.
def stationary_s(s0, i0, r0):
pass
np.logspace(-1, 1, 10)
pass
np.logspace(-1, 3, 10)
.¶pass
Herd immunity is the idea that a disease cannot grow into an epidemic when many people in a population are immune and not susceptible. The condition for herd immunity is obtained when
$$ \frac{dI(0)}{dt} < 0 $$Since
$$ \frac{dI(0)}{dt} = (\beta S(0) - \nu) I(0) $$we see that
$$ S(0) < 1/R_0 $$is the condition for herd immunity.
your answer: