{
"cells": [
{
"cell_type": "markdown",
"id": "d0b08398",
"metadata": {},
"source": [
"# Modeling COVID 19"
]
},
{
"cell_type": "markdown",
"id": "4b9f4f8a",
"metadata": {},
"source": [
"## Contents\n",
"\n",
"- [Modeling COVID 19](#Modeling-COVID-19) \n",
" - [Overview](#Overview) \n",
" - [The SIR Model](#The-SIR-Model) \n",
" - [Implementation](#Implementation) \n",
" - [Experiments](#Experiments) \n",
" - [Ending Lockdown](#Ending-Lockdown) "
]
},
{
"cell_type": "markdown",
"id": "6bc9d150",
"metadata": {},
"source": [
"## Overview\n",
"\n",
"This is a Python version of the code for analyzing the COVID-19 pandemic\n",
"provided by [Andrew Atkeson](https://sites.google.com/site/andyatkeson/).\n",
"\n",
"See, in particular\n",
"\n",
"- [NBER Working Paper No. 26867](https://www.nber.org/papers/w26867) \n",
"- [COVID-19 Working papers and code](https://sites.google.com/site/andyatkeson/home?authuser=0) \n",
"\n",
"\n",
"The purpose of his notes is to introduce economists to quantitative modeling\n",
"of infectious disease dynamics.\n",
"\n",
"Dynamics are modeled using a standard SIR (Susceptible-Infected-Removed) model\n",
"of disease spread.\n",
"\n",
"The model dynamics are represented by a system of ordinary differential\n",
"equations.\n",
"\n",
"The main objective is to study the impact of suppression through social\n",
"distancing on the spread of the infection.\n",
"\n",
"The focus is on US outcomes but the parameters can be adjusted to study\n",
"other countries.\n",
"\n",
"We will use the following standard imports:"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "acdabba3",
"metadata": {
"hide-output": false
},
"outputs": [],
"source": [
"import matplotlib.pyplot as plt\n",
"plt.rcParams[\"figure.figsize\"] = (11, 5) #set default figure size\n",
"import numpy as np\n",
"from numpy import exp"
]
},
{
"cell_type": "markdown",
"id": "5f45273f",
"metadata": {},
"source": [
"We will also use SciPy’s numerical routine odeint for solving differential\n",
"equations."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "d88851a9",
"metadata": {
"hide-output": false
},
"outputs": [],
"source": [
"from scipy.integrate import odeint"
]
},
{
"cell_type": "markdown",
"id": "d3479d2b",
"metadata": {},
"source": [
"This routine calls into compiled code from the FORTRAN library odepack."
]
},
{
"cell_type": "markdown",
"id": "09d09294",
"metadata": {},
"source": [
"## The SIR Model\n",
"\n",
"In the version of the SIR model we will analyze there are four states.\n",
"\n",
"All individuals in the population are assumed to be in one of these four states.\n",
"\n",
"The states are: susceptible (S), exposed (E), infected (I) and removed ®.\n",
"\n",
"Comments:\n",
"\n",
"- Those in state R have been infected and either recovered or died. \n",
"- Those who have recovered are assumed to have acquired immunity. \n",
"- Those in the exposed group are not yet infectious. "
]
},
{
"cell_type": "markdown",
"id": "3dd38a3d",
"metadata": {},
"source": [
"### Time Path\n",
"\n",
"The flow across states follows the path $ S \\to E \\to I \\to R $.\n",
"\n",
"All individuals in the population are eventually infected when\n",
"the transmission rate is positive and $ i(0) > 0 $.\n",
"\n",
"The interest is primarily in\n",
"\n",
"- the number of infections at a given time (which determines whether or not the health care system is overwhelmed) and \n",
"- how long the caseload can be deferred (hopefully until a vaccine arrives) \n",
"\n",
"\n",
"Using lower case letters for the fraction of the population in each state, the\n",
"dynamics are\n",
"\n",
"\n",
"\n",
"$$\n",
"\\begin{aligned}\n",
" \\dot s(t) & = - \\beta(t) \\, s(t) \\, i(t)\n",
" \\\\\n",
" \\dot e(t) & = \\beta(t) \\, s(t) \\, i(t) - σ e(t)\n",
" \\\\\n",
" \\dot i(t) & = σ e(t) - γ i(t)\n",
"\\end{aligned} \\tag{1.1}\n",
"$$\n",
"\n",
"In these equations,\n",
"\n",
"- $ \\beta(t) $ is called the *transmission rate* (the rate at which individuals bump into others and expose them to the virus). \n",
"- $ \\sigma $ is called the *infection rate* (the rate at which those who are exposed become infected) \n",
"- $ \\gamma $ is called the *recovery rate* (the rate at which infected people recover or die). \n",
"- the dot symbol $ \\dot y $ represents the time derivative $ dy/dt $. \n",
"\n",
"\n",
"We do not need to model the fraction $ r $ of the population in state $ R $ separately because the states form a partition.\n",
"\n",
"In particular, the “removed” fraction of the population is $ r = 1 - s - e - i $.\n",
"\n",
"We will also track $ c = i + r $, which is the cumulative caseload\n",
"(i.e., all those who have or have had the infection).\n",
"\n",
"The system [(1.1)](#equation-sir-system) can be written in vector form as\n",
"\n",
"\n",
"\n",
"$$\n",
"\\dot x = F(x, t), \\qquad x := (s, e, i) \\tag{1.2}\n",
"$$\n",
"\n",
"for suitable definition of $ F $ (see the code below)."
]
},