04_Exponential_growth.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is my fourth post on modeling the Coronavirus epidemic.  I recommend starting with [the first post](http://www.davidketcheson.info/2020/03/17/SIR_model.html) and reading them in order.\n",
    "\n",
    "So far, we've [learned about the SIR model](http://www.davidketcheson.info/2020/03/17/SIR_model.html) and [used available data](http://www.davidketcheson.info/2020/03/19/SIR_Estimating_parameters.html) combined with the model to [predict the epidemic](http://www.davidketcheson.info/2020/03/19/SIR_predictions.html).  We're now going to detour into some additional mathematical ideas that will help us get further understanding."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Exponential growth\n",
    "\n",
    "In the second post, we saw that the initial spread of a disease follows a differential equation of the form\n",
    "\n",
    "$$\n",
    "\\frac{dI}{dt} = \\beta I(t).\n",
    "$$\n",
    "\n",
    "Once you're comfortable with the idea that $dI/dt$ just means \"the rate of change of $I$\", you realize that this is one of the simplest equations imaginable.  So it shouldn't be surprising that this equation comes up a lot in the real world.  Its solution is\n",
    "\n",
    "$$\n",
    "I(t) = e^{\\beta t} I(0).\n",
    "$$\n",
    "\n",
    "Here $e\\approx 2.72$ is **Euler's number**.  This equation tells us that the number of infected grows very quickly.  In the case of Coronavirus, the number $I$ can double in about 3 days ([you can see more estimates of the doubling time here](https://ourworldindata.org/coronavirus#growth-country-by-country-view)).  We refer to this kind of growth (where a given quantity doubles over a certain time interval) as **exponential growth**.  \n",
    "\n",
    "Financial advisors love to talk about the power of exponential growth because that's also how compound interest works.  Of course, with most investments it takes a lot more than 3 days to double your money, but the principle is the same.  In fact, exponentially growing functions are all around us, and learning how they behave can help you to understand a lot of things."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's look again at our first prediction from the last post:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib\n",
    "import matplotlib.dates as dates\n",
    "from scipy.integrate import solve_ivp\n",
    "import pandas as pd\n",
    "plt.style.use('seaborn-poster')\n",
    "matplotlib.rcParams['figure.figsize'] = (10., 6.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def SIR(beta=0.25, gamma=0.05, N = 7770000000, confirmed_fraction=0.15,T=400,plot='SIR',\n",
    "        print_params=True,ymax=None):\n",
    "    \"\"\" Model the current outbreak using the SIR model.\"\"\"\n",
    "\n",
    "    du = np.zeros(3)\n",
    "    u0 = np.zeros(3)\n",
    "    \n",
    "    def f(t,u):\n",
    "        du[0] = -beta*u[1]*u[0]/N\n",
    "        du[1] = beta*u[1]*u[0]/N - gamma*u[1]\n",
    "        du[2] = gamma*u[1]\n",
    "        return du\n",
    "\n",
    "    # Initial values\n",
    "    total_cases = 199258./confirmed_fraction\n",
    "    u0[2] = 81972./confirmed_fraction+7956  # Initial recovered\n",
    "    u0[1] = total_cases-u0[2] # Initial infected\n",
    "    u0[0] = N - u0[1] - u0[2]\n",
    "\n",
    "    #dt = 0.01\n",
    "    times = np.arange(0,T)\n",
    "    solution = solve_ivp(f,[0,T],u0,t_eval=times)\n",
    "    S = solution.y[0,:]\n",
    "    I = solution.y[1,:]\n",
    "    R = solution.y[2,:]\n",
    "\n",
    "    today = '19/03/2020'\n",
    "    start = dates.datestr2num(today)\n",
    "    mydates = np.arange(T)+start\n",
    "    \n",
    "    fig = plt.figure(figsize=(12,8)); leg=[]\n",
    "    if 'S' in plot:\n",
    "        leg.append('Susceptible')\n",
    "        plt.plot_date(mydates,S/1.e9,'-b',lw=3)\n",
    "    if 'I' in plot:\n",
    "        leg.append('Infected')\n",
    "        plt.plot_date(mydates,I/1.e9,'-',lw=3,color='brown')\n",
    "    if 'R' in plot:\n",
    "        leg.append('Recovered')\n",
    "        plt.plot_date(mydates,R/1.e9,'-g',lw=3)\n",
    "    ax = plt.gca()\n",
    "    ax.xaxis.set_major_locator(dates.MonthLocator())\n",
    "    ax.xaxis.set_major_formatter(dates.DateFormatter('%b'))\n",
    "    fig.autofmt_xdate()\n",
    "    plt.legend(leg)\n",
    "    plt.ylabel('Individuals (in billions)')\n",
    "    plt.xlim(737503,737503+T)\n",
    "    if ymax: plt.ylim(-ymax/10,ymax)\n",
    "    if print_params:\n",
    "        plt.text(737503+3*T/4, 6,r'$\\beta$ = {}, $\\gamma$ = {}'.format(beta,gamma),fontsize=15,bbox=dict(facecolor='white', alpha=1.0))\n",
    "    return fig\n",
    "\n",
    "SIR(beta=0.25, gamma=0.05, confirmed_fraction=0.15, plot='I',print_params=False);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Right now, we are at the very left edge of this plot.  How many people are infected at the start of the plot?  It looks like zero, but we know there are over 100 thousand current cases.  It just looks like zero because the scale here is in billions, and 100 thousand is much too small to see on that scale!  This is a common problem with exponentially growing functions.  We can try to solve the problem by zooming in on the left end of the graph, showing results for just the next 30 days:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "SIR(T=30,plot='I',print_params=False);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It's a little better, but even on this scale the current number of infections is too small to be seen.  We can try again to fix it by just changing the vertical scale:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "SIR(T=30,plot='I',print_params=False,ymax=0.01);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can see that the starting value is not zero, but we can't see the right part of the graph at all!  Let's find a better solution."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Logarithmic scaling\n",
    "In the plots above, the scale of the $y$-axis is **linear**.  That means that equal distances in $y$ represent equal changes in the value of the function.  The problem with using this for an exponentially growing function is that the changes in the function at early times are tiny compared to the later growth.\n",
    "\n",
    "Instead, we can visualize the growth using a **logarithmic** scale:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def SIR_semilog(beta=0.25, gamma=0.05, N = 7770000000, confirmed_fraction=0.15,T=400,plot='SIR',\n",
    "        print_params=True,ymax=None):\n",
    "    \"\"\" Model the current outbreak using the SIR model.\"\"\"\n",
    "\n",
    "    du = np.zeros(3)\n",
    "    u0 = np.zeros(3)\n",
    "    \n",
    "    def f(t,u):\n",
    "        du[0] = -beta*u[1]*u[0]/N\n",
    "        du[1] = beta*u[1]*u[0]/N - gamma*u[1]\n",
    "        du[2] = gamma*u[1]\n",
    "        return du\n",
    "\n",
    "    # Initial values\n",
    "    total_cases = 199258./confirmed_fraction\n",
    "    u0[2] = 81972./confirmed_fraction+7956  # Initial recovered\n",
    "    u0[1] = total_cases-u0[2] # Initial infected\n",
    "    u0[0] = N - u0[1] - u0[2]\n",
    "\n",
    "    #dt = 0.01\n",
    "    times = np.arange(0,T)\n",
    "    solution = solve_ivp(f,[0,T],u0,t_eval=times)\n",
    "    S = solution.y[0,:]\n",
    "    I = solution.y[1,:]\n",
    "    R = solution.y[2,:]\n",
    "\n",
    "    today = '19/03/2020'\n",
    "    start = dates.datestr2num(today)\n",
    "    mydates = np.arange(T)+start\n",
    "    \n",
    "    fig = plt.figure(figsize=(12,8)); leg=[]\n",
    "    if 'S' in plot:\n",
    "        leg.append('Susceptible')\n",
    "        plt.semilogy(mydates,S/1.e9,'-b',lw=3)\n",
    "    if 'I' in plot:\n",
    "        leg.append('Infected')\n",
    "        plt.semilogy(mydates,I/1.e9,'-',lw=3,color='brown')\n",
    "    if 'R' in plot:\n",
    "        leg.append('Recovered')\n",
    "        plt.semilogy(mydates,R/1.e9,'-g',lw=3)\n",
    "    ax = plt.gca()\n",
    "    ax.xaxis.set_major_locator(dates.MonthLocator())\n",
    "    ax.xaxis.set_major_formatter(dates.DateFormatter('%b'))\n",
    "    fig.autofmt_xdate()\n",
    "    plt.legend(leg)\n",
    "    plt.ylabel('Individuals (in billions)')\n",
    "    plt.xlim(737503,737503+T)\n",
    "    if ymax: plt.ylim(-ymax/10,ymax)\n",
    "    if print_params:\n",
    "        plt.text(737503+3*T/4, 6,r'$\\beta$ = {}, $\\gamma$ = {}'.format(beta,gamma),fontsize=15,bbox=dict(facecolor='white', alpha=1.0))\n",
    "    return fig"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "SIR_semilog(plot='I',T=40,print_params=False);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Look carefully at the $y$ axis here.  As you can see, each of the evenly space labels on the axis represents a value **10 times greater** than the value below it.  In other words, equal distances on this axis represent equal **ratios**.  The great thing about this scaling is that we can easily see how the function varies over the whole graph.  It might seem strange that it looks almost like a straight line, but that's exactly how an exponentially growing function should look on this kind of plot.  Remember, equal distances in $y$ represent equal ratios, and we said that this kind of function doubles over each interval of some fixed size."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With this plot, it's easy to answer questions like \"when do we expect to have more than 1 million cases?\"  We couldn't possibly answer that by looking at the previous plots.\n",
    "\n",
    "We can also see more easily how rapidly the epidemic goes from something small to a true global crisis.  **At present there are less than a million cases, but before the end of April we would have (without mitigation) over 1 billion**.  Never in history have that many human beings been ill at the same time with a single disease."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In technical terms, we say this is a **semi-log plot**, because the $y$ axis is logarithmic while the $x$ axis is linear.  If both the $x$ and $y$ axes were logarithmic, we would call it a **log-log plot**.\n",
    "\n",
    "Understanding logarithmic plots is a bit of a superpower, because it allows you to study data by looking at plots like the last one above, which let you simultaneously see the parts where the function is small and big."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# The model predictions on a log scale"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here's a semilog plot of our basic model for the epidemic over the next year:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "SIR_semilog(plot='SIR',T=365,print_params=False);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Remember, this is exactly the same data that's in the first plot near the top of this post.  We're just looking at it in a different way.  But with this new plot we can much more easily see how many susceptibles remain at the end of a year: about a tenth of a billion, or 100 million.  \n",
    "\n",
    "We can also see that there are still about 1000 infected people in this model at the end of a year.  Notice that after the infection peaks, it declines in a way that also looks like a straight (downward-trending) line on this plot; that means that the decrease is also **exponential**.  In other words, after we pass the peak of infection, the number of infected individuals will consistently reduce by a factor of two over a certain time interval.  Looking back at the model, we see that the rate of decrease is determined by $\\gamma$.  In the late stages of the epidemic, because the fraction of susceptible people is small, we have approximately:\n",
    "\n",
    "$$\n",
    "\\frac{dI}{dt} = -\\gamma I(t)\n",
    "$$\n",
    "\n",
    "whose solution is\n",
    "\n",
    "$$\n",
    "I(t+\\tau) = e^{-\\gamma \\tau} I(t),\n",
    "$$\n",
    "\n",
    "which, for our estimate $\\gamma \\approx 0.05$, implies that the number of infections is reduced by half about every 14 days."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Assessing the effectiveness of mitigation\n",
    "Armed with our new superpower, let's look at our data from specific countries with this logarithmic scaling.  We'll start again with Italy."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "url = \"https://raw.githubusercontent.com/CSSEGISandData/COVID-19/master/csse_covid_19_data/csse_covid_19_time_series/time_series_19-covid-Confirmed.csv\"\n",
    "data = pd.read_csv(url)\n",
    "today = '3/19/20'  # Update this to include more data\n",
    "days = pd.date_range(start='1/22/20',end=today)\n",
    "dd = np.arange(len(days))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "row=data['Country/Region'].tolist().index('Italy')\n",
    "total_cases = [data[day.strftime('%-m/%-d/%y')][row] for day in days]\n",
    "plt.semilogy(total_cases[:],'.k');\n",
    "plt.xlabel('Days from January 22nd'); plt.ylabel('Total cases');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Several things become clearer with this scaling.  Notice that the plot is **not** a straight line, but we can make good guesses as to why.  Before day 30 (February 21st), there were only one or two known cases, and after day 30 there was a very abrupt increase.  It seems likely that the virus was spreading before day 30 and the new cases were only detected later.  This matches with [statements from experts in Italy](https://en.wikipedia.org/wiki/2020_coronavirus_pandemic_in_Italy#cite_note-250).\n",
    "\n",
    "Recall our full model equation for $I(t)$:\n",
    "\n",
    "$$\n",
    "\\frac{dI}{dt} = (\\beta \\frac{S}{N}-\\gamma) I \\\\\n",
    "$$\n",
    "\n",
    "Since only a tiny fraction of the whole Italian population is infected, we have $S/N\\approx 1$, so the slope of our exponential growth line on a semilog plot should be $\\beta-\\gamma \\approx 0.2$.  Let's see how a line with that slope matches the data:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "beta=0.25; gamma=0.05;\n",
    "plt.semilogy(total_cases[:],'.k');\n",
    "plt.xlabel('Days from January 22nd'); plt.ylabel('Total cases');\n",
    "beta=0.25; gamma=0.05;\n",
    "plt.semilogy(np.arange(7,58),np.exp((beta-gamma)*dd[7:58])/1.5);\n",
    "plt.legend(['Infected cases',r'Exponential growth with rate $\\beta-\\gamma = 0.2$']);\n",
    "plt.title('Italy');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that we are doing essentially the same thing that we did when trying to determine $\\beta$ in the second post of this series; but now we are looking at the results on a log scale to reveal more detail.\n",
    "\n",
    "It seems plausible (not utterly convincing) that the virus has been spreading at this expected rate in Italy for about 50 days now.  However, notice that in the last week the slope seems to have decreased.  Since Italy is now making great efforts to detect all new cases, it seems most likely that this is the result of mitigation.  It's probably too soon to try to assess the effectiveness of that mitigation, but let's make an attempt anyway:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.semilogy(total_cases[:],'.k');\n",
    "plt.xlabel('Days from January 22nd'); plt.ylabel('Total cases');\n",
    "beta=0.25; gamma=0.05;\n",
    "plt.semilogy(np.arange(7,50),np.exp((beta-gamma)*dd[7:50])/1.5);\n",
    "plt.semilogy(np.arange(51,58),np.exp(0.14*dd[51:58])*14);\n",
    "plt.legend(['Infected cases',\n",
    "            r'Exponential growth with rate $\\beta-\\gamma = 0.2$',\n",
    "            r'Exponential growth with rate $0.14$']);\n",
    "plt.title('Italy');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We see that the slope over the last several days is about 0.14, corresponding to a mitigation factor $q \\approx 0.7$ in the model I introduced in the last post."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "row=data['Country/Region'].tolist().index('Spain')\n",
    "total_cases = [data[day.strftime('%-m/%-d/%y')][row] for day in days]\n",
    "\n",
    "plt.semilogy(total_cases[:],'.k');\n",
    "plt.xlabel('Days from January 22nd'); plt.ylabel('Total cases');\n",
    "beta=0.25; gamma=0.05;\n",
    "plt.semilogy(np.arange(8,58),np.exp((beta-gamma)*dd[8:58])/5);\n",
    "plt.legend(['Infected cases',r'Exponential growth with rate $\\beta-\\gamma = 0.2$']);\n",
    "plt.title('Spain');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For Spain, we see a similar pattern, but no sign of any impact from mitigation yet."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "row=data['Country/Region'].tolist().index('Korea, South')\n",
    "total_cases = [data[day.strftime('%-m/%-d/%y')][row] for day in days]\n",
    "\n",
    "plt.semilogy(total_cases[:],'.k');\n",
    "plt.xlabel('Days from January 22nd'); plt.ylabel('Total cases');\n",
    "beta=0.25; gamma=0.05;\n",
    "plt.semilogy(np.arange(0,58),np.exp((beta-gamma)*dd[0:58])/0.7);\n",
    "plt.legend(['Infected cases',r'Exponential growth with rate $\\beta-\\gamma = 0.2$']);\n",
    "plt.title('Korea, South');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For South Korea, up to about day 40 we see the same pattern (with again a similar plateau in late February followed by a rapid rise when testing increased).  But afterward we see completely different behavior, as the curve flattens.  Given that South Korea has perhaps the most agressive testing strategy in the world, it seems unlikely that this can be attributed to new infections going undetected.  Instead, the evidence suggests that mitigation strategies have been very successful.  We can measure this success by looking at how much the slope of the curve has changed:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "row=data['Country/Region'].tolist().index('Korea, South')\n",
    "total_cases = [data[day.strftime('%-m/%-d/%y')][row] for day in days]\n",
    "\n",
    "plt.semilogy(total_cases[:],'.k');\n",
    "plt.xlabel('Days from January 22nd'); plt.ylabel('Total cases');\n",
    "beta=0.25; gamma=0.05;\n",
    "plt.semilogy(np.arange(0,42),np.exp((beta-gamma)*dd[0:42])/0.7);\n",
    "plt.semilogy(np.arange(44,58),np.exp(0.02*dd[44:58])*2800);\n",
    "plt.legend(['Infected cases',\n",
    "            r'Exponential growth with rate $\\beta-\\gamma = 0.2$',\n",
    "            r'Exponential growth with rate $0.02$']);\n",
    "plt.title('Korea, South');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "An approximate fit to the data from recent days suggests that the growth rate has been cut to approximately $0.02$; this would correspond to a value of $q$ (from my previous post) of about 1/4, meaning each infected person on average transmits the disease to only 1/4 as many people as they naturally would.  But again it is probably too early to estimate this number with confidence."
   ]
  }
 ],
 "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
}