SimpleEpidemicModel.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# The simplest epidemic model\n",
    "\n",
    "### The deterministic SIR model, following Kermack & McKendrick (1927)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from matplotlib import pyplot as plt\n",
    "from ipywidgets import interactive,HBox,VBox\n",
    "import ipywidgets as widgets\n",
    "from scipy.integrate import odeint"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# SIR model in non-dimensional form \n",
    "def SIR(x,t,Ro,tau,d,t_start_qua,t_stop_qua,Ro_qua):\n",
    "    if (t >= t_start_qua/tau) & (t <= t_stop_qua/tau):\n",
    "        Ro_t = Ro_qua\n",
    "    else:\n",
    "        Ro_t = Ro\n",
    "    S,I,R,P = x\n",
    "    xdot = np.array([-Ro_t*S*I, Ro_t*S*I - I, (1.-d)*I, -d*I])\n",
    "    return xdot\n",
    "\n",
    "def solve(N,L,x0,Ro,tau,d,t_start_qua,t_stop_qua,Ro_qua):\n",
    "    t = np.linspace(0,L,N)/tau\n",
    "    x = odeint(SIR, x0, t, args=(Ro,tau,d,t_start_qua,t_stop_qua,Ro_qua))\n",
    "    # output\n",
    "    x = x*100. # convert to %\n",
    "    return x.T"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "b1134eea4c4d4452940fefde23b92b86",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "092d198f693444dfa2491f1e5b3be7a7",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(VBox(children=(FloatSlider(value=2.5, description='$R_0$', max=6.0, min=0.5, step=0.25), IntSli…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%matplotlib widget\n",
    "\n",
    "# Initial parameter setting\n",
    "Ro = 2.5\n",
    "tau = 14. # duration of infection (days)\n",
    "d = 0.   # death rate (% infected who die)\n",
    "N = 120   # no. timesteps\n",
    "L = 400.  # length run (days)\n",
    "t_start_qua = 0. # time quarantine starts (days)\n",
    "t_stop_qua = 0.1 # time quarantine ends (days)\n",
    "Ro_qua = 2.5      # Ro during quarantine\n",
    "\n",
    "# Initial conditions\n",
    "x0 = np.array([0.999, 0.001, 0., 1.])\n",
    "\n",
    "# Solve with initial parameters and plot\n",
    "S,I,R,P = solve(N, L, x0, Ro, tau, d, t_start_qua, t_stop_qua, Ro_qua)\n",
    "fig,ax = plt.subplots(1,1,figsize=[8,3])\n",
    "time = np.linspace(0,L,N)\n",
    "l1,=ax.plot(time,S)\n",
    "l2,=ax.plot(time,I)\n",
    "l3,=ax.plot(time,R)\n",
    "#l4,=ax.plot(time,100-P,'k-')\n",
    "rect=ax.axvspan(t_start_qua,t_stop_qua,color='0.5',alpha=0.5,linewidth=0,zorder=1)\n",
    "ax.set_xlabel('time (days)',fontsize=12)\n",
    "ax.set_ylabel('fraction of population (%)',fontsize=12)\n",
    "ax.legend(['Susceptible','Infected','Removed'],loc=(1.1,0.7))\n",
    "ax.set_xlim(0,L)\n",
    "ax.set_xticks(np.arange(0,400,30))\n",
    "ax.set_yticks(np.arange(0,110,10))\n",
    "ax.yaxis.set_ticks_position('both')\n",
    "ax.tick_params(labeltop=False, labelright=True)\n",
    "ax.grid(axis='y',linestyle=':',linewidth=0.5)\n",
    "plt.tight_layout()\n",
    "fig.canvas.header_visible=False\n",
    "\n",
    "# Interactive plot\n",
    "Ro_widget = widgets.FloatSlider(min=0.5, max=6, step=0.25, value=Ro,\n",
    "                                continuous_update=True, description=r'$R_0$')\n",
    "tau_widget = widgets.IntSlider(min=1, max=30, step=1, value=tau,\n",
    "                                continuous_update=True, description=r'$\\tau$ (days)')\n",
    "#d_widget = widgets.IntSlider(min=0, max=10, step=1, value=1,\n",
    "#                             continuous_update=True, description='d (%)')\n",
    "t_start_qua_widget = widgets.IntSlider(min=0, max=300, step=5, value=t_start_qua,\n",
    "                                       continuous_update=True, description='$t_\\mathrm{start\\ qua}$ (day)')\n",
    "t_stop_qua_widget = widgets.IntSlider(min=0, max=300, step=5, value=t_stop_qua,\n",
    "                                       continuous_update=True, description='$t_\\mathrm{end\\ qua}$ (day)')\n",
    "Ro_qua_widget = widgets.FloatSlider(min=0., max=3., step=0.25, value=Ro_qua,\n",
    "                                continuous_update=True, description='$R_{0q}$')\n",
    "def make_plot(Ro, tau, t_start_qua, t_stop_qua, Ro_qua):\n",
    "    S,I,R,P = solve(N, L, x0, Ro, tau, d/100., t_start_qua, t_stop_qua, Ro_qua)\n",
    "    l1.set_ydata(S)\n",
    "    l2.set_ydata(I)\n",
    "    l3.set_ydata(R)\n",
    "    #l4.set_ydata(100-P)\n",
    "    xy = rect.get_xy()\n",
    "    if t_start_qua < t_stop_qua:\n",
    "        xy[:,0] =[t_start_qua,t_start_qua,t_stop_qua,t_stop_qua,t_start_qua]\n",
    "    else:\n",
    "        xy[:,0] =[0,0,0.1,0.1,0]\n",
    "    rect.set_xy(xy)\n",
    "\n",
    "w=interactive(make_plot,\n",
    "          Ro          = Ro_widget, \n",
    "          tau          = tau_widget, \n",
    "          #d           = d_widget, \n",
    "          t_start_qua = t_start_qua_widget,\n",
    "          t_stop_qua  = t_stop_qua_widget,\n",
    "          Ro_qua      = Ro_qua_widget)\n",
    "\n",
    "items = w.children\n",
    "left_box = VBox([items[0], items[1]])\n",
    "right_box = VBox([items[2], items[3], items[4]])\n",
    "HBox([left_box, right_box])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### The short story\n",
    "The interactive app above solves the deterministic SIR model with paremeter settings adjustable via the sliders. To understand what the parameters mean, read on.\n",
    "\n",
    "### The long story\n",
    "In this time of pandemic, the media are full of terms such as \"reproduction number\", \"herd immunity\", \"non-pharmacological intervention\" and \"second wave infections\" after quarantines are lifted. What do these terms mean? How can we have an informed debate without really understanding what they mean?\n",
    "\n",
    "With a bit of math, it's easy to get a good sense for these concepts. The foundations of mathematical epidemiology were laid in the 1920's, especially through the work of Kermack and McKendrick [(see Anderson 1991 for a historical account)](https://link.springer.com/article/10.1007/BF02464422). Their basic insights are still valid today, though very much complexity has been added to the modelling. What follows is a re-derivation of their original work.\n",
    "\n",
    "### The model\n",
    "Suppose a novel virus is introduced to a population consisting of $P$ individuals. As the virus spreads, the population will be split into three classes of individuals: \n",
    "- *Susceptible*: those who have not yet been exposed to the virus;\n",
    "- *Infected*: those who have been exposed, have been infected and can transmit the virus;\n",
    "- *Removed*: those who have been infected, have gone through the disease and have recovered and become immune or have died from the disease.\n",
    "\n",
    "At any given time, there are $S$, $I$ and $R$ individuals in each class (hence the name SIR model). For the virus to spread, an infected individual must make contact with a susceptible individual. Make the following assumptions:\n",
    "1. each individual encounters a fixed number $n$ other individuals each day;\n",
    "2. of these $n$, a fixed fraction $S/P$ is susceptible;\n",
    "3. in each encounter between an infected and a susceptible individual, virus transmission occurs with fixed probability $p$;\n",
    "4. an infected individual remains infected and can transmit the virus for a time $\\tau$ days, after which they recover and become immune (or die);\n",
    "5. none of the population changes its behavior in response to the virus.\n",
    "\n",
    "With these assumptions, the number of new individuals becoming infected each day is $p\\,n\\,\\frac{S}{P}\\,I,$ and the number of people recovering is $I/\\tau$. We can then write a system of equations for the rate of change of each population class:\n",
    "\n",
    "\\begin{equation}\n",
    "\\frac{dS}{dt} = - p\\,n\\,\\frac{S}{P}\\,I, \\ \\ \\ \\ \\ \\ \\frac{dI}{dt} =  p\\,n\\,\\frac{S}{P}\\,I \\ - \\ \\tau^{-1}\\, I, \\ \\ \\ \\ \\ \\ \\frac{dR}{dt} = \\tau^{-1}\\, I \n",
    "\\end{equation}\n",
    "\n",
    "We then adimensionalize the equations by scaling $S$, $I$ and $R$ by $P$ and time $t$ by $\\tau$, to give:\n",
    "\n",
    "\\begin{equation}\n",
    "\\frac{dS}{dt} = - R_o \\, S \\, I, \\ \\ \\ \\ \\ \\ \\frac{dI}{dt} =  R_o \\, S \\, I \\ - \\ I, \\ \\ \\ \\ \\ \\ \\frac{dR}{dt} = I \n",
    "\\end{equation}\n",
    "\n",
    "where\n",
    "\n",
    "\\begin{equation}\n",
    "R_0 \\ = \\ p\\, n \\, \\tau.\n",
    "\\end{equation}\n",
    "\n",
    "Note that in this adimensional form, $S$, $I$ and $R$ refer to the *fraction* of the population in each class.\n",
    "\n",
    "### Reproduction number\n",
    "$R_0$ is known as the *reproduction number*. It can be interpreted as the mean number of people that an infected individual will in turn infect. It is the fundamental non-dimensional parameter controlling the spread of an epidemic. From the second equation set above, it can be seen that in the early stages of the epidemic, when $S\\sim1$ and $I\\sim0$, $I$ grows as\n",
    "\n",
    "$$ I \\ \\sim \\ \\exp{\\left( \\frac{R_o-1}{\\tau} \\, t \\right)}. $$ \n",
    "\n",
    "Thus, if $R_0>1$ the virus goes viral and spreads exponentially; while if $R_0<1$ the outbreak quickly subsides. For the ongoing COVID-19 epidemic, estimates of $R_0$ are in the ballpark of 2-3 [(Liu et al. 2020)](https://academic.oup.com/jtm/article/27/2/taaa021/5735319). Taking $R_0$ = 2.5 and $\\tau$=8 days gives a growth rate of about 30% per day, roughly what is seen in the data (https://www.ft.com/coronavirus-latest).\n",
    "\n",
    "In the interactive app, you can see that setting $R_0$ higher makes both the growth and the peak of infections higher. Making $\\tau$ smaller postpones the peak of infections but not its amplitude.\n",
    "\n",
    "You can also see that once the main infection peak has passed, the model settles into a stable steady state with a residual susceptible population--the virus stops short of infecting the whole population, though the size of the steady-state susceptible population gets smaller as $R_0$ increases.\n",
    "\n",
    "### Herd immunity\n",
    "Now suppose the virus is introduced to a population where a certain fraction is *already immune* to it (i.e. at the initial time $R>0$ and $S<1$, with $I\\sim0$). From the second equation set above, it can be seen that the epidemic will not spread if\n",
    "$ R_0\\, S -1 < 0$, which requires $S < 1/R_0$. This leads to the concept of *herd immunity*: an epidemic will not spread if there is initially a fraction \n",
    "\n",
    "$$ 1 \\ - \\ \\frac{1}{R_0}$$ \n",
    "\n",
    "of the population that is already immune. For $R_0$=2.5, this fraction is 60%, a number often repeated in the media.\n",
    "\n",
    "In the interactive app--which assumes no previous immunity--you can see that once an epidemic starts it does not stop at the immunity threshold; $R$ overshoots to a higher value.\n",
    "\n",
    "### Non-pharmacological intervention and second wave infections\n",
    "With no vaccine or anti-viral drugs to fight the disease, the only way to control the spread of the epidemic is by changing the behaviour of individuals so as to reduce $R_0$. Recall that $R_0 = p\\,n\\,\\tau$: we can reduce the probability of transmission $p$ by washing our hands, coughing into elbows etc., and we can reduce $n$ by avoiding contact with other people. *Non-pharmacological intervention* means that authorities mandate such actions, by imposing social distancing norms, quarantines, curfews etc. \n",
    "\n",
    "In the interactive app, you can play with a \"quarantine\" period starting at time $t_\\mathrm{start\\ qua}$ and ending at $t_\\mathrm{end\\ qua}$, during which time $R_0$ takes on the value $R_{0q}$ (whose value is presumably less than that of the basic, non-intervention $R_0$). You can see that the quarantine will temporarily suppress the growth of infections, but infections shoot back up in a *second wave* if the herd immunity threshold has not been reached.\n",
    "\n",
    "### More complex models\n",
    "The minimalist SIR model described here can be extended in many ways to relax the highly unrealistic assumptions made. For example, the infected class $I$ can be segmented into a set of non-symptomatic carriers and people with symptoms of various severity. For a very nice interactive example (with much better graphics than mine) see Gabriel Goh's [Epidemic Calculator](https://gabgoh.github.io/COVID/index.html).\n",
    "\n",
    "Seasonality driven by weather-related variations in virus transmission rates can also be incorporated, as explained [here](https://smw.ch/article/doi/smw.2020.20224). An interesting result is that the model's long-term steady state shows damped fluctuations in response to perturbations, which can resonate with the seasonal forcing.\n",
    "\n",
    "\n",
    "A much more flexible approach is to make the model stochastic, where all parameters take random values drawn from certain probability distributions. See my colleague Tom Britton's intro talk on [stochastic epidemic modelling](https://www.youtube.com/watch?v=gSqIwXl6IjQ).\n",
    "\n",
    "The most sophisticated models used for actual operational predictions and decision-making use agent-based approaches with full geographical resolution, see for example the [Imperial College report](https://www.imperial.ac.uk/media/imperial-college/medicine/sph/ide/gida-fellowships/Imperial-College-COVID19-NPI-modelling-16-03-2020.pdf) led by Niel Ferguson that has attracted much attention.\n",
    "\n",
    "I am not an epidemiologist but a climate scientist, and the analogy with climate modelling is interesting. The simple SIR model presented here is analogous to the basic zero-dimensional energy balance model we teach in undergraduate climate classes: good for physical insight, and gets some of the basic numbers in the right ballpark. The Imperial College-style models are the counterpart to full-complexity climate models, which are hard or impossible to understand in full but provide more spatial and process detail. \n",
    "\n",
    "### Acknowledgements\n",
    "Discussion and inspiration from Ken Haste Andersen, see his epidemic simulator [here](http://oceanlife.dtuaqua.dk/Corona). \n",
    "\n",
    "This is a [Jupyter notebook](https://jupyter.org) served through [Voilà](https://voila.readthedocs.io). Code is [here](https://github.com/rodrigo-caballero/voila).\n",
    "\n",
    "[Rodrigo Caballero](http://climdyn.misu.su.se), March 2020. "
   ]
  }
 ],
 "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.7.6"
  },
  "latex_envs": {
   "LaTeX_envs_menu_present": true,
   "autoclose": false,
   "autocomplete": true,
   "bibliofile": "biblio.bib",
   "cite_by": "apalike",
   "current_citInitial": 1,
   "eqLabelWithNumbers": true,
   "eqNumInitial": 1,
   "hotkeys": {
    "equation": "Ctrl-E",
    "itemize": "Ctrl-I"
   },
   "labels_anchors": false,
   "latex_user_defs": false,
   "report_style_numbering": false,
   "user_envs_cfg": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}