{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Simulation of Epidemological SIR-Model\n", "\n", "`English version last change on 2020-08-13; German version created on 2020-04-02; author: Olaf Behrendt`\n", "\n", "**Contents**\n", "\n", "* [Preface](#Preface)\n", "* [Introduction](#Introduction)\n", " * [An Illustrative First Example](#An-Illustrative-First-Example)\n", " * [Motivation of State Transition Rates](#motivation-of-state-transition-rates)\n", "* [Dynamics of the SIR-model](#Dynamics-of-the-SIR-model)\n", " * [System of Differential Equations](#System-of-Differential-Equations)\n", " * [Verbalization of SIR system of ODEs](#Verbalization-of-SIR-system-of-ODEs)\n", " * [Euler Method](#Euler-Method)\n", "* [Implementation of Euler's Method](#Implementation-of-Euler's-Method)\n", " * [Preparations](#Preparations)\n", " * [Calculate Time Step](#Calculate-Time-Step)\n", " * [Main Function and Plotting](#Main-Function-and-Plotting)\n", "* [Simulations](#Simulations)\n", " * [Parameters for SARS-CoV-2](#Parameters-for-SARS-CoV-2)\n", " * [Running the Simulations](#Running-the-Simulations)\n", " * [Summary of Scenario Simulations](Summary-of-Scenario-Simulations)\n", "* [Glossary](#Glossary)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Preface\n", "\n", "In the following presentation we introduce a classical model for the development of an infectious deseases, the SIR-model. It was first published in 1927 in the proceedings of the Royal Society by Scottish biochemist [William Ogilvy Kermack](https://en.wikipedia.org/wiki/William_Ogilvy_Kermack) and [Anderson Gray McKendrick](https://en.wikipedia.org/wiki/Anderson_Gray_McKendrick). McKendrick was a physican and did not have a mathematical degree, nervertheless he was a brilliant mathematician and pioneer in mathematical epidemiology. Apart from beeing a biochemist William Kermak also had a degree in mathematics. He survived the active service in the first world war, but unfortunately was blinded by a labratory experiment at the age of 26. \n", "\n", "I am trying to present the SIR-model in an most accessible way, so that anyone having had basic calculus at secondary school should have no difficulties. Specifically we some general familarity with simple equations, functions and elementary probability. But even if the technical hurdles are quite low, the reader probably has to spent some amount of time and concentration to follow the presentation. Anyways, if the reader has issues understanding the text, in doubt, the author has failed to clearly and consistently express his thoughts.\n", "\n", "An important part and probably the most \"enjoyable\" part of this notebook are the simulations using Python. The reader can change parameters and re-run simulation, so that she gets a better intuitive understanding of the quantitative nature of the model. The adventurous reader can even change or extend the SIR-model to his likings (e.g. add an compartmen/equation for the SEIR-model).\n", "\n", "## Introduction\n", "\n", "The SIR-model (**s**usceptible-**i**nfected-**r**emoved model) is a mathematical model describing the temporal devolpment of an infectious desease with immunization. An individual of the considered population of size $N$ resides in exactly one of three groups or compartments at a certain point of time:\n", "\n", "* $S(t)$ is the number of **susceptible** individuals at time point $t$. An individual is susceptible if she can be infected.\n", "* $I(t)$ is the number of **infected** individuals at time point $t$. \n", "* $R(t)$ is the number of **recovered** individuals at time point $t$. We assume that recovered individuals are **immune**, that is an recovered individual will remain in state _recovered_ forever.\n", "\n", "It is clear that the sum $S(t)+I(t)+R(t)=N$ is constant for all $t$ in the considered time range. Also note that we regard deceased individuals as _recovered_, which might seem a bit strange at first glance.\n", "\n", "### An Illustrative First Example\n", "\n", "At the start of the simulation $t=0$ we assume ten infected individuals in a population of size $N=1000$ (e.g. a village). So $S(0)=990$, $I(0)=10$, $R(0)=0$.\n", "Further we need some assumptions about the _contact rate_ $\\beta$ and _recovery rate_ $\\gamma$, two important terms we will shortly introduce in more detail.\n", "Please note that we cut of the digits after the decimal point, so we might observe some rounding errors.\n", "\n", "\n", "\n", "$$\n", "\\begin{array}{r|r|r|r|r|r}\n", "\\text{day } t & S(t) & I(t) & R(t) & \\text{new infects}\\\\\n", "\\hline\n", " 0 & 990 & 10 & 0 & 0\\\\\n", " 1 & 980 & 17 & 3 & 10\\\\\n", " 2 & 964 & 27 & 9 & 16\\\\\n", " 3 & 938 & 44 & 18 & 26\\\\\n", " 4 & 896 & 71 & 33 & 42\\\\\n", " 5 & 832 & 112 & 57 & 64\\\\\n", " 6 & 739 & 167 & 94 & 93\\\\\n", " 7 & 616 & 235 & 149 & 124\\\\\n", " 8 & 471 & 301 & 228 & 145\\\\\n", " 9 & 329 & 343 & 328 & 142\\\\\n", "10 & 216 & 341 & 442 & 113\\\\\n", "11 & 143 & 301 & 556 & 74\\\\\n", "12 & 100 & 244 & 657 & 43\\\\\n", "13 & 75 & 187 & 738 & 24\\\\\n", "14 & 61 & 139 & 800 & 14\\\\\n", "15 & 53 & 101 & 846 & 8\\\\\n", "16 & 47 & 73 & 880 & 5\\\\\n", "17 & 44 & 52 & 904 & 3\\\\\n", "18 & 42 & 37 & 921 & 2\\\\\n", "19 & 40 & 26 & 934 & 2\\\\\n", "20 & 39 & 18 & 942 & 1\\\\\n", "21 & 38 & 13 & 949 & 1\\\\\n", "22 & 38 & 9 & 953 & 0\\\\\n", "23 & 38 & 6 & 956 & 0\\\\\n", "\\vdots & \\vdots & \\vdots & \\vdots & \\vdots\\\\\n", "\\end{array}\n", "$$\n", "\n", "For more thourough background on compartemental models in general, please refer to [The Mathematics of Infectious Diseases, H.W. Hethcote, 2000 in SIAM Rev., 42(4), 599–653.](https://epubs.siam.org/doi/abs/10.1137/s0036144500371907)." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Assumptions and Definitions\n", "\n", "Let us now introduce a more exact description of the SIR-model together with some adequate assumptions in order to better understand its motivation and dynamics.\n", "\n", "* Following state transitions are possibly: $S\\rightarrow I\\rightarrow R$. In other words, if a susceptible individual in compartment $S$ is infected he moves into compartment $I$. Then, after some time, he recovers (or deceases) and finally moves into compartment $R$. Especially we assume that immunity does not wear off with time, that is there is no transition from state _recovered_ to state _susceptible_.\n", "* $S(t)+I(t)+R(t)=N$ (or all $t\\geq 0$). That means we do not consider births or deaths or any other demographic changes or travel etc..\n", "* If an individual is infectied he is instantly infectious. This is a simplification since usually an infected individual only becomes infectious after the [_pre-infectious period_](#latent-period).\n", "* Usually the unit of a time period is one day.\n", "* $\\pmb\\beta$ (beta) **contact rate** or spreading rate or **infection rate** is the average number of contacts of an infected individual per unit of time. E.g., a contact rate of $\\beta=1.0$ means that an infected individual on average infects one other person per day. The contact rate is implicitly composed of several factors, like\n", " * The (biological) infectivity of the pathogene\n", " * The number of contacts of an infected person with susceptible persons. Let us assume the infectivity is to transfer the desease is $0.25\\%$. Then the probability of infection when the infected individual meets two susceptible persons a day is $1-(0,75)^2\\approx 44\\%$, when meeting five susceptible persons a day it is about $76\\%$ and so on. Also note that if a large proportion of the population is immune, the contact rate will dramatically decrease, which is called [herd immunity](https://en.wikipedia.org/wiki/Herd_immunity).\n", " * There are a lot more factors effecting the contact rate, like duration of contact, proximity or locality (in the open or in a closed room), wearing face masks, ventilation, temperature, huminidy, etc..\n", " \n", " \n", "* $\\pmb\\gamma$ (gamma) **recovery rate** is the rate at which an infected person recovers. More precisely formulated, we can assume that the the time periods for a change from compartment $I$ to compartment $R$ are [exponetially distributed](https://en.wikipedia.org/wiki/Exponential_distribution). This means that the mean duration in compartment $I$ is equal to $1/\\gamma$ and that the probability after $t$ units of time is equal to $\\exp(-\\gamma t)$. For example if we assume $\\gamma=0.2$, then we have a mean duration period in compartment $I$ of $D_I=1/\\gamma=1/0.2=5$ days. The probabilty that a person is still infectious after one day is $exp(-0.2)\\approx 0.82$ and that person is still infectious after five days is $exp(-0.2\\cdot5)\\approx 0.37$ and so on. The technicalities are not so important for the further understanding. The take-away point here is that \n", "$$1/\\gamma = D_I$$\n", "constitutes the **mean infectious period**. Also note that the SIR-model does not explicitly model a [pre-infectious period](#pre-infectious-period) (but it can quite easily be extended into a SEIR-model).\n", "* The **[basic reproduction number](https://en.wikipedia.org/wiki/Basic_reproduction_number)** $R_0$ is defined as the average number of secondary infections produced by\n", "a single infected individual introduced into a completely susceptible population. If infection rate $\\beta$ and recovery rate $\\gamma$ are equal, then the overal number of infected persons is constant (on average). If $\\beta > \\gamma$ the average number of infections rises and for $\\beta < \\gamma$ the average number of infections drops. For the SIR-model it holds that \n", "\n", " $$\n", " R_0=\\frac{\\beta}{\\gamma}.\n", " $$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Motivation of State Transition Rates\n", "\n", "Let us denote the population size by $N$. \n", "Then the total number of transitions $S\\rightarrow I$, or in other words the number of individuals who get infected during a time period is given by\n", "\n", "$$S\\left[\\beta\\left(\\frac{I}{N}\\right)\\right] = \\frac{\\beta IS}{N}.$$\n", "\n", "But how can we motivate the left side of the above equation (the right side is just a simple reformulation)? Let's break down the expression on the left hand side of the above equation:\n", "\n", "* Since $N$ is the total population size we can interprete $\\frac{I}{N}$ as the [classical probability](https://en.wikipedia.org/wiki/Probability_interpretations#Classical_definition) (sometimes also called _Laplace probability_) for the random event _\"contact with an infected individual\"_. \n", "* It follows that $\\beta\\frac{I}{N}$ is the probability to meet an infected individual **and** to get infected.\n", "* But only susceptible individuals can be infected, so that the average number of newly infected individuals is given by $S\\left[\\beta\\left(\\frac{I}{N}\\right)\\right] = \\frac{\\beta IS}{N}$. Voilà.\n", " \n", "To sum up, per unit time period a total of $\\beta IS\\over N$ individuals transition from compartment $S$ into compartment $I$, which we write as $S\\stackrel{\\beta\\cdot IS\\over N}\\longrightarrow I$.\n", "The average number of recovered individuals per unit time period is given by $\\gamma I$ so we can write all transitions in the SIR-model as\n", "\n", "$$S\\stackrel{\\beta IS\\over N}\\longrightarrow I \\stackrel{\\gamma I}\\longrightarrow R.$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Dynamics of the SIR-model\n", "\n", "### System of Differential Equations\n", "\n", "We can model the [dynamics](#dynamical-system) of the SIR-model $-$ or more specifically the evolution of the numbers of individuals in the respective compartments as a function of time $-$ as a system of three [ordinary differential equations](#ordinary-differential-equation) (ODEs) as shown below. If you don't know what an ODE is, don't worry, just keep on reading, you should get an sufficient intuitive understanding.\n", "\n", "$$\n", "\\begin{aligned}\n", "{\\frac {dS(t)}{dt}} &= -{\\beta I(t)S(t)\\over N}\\\\\n", "{\\frac {dI(t)}{dt}} &= {\\beta I(t)S(t)\\over N}-\\gamma I(t)\\\\\n", "{\\frac {dR(t)}{dt}} &= \\gamma I(t)\n", "\\end{aligned}\n", "$$\n", "\n", "#### Verbalization of SIR system of ODEs\n", "\n", "We can regard $\\frac {df(t)}{dt}$ as the change in the value of the function $f$ in a (\"very short\") time period $dt$. So we can verbalize the three equations above in the following way:\n", "\n", "* The number of **susceptible** individuals $S(t)$ is reduced by the number of newly infected $\\beta I(t)S(t)\\over N$ during time period $dt$.\n", "* The number of **infected** $I(t)$ is increased by the number of newly infected $\\beta I(t)S(t)\\over N$ and decreased by the newly recovered individuals $\\gamma I(t)$ during time period $dt$.\n", "* The number of **recovered** $R(t)$ is increased by the number of newly recovered individuals $\\gamma I(t)$ during time period $dt$.\n", "\n", "Not that our requirenment $S(t)+I(t)+R(t)=N$ for all $t$ from section [Assumptions und Definitions](#Assumptions-and-Definitions) is fullfilled, since for each time step the total number of individuals does not change:\n", "\n", "$${dS(t)\\over dt}+{dI(t)\\over dt}+{dR(t)\\over dt}= -{\\beta I(t)S(t)\\over N} + \\left({\\beta I(t)S(t)\\over N}-\\gamma I(t)\\right) + \\gamma I(t) = 0.$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Euler Method\n", "\n", "How can we solve the SIR ODE system? It turns out we can use the [Euler method](https://en.wikipedia.org/wiki/Euler_method) (named after Swiss mathematician Leonhard Euler 1707-1783) to numerically solve the ODEs. \n", "\n", "Remember from school, that $\\frac {df(t)}{dt}=f'(t)$ is the [derivative](https://en.wikipedia.org/wiki/Derivative) of the function $f$. For example for a quadratic function $f(x)=ax^2+bx+c$ the derivative $f'$ is a straight line $\\left(f'(x)=2ax+b\\right)$. Locally we can approximate a function with the help of its first derivative: For a quadratic function, the first derivative is a straight line and we can approximate $f$ at any x-value $x_0$ by a straight line going through the coordinate $\\left(x_0,f(x_0)\\right)$ with gradient $f'(x_0)$. This is also called a (first-order) [Taylor approximation](https://en.wikipedia.org/wiki/Taylor%27s_theorem) (Brook Taylor was an English mathematician who lived 1685 - 1731). For a differentiable function $f$ the first Taylor approximation in point $x_0$ formally looks like this:\n", "\n", "$$\n", "f(x_0+h)\\;\\approx\\; f(x_0)+f'(x_0)h \\;\\overset{\\text{def}}=\\; f(x_0)+{df(x_0)\\over dt}h\n", "$$\n", "\n", "A picture may be better than a many words:\n", "\n", "
\n",
" \n",
"\n",
"