DiffEq_Lab.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "metropolitan-cabinet",
   "metadata": {},
   "source": [
    "# Differential Equations Lab: The SIR Epidemic Model\n",
    "\n",
    "\n",
    "## Today's Lab\n",
    "\n",
    "#### In today's lab, you'll learn:\n",
    "\n",
    "1. What differential equations are, including:\n",
    " - the differences between independent variables, dependent variables and derivatives in a differential equation\n",
    " - assumptions made modeling with differential equations  \n",
    " (From MAA website)  \n",
    "2. To develop the SIR Model for the spread of an infectious disease, including the concepts of contact number and herd immunity;\n",
    "3. To develop a version of Euler's Method for solving a system of differential equations  \n",
    " Pre-reqs:  \n",
    " The concept of derivative and the Chain Rule;  \n",
    " The concept of autonomous first-order differential equation;  \n",
    " At least surface understanding of Euler's method as a tool for generating an approximate numerical or graphical solution of a first-order equation;  \n",
    "Some experience with the selected CAS.\n",
    "\n",
    "Be aware that there are many different models that are used in mathematical epidemiology today (e.g. SI, SIS, SIR, SIRS, SEIR, SEAIR, SEIRD, and SIR with vital dynamics, along with others) and that this lab will focus specifically on the SIR model.\n",
    "\n",
    "\n",
    "This lab covers the SIR definition of the MAA website located here:   \n",
    "https://www.maa.org/press/periodicals/loci/joma/the-sir-model-for-spread-of-disease-introduction  \n",
    "You should check out the MAA website, but not right now. Instead, let's get started!  "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "continent-daisy",
   "metadata": {},
   "source": [
    "#### Background\n",
    "\n",
    "A mathematical description of the spread of a disease in a population is the SIR model, which divides a fixed population of $N$ individuals into three \"compartments\" which may vary as a function of time, $t$:  \n",
    "\n",
    "$S(t)$ are those susceptible but not yet infected with the disease;  \n",
    "$I(t)$ is the number of infectious individuals;  \n",
    "$R(t)$ are those individuals who have recovered from the disease and now have immunity to it.  \n",
    "\n",
    "The fraction of susceptible, infectious and recovered people in the population is:   \n",
    "$s(t) = \\frac{S(t)}{N}$  \n",
    "$i(t) = \\frac{I(t)}{N}$  \n",
    "$r(t) = \\frac{R(t)}{N}$  \n",
    "\n",
    "The SIR model describes the change in the population of each of these compartments in terms of two parameters, $\\beta$ and $\\gamma$, which are both positive constants. $\\beta$ describes the effective contact rate of the disease: an infected individual comes into contact with $\\beta N$ other individuals per unit time (of which the fraction that are susceptible to contracting the disease is $S/N$). $\\gamma$ is the mean recovery rate: that is, $1 / \\gamma$ is the mean period of time during which an infected individual can pass it on.\n",
    "\n",
    "<img src=\"SIR.png\"/>  \n",
    "\n",
    "The following assumptions apply to the SIR model:  \n",
    "\n",
    "No one is added to the susceptible group, since we are ignoring births and immigration. The only way an individual leaves the susceptible group is by becoming infected. We assume that the time-rate of change of  $S(t)$,  the number of susceptibles, depends on the number already susceptible, the number of individuals already infected, and the amount of contact between susceptibles and infecteds. In particular, suppose that each infected individual has a fixed number  $\\beta$  of contacts per day that are sufficient to spread the disease. Not all these contacts are with susceptible individuals. If we assume a homogeneous mixing of the population, the fraction of these contacts that are with susceptibles is  $s(t)$.  Thus, on average, each infected individual generates  $ \\beta s(t)$  new infected individuals per day. [With a large susceptible population and a relatively small infected population, we can ignore tricky counting situations such as a single susceptible encountering more than one infected in a given day.]\n",
    "\n",
    "We also assume that a fixed fraction  $\\gamma$  of the infected group will recover during any given day. This model considers deaths due to an epidemic as part of the \"recovered\" group. A \"recovered\" person can still feel miserable, and might even die later from the epidemic or complications.\n",
    "\n",
    "From these assumptions, we get the following equations\n",
    "\n",
    "$ \\frac{dS}{dt}=- \\beta SI, $  \n",
    "$ \\frac{dI}{dt}=\\beta SI -\\gamma I, $  \n",
    "$ \\frac{dR}{dt}=\\gamma I $  \n",
    "\n",
    "[//]: # (Should the above equations be divided by the population N? Probably not, gonna change that)\n",
    "\n",
    "The change in the susceptible, infected, and recovered fractions of the population is given by the following equations\n",
    "\n",
    "$ \\frac{ds}{dt}=- \\beta s(t) i(t), $  \n",
    "$ \\frac{di}{dt}=\\beta s(t) i(t) - \\gamma i(t), $  \n",
    "$ \\frac{dr}{dt}= \\gamma i(t) $  "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ceramic-patrick",
   "metadata": {},
   "source": [
    "#### Below is a graph of an SIR model, feel free to play around with the parameters."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "municipal-lecture",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAEHCAYAAACp9y31AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Z1A+gAAAACXBIWXMAAAsTAAALEwEAmpwYAABL7UlEQVR4nO3dd3xUdbr48c/UTHpCQkhIQgmBNDqE3qV3FFlw14Luxbarrru6etfrRa8rrP70rle8d3ddl9V1KYoFWIrSpEMoQXoIIZEUIKS36XN+f4wZiSSEkEwm5Xnzmte0U56TQ55z8j3f83xViqIoCCGEaDfUng5ACCFE85LEL4QQ7YwkfiGEaGck8QshRDsjiV8IIdoZSfxCCNHOuC3xP/zww4SFhdG7d+9av1cUhaeeeorY2Fj69u3L8ePH3RWKEEKIG7gt8T/00ENs3bq1zu+3bNlCeno66enp/OUvf+Hxxx93VyhCCCFu4LbEP2bMGDp06FDn9+vXr+eBBx5ApVIxbNgwSkpKuHLlirvCEUII8T2PtfHn5uYSHR3teh8VFUVubq6nwhFCiHZD6+kA3OXyZfjb3xo+n8EAHTo4H6Gh0KkThIdDUBCoVE0ephBCNDuPJf7IyEiys7Nd73NycoiMjKx12kuXLjV4+Q4HvPBCDGYzXLiQhdWqwmz+4WEyqamqUmE0qqmsVFFRoaaiQk1xsYraWpy8vBx07myjc2cbkZE2wsPtLfJAEBMTA9zZz6y1am/b3N62F2SbGzN/bTyW+GfPns2KFStYuHAhhw8fJjAwkIiIiCZbvlrtPHs3GCAkxHFb8ygKGI0qSkrUlJRoKCrScP26hoICDVVVajIz9WRm6gHw9nYQE2OlRw8rXbpY0WiaLHQhhHArtyX+RYsW8c0331BQUEBUVBSvvPIKVqsVgMcee4zp06ezefNmYmNj8fHxYeXKle4K5bapVODjo+DjY6dzZ3uN78rLVeTm6sjN1XL5spayMg1nznhx5owXvr4OEhPN9O5tISDg9g4yQgjhKW5L/KtXr77l9yqVivfee89dq29y/v4K8fEW4uMtKAoUFmrIyNCRlqanuFjDkSPeHD1qoFcvK0OHGgkOlgOAEKJlarMXd91JpYLQUDuhoXaGDDGRl6fl9Gk9Fy7oSUvTk56uIz7ewrBhRvz9ZbiDtkpRFGw2W7OvNy8vD8D1F3R7INt8a1qtFlUDLjpK4m8klQoiI50XfIcPN5GSYuDcOT1nz3qRkaFj1CgjSUmWFnkhWDSOzWbD398fg8HQoF+6xgoODgaguLi42dbpabLNtVMUBZPJRHl5OTqd7raXLbV6mlBAgIOJE6u4//4yune3YDar2bHDly++8KO8XDJ/W9TcSV+IG6lUKgwGQ4Pnk8TvBkFBDmbNqmTq1AoMBgfZ2TpWrw4gL0+6/rQ1kvSFp93J/0FJ/G6iUkFcnJX77y+jSxcrRqOazz7z58wZvadDE23MW2+9xfDhwxk1ahRjxozh6NGjHonj1KlTbNu2zfV+y5Yt/PGPfwTgySefZP369TfNs2/fPhYuXNhcIYrvSRu/m/n4KMyZU8Hevd6cOGFg+3ZfSkvVDB9uknZ/0WgpKSl89dVXfPPNN3h5eVFYWIjFYvFILKdOneLEiRNMmjQJgGnTpjFt2jSPxCJuTc74m4FaDWPHGrnrrkrUaoUjR7xJSWl4u5wQP3bt2jVCQkLw8vICICQkhIiICPr160dhYSEAqampzJo1C4D9+/czZswYxowZw9ixYykvLwfgnXfeYeTIkYwePZpXXnkFgMzMTObPn8/48eOZPn06Fy5cAJxn788++ywTJkwgOTmZr776CovFwrJly/jiiy8YM2YMn3/+OatWreL55593xbp79+4a8/xYZWUlv/jFL5g4cSJjx45l8+bN7vvBtXNyxt+Meve2oNcrbN3qy6FD3mg0CoMHmz0dlmgCy5e750D+wgumW34/fvx43nzzTZKTkxk7dizz5s1j5MiRdU6/YsUK3njjDYYNG0ZFRQUGg4Ft27axefNmtm3bho+Pj6sXya9+9SveeustevTowdGjR3nuuedczTWXL19m+/btZGZmMnv2bI4dO8aLL77IiRMneOONNwBYtWpVjXXXNs+N3n77bcaMGcOKFSsoLS11HQB8fX0b/HMTtyaJv5n16mXFbq9i2zYf9u/3Qa+Hvn0l+Ys74+fnx65duzh48CB79+7lkUce4eWXX65z+qFDh/LSSy9x7733MnPmTPz8/Ni9ezf33XcfPj4+gLMbYUVFBSkpKSxevNg1r9n8w//TuXPnolar6dGjB926dSM9Pb3eWOubZ9euXWzZsoUVK1YAYDKZyMnJIS4urkE/E1E/SfwekJBgwW6HHTt82b3bmw4d7ERFNf+NQKLp1Hdm7k4ajYZRo0YxatQoEhMTWbNmDVqtFofDeff4jQn7mWeeYfLkyWzbto1p06axbt26WpfpcDgIDAxkz549tX7/454kt9OzpL55FEXhww8/pGfPnvUuSzSOtPF7SO/eFgYPNuFwqNi82Vf6+Ys7kp6eTkZGhuv96dOniY6OpkuXLpw4cQKADRs2uL7PzMwkMTGRp59+mgEDBpCens64ceNYtWoVVVVVgPOGoYCAALp06cKXX34JOJPy6dOnXctZv349DoeDzMxMsrKyiI2Nxc/Pj4qKijpjrW2eG02YMIH3338fRXHe7X7y5MlG/WxE3STxe9Dw4Uaio51dPTdv9sNur38eIW5UWVnJE088wbBhwxg1ahRpaWn89re/5fnnn+ff//3fmTBhApobSsf+6U9/YsSIEYwaNQqdTsfEiROZOHEi06ZNY8KECa42doC//OUvfPzxx4wePZrhw4fXuNgaFRXFxIkTWbBgAW+//TYGg4HRo0eTlpbmurj7Y7XNc6Pf/OY3WK1WRo0axfDhw3n99dfd9FMTKqX68NqCNbYedUuu4V1VpWLNmgDKy9X062di3Dhjo5bXGra5qXlqm61WK2FhYc26TvB8+YInn3ySyZMnM2fOnGZbp6e32RMass35+fk3lWy4VT1+OeP3MB8fhenTK1CrFb791kBOjlx2EUK4lyT+FiA83E5ysvPi4I4dPrSjAoSiFXrvvfea9WxfND1J/C3E4MEmOnSwU1KiISXF29PhCCHaMEn8LYRWC3fdVYlKpXD8uBf5+VLQTQjhHpL4W5DOne307WvG4VCxY4cPLf+yuxCiNZLE38KMGGHE19dBfr6WCxduf2AFIYS4XZL4Wxi9HoYNc3bpPHDAW/r2i1uKjo6ud5qDBw8yfPhwxowZg9HYsO7CmzZt4vz5826JS3iOJP4WKDHRQocOdsrKNJw65eXpcEQr9+mnn/KrX/2KPXv24O3dsI4DmzdvJi0tzU2RCU+RxN8CqdXOJh+Aw4cNmKWGm6jHvn37mDVrFg8++CBDhw5lyZIlKIrCRx99xJdffsnrr7/OkiVLAPif//kf7rrrLkaNGsWyZctcy1izZg2jRo1i9OjRPPbYYxw+fJgtW7bwn//5n4wZM4bMzMw6SzV/9913TJ48mZEjR/L73//eIz8DcfvkbqEWKibGSkSEjStXtBw/bmD4cM8VARP1Myxf7pblml544banPXnyJAcOHCAiIoKpU6dy+PBhHnjgAQ4fPuy603bnzp1cunSJ7du3oygK9913HwcOHCA4OJi33nqLrVu3EhISQnFxMcHBwUybNq3GXbpz586ttVTziy++yMMPP8zChQv561//6pafhWg6kvhbKJUKRo2q4tNPA0hNNTBggBmDQbr5iLoNHDiQyMhIAPr06cPly5cZNmxYjWl27drFrl27GDt2LOCs9ZORkYHRaGTOnDmEhIQAP5QLuNGtSjUfPnyYDz/8EIAFCxa4BnMRLZMk/hasc2c7XbpYuXxZx8mTXgwZImf9LVVDzszdpXoULnCWarbZbi71rSgKv/rVr3jooYdqfP6Xv/yl3uU3tFSzaLmkjb+FGzzYmexPnPCSUg6i0SZMmMDHH3/sKp+cl5fH9evXGT16NOvXr6eoqAj4oTDYjaWWb1WqeejQoa6KnHXV+BcthyT+Fi4qykZYmA2jUc25c3pPhyNauQkTJjB//nymTJnCyJEjeeihh6ioqCAhIYFnn32WmTNnMnr0aF566SUA5s2bx7vvvsvYsWPJzMyss1TzsmXL+Otf/8rIkSO5cuWKJzdR3AYpy9wKXLigY8sWPwID7TzwQBnqWxyu28o2N4SUZW77ZJtvTcoyt0GxsVYCAuyUlmrIyJC7eYUQjSOJvxVQq2HQIGfviaNHDVLDRwjRKJL4W4mEBDMGg7OGz9WrUrlTCHHnJPG3EjodJCVZAKSMgxCiUSTxtyK9ezubey5c0GMySZ9pIcSdkcTfigQFOejSxYrdrpKunUKIOyaJv5Xp08d51n/6tJdc5BWEhoYyZswYRowYwaJFiygtLfV0SA22fPly3n333QbPt2XLFv74xz8CN5ePnjVrFqmpqbec/8qVKzz44IMNXq87pKam8sL3d3/v27ePw4cPu3V9bk38W7duJS4ujtjYWJbXUsTq8uXLjB8/ngEDBtC3b1/XzSCibt27W/HxcVBUpCEvTyputHfe3t7s2bPHVWitpRRIUxQFh8Ph1nVMmzaNZ555Briz8tERERGu+kKeNmDAAFeO3LdvHykpKW5dn9sSv91u58knn2TLli2cPXuW1atXc/bs2RrTvPbaayxYsIDU1FTWrFnDE0884a5w2gyNBpKSnGf9cpFX3Cg5Odl112xd5ZPz8/O5//77GT16NKNHj3adWb733nuMGDGCESNG8H//938AvPLKKzUOJDeemddW2vny5csMGTKExx9/nBEjRpCbm1tnCei33nqL5ORkpk2bxsWLF2/aFrvdTv/+/VEUhdLSUjQajatG0IwZM8jIyGDVqlU8//zztZaPBli/fj0TJ04kOTmZgwcP3rSOy5cvM2LECABWrVrFz372M+bNm0e/fv14//33ee+99xg7diyTJk1y3UQ1a9Ys/v3f/50JEyYwdOhQjh8/zgMPPMDgwYNd5ahvXC7Au+++60rqs2bNYunSpTfFtW/fPhYuXMjly5f5+9//zp/+9Cf69+/P3r176d+/P9bv67WUlZXVeH+n3HbKmJKSQmxsrOvusYULF7J+/XoSExNd06hUKsrKygAoLS2lc+fOtS7rVneg3Y7Gzt/ShIRARgYUFkJ4OPj43DxNW9vm29Hc25yXl+e6u3LpN0vdso6l4+pebvW6g4ODsdvtHDx4kEceeYTg4GDmz5/Pn/70J3r27Mnhw4d58cUX2blzJ4899hgTJ07kmWeewW63U1FRwcWLF1m7di1Hjx5FURSGDh3KtGnTePDBB3nmmWd47rnnANi4cSNfffUVR44cITc3l2PHjqEoCrNnz+bUqVN06dKFjIwM/vGPfzBs2DC+/vrrWqfz9fVl/fr1nDx5EpvNxsCBAxk+fPhNFUETEhK4cuUKmZmZDBw4kL179zJ06FCuXLnC4MGDOX36NF5eXkydOpU5c+Ywc+ZM5s+fD4BWq0Wr1XLs2DE2b97M22+/zfbt22ssv/qAEhwcjK+vL2lpaaSmpmIymYiNjeUPf/gDJ0+e5Fe/+hUbNmzgmWeeQavVEhAQQGpqKu+88w73338/x44do0OHDvTo0YMXX3yRwMBA13IBfHx8cDgcBAcH1xmXv78/Op2Ofv368fjjj+Pn58dvfvMbwFlm48CBA8ydO5dPP/2U+fPn33THuNForDN/1sZtiT83N7fG8GtRUVE3tVstXbqUyZMn8+6771JZWXnTjhG1CwyE2FhIT4fTp2HIEE9HJDzFaDTSv39/cnNzSUhIYNKkSVRUVHDgwAHuvfde13TV5ZN37tzJRx99BDgreAYGBrJv3z7mzZuHr68vAHfffTd79+7lqaeeIj8/31XILTg4mOjoaN555x2+/vprBgwYADjLNaenp9OlSxe6du3qKgX99ddf1zpdeXk58+bNw+f7M5bZs2fXum2jR49mz549ZGZm8uKLL/L+++8zduxYkpOTb+tnc/fddwMwaNAgsrKy6p1+/Pjx+Pv74+/vT2BgILNmzQKcJa5Pnjzpmq463j59+pCUlERERATgPPHIzs4mKCioSeP6+c9/zhtvvMHcuXNZuXIl77//fr3z1MejjcSrV6/moYce4te//jUHDx7k/vvv5/Tp06h/VIymvdfqqU2HDjqKi/3Yvt1GaGi56/O2vM118WStnuomgKf7Pe2WddRWp+XGGi7e3t7s2rWLqqoq5s+fz5tvvsmiRYsIDAxk165dNy1LURSKi4trlHCuqqrCaDS61mU0GqmqqqK4uJiZM2fy0UcfkZ+fz6xZsyguLsZoNPL000/fVNr58uXLGAyGGsupbbr/+7//q7E+k8nkWt+NBgwYwN/+9jeuXr3Kq6++yptvvsmWLVsYNGgQxcXFVFZWYjabKS4uxmKxUFFR4VqGzWZzfVdeXo7FYrlp+aWlpdjtdteyfvzzro7RaDRSWVlJcXFxjeVWVlaiVqtd8zgcDoqKitBqtTX+b1TP9+P5b4yrvLzcNY/RaKyRAxMTE8nIyGDjxo1YLBYiIyNv2pby8vKb/v83qlaPw+EgNTWVTZs2sXPnTvLz8+ubBYDIyEiys7Nd73NyclyDRFT74IMPWLBgAQDDhw/HZDJRUFBwW8tv72JirOj1CteuaSkqks5Z7Z2Pjw/Lly/nvffew8fHp87yyWPGjOFvf/sb4GxHLysrc1XZrKqqorKykk2bNjF8+HDAWZ3z888/Z8OGDa5RuOoq7fxjdU03YsQINm/ejNFopLy8nK+++qrWbRo4cCApKSmo1WoMBgP9+/fn73//e43282o3lo/2tLCwMAoKCigqKsJsNte5fXWpbVt+8pOfsGTJEu67774mibHOjJGRkcGSJUuIjY3lhRdeYPXq1fzv//4vEydOZNiwYaxcufKWV+2Tk5NJT08nMzMTi8XCmjVrbvqTrkuXLuzYsQOAc+fOYTKZ6NixY5NsWFun00FsrPNO3vPnpU+/gL59+5KUlMRnn312y/LJ+/btY+TIkYwfP57z58/Tr18/Fi1axMSJE5k0aRL3338/ffv2BZzt7BUVFURERBAeHg7UXdr5x+qarl+/fsybN48xY8awYMECV1PQj3l5eREZGcngwYMBZ9NPRUVFjeuE1X5cPtqTdDodzz33HBMnTuTuu++mV69eDZp/6tSpbNq0yXVxF+Dee++lpKSEe+65p0lirLMs86JFi3j88ccZPXr0TSPr5Ofns2rVKoKDg2/ZD3bz5s2ui0gPP/wwv/vd73j55ZcZPHgws2fP5uzZs/zbv/0bFRUVqFQq3njjDSZPnnzTcqSpp3Y5OVo++8yfgAA7Dz1UhkrV9re5NlKWue1r79u8fv16tmzZwp/+9Kdap21oWWapx9+KKQqsXBlIebmae+4pJyrK1ua3uTaS+Nu+9rzNS5YsYfv27axdu5bY2Nhap21o4r+ti7sHDhwgKyurxhieDzzwwO3MKtxIpYK4ODNHj3pz/ryeqKibx1gVQrRuf/jDH5p8mfUm/vvvv5+MjAz69++PRuMsB6xSqSTxtxDx8RaOHvUmPV3P2LFVng5HCNEK1Jv4jx49ytmzZ29q5xctQ0iIg7AwG/n5WrKydMTFeTqi9kVRFPndEB51J6319fYD7N27N1evXr2jgETz6NXL2bvnwgXp3dPcTCbTHf3iCdEUFEXBZDI1eL56z/gLCgpITExkyJAhNW762LBhQ4NXJtyjVy8L+/b5kJWlw2IBveT/ZqHVaikvL6e8vLz+iZuQ0WgEaPb1epJs861ptQ27F7feqZcuXdqgBYrm5++vEBFh48oVLWlp0KePpyNqH1Qq1U09KZpDdU2W9tRzS7a5adXb1DN27Fji4+NdZzYJCQmMHTu2yQMRjVPd3HPmjIcDEUK0ePWe8X/yySc899xzjBs3DkVR+OUvf8mbb77pqoInWobYWAt79niTng530OQnxE2qr10oKLW+rzHtDZ/Vdc2jtvlul8nm/E9ttpvbzTUVo9XZ1ONQHKhVTVuWpd7E//vf/54jR464blS5fv06EydOlMTfwvj5KXTubKOqCtLSwN/f0xEJRVGwOCxY7BbMdjNWh9X1sDlsrmebYsPmsGFX7NgcNhyKA5vifLYrduwOOw7FgQPn+9DSUByKg2v513AoDhQU5/ffv1YUBdc/RcGBAxRu+q46xttJ7J4WnNsOb+DKc27zzLCZdDB0aNJl15v4HQ5HjbsTQ0JC3D6yjrgzvXpZOXHCWar5+xpbogkpioLZbqbCWkG5tZwqWxVV1ioqbZWYbCaMdiNGmxGz3YzZbsZit7gliVbpnPdrFBubJwmqcHZXra3bavV3rve3MU19n9e2PIPWAICXpmkGH7qddXuat9YbqP1n2lj1Jv6pU6cyZcoUFi1aBMDatWuZPn16kwciGi821sK33zoHaRkwQIXB0PLO3Fo6RVGosFZQbC52PcosZZRaSimzlGFzNOzuaJ1ah16jR6/Wo9fo0al1rodWrUWr1qJT6dCoNWhVWjRqDRrVDw+1Su36TIUKjUpDt67d0Kg1ZF/ORoUKtUqNSvX9Mz88q1QqXP9UPzzDD4mvugmhtu9a0v0JUoqkadWb+N98800+++wz9u/fDzjrRsybN6/JAxGN5+Oj0L07XLoEGRk6kpIsng6pRVMUhRJLCVcrr3LNeI0CYwEFpgLMdnOd8+jUOvz1/vhqffHV+eKr9cVb642P1geD1oC31huDxoCXxgsvjVeTt80CdA3qCoDFR/avuDO31fnznnvuabJyoMK9kpKcif/iRb0k/h9RFIVCUyHZFdlkV2RzpfIKJvvNV8INWgMdvDoQ7BVMsFcwgV6BBOoDCdAHNFlTgxCeVGfiHzVqFPv27cPf37/Gn3zVt6hXj5UrWpb4ePjXv+DyZS1mswovr/bd3GNz2MiuyCazLJNLZZeotFbW+N5H60O4TzjhPuGEeofS0bsjvlrfFtXMIURTqzPx79u3D2hfd8q1Bb6+0LUrFBWpyMzUER/f/s76FUUhtzKXc8XnuFhyEYvjh5+Br86XaL9oov2iifSNJEAfIEletDu3VZ3zH//4R72fiZYjMRFSUyE9vX0lfpPNxOmi05wqPEWZ5Ye/SDt6dyQmIIaYgBg6eneURC/avXoT/5kf3Qpqs9k4duyY2wISjRcf73z+7rv2Ubun2FjM/uz97Dq3y9Xrxl/nT1xwHAnBCU3eB1qI1q7OxL9s2TJef/11jEYjAQEBgPNPaL1ez5IlS5otQNFwAQG4avdkZeno1cvq6ZDcosxSRsq1FK5mX3Xe9OSw0dW/K/1C+9HNv5uc2QtRhzoT/4svvuh6LFu2rDljEk0gNtbClStaLl7Ut7nEb7FbOJJ/hOPXj+NQHHQI7kC/Tv2IDouWs3shbkO9TT3Lli2juLiY9PT0GnWfx4wZ49bAROPExlrZuxeysnRYreCBIpJNTlEULpRcYO+Vva7eOXFBcfxkyE8I8QlpVzf3CNEY9Sb+v/71r7zzzjvk5OTQv39/Dh06xPDhw9m5c2dzxCfuUEDADyNzXb6so0eP1n3WX2mtZEfODjLLMgHo5NOJcZHjCPcJJ8QnxMPRCdG61Htb4TvvvMORI0fo2rUru3btIjU1laCgoGYITTRWz57OHj0XL7bu0/2M0gz+eeGfZJZl4qXx4q6ou/hJ7E8I9wn3dGhCtEr1nvEbDAYMBmeBJLPZTHx8PGlpaW4PTDRejx5W9u+HzEwddjtoNJ6OqGHsDjt7r+zl24JvAYj2i2ZS9CT89VJ6VIjGqDfxR0VFUVJSwty5c5k0aRLBwcF07dq1OWITjRQc7CAkxE5hoYbsbC3dujWswJgnGW1GNn+3mZyKHNQqNaMiRtE/tL/01BGiCdSb+L/44gvAOQTj+PHjKS0tZerUqW4PTDSN2FgLhYXeZGToW03iLzAWsDFrI2WWMny0PszoNoPOvp09HZYQbUadib+oqOimz/p8P5hrRUUFHTpIt7nWIDbWyuHD3mRk6Bg/HtRNXyyySV2tusr6S+sx2U2EeYcxs9tMadoRoonVmfgHDRqESqWqdZgzlUolXedaiZAQO4GBdkpLNeTlaYmKarln/TkVOWzI3IDVYaV7QHemdZ2GTt26L0wL0RLVmfgzMzObMw7hJioV9Oxp5ehRDRcv6lps4s8qy+JfWf/CrtjpFdSLydGT0ahb2dVoIVqJetv49+zZU+vncgNX69Gjh4WjRw1kZOgZO9ZIS7s+mleZx6bvNmFX7PTu0JvxUePdMoCJEMLptkbgqmYymUhJSWHQoEFyA1cr0qmTHX9/B+Xlaq5e1RARYfd0SC4FxgI2Zm7E5rCR2CGRCVETpOeOEG5Wb+LfuHFjjffZ2dk888wz7opHuIFK5TzrP3HCwMWLeiIijJ4OCXAWWVuf6byQ2z2gO3dF3SVJX4hm0OC/p6Oiojh37pw7YhFudONdvLVcr292FruFjZkbqbBWEOkbybSu06R5R4hmUu8Z/y9/+UvXWZjD4eDEiRMMHDjQ7YGJphURYcfX10FZmYb8fA2dOnmuuUdRFHbk7KDAVECwVzCzus+S3jtCNKN6E//gwYN/mFirZdGiRYwcOdKtQYmmV93cc/Kks7mnUyfPNfccv36cCyUX0Kv1zOw2UwYwF6KZ1fu39YMPPsiiRYsYMGAAffv2JTk5+bYXvnXrVuLi4oiNjWX58uW1TvPJJ5+QmJhIUlIS99133+1HLhqsZ09nhc70dM8191wuv8z+K/sBmNxlstTPF8ID6j3j37x5M48++ig9evRAURQyMzP585//zLRp0245n91u58knn2Tbtm1ERUWRnJzM7NmzSUxMdE2Tnp7OsmXL2L9/P8HBweTn5zd+i0SdOne24e3toLRUw/XrGsLCmre5p8paxVeXv0JBITksmR6BPZp1/UIIp3oT/7PPPsuuXbuIjY0FICMjgxkzZtSb+FNSUoiNjSUmJgaAhQsXsn79+hqJ//333+fJJ58kODgYgLCwsDveEFE/tdpZwuHUKS/S03XNmvir2/WrbFVE+kYyLHxYs61bCFFTvYnf39/flfQBYmJi8Pevv3ZKbm4u0dHRrvdRUVEcPny4xjQXLlwAYOTIkdjtdpYuXVprAbjqg8edauz8rVFd26xSQU4OlJQE0707zXYz14mrJyjWFBMeGs6S5CUEGYKafB3tbT+3t+0F2eamclsXd6dPn86CBQtQqVR8+umnJCcn8/nnnwNw99133/HKbTYb6enpfPPNN+Tk5DBmzBhOnTolA724Udeu4OMDhYWQnw+dOrl/ncXGYrakbwFges/pbkn6QojbV2/iN5lMdOrUid27dwPQsWNHjEYjGzduRKVS1Zn4IyMjyc7Odr3PyckhMjKyxjRRUVEMHToUnU5H9+7d6dWrF+np6TddQL7TgnDVR8r2VFDudrY5ONiH3Fwvtm0zMmKEqc7pmoKiKHx26TOuVlwlNjAWv0q/Jt8f7W0/t7ftBdnmxsxfm3oT/8qVK+9opcnJyaSnp5OZmUlkZCRr1qxh1apVNaaZO3cuq1evZvHixRQUFHDhwoV2+adcc+vZ08Lp015cuKBn+HCTW5t7zhefJ7ciF2+tt5RjEKKFqLc7Z05ODvPmzSMsLIywsDDuuececnJy6l2wVqtlxYoVTJkyhYSEBBYsWEBSUhIvv/wyGzZsAGDKlCmEhISQmJjI+PHjefPNNwkJkYGz3S0qyoaPj7N3T36++ypgmmwm9l7ZC8CoiFF4a73dti4hxO2r94x/8eLF3HfffXz66acAfPzxxyxevJht27bVu/Dp06czffr0Gp+9+uqrrtcqlYq3336bt99+u6Fxi0ZQq51n/d9+ayAtzX03cx24egCjzUikbyQJwQluWYcQouHqPeO/fv06ixcvRqvVotVqeeihh7h+/XpzxCbcqFcvZ+2e9HS9W27mulp1ldOFp1Gr1IyLHCdNPEK0IPUm/pCQED7++GPsdjt2u52PP/5YmmPagIgIOwEBdioq1OTl1fuHX4MoisKunF0oKAzsOJBQ79AmXb4QonHqTfx/+9vf+OSTTwgPDyc8PJx169bd8QVf0XKoVNCrl7OEw4ULTVsg7ULJBfKN+fjp/EgOu/0SH0KI5lHvqV7Xrl1dF2NF29Krl3NkrgsX9IwZY0TTBNd5bQ4bB64eAGBY+DD0Gn3jFyqEaFL1nvFfunSJWbNm0bFjR8LCwpgzZ0676kvbloWG2gkOtmMyqcnObprmnpOFJymzlBFiCJELukK0UPUm/vvuu48FCxZw5coV8vLyuPfee1m0aFFzxCbczNnc47zIm5bW+DNzs93MkWtHABgZMVIGVhGihar3N7Oqqor777/f1avnZz/7GSaTe+/2FM0nPr56ZC49ZnPjlnU0/ygmu4lIv0i6+XdrfHBCCLeoN/FPmzaN5cuXk5WVxXfffccbb7zB9OnTKSoqoqioqDliFG4UFOQgMtKKzabi4sU7P+uvtFZyouAE4LxZS7pvCtFy1duw+8knnwDw5z//ucbna9asQaVSSXt/GxAfbyE3V8e5c3qSkix3tIzU66nYHDZiAmII9wlv4giFEE2p3sSfmZnZHHEID+rZ08Lu3T7k5uooLVUTGOho0PxGm5GThScBGNJpiDtCFEI0Ibn6JvDygthY55n+uXMNb+45fv04VoeVbv7d6OTTDHWehRCNIolfAD9c5D13rmElHEw2EycLnGf7yZ3kZi0hWoM6E//+/c4Bsc2N7eohWoXoaBt+fg7KyjQNKuFwouAEFoeFaL9oOvt2dmOEQoimUmfif+qppwAYPnx4swUjPEethvh450H+zJnba+6x2C18W/AtIG37QrQmdZ7a6XQ6lixZQm5urusgcKP/+Z//cWtgovklJlo4etSb9HRnCQeD4dZtPmeLzmKym4jwiSDKL6qZohRCNFadif9f//oX27dv56uvvmLQoEHNGZPwkOBgB9HRVrKznV07Bwyou5nPoThc/fYHhg1spgiFEE2hzsQfGhrKwoULSUhIoF+/fs0Zk/Cgvn3NZGfrOHXKi/79zXUOy5hZlkmppZQAfQAxATJcphCtyW3V47+ToRdF6xQTY8XX10FxsYacnLov8qZeTwWgf2h/qckjRCtT72/s4sWLmT17Nnl5eeTl5TFr1iwWL17cHLEJD1CroXdvZxPPyZNetU5zreoauZW56NV6EjskNmd4QogmUG/iz8/Pl6EX25mkJDNqtcKlSzoqKm5u66lu2+8d0hsvTe0HByFEy1Vv4g8NDZWhF9sZf3+F7t2tOBwqzpypmdgrrBVcKLmAChX9QuXajxCtUYOGXoyIiJChF9uJvn1/aO6x2X74/EzRGRyKgx6BPQjQB3goOiFEY8jQi6JW0dE2QkNtFBRoOX9eT+/eFhyKg9OFpwHoE9LHwxEKIe6UdMcQtVKpYNAg51n/8eMGFMXZhbPCWkGQVxDRftEejlAIcack8Ys69expwd/f2bUzM1PHqcJTgPNsXwZaEaL1ksQv6qTRwIABzmE29x418l35d2hUGhKDpQunEK3ZbSf+Q4cOMXXqVMaNG8eXX37pxpBES5KUZMbLy8HJ66cpL1fRK6gXBq3B02EJIRqhzou7V69eJTz8hyH03n77bb744gsURWHo0KHMnTu3OeITHqbXQ2LvSg6cOYUxT0ufAXJRV4jWrs7E/9hjjzFw4ECef/55DAYDQUFBrFu3DrVaTUCAdONrTwJizmI/X4WlqDOURYJvw4ZmFEK0LHU29Xz55ZcMGDCAmTNn8tFHH/HHP/4Rs9lMYWGhNPW0M5mmM4R3stHR1peDB308HY4QopFu2cY/a9YsvvrqK0pLS5k3bx69evXiqaeeomPHjs0Vn/Cwcks52eXZREdCZ2082dk6srNvf4QuIUTLU2fi37BhA+PHj2fq1Kn07t2btWvXsn79ehYuXEhGRkZzxig86GzxWRQU4kJiGDrQ+d/lwAHvBo3LK4RoWeo8dXvppZdISUnBaDQyZcoUUlJSeOutt0hPT+d3v/sda9asac44hQcoisLZorMAJHZIJCLCxIkTXly9qiUrS0f37lYPRyiEuBN1Jv7AwEA+//xzqqqqCAsLc33es2dPSfrtRE5FDmWWMvx1/kT7RaNWQXKyiT17fNi/35suXaxoNJ6OUgjRUHU29XzxxRcUFhZis9lYtWpVc8YkWoizxc6z/YQOCa7BVvr0MRMQYKewUFNnvX4hRMtWZ+I3GAz88pe/5LHHHquz+2ZFRcUtF75161bi4uKIjY1l+fLldU732WefoVKpOHr06G2GLdzNbDeTXpIOUONOXa0Wxo0zAnDokHet9fqFEC1bnYl/zpw5/PrXv2bPnj1UVla6Pr906RIffPABU6ZMYevWrXUu2G638+STT7JlyxbOnj3L6tWrOXv27E3TlZeX88477zB06NBGbopoShdKLmBX7ET5RRHoFVjju+7drXTvbsFiUbFvn3TvFKK1qbONf8eOHWzevJk///nP7N+/n+LiYrRaLXFxccyYMYMPP/ywxp29P5aSkkJsbCwxMc6BuBcuXMj69etJTKxZ5+U//uM/+O1vf8ubb75Z57Kql3GnGjt/a9TYbf6m9BuCg4OZFD+JmPCbl/XQQ7BiBeTnOyt5du/eqNU1ifa2n9vb9oJsc1O5ZYfs6dOnM3369DtacG5uLtHRP5TujYqK4vDhwzWmOX78ONnZ2cyYMeOWiV80ryJjEZdLL6NT60gITah1mqAgGDMGdu6ETZvg0UdBp2veOIUQd8Zjd+I4HA6effZZ/v73v9c77aVLl+5oHdVHyjudvzVqim0+dPUQxcXFxAfHk3s5t87pnH/wBZCeruHjj02MHWu843U2Rnvbz+1te0G2uTHz18ZtZZkjIyPJzs52vc/JySEyMtL1vry8nNOnTzNu3Di6devGoUOHmD17tlzg9TBFUThffB6A+OD4W06r1cKUKZWo1QonThj47ju5o1eI1sBtiT85OZn09HQyMzOxWCysWbOG2bNnu74PDAykoKCArKwssrKyGDZsGBs2bGDw4MHuCknchitVVyi1lOKr872tUbY6dbIzbJizZv+2bb6YTNLLR4iW7paJ3263Ex9/67O+umi1WlasWMGUKVNISEhgwYIFJCUl8fLLL8sYvi1Y9dl+XFCcq+9+fQYNMtG5s43KSjU7dvhIOQchWrhb/m2u0WiIi4vj8uXLdOnSpcELr+3i8KuvvlrrtN98802Dly+als1h40LJBQASgmu/qFsbtRomT67kn/8M4OJFPceO2Rg82OyuMIUQjVRvo2xxcTFJSUkMGTIEX19f1+dy1t72ZJVnYbab6ejdkVDv0AbNGxjoYPLkSjZt8uPAAW9CQhxSy0eIFqrexP9f//VfzRGHaAFubOa5E7GxVoYNM3LokDdbt/qyYEEZISEyaIsQLU29jbhjx46lW7duWK1Wxo4dS3JyMgMHDmyO2EQzMtlMZJZlokJFXPCdJX6AIUNM9OzpvKt340Y/KivlYq8QLU29if/9999n/vz5PProo4DzxiwZb7ftuVh6EYfiIMovCj+d3x0vR6WCiRMr6djRRmmphi++8MdolOQvREtSb+J/77332L9/v6tQW8+ePcnPz3d7YKJ5nS/5vpmnEWf71fR6mDu3gg4dnFU8v/jCD7NZkr8QLUW9id/Lywu9Xu96b7PZUKnkl7gtKbeUk1uRi0alITYwtkmW6eOjMG9eOYGBdq5f1/Lll37Sx1+IFuK22vhff/11jEYj27Zt495772XWrFnNEZtoJmklaQDEBMbgpWm6Gvt+fgr33FNOQICdq1e1fPKJP2VlbrtnUAhxm+r9LVy+fDkdO3akT58+/PnPf2b69Om89tprzRGbaCZpxc7Ef6e9eW7F319h/vxyQkNtFBdrWLvWn2vXZNguITyp3u6carWaBx98kKFDh6JSqYiLi5OmnjakwFhAgakAg8ZAN/9ublmHM/lXsGmTL9nZOtat82fChCoSEixuWZ8Q4tbqPePftGkTPXr04KmnnuIXv/gFsbGxbNmypTliE82gupmnZ1BPNGr3nYl7eSnMmVNBQoIZm03F11/78vXXPlgk9wvR7Oo94//1r3/Nrl27iI11XvTLyMhgxowZTJs2ze3BCfdSFMWtzTw/ptHApElVdO5sY/duH86d8+LqVS0TJ1bSubPd7esXQjjVe8bv7+/vSvrgrPHs7+/v1qBE88irzKPcWo6/zp/Ovp2bZZ0qFfTubeEnPymjQwc7xcUa1q3zZ+dOH+nyKUQzqfOM//PPPwdg8ODBTJ8+nQULFqBSqfj0009JTk5utgCF+7hKNAQ3/3Wb0FAHCxeWkZJi4PhxA6dOeZGRoWPoUCNJSRY0cv1XCLepM/Fv3LjR9bpTp07s3r0bgI4dO2I0emakJdF07A47F0svAs3TzFMbnQ5GjjQRF2dh505frlzRsmuXL8ePGxg+3EjPnlbU0vtTiCZXZ+JfuXJlc8YhmllWeRYmu4lQQ2iDK3E2tdBQB/feW056uo6DB70pKdGwdasfBw7Y6d/fTGKiGa+mu71AiHav3ou7mZmZvPvuu2RlZWGz2VyfS1nm1q26N09TlGhoCioV9OplJTbWytmzeo4cMVBWpmHPHh8OHfImLs5CfLyZiAg70ptYiMapN/HPnTuXRx55hFmzZqGWv7vbBLPdzKVS5wDOnmrmqYta7bz4m5ho4dIlHSdOeJGbq+PUKS9OnfIiKMhOz54WevSwEhYmBwEh7kS9id9gMPDUU081RyyimWSUZmBX7ET6RuKvb5k9tNRqZ33/2FgrBQUazp3Tc/68npISDUeOeHPkiDf+/g66drVSVQXdu3s6YiFaj3oT/9NPP80rr7zC5MmT8bqhoVVq8rde54rPAS2nmac+oaF2Ro82MnKkkexsLZcu6cjI0FNerub0aS9yc6unDCAiwkZ4uI1Onex06GBHW+//cCHan3p/LU6dOsU//vEPdu7c6WrqUalU7Ny50+3BiaZ3YyXOXkG9PB1Og6jV0LWrja5dbYwbZyQ/X0N2thaTKZjLl+H6dQ3FxRrOnvX6fnqF4GAHISF2goOdj6AgB4GBDry8FGkmEu1WvYn/008/5dKlSzVKM4vWK60kDQWF7gHdm7QSZ3NTqaBTJzudOtmJiQGbDVJSyrh6VcvVq1ry8zWUlKgpLNRQWHjzTQF6vUJAgB0/PwU/Pwe+vg58fBz4+ip4ezvw9lbw9lbkACHapHoTf+/evSkpKSEsLKw54hFupCiK66athOAED0fTtLRaCA+3Ex5uB8wAWK1QWOj8K6CoSE1xsYbSUjWlpRosFhUFBVoKCm69XJXKmfz1+h+eq1/rdKDVKuh0zodW63yv1YJGo6DVKmg0fP9wvlarFdRq52cqVfVnztdqNa6HHGyEO9Wb+EtKSoiPjyc5OblGG79052x9rhuvU2gqxKA10NW/q6fDcTud7saDwQ8UBUwmFeXlaioqqh8qqqrUVFU5n00mFSaTCrO5+nXzx69WK98/Ow8EarXzr4+QEOf7oqJAVCpcD3AeQKrV/Lyu55rT38qPv7/dg1NDD2K1Td+xo/P5+vWW2RnBHaq3edAgNQEBjiZddr2J/5VXXmnSFQrPqR5esVdQL7dW4mzpVCq+b8qxExZ26+JwDgeYzSrXw2L54WG1qrBawWpVYbM5X9tsKux2FTYbrmeHo/rZ+ZndDoqi+v6980DkcDg/d65T9aNnV+QAVFXx/XP76V5dfeAtLm4/V+urt7lfv6Zfdr0/xbFjxzb9WkWzcygOVyXO+KB4D0fTeqjVuNr7m5PzYFDzGZwHi27dglEUuHSp5PvPq+dRoSi4Pvvx6xufqw8iNT+r+bq+7+qa7nY0dPquXYMB+O67sobN2Ip16eLcZqOxac/24TYSv7+/v6uAl8ViwWq14uvrS1lZ+9kBbcHl8stU2aoI8goi3Cfc0+GIeqhU1FKozpktfX2d7/z8fpw9m/fg1Jyio53PVmv7Kd/dpYvz+dKlpl92vYm/vLzc9VpRFNavX8+hQ4eaPhLhVtUXdeOD42UENSHauQY1EqpUKubOnctXX33lrniEG5hsJi6WXkSFqs315hFCNFy9Z/zVdfkBHA4HR48exWAwuDUo0bQulFzArtiJ9osmQB/g6XCEEB5Wb+K/sS6/VqulW7durF+/3q1BiaZ1pugMAEkdkjwciRCiJag38Utd/tatwFhAvjEfL40XMYExng5HCNEC1Jn4X3311TpnUqlU/Md//IdbAhJN62zxWcBZflmn1nk4GiFES1Bn4vet7jN2g8rKSj744AMKCwsl8bcCdof9hxINHeSirhDCqc7E/+tf/9r1ury8nHfeeYeVK1eycOHCGt+JliuzPBOjzUiIIYRO3p08HY4QooW4ZXfOoqIiXnrpJfr27YvNZuP48eP84Q9/uO2CbVu3biUuLo7Y2FiWL19+0/dvv/02iYmJ9O3bl7vuuovvvvvuzrZC1OpMofOibmKHROm7L4RwqTPxP/fccyQnJ+Pv78+pU6dYunQpwcHBt71gu93Ok08+yZYtWzh79iyrV6/m7NmzNaYZMGAAR48e5eTJk8yfP5/nn3/+zrdE1FBqLuW78u/QqDTSd18IUUOdif+tt94iLy+P1157jc6dOxMQEEBAQAD+/v4EBNTfFzwlJYXY2FhiYmLQ6/UsXLjwpm6g48ePx8fHB4Bhw4aRk5PTyM0R1U4XnUZBoWdQT7y13p4ORwjRgtTZxu9wNK4wUG5uLtHVBTaAqKgoDh8+XOf0H3zwAdOmTav1u5iYxnVDbOz8rY3NYeMKVwgODmbmgJlEB0bXP1Mb0N72c3vbXpBtbiotosbpxx9/zNGjR9m9e7enQ2kTzl0/R5W1inC/cKICojwdjhCihXFb4o+MjCQ7O9v1Picnh8jIyJum2759O7///e/ZvXt3jYFebnTpDsvTVR8p73T+1igmJoYjeUcoLi6mv29/MjMzPR2S27W3/dzethdkmxszf23cNpJDcnIy6enpZGZmYrFYWLNmDbNnz64xTWpqKo8++igbNmyQoR2byLWKa1wuvYxerZe6+0KIWrkt8Wu1WlasWMGUKVNISEhgwYIFJCUl8fLLL7uGbXzuueeoqKjg3nvvpX///jcdGETDHc07CjjLL+s1eg9HI4Roidzaxj99+nSmT59e47MbS0Fs377dnatvd4w2IyeungCgT0gfzwYjhGix2s+gne3AqcJTWB1WYjvEEuod6ulwhBAtlCT+NsLmsHGi4AQAI6NHejYYIUSLJom/jThffB6jzUiEXwTdgrp5OhwhRAsmib8NUBSF49ePAzAieoTU5RFC3JIk/jYgsyyTYnMx/jp/EjsmejocIUQLJ4m/lVMUhWPXjwHQv2N/NGqNhyMSQrR0kvhbue/KvyOvMg+D1iBj6gohbosk/lZMURQOXj0IwKCOg/DS1F7yQgghbiSJvxXLLMsk35iPj9aHviF9PR2OEKKVkMTfSt14tp8clizlGYQQt00SfyuVXppOgakAP50fvUN6ezocIUQr0iLq8YuGsTvsP5ztd0pGq27gbrTb0RQUoK6qQmU0ojKbcfj44PD3x+Hvj+LnB3IvgBBtliT+Vii1IJUScwlBXkEkBd9eTx6V2YzuwgV0WVnosrNRWa11TmsPCMDarRvW7t2xRUeDRrqICtGWSOJvZcot5aRcSwFgXOS4evvtqyorMaSm4nXqFCqLxfW5PTgYR2AgDoMBRa9HXVWFurwcdWkpmrIyNCdPYjh5EoePD+a+fTH36YPy/fjIQojWTRJ/K7P3yl6sDis9AnvQ1b9r3RNarRhSUjCkpqKy250fRUVh6dULa7duKP7+tc+nKGiuXUOXmYn+4kU0RUV4HzqE4cgRLImJGIcMcTYFCSFaLUn8rcjl8sukl6SjVWsZ23ls3RNeukTAP/+JprQUAEtMDKbkZOzh4fWvRKXCHh6OPTwc07BhaHNy8EpNRZeVhdepU+jPncPUvz/mwYNR6hgqUwjRsknibyWsDivf5H4DwJCwIfjrazljt1ph40Y4dgxNaSn2kBAqJ0zA3rnzna1UpcIWHY0tOhp1YSHeBw+iz8jA++hRvM6cwThyJJbERLkQLEQrI4m/ldh3ZR/F5mKCvYIZ0HHATd+rCwvx27IFHA7QaDAOH45p0KAmuzDrCAmhcuZMTFeu4L1/P7rcXHy3b8fr1Cmqxo/H3qlTk6xHCOF+kvhbgcyyTE4WnEStUjO1y9Sbum/qz53DZ+dOVDYbxMbCvfdiqqpySyz2iAgq7rkH3YUL+Ozdi/baNfzXrsXcty/G4cNBmn+EaPEk8bdwVdYqtmVvA2B4+HDCfMJ++NJux3vPHgwnTwJgjo+HRx4BvR4uXXJfUCoV1rg4Srt1wzslBa8TJzB8+y36ixepGjsWa2ysNP8I0YJJ4m/BFEVhe852jDYjkX6RDOw40PWdqrISv82b0ebloajVVI0bh6VPH2fSby5eXhhHj8YSH4/Pzp1or17Fb/NmrN26UTVuHI7AwOaLRQhx26RkQwt24OoBMssy8dJ4MSV6CmqVc3dp8vIIWL0abV4eDj8/yufPdyZ9D7F37Ej5vfdSNX48Di8vdFlZBHz8MYYjR+D7rqRCiJZDzvhbqNOFpzmafxS1Ss30rtNdvXj0p07h8803qBwOrJGRVE6bhuLr6+FoAbUac9++WHr0wGfvXvRpaXgfOID+7Fmqxo3D1vUW9xwIIZqVJP4W6HL5ZXbl7gJgfOR4uvh3AasVn9278TpzBgBTv34YR49uceUUFF9fKqdOxZyQgM/u3WiKi/H/8kssPXpgHD1amn+EaAEk8bcweZV5bMrahENxMKjjIHqH9EZdXIzv5s1oCwpQNBqq7roLS0KCp0O9JVvXrpT99Kd4pabinZKCPiMDXWYm5v79MQ0ZIjd/CeFBkvhbkOzybDZmbcTqsNIrqBcjI0aiS0vDd8cOVFYr9sBAKmfMwN6xo6dDvT0aDebBg7HExeF94ABe589jOH4c/dmzmAYNwtyvH+h0no5SiHZHEn8LkVmWyaasTdgVO/HB8UzuOBq/rVvRX7gAgKVnTyrvuqtV9pNX/P2pmjIFc//+eO/diy43F5/9+zGkpmJKTsaclCQHACGakSR+D1MUhdSCVPZf2Y9DcdA3pC8TLV3x++cq1JWVKFotVaNHO3vttPK+8fZOnai45x60332H98GDaPPz8dm9G0NKCub+/TH37YtiMHg6TCHaPEn8HmS2m9mWvY2M0gwAkv2TmHiiEq+MjQDYOnemctIkHEFBHoyyialU2Lp1o7xrV3SXLmFISUGbn4/3wYMYjh7FnJCAuU8fHKGhno5UiDZLEr+HZJdnsyNnB6WWUvSKhlkVUfTel4bKZkPR6TAOHYp5wABQt9FbLVQqrD16YI2JQZudjeHYMXSXL2P4fhwAa2QklqQkLD16NO9NaUK0A5L4m1mltZK9eXtJK0kDu52IIgvzL3kTYvwOcLblV40eXXe9/LZGpcLWpQsVXbqgLijA69QpvM6dQ5eb67wWoNNh6dEDS69eztHAtPJfVojGkt+iZmK0GUm9nsq3Bd9iMVXgdb2Q0Xk6hlV2QAvYOnXCOGIEti5dPB2qxzhCQzGOH49xxAj0aWl4nT+P9soVvM6fx+v8eRS9Hmv37li6d8fWpQuKt7enQxaiVZLE72ZFpiJOFZ7idMEpHMUFaAoKSCpQmGTpQpBiwNaxIxXDh2Pt1q3VX7xtMl5eWPr2xdK3L+riYvQXLqC7eBFtQQH6tDT0aWkoKhX2sDCs0dHOUtTR0Z6OWohWQxK/G1RaK7lUdolz+ae4djUNTUkJmuIi4s3+DLN1JlIJwBoTQ3m/ftiioiTh34IjOBjT0KGYhg5FXVKCLiMD3Xffoc3NRXvtGtpr1yAjA1Qq/DUa7GFh2Dp1wh4Whr1DB+kmKkQtJPE3AavDyrWqa+QWZnA59xT5hZmoKypQV5Tj51ATb+/AQFsCHYOjscTFUZqQ0H7a8JuQIygI86BBmAcNAosFXW4u2txcsFggLw/t9etor193lbVQVCocQUHYO3RwDi4fHIw9KAhHYKBz4Hg54Ip2yq2Jf+vWrTz99NPY7XZ+/vOf88ILL9T43mw288ADD3Ds2DFCQkJYu3Yt3bp1c2dIjaIoChWmUsqKcykqzKawJIfCsjwKKq6gVFWisloBMChqujoCiHfE0iM0HlW3Hlh69KBMuig2ne/b+63du0NMDFgslKWkoL16Fe21a2gKClCXlKApLkZTXHzT7IpWi8PfH4efn/PZ1xfF1xeHj4/ztcGA4u3tvK9ADhCijXFb4rfb7Tz55JNs27aNqKgokpOTmT17NomJia5pPvjgA4KDg7l48SJr1qzht7/9LWvXrnVXSKAoYLOhcjiczzYbDosZq8WIxVSB2VKJ2VyJyVSOyVROlbmcSnMZVaYyyizllJlLsNstNy1WragIU3zorAqjS3APIjsnouocjS0yEnMrvNO2VdLrsXfujL1zZ8zVn9lsaIqKUBcXoykqQlNSgrqkBHVZGWqTqc6Dwo0UlQpFr0fx8vrhode7Huj1KFotik6HotVC9bNWi6LROF9rNKDRoFQ/q9XOz9Rq5+vqhxxgRDNRKYqiuGPBBw8eZOnSpXz11VcALFu2DIAXX3zRNc2UKVNYunQpw4cPx2azER4ezvXr11E1wS9AQeZptq15HZvDht1hx+awYVPsWLFjw4EFO1Yc2HA0aLm+Kj3BhmDC/MPpFBRFp9CudI5KQB8RBUFB8svbWphMUFoKZWU/PCoroaLC+TAaoarK+dxcVCrnAUGlqnkw+PFz9ePGg8WNn//4s+rXt3q+k9e1xX8723in2uvv1syZEBDQpIt02xl/bm4u0Tf0tIiKiuLw4cN1TqPVagkMDKSwsJDQJmgSCe3em0Uvrmr0ckQbZTA4HzJIvGiH2uhtoUIIIeritsQfGRlJdna2631OTg6RkZF1TmOz2SgtLSUkJMRdIQkhhMCNiT85OZn09HQyMzOxWCysWbOG2bNn15hm9uzZfPjhhwCsW7eOCRMmNEn7vhBCiLq5LfFrtVpWrFjBlClTSEhIYMGCBSQlJfHyyy+zYcMGAB555BEKCwuJjY3l7bffZvny5U22/q1btxIXF0dsbGyTLrclyc7OZvz48SQmJpKUlMQ777wDQFFREZMmTaJnz55MmjSJ4np6rrRGdrudAQMGMHPmTAAyMzMZOnQosbGx/OQnP8Fiubn3VWtWUlLC/PnziY+PJyEhgYMHD7b5/fzf//3fJCUl0bt3bxYtWoTJZGpz+/nhhx8mLCyM3r17uz6ra78qisJTTz1FbGwsffv25fjx43e+YqUNstlsSkxMjJKRkaGYzWalb9++ypkzZzwdVpPLy8tTjh07piiKopSVlSk9e/ZUzpw5ozz33HPKsmXLFEVRlGXLlinPP/+8J8N0i7feektZtGiRMmPGDEVRFOXee+9VVq9erSiKojz66KPK//7v/3oyvCb3wAMPKO+//76iKIpiNpuV4uLiNr2fc3JylG7duilVVVWKojj378qVK9vcft69e7dy7NgxJSkpyfVZXft106ZNytSpUxWHw6EcPHhQGTJkyB2vt00m/gMHDiiTJ092vX/99deV119/3YMRNY/Zs2crX3/9tdKrVy8lLy9PURTnwaFXr14ejqxpZWdnKxMmTFB27NihzJgxQ3E4HEpISIhitVoVRbl5/7d2JSUlSrdu3RSHw1Hj87a8n3NycpSoqCilsLBQsVqtyowZM5StW7e2yf2cmZlZI/HXtV+XLFmirFq1qtbpGqpN9uqprStpbm6uByNyv6ysLFJTUxk6dCjXrl0jIiICgPDwcK5du+bh6JrWM888wxtvvIH6+7EKCgsLCQoKQvt9yea2tr8zMzPp2LEjixcvZsCAAfz85z+nsrKyTe/nyMhIfvOb39ClSxciIiIIDAxk0KBBbXo/V6trvzZlXmuTib+9qaio4J577uGPf/wjAT+60UOlUrWpC+b/+te/CAsLY9CgQZ4OpdnYbDaOHz/O448/TmpqKr6+vjddt2pr+7m4uJj169eTmZlJXl4elZWVbN261dNhNTt37dc2mfhvpytpW2G1Wrnnnnv46U9/yt133w1Ap06duHLlCgBXrlwhLCzMkyE2qf3797Nhwwa6devGwoUL2blzJ08//TQlJSXYbDag7e3vqKgooqKiGDp0KADz58/n+PHjbXo/b9++ne7du9OxY0d0Oh133303+/fvb9P7uVpd+7Up81qbTPy305W0LVAUhUceeYSEhASeffZZ1+c3dpP98MMPmTNnjqdCbHLLli0jJyeHrKws1qxZw4QJE/jnP//J+PHjWbduHdD2tjk8PJzo6GjS0tIA2LFjB4mJiW16P3fp0oVDhw5RVVWFoiiubW7L+7laXft19uzZfPTRRyiKwqFDhwgMDHQ1CTXYnV6QaOk2bdqk9OzZU4mJiVFee+01T4fjFnv37lUApU+fPkq/fv2Ufv36KZs2bVIKCgqUCRMmKLGxscpdd92lFBYWejpUt9i1a5erV09GRoaSnJys9OjRQ5k/f75iMpk8HF3TSk1NVQYNGqT06dNHmTNnjlJUVNTm9/PLL7+sxMXFKUlJScrPfvYzxWQytbn9vHDhQiU8PFzRarVKZGSk8te//rXO/epwOJQnnnhCiYmJUXr37q0cOXLkjtfrtiJtQgghWqY22dQjhBCibpL4hRCinZHEL4QQ7YwkfiGEaGck8QshRDsjiV+0OYWFhfTv35/+/fsTHh5OZGQk/fv3x8/PjyeeeMIt67xy5QqTJ0++6fOHHnrI1e9ciJbCbUMvCuEpISEhnDhxAoClS5fi5+fHb37zG7euc+vWrUyZMsWt6xCiqcgZv2g3vvnmG1f9/qVLl/Lggw8yevRounbtyueff87zzz9Pnz59mDp1KlarFYBjx44xduxYBg0axJQpU1y30v/Y1q1bmTZtGoqi8Itf/IK4uDgmTpxIfn6+a5pXX32V5ORkevfuzZIlS1AUhYyMDAYOHOiaJj093fX+hRdeIDExkb59+7r9wCXaF0n8ot3KyMhg586dbNiwgZ/97GeMHz+eU6dO4e3tzaZNm7Barfzyl79k3bp1HDt2jIcffpjf/e53Ny3HbreTlpZGYmIiX3zxBWlpaZw9e5aPPvqIAwcOuKb7xS9+wZEjRzh9+jRGo5F//etf9OjRg8DAQNdfKCtXrmTx4sUUFhbyxRdfcObMGU6ePMlLL73UXD8W0Q5I4hft1rRp09DpdPTp0we73c7UqVMB6NOnD1lZWaSlpXH69GkmTZpE//79ee2118jJyblpOYcPH3YVUNuzZw+LFi1Co9HQuXNnJkyY4Jpu165dDB06lD59+rBz507OnDkDwM9//nNWrlyJ3W5n7dq13HfffQQGBmIwGHjkkUf4/PPP8fHxaYafiGgvpI1ftFteXl4AqNVqdDqdq/ytWq3GZrOhKApJSUkcPHjwlsvZsmWL66BRF5PJxBNPPMHRo0eJjo5m6dKlmEwmAO655x5eeeUVJkyYwKBBgwgJCQEgJSWFHTt2sG7dOlasWMHOnTsbu8lCAHLGL0Sd4uLiuH79uivxW61W11n6jXbs2MHEiRMBGDNmDGvXrsVut3PlyhV27doF4EryoaGhVFRU1OjpYzAYmDJlCo8//jiLFy8GnGMslJaWMn36dP77v/+bb7/91q3bKtoXOeMXog56vZ5169bx1FNPUVpais1m45lnniEpKck1zfXr1zEYDPj7+wMwb948du7cSWJiIl26dGH48OEABAUF8W//9m/07t2b8PBwkpOTa6zrpz/9KV988YWrS2h5eTlz5szBZDKhKApvv/12M221aA+kOqcQjfDxxx+Tk5PDCy+80Kjl/L//9/8oLS3lv/7rv5ooMiHqJolfCA+bN2+eq4dRaGiop8MR7YAkfiGEaGfk4q4QQrQzkviFEKKdkcQvhBDtjCR+IYRoZyTxCyFEOyOJXwgh2pn/Dzq6AoCs7nWvAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np\n",
    "from scipy.integrate import odeint\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# Total population, N.\n",
    "N = 1000\n",
    "# Initial number of infected and recovered individuals, I0 and R0.\n",
    "I0, R0 = 1, 0\n",
    "# Everyone else, S0, is susceptible to infection initially.\n",
    "S0 = N - I0 - R0\n",
    "# Contact rate, beta, and mean recovery rate, gamma, (in 1/days).\n",
    "beta, gamma = 3/5, 1/3\n",
    "# A grid of time points (in days)\n",
    "t = np.linspace(0, 100, 100)\n",
    "\n",
    "# The SIR model differential equations.\n",
    "def deriv(y, t, N, beta, gamma):\n",
    "    S, I, R = y\n",
    "    dSdt = -beta * S * I / N\n",
    "    dIdt = beta * S * I / N - gamma * I\n",
    "    dRdt = gamma * I\n",
    "    return dSdt, dIdt, dRdt\n",
    "\n",
    "# Initial conditions vector\n",
    "y0 = S0, I0, R0\n",
    "# Integrate the SIR equations over the time grid, t.\n",
    "ret = odeint(deriv, y0, t, args=(N, beta, gamma))\n",
    "S, I, R = ret.T\n",
    "\n",
    "# Plot the data on three separate curves for S(t), I(t) and R(t)\n",
    "fig = plt.figure(facecolor='w')\n",
    "ax = fig.add_subplot(111, facecolor='#dddddd', axisbelow=True)\n",
    "ax.plot(t, S/N, 'b', alpha=0.5, lw=2, label='Susceptible')\n",
    "ax.plot(t, I/N, 'r', alpha=0.5, lw=2, label='Infected')\n",
    "ax.plot(t, R/N, 'g', alpha=0.5, lw=2, label='Recovered with immunity')\n",
    "ax.set_xlabel('Time /days')\n",
    "ax.set_ylabel('Number (% of population)')\n",
    "ax.set_ylim(0,1)\n",
    "ax.yaxis.set_tick_params(length=0)\n",
    "ax.xaxis.set_tick_params(length=0)\n",
    "ax.grid(b=True, which='major', c='w', lw=2, ls='-')\n",
    "legend = ax.legend()\n",
    "legend.get_frame().set_alpha(0.5)\n",
    "for spine in ('top', 'right', 'bottom', 'left'):\n",
    "    ax.spines[spine].set_visible(False)\n",
    "plt.show()\n",
    "\n",
    "#Credit to this website for this model: https://scipython.com/book/chapter-8-scipy/additional-examples/the-sir-epidemic-model/"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "altered-seller",
   "metadata": {},
   "source": [
    "#### Contact Number\n",
    "\n",
    "We assumed that the parameters  $\\beta$  and  $\\gamma$  could be estimated somehow, and therefore it would be possible to generate numerical solutions of the differential equations. In fact, as we have seen, the fraction  $\\gamma$  of infecteds recovering in a given day can be estimated from observation of infected individuals. Specifically,  $\\gamma$  is roughly the reciprocal of the number of days an individual is sick enough to infect others. For many contagious diseases, the infectious time is approximately the same for most infecteds and is known by observation.\n",
    "\n",
    "There is no direct way to observe  $\\beta$,  but there is an indirect way. Consider the ratio of  $\\beta$  to  $\\gamma$:\n",
    "\n",
    "$\\frac{\\beta}{\\gamma} = \\beta \\frac{1}{\\gamma} =$ the number of close contacts per day per infected times the number of days infected = the number of close contacts per infected individual.\n",
    "\n",
    "We call this the contact number, $c$, where $c = \\frac{\\beta}{\\gamma}$\n",
    "\n",
    "From our previous equations, we have:  \n",
    "$ \\frac{ds}{dt}=- \\beta s(t) i(t), $  \n",
    "$ \\frac{di}{dt}=\\beta s(t) i(t) - \\gamma i(t), $  \n",
    "\n",
    "Dividing the second equation by the first, we have:  \n",
    "$\\frac{di}{ds}=\\ -1+\\frac{c}{s}$  \n",
    "(Explain how we arrived at this step using the chain rule)\n",
    "\n",
    "Important conclusions about this equation:  \n",
    "- The only parameter that appears is  $c$,  the one we are trying to determine.  \n",
    "- The equation is independent of time. That is, whatever we learn about the relationship between  $i$  and  $s$  must be true for the entire duration of the epidemic.  \n",
    "- The right-hand side is an explicit function of  $s$,  which is now the independent variable.  \n",
    "\n",
    "Show that $i(s)$ must have the form:  \n",
    "\n",
    "$i = -s + \\frac{1}{c}ln(s) + q,$ where $q$ is a constant.\n",
    "\n",
    "And explain why $i + s - \\frac{1}{c}ln(s)$ must be independent of time.\n",
    "\n",
    "There are two times when we know (or can estimate) the values of  $i$   and  $s$: at $t = 0$  and  $t = \\infty$.  For our considerations,  $i(0)$  is approximately  $0$  and  $s(0)$  is approximately  $1$.  A long time after the onset of the epidemic, we have  $i(\\infty)$  approximately  $0$   again, and  $s(\\infty)$  has settled to its steady state value. If there has been good reporting of the numbers who have contracted the disease, then the steady state is observable as the fraction of the population that did not get the disease.\n",
    "\n",
    "For such an epidemic, explain why\n",
    "\n",
    "$c=\\frac{\\ln\\left(s_{\\infty}\\right)}{s_{\\infty}-1}$\n",
    "\n",
    "[Hint: Use the fact that the quantity $ i + s - \\frac{1}{c}ln(s) $ is the same at  $t = 0$  and at  $t = \\infty$.]\n",
    "Use one of your numerical solutions to estimate the value of  $s(\\infty)$. Compare your calculated value with the one you get by direct calculation from the definition,  $c = \\frac{\\beta}{\\gamma}$."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "organized-slovak",
   "metadata": {},
   "source": [
    "#### Herd Immunity\n",
    "\n",
    "Each strain of flu is a disease that confers future immunity on its sufferers. For such a disease, if almost everyone has had it, then those who have not had it are protected from getting it -- there are not enough susceptibles left in the population to allow an epidemic to get under way. This group protection is called herd immunity.  \n",
    "\n",
    "Earlier you experimented with the relative sizes of  $\\beta$  and  $\\gamma$,  and you found that, if  $\\beta$  is small enough relative to  $\\gamma$,  then  no epidemic can develop. If the contact number  $c = \\frac{\\beta}{\\gamma}$  is small enough, then there will be no epidemic. But another way to prevent an epidemic is to reduce the initial susceptible population artificially by inoculation.  \n",
    "\n",
    "So, what fraction of the population must be inoculated to obtain herd immunity? Or, put another way, how small must  $s_0$  be to insure that an epidemic cannot get started? It depends on the contact number.\n",
    "\n",
    "Explain why keeping an epidemic from getting started is the same as keeping  $\\frac{di}{dt}$  negative from  $t = 0$   on.  \n",
    "Write the right-hand side of the infected-fraction differential equation:  \n",
    "\n",
    "$\\frac{di}{dt}=\\beta s\\left(t\\right)i\\left(t\\right)- \\gamma i\\left(t\\right)  $\n",
    "\n",
    "in factored form. Explain why one factor is always positive and why the sign of other factor depends on the size of  $s(t)$.\n",
    "Explain why  $s(t)$  is a decreasing function, and thus has its largest value at  $t = 0$. It follows that, if $\\frac{di}{dt}$ is negative at time  $0$,   then it stays negative.\n",
    "Show that\n",
    "\n",
    "$i'\\left(0\\right)=\\left(\\beta s_{0}-\\gamma\\right)i_{0}$\n",
    "\n",
    "Explain why, if  $s_0$  is less than  $1/c$,   then no epidemic can develop.  \n",
    "\n",
    "\n",
    "From 1912 to 1928, the contact number for measles in the U.S. was 12.8. If we assume that  $c$  is still 12.8 and that inoculation is 100% effective -- everyone inoculated obtains immunity from the disease -- what fraction of the population must be inoculated to prevent an epidemic?  \n",
    "\n",
    "$s_0 < 1/12.8$, so at least 92.1875% of the population should be vaccinated.\n",
    "\n",
    "Suppose the vaccine is only 95% effective. What fraction of the population would have to be inoculated to prevent a measles epidemic?\n",
    "\n",
    "$s_0 < 1/12.8$ and 1/20 of $s_0$ is vaccinated but does not receive immunity. Similar calculation, but I think you have to multiply the total amount of people needed to be vaccinated by 21/20 to compensate? (I'm going off of the idea that if a vaccine had 50% efficacy then you would need to vaccinate twice as many people to have an equivalent immunization as a 100% effective vaccine.) The number I got was 96.796875%\n",
    "\n",
    "[//]: # (let me know if the above problem was solved right, thanks!)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "employed-press",
   "metadata": {},
   "source": [
    "## Questions\n",
    "\n",
    "Question 1: Explain briefly the modeling steps that lead to the SIR model"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "moved-belief",
   "metadata": {},
   "source": [
    "Question 2: Given a population and disease combination for which the SIR model is appropriate, what are the possible outcomes when a trace of infection is introduced into the population? How can you tell whether there will be an epidemic?"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "gothic-commonwealth",
   "metadata": {},
   "source": [
    "Question 3: Does \"epidemic\" mean that almost everyone will get the disease? If so, what keeps the spread of disease going? If not, what causes the epidemic to end before everyone gets sick?"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "interim-serial",
   "metadata": {},
   "source": [
    "Question 4: How can it happen that a large percentage of a population may get sick during an epidemic even though only a small percentage is sick at any one time?"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "headed-reviewer",
   "metadata": {},
   "source": [
    "Question 5: Explain briefly the key idea for finding solutions of an SIR model without finding explicit solution formulas."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "outer-casting",
   "metadata": {},
   "source": [
    "Question 6: Describe briefly the meaning and significance of contact number."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "instant-point",
   "metadata": {},
   "source": [
    "Question 7: Describe briefly the meaning and significance of herd immunity. How can an inoculation program lead to herd immunity?"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "short-tribune",
   "metadata": {},
   "source": [
    "Question 8: The contact number for poliomyelitis in the U.S. in 1955 was 4.9. Explain why we have been able to eradicate this disease even though we cannot eradicate measles. Give a careful argument -- \"smaller contact number\" is an observation, not an explanation."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "configured-grain",
   "metadata": {},
   "source": [
    "#### Answers"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "possible-peace",
   "metadata": {},
   "source": [
    "1: See the background section, but more importantly, how exactly should this be graded? More open-ended or more concrete?  \n",
    "2: Possible Outcomes: epidemic or no epidemic, but also diseases can be endemic vs. non-endemic, depending on the contact number.  \n",
    "3: More or less, it means \"the rapid spread of disease to a large number of people in a given population within a short period of time.\" But to answer the question more fully, the spread of disease goes until people have developed a resistance/immunity to the disease, at which point the disease becomes either endemic (constant, like the flu/measles) or non-endemic, and dies out (e.g. polio.)  \n",
    "4: Essentially, because diseases spread, and it takes time for the disease to infect a host and then others. There is a flu season, and also we have first hand experience with epidemics because of the COVID-19 pandemic.  \n",
    "5: We use the euler method for numerically solving our ODE, as it is difficult to know the parameters of an epidemic a priori.\n",
    "Also, I'm not sure what to use for the ODE and euler method, there are files for a bunch of them on the introduction, but there are also some other options, e.g. a lab going over all the same stuff but using a link for math.studio? I would use matlab, but the above python graph could also work, maybe?\n",
    "\n",
    "6: Contact number: \"it measures the relative contagiousness of the disease, because it tells us indirectly how many of the contacts are close enough to actually spread the disease.\"  \n",
    "7: Herd Immunity: Reducing $s_0$ enough so that no epidemic occurs, $s_0 < 0$ for no epidemic.  \n",
    "8: Might have to look into the actual numbers. I know historically that Salk invented a polio vaccine in the 1950s, so my guess is that we were able to inoculate enough people so that $s_0 < 1/4.9,$ or $s_0 <$ 20.4% of the population. However, with Measles, not only does it have a larger contact number, but I think it occurs mostly in children, especially if they have not been vaccinated. In other words, there is not herd immunity for Measles like there is for Polio. (At least in the US, other places still have these diseases endemically, which is why you might need booster shots when travelling abroad.)\n",
    "\n",
    "I feel like I've more or less done all I can do for now, now I just need to clean up this notebook so it is presentable."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "likely-semester",
   "metadata": {},
   "source": [
    "#### Sources\n",
    "- Smith D. & Moore L. (2004, December)  \n",
    "The SIR Model for Spread of Disease - Introduction. JOMA.  \n",
    "Retrieved May 18, 2021, from https://www.maa.org/press/periodicals/loci/joma/the-sir-model-for-spread-of-disease-introduction \n",
    "\n",
    "- https://scipython.com/book/chapter-8-scipy/additional-examples/the-sir-epidemic-model/"
   ]
  }
 ],
 "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.9.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}