
{
"cell_type": "markdown",
"id": "cc58ea1c",
"metadata": {},
"source": [
"### Parameters\n",
"\n",
"Both $ \\sigma $ and $ \\gamma $ are thought of as fixed, biologically determined parameters.\n",
"\n",
"As in Atkeson’s note, we set\n",
"\n",
"- $ \\sigma = 1/5.2 $ to reflect an average incubation period of 5.2 days. \n",
"- $ \\gamma = 1/18 $ to match an average illness duration of 18 days. \n",
"\n",
"\n",
"The transmission rate is modeled as\n",
"\n",
"- $ \\beta(t) := R(t) \\gamma $ where $ R(t) $ is the *effective reproduction number* at time $ t $. \n",
"\n",
"\n",
"(The notation is slightly confusing, since $ R(t) $ is different to\n",
"$ R $, the symbol that represents the removed state.)"
]
},
{
"cell_type": "markdown",
"id": "44d35746",
"metadata": {},
"source": [
"## Implementation\n",
"\n",
"First we set the population size to match the US."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "b42ea9ee",
"metadata": {
"hide-output": false
},
"outputs": [],
"source": [
"pop_size = 3.3e8"
]
},
{
"cell_type": "markdown",
"id": "5d6b1d66",
"metadata": {},
"source": [
"Next we fix parameters as described above."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "c2d9b771",
"metadata": {
"hide-output": false
},
"outputs": [],
"source": [
"γ = 1 / 18\n",
"σ = 1 / 5.2"
]
},
{
"cell_type": "markdown",
"id": "5b98b21c",
"metadata": {},
"source": [
"Now we construct a function that represents $ F $ in [(1.2)](#equation-dfcv)"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "4b496ce2",
"metadata": {
"hide-output": false
},
"outputs": [],
"source": [
"def F(x, t, R0=1.6):\n",
" \"\"\"\n",
" Time derivative of the state vector.\n",
"\n",
" * x is the state vector (array_like)\n",
" * t is time (scalar)\n",
" * R0 is the effective transmission rate, defaulting to a constant\n",
"\n",
" \"\"\"\n",
" s, e, i = x\n",
"\n",
" # New exposure of susceptibles\n",
" β = R0(t) * γ if callable(R0) else R0 * γ\n",
" ne = β * s * i\n",
"\n",
" # Time derivatives\n",
" ds = - ne\n",
" de = ne - σ * e\n",
" di = σ * e - γ * i\n",
"\n",
" return ds, de, di"
]
},
{
"cell_type": "markdown",
"id": "7ada4dac",
"metadata": {},
"source": [
"Note that `R0` can be either constant or a given function of time.\n",
"\n",
"The initial conditions are set to"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "e41693f3",
"metadata": {
"hide-output": false
},
"outputs": [],
"source": [
"# initial conditions of s, e, i\n",
"i_0 = 1e-7\n",
"e_0 = 4 * i_0\n",
"s_0 = 1 - i_0 - e_0"
]
},
{
"cell_type": "markdown",
"id": "43ae689e",
"metadata": {},
"source": [
"In vector form the initial condition is"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "6cfc620f",
"metadata": {
"hide-output": false
},
"outputs": [],
"source": [
"x_0 = s_0, e_0, i_0"
]
},
{
"cell_type": "markdown",
"id": "ca4a5515",
"metadata": {},
"source": [
"We solve for the time path numerically using odeint, at a sequence of dates\n",
"`t_vec`."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "434330f0",
"metadata": {
"hide-output": false
},
"outputs": [],
"source": [
"def solve_path(R0, t_vec, x_init=x_0):\n",
" \"\"\"\n",
" Solve for i(t) and c(t) via numerical integration,\n",
" given the time path for R0.\n",
"\n",
" \"\"\"\n",
" G = lambda x, t: F(x, t, R0)\n",
" s_path, e_path, i_path = odeint(G, x_init, t_vec).transpose()\n",
"\n",
" c_path = 1 - s_path - e_path # cumulative cases\n",
" return i_path, c_path"
]
},