compartmental_models/SEIRD Simulator in Python.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Simulation of SEIRD Model\n",
    "\n",
    "SEIRD is a type of compartmental models which are used in modelling of infectious disease using differential equations. These types of models divide the population into groups or compartments and the dynamics of these groups are expressed with the help of a system of differential equations.\n",
    "\n",
    "These system of equations are parametrized which capture the mechanistic nature of the disease. For simulation, you select values of these parameters and the resulting curves simulate the behaviour by solving the set of equations. Finally the results are plotted in a graph to visually understand the effect of the parameters.\n",
    "\n",
    "## SEIRD Model\n",
    "\n",
    "For completeness the SEIR model is produced below:\n",
    "\n",
    "<img src=\"images/seird_model.png\">\n",
    "\n",
    "$\\displaystyle \\frac{dS}{dt} = -\\frac{\\beta S I}{N}$<br><br>\n",
    "$\\displaystyle \\frac{dE}{dt} = \\frac{\\beta S I}{N} - \\sigma E$<br><br>\n",
    "$\\displaystyle \\frac{dI}{dt} = \\sigma E -  \\gamma I - \\mu I$<br><br>\n",
    "$\\displaystyle \\frac{dR}{dt} = \\gamma I$<br><br>\n",
    "$\\displaystyle \\frac{dD}{dt} = \\mu I$<br><br>\n",
    "$N = S + E + I + R + D$<br><br>\n",
    "Where,<br><br>\n",
    "$\\beta$ is infection rate or the rate of spread<br><br>\n",
    "$\\sigma$ is the incubation rate or the rate of latent individuals becoming infectious (average duration of incubation is $1/\\sigma$)<br><br>\n",
    "$\\gamma$ is the recovery rate or mortality rate. If the duration of indection is D then $\\gamma$ = 1/D<br><br>\n",
    "$\\mu$ is the mortality rate due to the disease\n",
    "\n",
    "#### A supplementary post to this notebook can be read at [Simulating Compartmental Models in Epidemiology using Python & Jupyter Widgets](https://dilbertai.com/2020/04/12/simulating-compartmental-models-in-epidemiology-using-python-jupyter-widgets/)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import sys\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "from scipy.integrate import odeint\n",
    "import plotly.graph_objects as go\n",
    "import plotly.io as pio\n",
    "pio.renderers.default = \"notebook\"\n",
    "%matplotlib inline\n",
    "plt.style.use('ggplot')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Jupyter Specifics\n",
    "from IPython.display import HTML\n",
    "from ipywidgets.widgets import interact, IntSlider, FloatSlider, Layout\n",
    "\n",
    "style = {'description_width': '100px'}\n",
    "slider_layout = Layout(width='99%')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def ode_model(z, t, beta, sigma, gamma, mu):\n",
    "    \"\"\"\n",
    "    Reference https://www.idmod.org/docs/hiv/model-seir.html\n",
    "    \"\"\"\n",
    "    S, E, I, R, D = z\n",
    "    N = S + E + I + R + D\n",
    "    dSdt = -beta*S*I/N\n",
    "    dEdt = beta*S*I/N - sigma*E\n",
    "    dIdt = sigma*E - gamma*I - mu*I\n",
    "    dRdt = gamma*I\n",
    "    dDdt = mu*I\n",
    "    return [dSdt, dEdt, dIdt, dRdt, dDdt]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def ode_solver(t, initial_conditions, params):\n",
    "    initE, initI, initR, initN, initD = initial_conditions\n",
    "    beta, sigma, gamma, mu = params\n",
    "    initS = initN - (initE + initI + initR + initD)\n",
    "    res = odeint(ode_model, [initS, initE, initI, initR, initD], t, args=(beta, sigma, gamma, mu))\n",
    "    return res"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "# ref: https://www.medrxiv.org/content/10.1101/2020.04.01.20049825v1.full.pdf\n",
    "initN = 1380000000\n",
    "# S0 = 966000000\n",
    "initE = 1\n",
    "initI = 1\n",
    "initR = 0\n",
    "initD = 0\n",
    "sigma = 1/5.2\n",
    "gamma = 1/2.9\n",
    "mu = 0.034\n",
    "R0 = 4\n",
    "beta = R0 * gamma\n",
    "days = 180"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [],
   "source": [
    "def main(initE, initI, initR, initD, initN, beta, sigma, gamma, mu, days):\n",
    "    initial_conditions = [initE, initI, initR, initN, initD]\n",
    "    params = [beta, sigma, gamma, mu]\n",
    "    tspan = np.arange(0, days, 1)\n",
    "    sol = ode_solver(tspan, initial_conditions, params)\n",
    "    S, E, I, R, D = sol[:, 0], sol[:, 1], sol[:, 2], sol[:, 3], sol[:, 4]\n",
    "    \n",
    "    # Create traces\n",
    "    fig = go.Figure()\n",
    "    fig.add_trace(go.Scatter(x=tspan, y=S, mode='lines+markers', name='Susceptible'))\n",
    "    fig.add_trace(go.Scatter(x=tspan, y=E, mode='lines+markers', name='Exposed'))\n",
    "    fig.add_trace(go.Scatter(x=tspan, y=I, mode='lines+markers', name='Infected'))\n",
    "    fig.add_trace(go.Scatter(x=tspan, y=R, mode='lines+markers',name='Recovered'))\n",
    "    fig.add_trace(go.Scatter(x=tspan, y=D, mode='lines+markers',name='Death'))\n",
    "    \n",
    "    if days <= 30:\n",
    "        step = 1\n",
    "    elif days <= 90:\n",
    "        step = 7\n",
    "    else:\n",
    "        step = 30\n",
    "    \n",
    "    # Edit the layout\n",
    "    fig.update_layout(title='Simulation of SEIRD Model',\n",
    "                       xaxis_title='Day',\n",
    "                       yaxis_title='Counts',\n",
    "                       title_x=0.5,\n",
    "                      width=900, height=600\n",
    "                     )\n",
    "    fig.update_xaxes(tickangle=-90, tickformat = None, tickmode='array', tickvals=np.arange(0, days + 1, step))\n",
    "    if not os.path.exists(\"images\"):\n",
    "        os.mkdir(\"images\")\n",
    "    fig.write_image(\"images/seird_simulation.png\")\n",
    "    fig.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "ed764ea271e643938faad7a3d9421b66",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "interactive(children=(IntSlider(value=1, description='initE', layout=Layout(width='99%'), max=100000, style=Sl…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "interact(main, \n",
    "         initE=IntSlider(min=0, max=100000, step=1, value=initE, description='initE', style=style, layout=slider_layout),\n",
    "         initI=IntSlider(min=0, max=100000, step=10, value=initI, description='initI', style=style, layout=slider_layout),\n",
    "         initR=IntSlider(min=0, max=100000, step=10, value=initR, description='initR', style=style, layout=slider_layout),\n",
    "         initD=IntSlider(min=0, max=100000, step=10, value=initD, description='initD', style=style, layout=slider_layout),\n",
    "         initN=IntSlider(min=0, max=1380000000, step=1000, value=initN, description='initN', style=style, layout=slider_layout),\n",
    "         beta=FloatSlider(min=0, max=4, step=0.01, value=beta, description='Infection rate', style=style, layout=slider_layout),\n",
    "         sigma=FloatSlider(min=0, max=4, step=0.01, value=sigma, description='Incubation rate', style=style, layout=slider_layout),\n",
    "         gamma=FloatSlider(min=0, max=4, step=0.01, value=gamma, description='Recovery rate', style=style, layout=slider_layout),\n",
    "         mu=FloatSlider(min=0, max=1, step=0.01, value=mu, description='Mortality rate', style=style, layout=slider_layout),\n",
    "         days=IntSlider(min=1, max=600, step=7, value=days, description='Days', style=style, layout=slider_layout)\n",
    "        );"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**References:**<br>\n",
    "1. SEIR and SEIRS Model https://www.idmod.org/docs/hiv/model-seir.html<br>\n",
    "2. Compartmental models in epidemiology https://en.wikipedia.org/wiki/Compartmental_models_in_epidemiology#The_SEIR_model<br>\n",
    "3. Solve Differential Equations in Python https://www.youtube.com/watch?v=VV3BnroVjZo<br>\n",
    "4. Computational Statistics in Python https://people.duke.edu/~ccc14/sta-663/CalibratingODEs.html<br>\n",
    "5. Ordinary Differential Equations (ODE) with Python and Jupyter https://elc.github.io/posts/ordinary-differential-equations-with-python/<br>\n",
    "6. SEIRS+ Model https://github.com/ryansmcgee/seirsplus<br>\n",
    "7. Stack Overflow https://stackoverflow.com/questions/40753159/why-is-scipy-minimize-ignoring-my-constraints<br>\n",
    "8. Lotka–Volterra equations https://en.wikipedia.org/wiki/Lotka%E2%80%93Volterra_equations<br>\n",
    "9. SEIR and Regression Model based COVID-19 outbreak predictions in India https://www.medrxiv.org/content/10.1101/2020.04.01.20049825v1.full.pdf<br>\n",
    "\n",
    "A simulator built with RShiny which provides many more parameters https://alhill.shinyapps.io/COVID19seir/"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}