None
This Jupyter/Python notebook presents models for the outbreak of an infectious disease into a susceptible population using standard epidemiological models. Model parameters are taken from a rapidly evolving scientific literature documenting the global COVID-19 outbreak. A control policy based on 'social distancing' is included in the model.
The notebook is organized as follows:
The executable Python code in this notebook can be edited and run Google's cloud servers by clicking on the "Open in Colab" button located in the header. Use the interactive sliders to adjust model parameters and perform 'what if' simulations.
COVID-19 is caused by the human coronavirus SARS-CoV-2. First identified in the 1960's, there are currently four human coronaviruses endemic to populations around the world:
These four common coronaviruses cause an upper respiratory disease that can progress to pneumonia. These endemic viruses cause about a quarter of all common colds. Most people will suffer from at least one during their lifetimes.
In recent decades, three additional coronaviruses that normally infect animals have evolved to infect humans. These include:
The last of these, now called SARS-CoV-2, first appeared in December, 2019, at a seafood market in Wuhan (Hubei, China). The rapid spread of SARS-CoV-2 in Wuhan, and subsequent appearance in other locations around the globe, has resulted in declaration of a global health emergency by the World Health Organization (WHO). Most countries are mobilizing to track the virus and control new outbreaks. At this stage, it is too early to know if efforts to contain the mitigate transmission of the virus will be successful in preventing COVID-19 from becoming a pandemic, or later an endemic disease with a global footprint.
The latest status on the global outbreak of COVID-19 can be found at the following links:
The purpose of this notebook is to demonstrate the modeling of an infectious epidemic using the latest available data for COVID-19, and to provide a framework for evaluating the performance of 'social distancing' and other mitigation strategies. The models and data used in this notebook have been extracted from a rapidly emerging and changing literature. Recent papers on COVID-19 can be found at the following links.
The SIR model is deterministic compartment model for the spread of an infectious disease that describes key phenomena encountered in epidemiology. In the SIR model, a population is broken into three non-overlapping groups corresponding to stages of the disease:
Neglecting demographic processes of birth and death from other causes, and assuming a negligible death rate due to infectious disease at issue, the progression of an epidemic can be modeled by rate processesl
$$\text{Susceptible} \xrightarrow{\frac{\beta S I}{N}} \text{Infectious} \xrightarrow{\gamma I} \text{Recovered} $$The rate processes are modeled as follows.
A model for the spread of an infectious disease in a uniform population is given by the deterministic SIR equations
\begin{align*} \frac{dS}{dt} & = -\frac{\beta S I}{N} \\ \frac{dI}{dt} & = \frac{\beta S I}{N} - \gamma I \\ \frac{dR}{dt} & = \gamma I \end{align*}The model becomes more generic by working with population fractions rather than raw population counts. To this end, define
\begin{align} s = \frac{S}{N} \qquad i = \frac{I}{N} \qquad r = \frac{R}{N} \end{align}After substitution, this results in a system of four equations.
\begin{align*} \frac{ds}{dt} & = -\beta s i \\ \frac{di}{dt} & = \beta s i - \gamma i \\ \frac{dr}{dt} & = \gamma i \end{align*}where $s + i + r = 1$ is an invariant.
The SIR model describes key epidemiological phenemena. Here is a brief synposis of the relevant results.
The following Python code implements a simulation of the SIR model. The parameter values were selected from the recent survey by Boldog, et al. (2020).
import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt
# parameter values
R0 = 2.4
t_infective = 5.1 + 3.3
# initial number of infected and recovered individuals
i_initial = 1/20000
r_initial = 0.00
s_initial = 1 - i_initial - r_initial
gamma = 1/t_infective
beta = R0*gamma
# SIR model differential equations.
def deriv(x, t, beta, gamma):
s, i, r = x
dsdt = -beta * s * i
didt = beta * s * i - gamma * i
drdt = gamma * i
return [dsdt, didt, drdt]
t = np.linspace(0, 180, 2000)
x_initial = s_initial, i_initial, r_initial
soln = odeint(deriv, x_initial, t, args=(beta, gamma))
s, i, r = soln.T
e = None
def plotdata(t, s, i, e=None):
# plot the data
fig = plt.figure(figsize=(12,6))
ax = [fig.add_subplot(221, axisbelow=True),
fig.add_subplot(223),
fig.add_subplot(122)]
ax[0].plot(t, s, lw=3, label='Fraction Susceptible')
ax[0].plot(t, i, lw=3, label='Fraction Infective')
ax[0].plot(t, r, lw=3, label='Recovered')
ax[0].set_title('Susceptible and Recovered Populations')
ax[0].set_xlabel('Time /days')
ax[0].set_ylabel('Fraction')
ax[1].plot(t, i, lw=3, label='Infective')
ax[1].set_title('Infectious Population')
if e is not None: ax[1].plot(t, e, lw=3, label='Exposed')
ax[1].set_ylim(0, 0.3)
ax[1].set_xlabel('Time /days')
ax[1].set_ylabel('Fraction')
ax[2].plot(s, i, lw=3, label='s, i trajectory')
ax[2].plot([1/R0, 1/R0], [0, 1], '--', lw=3, label='di/dt = 0')
ax[2].plot(s[0], i[0], '.', ms=20, label='Initial Condition')
ax[2].plot(s[-1], i[-1], '.', ms=20, label='Final Condition')
ax[2].set_title('State Trajectory')
ax[2].set_aspect('equal')
ax[2].set_ylim(0, 1.05)
ax[2].set_xlim(0, 1.05)
ax[2].set_xlabel('Susceptible')
ax[2].set_ylabel('Infectious')
for a in ax:
a.grid(True)
a.legend()
plt.tight_layout()
plotdata(t, s, i)
Given an outbreak in a susceptible population, the final state is reached when $i$ returns to a near zero value. A formula for the final state can be found by taking the ratio
\begin{align} \frac{di}{ds} & = \frac{\frac{di}{dt}}{\frac{ds}{dt}} = -1 + \frac{1}{s R_0} & \\ \end{align}Integrating,
\begin{align} \int_{i_0}^{i_f} di & = \int_{s_0}^{s_f} (-1 + \frac{1}{s R_0}) ds & \\ \end{align}\begin{align} i_f - i_0 & = -(s_f - s_0) + \frac{1}{R_0} \ln\frac{s_f}{s_0} & \\ \end{align}The starting point of an outbreak begins with a very small value $i_0 \approx 0$ and ends with a very small value $i_f \approx 0$. Setting $i_0 = i_f = 0$ gives
\begin{align} s_f - \frac{1}{R_0}\ln s_f & = s_0 - \frac{1}{R_0} \ln{s_0} & \\ \end{align}For the special case of an initially susceptible population, $s_0 = 1$ which gives
\begin{align} s_f - \frac{1}{R_0} \ln s_f & = 1 \end{align}import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import brentq, fminbound
def f(s, R0):
return s - np.log(s)/R0
def g(s0, R0):
rhs = f(s0, R0)
smin = fminbound(lambda s: f(s, R0), 0.0001, 1)
return brentq(lambda s: f(s, R0) - f(s0, R0), 0.001, smin)
R0 = np.linspace(1.01, 4.0, 100)
sf = [g(0.8, R0) for R0 in R0]
plt.plot(R0, sf)
[<matplotlib.lines.Line2D at 0x11b6c2050>]
The SEIR model extends the SIR model by adding an additional population compartment containing those individuals who have been exposed to the virus but not yet infective.
The compartment model can be diagrammed as follows.
$$\text{Susceptible} \xrightarrow{\frac{\beta S I}{N}} \text{Exposed} \xrightarrow{\alpha E} \text{Infectious} \xrightarrow{\gamma I} \text{Recovered} $$The rate processes are modeled as follows.
An elementary model for the spread of an infectious disease in a uniform population is given by the deterministic SEIR equations}
After substitution, this results in a system of four equations.
\begin{align*} \frac{ds}{dt} & = -\beta s i \\ \frac{de}{dt} & = \beta s i - \alpha e \\ \frac{di}{dt} & = \alpha e - \gamma i \\ \frac{dr}{dt} & = \gamma i \end{align*}where $s + e + i + r = 1$ is an invariant.
import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt
# parameter values
R0 = 4
t_incubation = 5.1
t_infective = 3.3
# initial number of infected and recovered individuals
e_initial = 1/20000
i_initial = 0.00
r_initial = 0.00
s_initial = 1 - e_initial - i_initial - r_initial
alpha = 1/t_incubation
gamma = 1/t_infective
beta = R0*gamma
# SEIR model differential equations.
def deriv(x, t, alpha, beta, gamma):
s, e, i, r = x
dsdt = -beta * s * i
dedt = beta * s * i - alpha * e
didt = alpha * e - gamma * i
drdt = gamma * i
return [dsdt, dedt, didt, drdt]
t = np.linspace(0, 160, 160)
x_initial = s_initial, e_initial, i_initial, r_initial
soln = odeint(deriv, x_initial, t, args=(alpha, beta, gamma))
s, e, i, r = soln.T
plotdata(t, s, i, e)
The addition of an exposed population compartment slows the outbreak, but doesn't appear to reduce the number of people ultimately infected by the disease.
What are the campus policy implications of these results?
The lack of a vaccine reduces the options for controlling the COVID-19 outbreak. Current efforts are focused on 'social distancing' designed to reduce transmission of the virus from individuals in the infective state to susceptible individuals.
For the purposes of modeling, we introduce a control parameter $u$ indicating the effectiveness of these efforts. $u=0$ corresponds to no controls, $u=1$ corresponds to perfect isolation of infective individuals. The purpose of this model is to explore how a social distancing stragtegy affects the outcome of an epidemic.
The compartment model can be diagrammed as follows.
$$\text{Susceptible} \xrightarrow{(1-u)\frac{\beta S I}{N}} \text{Exposed} \xrightarrow{\alpha E} \text{Infectious} \xrightarrow{\gamma I} \text{Recovered} $$The rate processes are modeled as follows.
After substitution, this results in a system of four equations.
\begin{align*} \frac{ds}{dt} & = -(1-u)\beta s i \\ \frac{de}{dt} & = (1-u)\beta s i - \alpha e \\ \frac{di}{dt} & = \alpha e - \gamma i \\ \frac{dr}{dt} & = \gamma i \end{align*}where $s + e + i + r = 1$ is an invariant.
import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt
# parameter values
u = 0.2
R0 = 2.4
t_incubation = 5.1
t_infective = 3.3
# initial number of infected and recovered individuals
e_initial = 1/20000
i_initial = 0.00
r_initial = 0.00
s_initial = 1 - e_initial - i_initial - r_initial
alpha = 1/t_incubation
gamma = 1/t_infective
beta = R0*gamma
# SEIR model differential equations.
def deriv(x, t, alpha, beta, gamma):
s, e, i, r = x
dsdt = -(1-u)*beta * s * i
dedt = (1-u)*beta * s * i - alpha * e
didt = alpha * e - gamma * i
drdt = gamma * i
return [dsdt, dedt, didt, drdt]
t = np.linspace(0, 160, 160)
x_initial = s_initial, e_initial, i_initial, r_initial
soln = odeint(deriv, x_initial, t, args=(alpha, beta, gamma))
s, e, i, r = soln.T
plotdata(t, s, e, i)
Social distancing has several beneficial effects:
The basic strategy is to slow transmission through 'social distancing' with the following goals:
Boldog, et al.
from scipy.integrate import odeint
import numpy as np
import matplotlib.pyplot as plt
u = 0.3
mu = 0
alpha = 1/5.1 # incubation period
R0 = 2.4
gamma = 3.3
beta = R0*gamma
N = 331000000
def SEIR(x, t):
S, E1, E2, I1, I2, I3, R = x
dS = -(1-u)*beta*S*(I1 + I2 + I3)/N
dE1 = -dS - 2*alpha*E1
dE2 = 2*alpha*E1 - 2*alpha*E2
dI1 = 2*alpha*E2 - 3*gamma*I1 - mu*I1
dI2 = 3*gamma*I1 - 3*gamma*I2 - mu*I2
dI3 = 3*gamma*I2 - 3*gamma*I3 - mu*I3
dR = 3*gamma*I3
return [dS, dE1, dE2, dI1, dI2, dI3, dR]
IC = [N, 1, 0, 0, 0, 0, 0]
t = np.linspace(0, 180, 1000)
soln = odeint(SEIR, IC, t)
s = soln[:, 0]
e = soln[:, 1] + soln[:, 2]
i = soln[:, 3] + soln[:, 4] + soln[:, 5]
r = soln[:, 6]
plotdata(t, s, i, e)
Ziff, Robert M., and Anna L. Ziff. "Fractal kinetics of COVID-19 pandemic." medRxiv (2020). https://www.medrxiv.org/content/10.1101/2020.02.16.20023820v1
Peng, Liangrong, et al. "Epidemic analysis of COVID-19 in China by dynamical modeling." arXiv preprint arXiv:2002.06563 (2020). https://arxiv.org/abs/2002.06563
Wu, Joseph T., Kathy Leung, and Gabriel M. Leung. "Nowcasting and forecasting the potential domestic and international spread of the 2019-nCoV outbreak originating in Wuhan, China: a modelling study." The Lancet (2020). https://www.thelancet.com/journals/lancet/article/PIIS0140-6736(20)30260-9/fulltext
Tuite, Ashleigh R., and David N. Fisman. "Reporting, Epidemic Growth, and Reproduction Numbers for the 2019 Novel Coronavirus (2019-nCoV) Epidemic." Annals of Internal Medicine (2020). https://annals.org/aim/fullarticle/2760912