
{
"cell_type": "markdown",
"id": "427d58f3",
"metadata": {},
"source": [
"## Experiments\n",
"\n",
"Let’s run some experiments using this code.\n",
"\n",
"The time period we investigate will be 550 days, or around 18 months:"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "d444168d",
"metadata": {
"hide-output": false
},
"outputs": [],
"source": [
"t_length = 550\n",
"grid_size = 1000\n",
"t_vec = np.linspace(0, t_length, grid_size)"
]
},
{
"cell_type": "markdown",
"id": "55f78431",
"metadata": {},
"source": [
"### Experiment 1: Constant R0 Case\n",
"\n",
"Let’s start with the case where `R0` is constant.\n",
"\n",
"We calculate the time path of infected people under different assumptions for `R0`:"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "8bc1d0c2",
"metadata": {
"hide-output": false
},
"outputs": [],
"source": [
"R0_vals = np.linspace(1.6, 3.0, 6)\n",
"labels = [f'$R0 = {r:.2f}$' for r in R0_vals]\n",
"i_paths, c_paths = [], []\n",
"\n",
"for r in R0_vals:\n",
" i_path, c_path = solve_path(r, t_vec)\n",
" i_paths.append(i_path)\n",
" c_paths.append(c_path)"
]
},
{
"cell_type": "markdown",
"id": "10964542",
"metadata": {},
"source": [
"Here’s some code to plot the time paths."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "989c2264",
"metadata": {
"hide-output": false
},
"outputs": [],
"source": [
"def plot_paths(paths, labels, times=t_vec):\n",
"\n",
" fig, ax = plt.subplots()\n",
"\n",
" for path, label in zip(paths, labels):\n",
" ax.plot(times, path, label=label)\n",
"\n",
" ax.legend(loc='upper left')\n",
"\n",
" plt.show()"
]
},
{
"cell_type": "markdown",
"id": "8340a09d",
"metadata": {},
"source": [
"Let’s plot current cases as a fraction of the population."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "4b04bde0",
"metadata": {
"hide-output": false
},
"outputs": [],
"source": [
"plot_paths(i_paths, labels)"
]
},
{
"cell_type": "markdown",
"id": "d4165872",
"metadata": {},
"source": [
"As expected, lower effective transmission rates defer the peak of infections.\n",
"\n",
"They also lead to a lower peak in current cases.\n",
"\n",
"Here are cumulative cases, as a fraction of population:"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "b0ff4b35",
"metadata": {
"hide-output": false
},
"outputs": [],
"source": [
"plot_paths(c_paths, labels)"
]
},
{
"cell_type": "markdown",
"id": "a8d60074",
"metadata": {},
"source": [
"### Experiment 2: Changing Mitigation\n",
"\n",
"Let’s look at a scenario where mitigation (e.g., social distancing) is\n",
"successively imposed.\n",
"\n",
"Here’s a specification for `R0` as a function of time."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "bde77f97",
"metadata": {
"hide-output": false
},
"outputs": [],
"source": [
"def R0_mitigating(t, r0=3, η=1, r_bar=1.6):\n",
" R0 = r0 * exp(- η * t) + (1 - exp(- η * t)) * r_bar\n",
" return R0"
]
},
{
"cell_type": "markdown",
"id": "d4c6575e",
"metadata": {},
"source": [
"The idea is that `R0` starts off at 3 and falls to 1.6.\n",
"\n",
"This is due to progressive adoption of stricter mitigation measures.\n",
"\n",
"The parameter `η` controls the rate, or the speed at which restrictions are\n",
"imposed.\n",
"\n",
"We consider several different rates:"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "52af61a9",
"metadata": {
"hide-output": false
},
"outputs": [],
"source": [
"η_vals = 1/5, 1/10, 1/20, 1/50, 1/100\n",
"labels = [fr'$\\eta = {η:.2f}$' for η in η_vals]"
]
},
{
"cell_type": "markdown",
"id": "a1358c89",
"metadata": {},
"source": [
"This is what the time path of `R0` looks like at these alternative rates:"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "0a53f3ec",
"metadata": {
"hide-output": false
},
"outputs": [],
"source": [
"fig, ax = plt.subplots()\n",
"\n",
"for η, label in zip(η_vals, labels):\n",
" ax.plot(t_vec, R0_mitigating(t_vec, η=η), label=label)\n",
"\n",
"ax.legend()\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"id": "7fbec267",
"metadata": {},
"source": [
"Let’s calculate the time path of infected people:"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "369ee547",
"metadata": {
"hide-output": false
},
"outputs": [],
"source": [
"i_paths, c_paths = [], []\n",
"\n",
"for η in η_vals:\n",
" R0 = lambda t: R0_mitigating(t, η=η)\n",
" i_path, c_path = solve_path(R0, t_vec)\n",
" i_paths.append(i_path)\n",
" c_paths.append(c_path)"
]
},