{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "Welcome to the third post in my series on modeling the Coronavirus pandemic. In the [first](http://www.davidketcheson.info/2020/03/17/SIR_model.html) and [second](http://www.davidketcheson.info/2020/03/19/SIR_Estimating_parameters.html) posts, we introduced the SIR model and used the available data to estimate the model parameters, in the absence of any mitigation (i.e. with no social distancing, quarantine, etc.). In this post we'll put model and parameters together to see what they predict for the current pandemic. We'll do this first assuming no mitigation, and then we'll make a very rough attempt to understand how mitigation may impact the predictions." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Remember, I am not an epidemiologist and I make no guarantees as to the accuracy of these predictions. I don't recommend that you make precise plans based on the details of these predictions. My goal is merely to show that with a bit of mathematics and common sense, we can have a reasonable idea of what the future holds." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Predictions in the absence of mitigation\n", "\n", "First, let's look at what we expect to happen without any mitigation. In other words, this is a scenario in which schools and workplaces remain open, people continue to shake hands or kiss cheeks when greeting one another, and so forth. **Note that the predictions here are NOT what we expect will actually happen, because we have intentionally ignored the effects of containment measures.**\n", "\n", "Here are the basic assumptions leading to the predictions below:\n", "\n", " - The dynamics of the COVID-19 epidemic follow the SIR model, with parameters $\\beta \\approx 0.25$ and $\\gamma \\approx 0.05$.\n", " - No containment measures (like quarantine, closures, and social distancing) are implemented.\n", " - The number of confirmed cases at present is about 15 percent of the actual cases (this is a very rough guess based on some expert opinions)." ] }, { "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", "from ipywidgets import interact, widgets\n", "import matplotlib.dates as dates\n", "from scipy.integrate import solve_ivp\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):\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", " T = 400\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(400)+start\n", " \n", " fig = plt.figure(figsize=(12,8))\n", " plt.plot_date(mydates,S/1.e9,'-',lw=3)\n", " plt.plot_date(mydates,I/1.e9,'-',lw=3)\n", " plt.plot_date(mydates,R/1.e9,'-',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(['Susceptible','Infected','Recovered'])\n", " plt.ylabel('Individuals (in billions)')\n", " plt.xlim(737503,737903)\n", " plt.text(737803, 6,r'$\\beta$ = {}, $\\gamma$ = {}'.format(beta,gamma),fontsize=15,bbox=dict(facecolor='white', alpha=1.0))" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "SIR(beta=0.25, gamma=0.05, confirmed_fraction=0.15);" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "There are several important things to notice here. \n", "\n", "## When does the maximum number of infected individuals occur, and how high is it?\n", "\n", "As we expect, there is an initial exponential growth of infection that eventually levels off and then decreases after most of the population has been infected. Recall from the first post that we expect the infection peak to occur when the fraction of susceptible individuals is $\\gamma/\\beta = 0.05/0.25 = 0.2$; i.e., when four-fifths of the world population has already been infected (including those who have recovered). With the current parameters, this peak occurs around the middle of May and results in almost 4 billion concurrent cases. Of course, the vast majority of those cases would not need medical care; based on our estimate of 15% reporting and less than 20% of reported cases being severe, perhaps only about 3% of all cases would require medical attention. Even so, that would mean a peak of 100 million cases requiring medical attention simultaneously, worldwide.\n", "\n", "## How many people catch the virus?\n", "\n", "In this scenario, almost everyone in the world is eventually infected; by day 400 only 54 million susceptible individuals remain.\n", "\n", "## How long does the epidemic last?\n", "\n", "After mid-May, the epidemic begins to trail off gradually, but it is not until about early November that the number of infected drops back to below 1 million. Thus (in this scenario) we should expect the severe epidemic to last at least several months. Of course, it would make sense for those who have already recovered from the virus -- which is most of the world population -- to go back to work/school in the early Fall, " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Variations on the average scenario\n", "\n", "In the last post, we saw that there is significant uncertainty in the value of $\\gamma$ and especially $\\beta$. How do these predictions change if we vary $\\beta$? We found values in the range $(0.2, 0.3)$, so let's look at what happens for each of the extremes of this interval:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "SIR(beta=0.2, gamma=0.05, confirmed_fraction=0.15)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "With a smaller value of $\\beta=0.2$, the virus spreads more slowly; the peak occurs near June 1st with just over 3 billion infected. Again, almost the entire world catches the virus (eventually 150 million susceptibles remain), and the number of cases does not drop below 1 million until November." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "SIR(beta=0.3, gamma=0.05, confirmed_fraction=0.15)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "With a larger values of $\\beta=0.3$, the virus spreads more quickly; the peak occurs in early May with just over 4 billion infected. The epidemic ends a bit sooner, with the number of cases dropping below 1 million some time in October. Only 19 million people remain susceptible after a year." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "As we can see, even with these fairly large changes in $\\beta$, the general picture remains the same. We can expect the epidemic to peak in May or June and last into the fall." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "interact(SIR,beta=widgets.FloatSlider(min=0.01,max=0.5,step=0.01,value=0.25,description=r'$\\beta$'),\n", " gamma=widgets.FloatSlider(min=0.01,max=0.5,step=0.01,value=0.05,description=r'$\\gamma$'),\n", " N = widgets.fixed(7770000000),\n", " confirmed_fraction=widgets.FloatSlider(min=0.01,max=1.0,step=0.01,value=0.15,description=r'Fraction C'));" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# The effect of mitigation" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "As I write this, almost every country in the world is adopting measures to mitigate the spread of the virus. This includes closing schools and workplaces, quarantining infected individuals and their contacts, and encouraging people to stay home." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "For simplicity we can view all of these mitigation measures as having a single effect: increasing the average time between encounters, or in other words reducing $\\beta$. Remember that $\\beta$ is the average number of people per day that a given individual has close contact with. Since $\\beta$ is a constant in our model but the mitigation techniques and their effectiveness may vary over time, we can incorporate mitigation by adding a new factor $q(t)\\in[0,1]$ multiplying $\\beta$:\n", "\n", "\\begin{align}\n", "\\frac{dS}{dt} & = -q(t)\\beta I \\frac{S}{N} \\\\\n", "\\frac{dI}{dt} & = q(t)\\beta I \\frac{S}{N}-\\gamma I \\\\\n", "\\frac{dR}{dt} & = \\gamma I\n", "\\end{align}\n", "\n", "What is the meaning of $q(t)$? If there were an absolute quarantine, with no human contact at all, we would have $q=0$, whereas if no measures are implemented then we would have $q=1$ (corresponding to the predictions above). In the real world, $q$ will be somewhere between these extremes." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "What is the correct value of $q$? Frankly, I have no idea and I doubt that even the experts can say with confidence. But we can hypothesize some values and explore their impact. For simplicity, we'll assume that some value $q<1$ is achieved through mitigation measures starting now and lasting for the next $N_q$ days. After $N_q$ days, these measures are lifted so that society (and $\\beta$) returns to normal.\n", "\n", "If we could maintain mitigation measures forever (i.e. $N_q=\\infty$) then this mitigation would have exactly the same effect as reducing $\\beta$. As we saw above, smaller values of $\\beta$ lead to a smaller infection peak, but a longer epidemic.\n", "\n", "Of course, we do not expect mitigation to last forever; people must go back to work eventually. This can have some surprising effects." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def SIR_mitigated(beta=0.25, gamma=0.05, N = 7770000000, confirmed_fraction=0.15, q=0.5, Nq=180):\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", " if t