01_SIR_Model.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "This post is the first in a series in which we'll use a simple but effective mathematical\n",
    "model to predict the ongoing Coronavirus outbreak.  As I write, the number of official\n",
    "confirmed global cases is just under 200,000 and many schools and workplaces all over the world have closed in order to slow its spread.  If you're like me, you're wondering:\n",
    "\n",
    "  - Am I likely to catch this virus?\n",
    "  - How long will it be until my school or workplace opens up again?\n",
    "\n",
    "In television and movies, math is often used to predict things like the behavior of a specific\n",
    "individual; math is a terrible tool for that kind of thing (indeed, math tells us that such things\n",
    "are inherently unpredictable).  But math is quite good at predicting the average behavior of large\n",
    "groups, and a little math can go a long way in telling us what will happen next with COVID-19.  My goal here is to help you understand those predictions with nothing more than high school math.\n",
    "\n",
    "*Disclaimer*: I am not an epidemiologist and I make no guarantees about the predictions\n",
    "we'll arrive at here.  By reading this you agree not to sue me.  What I write here does\n",
    "not reflect the opinion of my employer or anyone else.  My goal is simply to show you how\n",
    "some relatively straightforward mathematics can give you insight into what appears to be\n",
    "a very uncertain future.\n",
    "\n",
    "Each post in this series will be written in a Jupyter notebook, which you can download and experiment with yourself if you are so inclined.  The notebook for this first post is [here](https://github.com/ketch/covid-blog-posts/blob/master/The%20SIR%20Model%20-%20blog.ipynb)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Modeling the spread of infectious disease"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "An infectious disease spreads from one individual to another.  Consider the following simple model:\n",
    "\n",
    "  - On day zero, a single individual is infected\n",
    "  - On each subsequent day, each infected individual passes the disease to one more individual\n",
    "  \n",
    "How quickly does the number of infected individuals grow?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "1, 2, 4, 8, 16, ..."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "On each day, the number of infected doubles!  How many days would it take for everyone on earth to be infected?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "This is not a reasonable model for any of the diseases we know of.  Of course, the rate of new infections per infected person (1 per day) was an arbitrary choice and real values are likely to be smaller.  What other effects are missing from this model?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "  - **Recovery and immunity**: eventually, an individual recovers and can no longer infect others\n",
    "  - **Spread**: The infection can only spread to people who don't yet have it.  If most of the individuals in contact with an infected person are already infected, that person is less likely to spread the disease to someone new"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "Those two factors are vital to understanding the true dynamics of epidemics.  Of course, there are many other important details we have left out; for instance:\n",
    "\n",
    "  - The disease may affect different individuals in different ways.\n",
    "  - Individuals are spread out geographically.\n",
    "  - Certain individuals are likely to infect many others, while others are less likely.  This depends on many factors including culture, personality, and lifestyle, as well as the mode of transmission of the disease.\n",
    "  - Individuals might take actions to avoid getting infected (e.g. washing hands, avoiding sick people) or to avoid spreading the disease (e.g. staying home when sick).\n",
    "  \n",
    "All of these effects (and many others) will influence the spread of the disease.  A model that tries to incorporate them all would be very complex."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## The SIR model\n",
    "\n",
    "One of the simplest but most relevant models is based on the idea that the population consists of three groups:\n",
    "\n",
    "  - **S(t)** Susceptible (those who have not yet been infected)\n",
    "  - **I(t)** Infected (those who can currently spread the disease)\n",
    "  - **R(t)** Recovered (those who are now immune and cannot contract or spread the disease)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "When we write $S(t)$, we mean that the number of susceptible individuals ($S$) is given as a function of time ($t$).  We can't write down this function exactly; instead it will be described by a *differential equation*.  Now, differential equations are a bit like [whale sharks](https://en.wikipedia.org/wiki/Whale_shark): they sound scary at first, but in reality they are simple and friendly creatures.\n",
    "\n",
    "A differential equation is just a description of how some quantity changes.  In this case, we will have three differential equations, describing the rate of change of $S$, $I$, and $R$.  The idea is that susceptible people can become infected and infected people can become recovered:\n",
    "\n",
    "$$ S \\to I \\to R$$\n",
    "\n",
    "To define differential equations for the three groups, we only need to determine the rate at which each of these transitions occurs."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Rate of infection\n",
    "\n",
    "In our first simple model, we assumed the rate of infection was proportional to the number of infected.  This is very reasonable, but for someone new to become infected we need both an infected individual **and** a susceptible one.  If we imagine that people encounter each other randomly at some rate $\\beta$, then the rate of new infections is just the number of infected multiplied by the probability of encountering a susceptible individual:\n",
    "\n",
    "$$\n",
    "\\frac{dI}{dt} = \\beta I \\frac{S}{N}.\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "When we write $\\frac{dI}{dt}$, we mean *the rate of change of $I$ over time*.  So this equation tells us how quickly the number of infected people will increase.\n",
    "Here $N=S+I+R$ is the total population, so $S/N$ is the probability that a randomly chosen individual is susceptible."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "Of course, since new infected people were previously susceptible, the number of susceptible individuals must decrease at the same rate:\n",
    "\n",
    "$$\n",
    "\\frac{dS}{dt} = -\\beta I \\frac{S}{N}.\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Rate of recovery\n",
    "\n",
    "The other transition is from infected to recovered.  A proper model for this should involve a time delay, since (for many diseases) new infected individuals typically become recovered after a certain interval of time.  For instance, with the flu or the new Coronavirus, the number of new recovered individuals might depend on how many became infected about one or two weeks ago.  Incorporating such an effect would lead to a more complicated model known as a **delay differential equation**."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "Instead, we will simply assume that over any time interval, a certain fraction of the infected become recovered.  Denoting the recovery rate by $\\gamma$, we have\n",
    "\n",
    "$$\n",
    "\\frac{dR}{dt} = \\gamma I.\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "The number of infected must decrease at the same rate, so we must modify our ODE for $I(t)$ to read\n",
    "\n",
    "$$\n",
    "\\frac{dI}{dt} = \\beta I \\frac{S}{N}-\\gamma I.\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### The full model\n",
    "\n",
    "\\begin{align}\n",
    "\\frac{dS}{dt} & = -\\beta I \\frac{S}{N} \\\\\n",
    "\\frac{dI}{dt} & = \\beta I \\frac{S}{N}-\\gamma I \\\\\n",
    "\\frac{dR}{dt} & = \\gamma I\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "Notice that if we add the 3 equations together, we get \n",
    "\n",
    "$$\n",
    "\\frac{dN}{dt} = 0.\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "What do $\\beta$ and $\\gamma$ really mean?  We can think of $\\beta$ as the number of others that one infected person encounters per unit time, and $\\gamma^{-1}$ as the typical time from infection to recovery.  So the number of new infections generated by one infected individual is, on average, $$\\beta/\\gamma = R_0,$$ the **basic reproduction number**."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### SIR dynamics\n",
    "\n",
    "Notice that $S(t)$ can only decrease and $R(t)$ can only increase, but $I(t)$ may increase or decrease.  A key question is, under what conditions will $I(t)$ increase?  This will tell us whether a small number of cases could become an epidemic."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "We can write\n",
    "\n",
    "$$\n",
    "\\frac{dI}{dt} = \\left(\\beta \\frac{S}{N}-\\gamma \\right) I\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "from which we see that $I(t)$ grows if \n",
    "$$\\beta S/N > \\gamma.$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "$$\n",
    "\\frac{dI}{dt} = \\left(\\beta \\frac{S}{N}-\\gamma \\right) I\n",
    "$$\n",
    "\n",
    "Initially in a population we have \n",
    "$$S/N \\approx 1,$$\n",
    "\n",
    "so an epidemic of some size can occur if $\\beta > \\gamma$.  As the epidemic grows, the ratio $S/N$ becomes smaller, so eventually the spread slows down.\n",
    "\n",
    "What fraction of the population must been infected before $I(t)$ will start to decrease?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "The epidemic will begin to subside when \n",
    "$$S/N = (\\beta/\\gamma)^{-1} = R_0^{-1}.$$\n",
    "\n",
    "This determines the infection peak.  After this point, there will still be new infections but the overall number of infected will decrease."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## An example\n",
    "So what does an epidemic look like, using the SIR model?  We can easily compute the solution using standard numerical methods; I've omitted the code here since I want to focus on the model, but feel free to look at [the original Jupyter notebook with code](https://github.com/ketch/covid-blog-posts/blob/master/The%20SIR%20Model%20-%20blog.ipynb) and modify the parameters yourself."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from ipywidgets import interact, widgets\n",
    "plt.style.use('seaborn-poster')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "beta = 1.0\n",
    "gamma = 0.1\n",
    "N = 1.\n",
    "\n",
    "# Initial values\n",
    "\n",
    "R = [0]\n",
    "I = [1.e-6]\n",
    "S = [N - I[0] - R[0]]\n",
    "\n",
    "def change(S,I,R):\n",
    "    S_change = -beta*I*S/N\n",
    "    I_change = beta*I*S/N - gamma*I\n",
    "    R_change = gamma*I\n",
    "    return S_change, I_change, R_change\n",
    "\n",
    "\n",
    "T = 100\n",
    "dt = 0.01\n",
    "times = np.linspace(0,T,int(T/dt)+1)\n",
    "\n",
    "for t in times[1:]:\n",
    "    # Just use Euler's method for this introductory example\n",
    "    S_change, I_change, R_change = change(S[-1],I[-1],R[-1])\n",
    "    I.append(I[-1]+dt*I_change)\n",
    "    S.append(S[-1]+dt*S_change)\n",
    "    R.append(R[-1]+dt*R_change)\n",
    "\n",
    "plt.figure(figsize=(12,6))\n",
    "plt.plot(times,S,lw=3);\n",
    "plt.plot(times,I,lw=3);\n",
    "plt.plot(times,R,lw=3);\n",
    "plt.legend(['Susceptible','Infected','Recovered'])\n",
    "plt.xlabel('Time');\n",
    "plt.xlim(0,100); plt.ylim(-.01,1.01);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here I've set $N=1$ so the numbers on the vertical axis represent a fraction of the current population.  Initially, only a tiny fraction is infected while the remainder is susceptible.  The plot above shows the typical behavior of an epidemic: an initial rapid exponential spread, until much of the population is infected or recovered, at which point the number of infections begins to decline.\n",
    "\n",
    "According to our analysis above, the number of infections should begin to decrease when the susceptible fraction $S/N$ is equal to $\\gamma/\\beta$.  Here I've taken $\\beta=1$ and $\\gamma=1/10$, so the infection peak should occur when the susceptible population has dropped to 1/10.  Let's check:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(12,6))\n",
    "plt.plot(times,S,lw=3);\n",
    "plt.plot(times,I,lw=3);\n",
    "plt.legend(['Susceptible','Infected'])\n",
    "plt.xlabel('Time');\n",
    "plt.plot([18.,18.],[0,1],'--k')\n",
    "plt.plot([0,30],[0.1,0.1],'--k');\n",
    "plt.xlim(0,100);\n",
    "plt.ylim(-.01,1.01);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The results are in perfect agreement with what we predicted.  In the [original notebook](https://github.com/ketch/covid-blog-posts/blob/master/The%20SIR%20Model%20-%20blog.ipynb) there is an interactive model where you can adjust the parameters and see the results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def SIR_model(beta=1.0,gamma=0.1):\n",
    "\n",
    "    R = [0]\n",
    "    I = [1.e-6]\n",
    "    S = [N - I[0] - R[0]]\n",
    "\n",
    "    def change(S,I,R):\n",
    "        S_change = -beta*I*S/N\n",
    "        I_change = beta*I*S/N - gamma*I\n",
    "        R_change = gamma*I\n",
    "        return S_change, I_change, R_change\n",
    "\n",
    "\n",
    "    T = 100\n",
    "    dt = 0.01\n",
    "    times = np.linspace(0,T,int(T/dt)+1)\n",
    "\n",
    "    for t in times[1:]:\n",
    "        S_change, I_change, R_change = change(S[-1],I[-1],R[-1])\n",
    "        I.append(I[-1]+dt*I_change)\n",
    "        S.append(S[-1]+dt*S_change)\n",
    "        R.append(R[-1]+dt*R_change)\n",
    "\n",
    "    plt.figure(figsize=(12,6))\n",
    "    plt.plot(times,S,lw=3);\n",
    "    plt.plot(times,I,lw=3);\n",
    "    plt.plot(times,R,lw=3);\n",
    "    plt.legend(['Susceptible','Infected','Recovered'])\n",
    "    plt.xlabel('Time');\n",
    "    plt.xlim(0,100); plt.ylim(-.01,1.01);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "interact(SIR_model,beta=widgets.FloatSlider(min=0.01,max=1.0,step=0.01,value=1.0,description=r'$\\beta$'),\n",
    "         gamma=widgets.FloatSlider(min=0.01,max=1.0,step=0.01,value=0.1,description=r'$\\gamma$'));"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the next post, we will look at real-world data to estimate the values of $\\beta$ and $\\gamma$."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.5"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}