
{
"cell_type": "markdown",
"id": "6d36ea27",
"metadata": {},
"source": [
"These are current cases under the different scenarios:"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "d5deb9c0",
"metadata": {
"hide-output": false
},
"outputs": [],
"source": [
"plot_paths(i_paths, labels)"
]
},
{
"cell_type": "markdown",
"id": "35759fd8",
"metadata": {},
"source": [
"Here are cumulative cases, as a fraction of population:"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "36932e35",
"metadata": {
"hide-output": false
},
"outputs": [],
"source": [
"plot_paths(c_paths, labels)"
]
},
{
"cell_type": "markdown",
"id": "9e720da9",
"metadata": {},
"source": [
"## Ending Lockdown\n",
"\n",
"The following replicates [additional results](https://drive.google.com/file/d/1uS7n-7zq5gfSgrL3S0HByExmpq4Bn3oh/view) by Andrew Atkeson on the timing of lifting lockdown.\n",
"\n",
"Consider these two mitigation scenarios:\n",
"\n",
"1. $ R_t = 0.5 $ for 30 days and then $ R_t = 2 $ for the remaining 17 months. This corresponds to lifting lockdown in 30 days. \n",
"1. $ R_t = 0.5 $ for 120 days and then $ R_t = 2 $ for the remaining 14 months. This corresponds to lifting lockdown in 4 months. \n",
"\n",
"\n",
"The parameters considered here start the model with 25,000 active infections\n",
"and 75,000 agents already exposed to the virus and thus soon to be contagious."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "fcf256ba",
"metadata": {
"hide-output": false
},
"outputs": [],
"source": [
"# initial conditions\n",
"i_0 = 25_000 / pop_size\n",
"e_0 = 75_000 / pop_size\n",
"s_0 = 1 - i_0 - e_0\n",
"x_0 = s_0, e_0, i_0"
]
},
{
"cell_type": "markdown",
"id": "40a68958",
"metadata": {},
"source": [
"Let’s calculate the paths:"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "680e6339",
"metadata": {
"hide-output": false
},
"outputs": [],
"source": [
"R0_paths = (lambda t: 0.5 if t < 30 else 2,\n",
" lambda t: 0.5 if t < 120 else 2)\n",
"\n",
"labels = [f'scenario {i}' for i in (1, 2)]\n",
"\n",
"i_paths, c_paths = [], []\n",
"\n",
"for R0 in R0_paths:\n",
" i_path, c_path = solve_path(R0, t_vec, x_init=x_0)\n",
" i_paths.append(i_path)\n",
" c_paths.append(c_path)"
]
},
{
"cell_type": "markdown",
"id": "f27917ba",
"metadata": {},
"source": [
"Here is the number of active infections:"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "efcbc9fb",
"metadata": {
"hide-output": false
},
"outputs": [],
"source": [
"plot_paths(i_paths, labels)"
]
},
{
"cell_type": "markdown",
"id": "a6fc19a9",
"metadata": {},
"source": [
"What kind of mortality can we expect under these scenarios?\n",
"\n",
"Suppose that 1% of cases result in death"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "2eb91e32",
"metadata": {
"hide-output": false
},
"outputs": [],
"source": [
"ν = 0.01"
]
},
{
"cell_type": "markdown",
"id": "46a991c6",
"metadata": {},
"source": [
"This is the cumulative number of deaths:"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "bb88cbbe",
"metadata": {
"hide-output": false
},
"outputs": [],
"source": [
"paths = [path * ν * pop_size for path in c_paths]\n",
"plot_paths(paths, labels)"
]
},
{
"cell_type": "markdown",
"id": "f7bf2cdb",
"metadata": {},
"source": [
"This is the daily death rate:"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "c7ba555d",
"metadata": {
"hide-output": false
},
"outputs": [],
"source": [
"paths = [path * ν * γ * pop_size for path in i_paths]\n",
"plot_paths(paths, labels)"
]
},
{
"cell_type": "markdown",
"id": "edddb263",
"metadata": {},
"source": [
"Pushing the peak of curve further into the future may reduce cumulative deaths\n",
"if a vaccine is found."
]
}
],
"metadata": {
"date": 1714442510.084284,
"filename": "sir_model.md",
"kernelspec": {
"display_name": "Python",
"language": "python3",
"name": "python3"
},
"title": "Modeling COVID 19"
},
"nbformat": 4,
"nbformat_minor": 5
}