{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Tipping point models, epidemics, and vaccination\n",
    "\n",
    "The SIS model from epidemiology as an example of a model exhibiting a tipping point. This model allows us to calculate the percentage of a population that needs to be vaccinated in order to get population level herd immunity.\n",
    "\n",
    "There's a nice animation [here](http://www.theguardian.com/society/ng-interactive/2015/feb/05/-sp-watch-how-measles-outbreak-spreads-when-kids-get-vaccinated).\n",
    "\n",
    "The discussion here is heavily based on the discussion of tipping points in Scott Page's very nice *coursera* course on [Model Thinking](https://www.coursera.org/learn/model-thinking/home/week/4). I essentially transcribed his algebra in SymPy.\n",
    "\n",
    "For more information on the models discussed here, hop over to [Wikipedia](https://en.wikipedia.org/wiki/Compartmental_models_in_epidemiology)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Diffusion model\n",
    "\n",
    "The first model we develop is a simple diffusion model. It has no tipping, but shows that the spread of infectious diseases is an accelating process, that becomes quicker as more people are infected.\n",
    "\n",
    "First some general definitions of symbols."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import sympy as sm\n",
    "sm.init_printing()\n",
    "from sympy.abc import N, c, t, tau\n",
    "Wt = sm.symbols('Wt')\n",
    "Wtt, Penc = sm.symbols('Wtt,Penc',cls=sm.Function)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let $W_t$ be the number of sick people in the population, whose overall size is $N$. Here is the formula for $P_{enc}$,  the probability that a healthy person and a sick person meet. These are the cases in which an infection can happen. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "Penc = (Wt/N)*((N - Wt)/N)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Using this we can write down the number of infected individuals in the next time period. It is $W_t$ plus the number of newly infected inviduals which depends on $P_{enc}$ as well as on $c$, the probability the people in the population interact, and $\\tau$, the transmission rate of the disease. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "Wtt=N*c*Penc*tau + Wt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAMoAAAAsBAMAAADWVjoEAAAAMFBMVEX///8AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAv3aB7AAAAD3RSTlMAdt27RKvvEFRmiTIimc3h244VAAAACXBIWXMAAA7EAAAOxAGVKw4bAAADE0lEQVRYCbVXPWhTURT+mrz60teXNCoI4mCm4h+kgwgOmgwORZAGB+tccdVURTcJ7VjBuDkoiSCE4tAugoNDHBwExSwuCja4OAhpbFVoaRvPuQ+ad09e20vJPcPtOfc73/1yz+097z1g3/bBjPloz7RDnbO4tjKJmZVyT647Bgzk/gAz/yY1UJAG6xoaETjrgEcLDRcIPKMnpOoUx1dp+KkDguQ2dbg38jcDFY+gRFvHn3M4nqsDVfa6JkkvulC05/4FUmvAJYKTTT1Hba020YBT0gFJuq/DEdEWcPR38GtjVQ33Mxw2Ym14RQ0ABCmVF3hPOAV/bh0DdeBK7kIa3sFjONLpdNJ08A1ObgytISlZkqQyZVI4HoVHZZ7mqceAfwIXvVfv336kMF6mwSlhBXxomukkDDU1NCJYxpy7Ghzvd4CKdv0N4iovVqI/VKyFPB+aZjoJiYyGRgS5YhpbTp6RU8BSkZ3XPOBAngYq1uexKjnAwxbbL3Z1EobbPLmbLb0DplRJXLo2N1TqvBqVCiHxjZ6y6yQDlYWnwOgzXtZp4wHdUbLTalQVG6f5zaqKQ4NGKsLJhLBI9/gYsJxmyGsmFunmJOAHBVCnXyMgm2c4bDrJ4PQrdeCqWiFZmMYXuN8wkFGxx4Vit6LC8KCTkKyGwSj/Mk2+VIB/vgxv/gmduNoa18E9x010sIenk5BaBEQH7aHsOCGa5455uEOQ6KA7J0tEdUs5GRVzt5QdlPMSQVGiKN25eL3r7+a5TUJlB2WCkQo/xUws+DWigzLRSAWGT2R12yA6qLmKyUa2c0QHtaQiOqglFdFB4bVaP261WoXtzfbFER00vJcReiz2xQDRQcMqfdlFsIjooJZURAe1pCI6qCWVoG7aaHb3Nco+ArdoSjrcySORvWmcb7quyMtRz0zZFnFmN+iVQCj3PfTu0VcBvcnYtSSyAL/JWLUkKmU0rErQ4uMYyaBqW6VGH4TBa7RNpQb8Vev/Ylys3Ceb2+C1nRIw8ZUcq8aflDH65LBrMVp+qG1XA3dvnySFwt4q/wHZ3ixJ4HZQCwAAAABJRU5ErkJggg==\n",
      "text/latex": [
       "$$Wt + \\frac{Wt c \\tau \\left(N - Wt\\right)}{N}$$"
      ],
      "text/plain": [
       "     Wt⋅c⋅τ⋅(N - Wt)\n",
       "Wt + ───────────────\n",
       "            N       "
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Wtt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To understand the dynamics of this process it is enough study how $P_{enc}$ changes as $W_t$ increases. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAHwAAAAsBAMAAABVvsF6AAAAMFBMVEX///8AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAv3aB7AAAAD3RSTlMAdt27RKvvEFRmiTIimc3h244VAAAACXBIWXMAAA7EAAAOxAGVKw4bAAACf0lEQVRIDd1Wv4vTYBh+mqbXNja9gqtDBB2Ug3YSHKQZHLrIdXMTFwen1h84esVbHITrH6Dcd1MX4bq4ehl1kAYcBAfJ5iBovB7KQY/6vLmCyZeKabuID+TlfZ/3efjefE2+BkjjbZqaxzwTMlcfu4XOxW6pfmmmKTbINo+AJz9uJmyatOBJd1MB+y7wivkGL1Q9hvyY4TOvOJLSYiC9kX9qV0ApFOKFhFbTA5RkMWjSPWlttYFvfZguYAdCRCMMNn2YPVZxaNJH0lsPkNutwWJqKIaKwwDfCGF1JYtBk1Zd9tYcvDlQOAfcaF6rcdd8MfjlY9iSxKFJLVEat6FGfpFG7Ig232cweziMBmL+G5q0HLBVDUvDrYbMjk8SjB4D5953r0sZhyYtOWzmf+Z4/5H0smjXXAbOPWooJsDjr4IvkmrSMyE5e1Lj/StmRT4qMzuHyU/k1hLQpJHdOu7BuOtRZ4bozoZvsTpRCS8LTWo65HJ8vqrfRWkFpeFs6wYs666QcWjSaOvKE6DgiMpu32OMfg6pd3kloUltxbbJYLlMULnaZ+RIxSuHzArCJaBJq5w1jeihTdNp5n6aIhO9MnM7Grmn1adl3ptLp8hikKKEkOMiC/60TMbD6nmWNf5nzXQl/PM7Y73/sMqML3FrFfsGDrwl/GenLkr1O90LGLlL2NHku1Pt0rkjYVGY2zzGovP//KJW0VsPT/hPK4m/jN1GHZAjN/Pxk1jFxm4fXNhU5jDRyFa0sO5AAa+fbi+zdQNYR/LR0JlOs62XVPmojKONT9JZKwU032UVp3Rmj19dH1N0VkK+bIzoLzurI6EzWJXDBLVA8aAjX5ntvzh+AVU/69IssXRFAAAAAElFTkSuQmCC\n",
      "text/latex": [
       "$$\\frac{Wt \\left(N - Wt\\right)}{N^{2}}$$"
      ],
      "text/plain": [
       "Wt⋅(N - Wt)\n",
       "───────────\n",
       "      2    \n",
       "     N     "
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Penc"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First, let's see what happens when $W_t$ is small and $N$ is large."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAADcAAAAqBAMAAADlvd4wAAAAMFBMVEX///8AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAv3aB7AAAAD3RSTlMAdrtEVN3vqxDNIomZZjLe39VDAAAACXBIWXMAAA7EAAAOxAGVKw4bAAABE0lEQVQ4EWNgYBD6b8DAEp/ewIAV6CswMHCiy7FkgNUyO/1iYOBG09a+/QdYhLv5DwODCZokAy9EkochnoFhB07J+xcYDuCQNGHgD2DYgENyBwP3Z2YDHJIHGLg+QRzLMRMMHoAUQh0ENFF/DbpGqCSzAwOD/WkcktwNDAysn3FIsgLFeT9gSDKCAqEpPxZIFqBL8pT9PIsuNnz5/3EDzGCjIBi4D53ArXsXQyVuySiG9wswZGFZJphhvQGGJAM8y/g3YEgiskwghhwDPMtwH8CUhGeZO5hyDDwMkCzDvIH5AYY0LMu8dHFqwJCEZZn8//8x5ICZEJZlMOVAmRBLloEoxJFlIJI4sgxEEkeWAUtizTIAYgtqLh3SdB8AAAAASUVORK5CYII=\n",
      "text/latex": [
       "$$\\frac{N - 1}{N^{2}}$$"
      ],
      "text/plain": [
       "N - 1\n",
       "─────\n",
       "   2 \n",
       "  N  "
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Penc.subs([(Wt,1)])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAABMAAAAqBAMAAACjC32aAAAALVBMVEX///8AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAADAOrOgAAAADnRSTlMAVO8Qq4lmdpnNu0TdIoKkHS0AAAAJcEhZcwAADsQAAA7EAZUrDhsAAACTSURBVBgZY2BABYzKML5JyGMYk4GN9kx2uBUcTi/b4DbTgPEOBh6Qb3jVuw0MjOdUDEAmrFvAwMAJZnFLvWVgYAWbymr9goFhI5jJwXCOgSEYysxLYGgAMzcy8B1gCAAzgxlYX3NvADMbGHhfQQwAya27ARbkFmBg2NcKZrIaMDDwvAYzeYAk2wMQ00rvGJB0AGIAFtkr3Lqms2kAAAAASUVORK5CYII=\n",
      "text/latex": [
       "$$\\frac{1}{N}$$"
      ],
      "text/plain": [
       "1\n",
       "─\n",
       "N"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Penc.subs([(Wt,1)]).subs([(N-1,N)])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Or you can do it in one step"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAABMAAAAqBAMAAACjC32aAAAALVBMVEX///8AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAADAOrOgAAAADnRSTlMAVO8Qq4lmdpnNu0TdIoKkHS0AAAAJcEhZcwAADsQAAA7EAZUrDhsAAACTSURBVBgZY2BABYzKML5JyGMYk4GN9kx2uBUcTi/b4DbTgPEOBh6Qb3jVuw0MjOdUDEAmrFvAwMAJZnFLvWVgYAWbymr9goFhI5jJwXCOgSEYysxLYGgAMzcy8B1gCAAzgxlYX3NvADMbGHhfQQwAya27ARbkFmBg2NcKZrIaMDDwvAYzeYAk2wMQ00rvGJB0AGIAFtkr3Lqms2kAAAAASUVORK5CYII=\n",
      "text/latex": [
       "$$\\frac{1}{N}$$"
      ],
      "text/plain": [
       "1\n",
       "─\n",
       "N"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Penc.subs([(Wt,1),(N-1,N)])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Thus, when $W_t$ is small, the probability of encountering someone who is sick is low. What happens when $W_t$ is large, i.e. $W_t\\approx N$?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAAoAAAAOBAMAAADkjZCYAAAAMFBMVEX///8AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAv3aB7AAAAD3RSTlMAiXZmMs1UEN0i77urRJlR0qN3AAAACXBIWXMAAA7EAAAOxAGVKw4bAAAAUUlEQVQIHWNgYFQWYWBgCGOomMDAvICBMYCB+wAD23cG/gMMvN8Y6h8w8H5imC/AwAIkHzCwfISKAGXZvjFwb2Bg/g7VxdDGUOXAwFCodIQBAG3HFgUteuAKAAAAAElFTkSuQmCC\n",
      "text/latex": [
       "$$0$$"
      ],
      "text/plain": [
       "0"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Penc.subs([(Wt,N)])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Zero?? How come?! The reason is simple (and symmetric): in this case the probability of encountering someone who is *healthy* is very low, hence the probability of an encouter between a sick person and a health person is low. To complete the picture, let's check what happens when $W_t=N/2$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAAsAAAArBAMAAABcLm8jAAAALVBMVEX///8AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAADAOrOgAAAADnRSTlMAVO8Qq4lmdpnNMt0iu/rAsXEAAAAJcEhZcwAADsQAAA7EAZUrDhsAAABfSURBVBgZY2CAASEDIIvZ1Q9EMTDkUUbVgbW36M1IAJtGFvEOBB4QqXUFRN0WMMWuAqaOTQFTAWCK7QCYMmcAUwFginEBmDJlAFMiHR2P2sE6INoZdMGc2HduYBpIAABFmBuw9gyMBgAAAABJRU5ErkJggg==\n",
      "text/latex": [
       "$$\\frac{1}{4}$$"
      ],
      "text/plain": [
       "1/4"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Penc.subs([(Wt,N/2)])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Ha!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For completeness let's put down some assumptions and graph $W_t$ over time. Assume N=100, $\\tau=0.5$ (which is very high). Start with only on person sick, i.e. $W_0=1$. Let $c=0.5$ too."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "Wupdate=Wtt.subs([(N,100),(tau,0.5),(c,0.5)])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "ws=[1]\n",
    "for t in range(0,49):\n",
    "    ws.append(Wupdate.subs([(Wt,ws[-1:][0])]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib notebook"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Populating the interactive namespace from numpy and matplotlib\n"
     ]
    }
   ],
   "source": [
    "%pylab inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEWCAYAAACJ0YulAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3dd3gc5bn38e+talmWi3DvNjbY4AY2xoGEnoReAyQEYgiBJC8QUji0FEIIoSRA4OQcEicUJySAAxxKKAEcmxYC7rjhLttyk1xkyVaztPf7x4yEMLK8lrRa7e7vc1177c7s7Mw9K+3e+5R5HnN3REREANLiHYCIiLQfSgoiIlJPSUFEROopKYiISD0lBRERqaekICIi9ZQUpF0ws9+b2U8bLH/XzLaY2S4zO8jMjjWzFeHyuS04zqtmNrl1oj6g4/7SzLaa2eZGnvuCmS07gH196r1pxRi/bmavt9b+JDGZrlOQWDOzAqAXUAPUAkuAPwNT3D3SyPaZQCkwyd0XhOumAy+6+4NtFXdrMbMBwHJgkLsXtXBfn3lvmrmfwcAaINPda1oSkyQXlRSkrZzl7nnAIOBu4CbgkX1s2wvoACxusG7QXsuJZBCwraUJIdTYeyPSapQUpE25+053fxG4GJhsZqMAzOzxsIrlEKCuKqXEzP5lZquAocBLYZVJtpkVmNkpdfs1s5+b2RPh4w5m9oSZbTOzEjObZWa9wudmmtm3wsdpZvYTM1trZkVm9mcz6xI+N9jM3Mwmm9m6sOrnx/s6LzPrEr6+ONzfT8L9nwK8AfQNY3+8kdeeYGaFDZYLzOwGM/vIzHaa2dPhOX3mvQm3H2Fmb5jZdjNbZmYXNdhXjpndF8a008zeNbMc4O0G+9llZp8zs8vN7N0Grz0mfO92hvfHNHhuppndYWbvmVmZmb1uZt2b+NNLglBSkLhw9w+BQuALe61fDhweLnZ195Pc/WBgHUFpo5O7V+1n95OBLsAA4CDgO0BFI9tdHt5OJEg6nYDf7bXN54FDgZOBn5nZyH0c87/DYw4Fjge+AVzh7m8CpwEbw9gv30/sdS4CTgWGAGOAyxt7b8wslyDp/A3oCXwN+F8zq9vuN8B44BggH7gRiADHNdhPJ3d/v+HBzSwfeBl4iOA9vB94ea82jEuAK8LjZgE3RHlu0o4pKUg8bST4omptewi+yIa5e627z3H30ka2+zpwv7uvdvddwC3AV80so8E2t7t7RVh/vwAYu/dOzCydoORzi7uXuXsBcB9wWQvO4SF33+ju24GXgHH72O5MoMDdH3P3GnefCzwLfMXM0oBvAte7+4bwvfh3FEkV4Axghbv/Jdzvk8DHwFkNtnnM3Ze7ewUwrYkYJYEoKUg89QO2x2C/fwH+CTxlZhvN7N6wgXZvfYG1DZbXAhkE9fZ1GvYWKicoTeytO8Ev5b331a8ZsR/IcSForzg6rCYrMbMSgmTXO4yrA7CqGcff+72Bz55TtDFKAlFSkLgws6MIvmDe3d+2+7Ab6NhguXfdA3ff4+63u/thBNUmZxJU5+xtI8GXap2BBD2kthxgLFsJSid772vDAe6nOdYDb7l71wa3Tu7+3TCuSuDgRl63v26He7830HbnJHGkpCBtysw6m9mZwFPAE+6+sJm7mk9Q1ZNpZhOArzQ4xolmNjqs1ikl+MKubWQfTwI/MLMhZtYJ+BXw9IF20XT3WoLqkzvNLM/MBgE/BJ5o1pkdmH8Ah5jZZeF7kWlmR5nZyLC776PA/WbW18zSwwblbKCYoG1h6D72+0q430vMLMPMLgYOC48nSUxJQdrKS2ZWRvDL9scEDZdXtGB/PyX4BbwDuJ2gobVOb+AZgoSwFHiLxr+gHyWoanqboM9+JXBdM+O5jqD0spqg9PO3cP8x5e5lwJeArxL8ut8M3ANkh5vcACwEZhFU1d0DpLl7OXAn8F5Y7TRpr/1uIyhh/QjYRtBAfaa7b431OUl86eI1ERGpp5KCiIjUU1IQEZF6SgoiIlJPSUFEROpl7H+T9qt79+4+ePDgeIchIpJQ5syZs9XdezT2XEInhcGDBzN79ux4hyEiklDMbO+r1eup+khEROopKYiISD0lBRERqaekICIi9ZQURESkXsySgpk9Gk5xuKjBuvxw2sAV4X23cL2Z2UNmtjKcgvDIWMUlIiL7FsuSwuME0wk2dDMw3d2HA9PDZQimKxwe3q4GHo5hXCIisg8xu07B3d82s8F7rT4HOCF8PBWYCdwUrv+zB0O2/sfMuppZH3ffFKv4RCQ1uTs1EWdPbYQ9NU51bYQ9tRFqap09kQi14XPBvVMbcWoiESIRqHWnNhKhNgK1kQgRh9qIE/HwFm7j7kSccH1wzEjEcahfdgen7vngsX/quU/WU/+4/iQ4eWQvxg7o2urvT1tfvNar7ove3TeZWc9wfT+CcfbrFIbrPpMUzOxqgtIEAwcOjG20IhI37k5pZQ07y/dQUlHNrsoayqpq2FVZw66q4FZWWUNFdQ3l1bWU76mlsrqW8upaKvbUUrmnluqaCFX1t1qqaiJU10TifWqtomfnDkmRFPbFGlnX6EQP7j4FmAIwYcIETQYhkkDcneJdVWzeWUlxWRVFZVUUlVZRvKuSotIqtu2uZkd5dZgI9lAbafojnpludMzKoGNWOjmZ6eRkpdMxK528Dhn0yMsmOyON7Ix0sjPTyM5IIytczko3MtLTyExPIyvdyExPC5eNjLQ00tOMzHQjPe2T5U/dzEhLg4y0NNIM0tKMNAvWm0F6uJxmYA3uDUgzw9L45LGBEd43fMwnr7HwG9Kssa/K1tXWSWFLXbWQmfUBisL1hcCABtv1J5hFSkQSTE1thIJt5awsKmPttnIKd1Swfkc567cHj6sa+aXerWMmPfKy6d4pm5G9O9O1YyZdO2bSrWMWXTtm0SUnk7wOGXTKzqi/79Qhg+yM9DicYXJr66TwIjAZuDu8f6HB+mvN7CngaGCn2hNE2r8tpZXMX1/C8s1lLC/axYotZawu3k117Sdf/F1yMunfLYfhPfM4aURP+nfrSJ8uHeiRl03Pzh3o3ilLX+7tSMySgpk9SdCo3N3MCoHbCJLBNDO7ElgHXBhu/gpwOrASKKdlc/eKSAxU1dSyaEMp89btYN76Euat3cHGnZX1z/fvlsMhvfI4/tAeHNIzj0N65TGoe0c6d8iMY9RyoGLZ++hr+3jq5Ea2deCaWMUiIgfO3VlZtIsZy4qYuayY2QU76ksA/brmcOSgbnxrYDfGDujKiN555Ga3lyZKaQn9FUWkXkV1Le+t3FqfCDaUVAAwoncek48ZxPhB+Rw5sCs9O3eIc6QSK0oKIikuEnFmFWznmTmFvLJwE7ura8nNSufYYd255sRhnHBoD/p2zYl3mNJGlBREUtT67eU8O7eQZ+cWsn57BblZ6Zw+ug/njOvHUUO6qfE3RSkpiKQQd+e9ldt4+K2VvLdyG2ZwzMEH8YNTDuHUUb3pmKWvhFSn/wCRFBCJOG8u3cL/zFzFgvUl9OqczY++eAjnj+9PP1UNSQNKCiJJrKY2wssLN/G/M1axbEsZA/M78qvzRnPB+H6qHpJGKSmIJKkZHxdx+0uLKdhWzvCenfjtxeM4c0wfMtI1jYrsm5KCSJLZUFLB7S8u5vUlWzi4Ry6/v3Q8XzqsF2lpsR83RxKfkoJIkqiuifDIu2t4aPoKHOfGUw/lW58fSlaGSgYSPSUFkSTw/qpt/PSFRaws2sWXDuvFz846jP7dOsY7LElASgoiCay6JsLdr37Mo++toX+3HB6ZPIGTR/aKd1iSwJQURBLUxpIKrv3bXOauK+HyYwZz06kjyMlSjyJpGSUFkQT01vJivv/UPKprIvzPJUdyxpg+8Q5JkoSSgkgCqY04D05fwX//awWH9Mzjfy89koN7dIp3WJJElBREEsS2XVVc/9R83l25lQuO7M8vzx2l6iJpdUoKIglgQ0kFl/3pAwpLKrjngtFcNGFAm8zXK6lHSUGknVuzdTdf/+N/KKus4a/fOpqjBufHOyRJYkoKIu3Y0k2lXPbIh0TcefLqSYzq1yXeIUmS2++ljmZ2rJnlho8vNbP7zWxQ7EMTSW1z1+3g4j+8T0aaMe3bSgjSNqK5/v1hoNzMxgI3AmuBP8c0KpEU9++VW7n0Tx/QLTeLv3/ncwzrmRfvkCRFRJMUatzdgXOAB939QUD/oSIx8q+Pt3D547Po3y2Hv3/7cwzI13AV0naiaVMoM7NbgEuB48wsHciMbVgiqWnO2u1894m5jOidx9QrJtItNyveIUmKiaakcDFQBVzp7puBfsCvYxqVSApaXbyLb02dTd+uOTyuhCBxst+SQpgI7m+wvA61KYi0qq27qrj8sVmkmfH4FUeRr4QgcbLPpGBmZYA39hTg7t45ZlGJpJCK6lqunDqborJKnrxqEoMOyo13SJLC9pkU3F2NySIxVhtxvvfUPD4qLOEPl47niIHd4h2SpLioL14zs55Ah7rlsBpJRJrJ3bn9pcW8sWQLt599OF86vHe8QxKJ6uK1s81sBbAGeAsoAF6NcVwiSe+P76zmz++v5aovDGHyMYPjHY4IEF3vozuAScBydx8CnAy8F9OoRJLcv1du5a5XP+aM0X245bSR8Q5HpF40SWGPu28D0swszd1nAONiHJdI0iouq+L6p+cztHsuv75wDGlpGu1U2o9o2hRKzKwT8DbwVzMrAmpiG5ZIcopEnB9Om09pxR7+cuVEOmZpTEppX6IpKZwDVAA/AF4DVgFnteSgZvYDM1tsZovM7Ekz62BmQ8zsAzNbYWZPm5k6akvSefitVbyzYiu3nXU4I3qrV7e0P/tNCu6+291r3b3G3ae6+0NhdVKzmFk/4HvABHcfBaQDXwXuAR5w9+HADuDK5h5DpD2aXbCd+99Yzplj+vC1iQPiHY5Io6LpfVRmZqXhrdLMas2stIXHzQByzCwD6AhsAk4Cngmfnwqc28JjiLQbO3ZX870n59G/Ww53nT9as6ZJuxXNMBefuojNzM4FJjb3gO6+wcx+A6wjqJZ6HZgDlLh7XVtFIcEYS59hZlcDVwMMHDiwuWGItBl357+eWUDxriqe++6x5HXQeJLSfkXTpvAp7v48wa/6ZjGzbgTtFEOAvkAucFpjh9rH8ae4+wR3n9CjR4/mhiHSZh59r4A3lxZx6+kjGd1fE+VI+7bfkoKZnd9gMQ2YwD6+sKN0CrDG3YvD/T8HHAN0NbOMsLTQH9jYgmOItAtLNpZy96tL+eJhvbhcF6hJAoimP1zDnkY1BFc0n9OCY64DJplZR4Lqo5OB2cAM4CvAU8Bk4IUWHEMk7mpqI9z07Ed0ycnk3gvGqB1BEkI0bQpXtOYB3f0DM3sGmEuQZOYBU4CXgafM7Jfhukda87gibe3R99awcMNOfnfJEZobQRJGU0Nn/zdNVBO5+/eae1B3vw24ba/Vq2lBA7ZIe7J2227uf2M5p4zsxRmj+8Q7HJGoNdXQPJugV1AH4EhgRXgbB9TGPjSRxOTu3PLcQjLT0vjluaNUbSQJpan5FKYCmNnlwInuvidc/j1BN1IRacS02ev596pt3HneKHp36bD/F4i0I9F0Se0LNLxWoVO4TkT2UlRayS9fXsrRQ/L52lG6jkYSTzS9j+4G5pnZjHD5eODnMYtIJIH97IXFVNdEuPsCjX4qiSma3kePmdmrwNHhqpvdfXNswxJJPK8t2sRrizdz06kjGNJd8yxLYtpn9ZGZjQjvjySoLlof3vqG60QktLNiDz99YTGH9+3MVV8YEu9wRJqtqZLCDwnGGLqvkeecFgx1IZJsHnhjOdt2VfHY5UeRkX7Ao8eItBtN9T66Orw/se3CEUk8y7eU8Zf/rOXrRw9iVD+NbSSJLZqhsxeY2S1mdnBbBCSSSNydO/6xhE7ZGfzwi4fEOxyRFoumnHs2wcVq08xslpndYGbqaycCTF9axDsrtvL9U4ZrKAtJCtHMvLbW3e919/HAJcAYYE3MIxNp56pqavnly0sY1rMTl04aFO9wRFpFVLOGm9lg4CLgYoJSw42xC0kkMUz9dwEF28qZ+s2JZKpxWZJENPMpfABkAn8HLnT31TGPSqSdKy6r4qHpKzl5RE+OP0STPUnyiKakMNndP455JCIJ5L7Xl1FVU8uPzxgZ71BEWlVTQ2df6u5PAKeb2el7P+/u98c0MpF2atGGnTw9ez1XHjuEoT06xTsckVbVVEmh7jr9vCa2EUkp7s4vXlpCfscsrjt5eLzDEWl1TV289ofw/va2C0ekfXt54SY+LNjOr84bTZeczHiHI9Lqmqo+eqipF7Zk5jWRRFRdE+He15YxonceFx81IN7hiMREU/3o5qCZ10TqPT17Peu2l3PTqSNI17DYkqQ085pIFMqra3ho+gomDs7nhEPVBVWSl2ZeE4nCY+8VUFxWxY2nHqo5lyWpaeY1kf0oKa/m92+t4pSRPZkwOD/e4YjElGZeE9mPh99axa6qGm748qHxDkUk5qIdsCUdKAZ2AIeY2XGxC0mk/di8s5LH3yvgvHH9GNG7c7zDEYm5aMY+uodgILzFQCRc7cDbMYxLpF14cPoKIu78QHMlSIqIpk3hXOBQd6+KdTAi7cnq4l1Mm72eyyYNYkB+x3iHI9Imoqk+Wk0wSqpISrnv9eVkZ6RxzYnD4h2KSJuJpqRQDsw3s+lAfWlBVzRLMvuosISXF27ieycPp0dedrzDEWkz0SSFF8ObSMr49T+X0a1jJld9YUi8QxFpU9F0SZ3aFoGItBezCrbzzoqt3Hr6CPI6qOZUUktTA+JNc/eLzGwhQW+jT3H3MTGNTCROHnhjOd07ZXPZpMHxDkWkzTVVUrg+vD+ztQ9qZl2BPwGjCBLON4FlwNPAYKAAuMjdd7T2sUWa8p/V2/j3qm385IyR5GSlxzsckTa3z95H7r4pvF/b2K2Fx30QeM3dRwBjgaXAzcB0dx8OTA+XRdqMu3P/G8vpkZfNpZMGxTsckbiI9ormVmNmnYHjgEcA3L3a3UuAc4C69oupBNdHiLSZ91dt48M12/l/JxxMh0yVEiQ1tXlSAIYSDJnxmJnNM7M/mVku0KtB6WQT0LOxF5vZ1WY228xmFxcXt13UktTcnQfeXE6vztl8beLAeIcjEjfxSAoZBJP2POzuRwC7OYCqInef4u4T3H1Cjx4a115ax7srtzKrYAfXnjhMpQRJaU31Pmq011GdFvQ+KgQK3f2DcPkZgqSwxcz6uPsmM+sDFDVz/yIHpK4toW+XDlykaTYlxTXV+6iu19E14f1fwvuvE1zl3CzuvtnM1pvZoe6+DDgZWBLeJhPM3zAZeKG5xxA5EG8tL2beuhLuPG8U2RkqJUhqa2o6zrUAZnasux/b4Kmbzew94BctOO51wF/NLItgbKUrCKqyppnZlcA64MIW7F8kKu7OA28sp1/XHC4cr1KCSDTDXOSa2efd/V0AMzsGyG3JQd19PjChkadObsl+RQ7Uvz4uYkHhTu4+fzRZGfFoYhNpX6JJClcCj5pZF4I2hp0EF5uJJDR357dvrmBAfg4XjO8f73BE2oVoxj6aA4wNry8wd98Z+7BEYm/60iIWbtjJvReMITNdpQQRiKJLqpn1MrNHgKfdfaeZHRbW+4skLHfnt9OXMyA/h/OO7BfvcETajWh+Hj0O/BPoGy4vB74fq4BE2sK/Pi5i0YZSrj1xmEoJIg1E82no7u7TCOdndvcaoDamUYnEkLvz4PQV9O+Ww/lHqi1BpKFoksJuMzuI8EI2M5tE0NgskpBmLCvio8KdKiWINCKa3kc/JJh57eDw+oQe6BoCSVDuzoNvqpQgsi/RJIXFwPHAoYARzHugn1eSkGYuK2ZB4U7u0nUJIo2K5lPxvrvXuPtid1/k7nuA92MdmEhrC3ocraBf1xwuUClBpFFNDYjXG+gH5JjZEQSlBIDOQMc2iE2kVb21vJgF60v41XkqJYjsS1PVR18GLgf6A/c3WF8K3BrDmERaXd3Vy/265vAVXb0ssk9NDYg3FZhqZhe4+7NtGJNIq3t7xVbmrw9GQlUpQWTfovl0jDezrnULZtbNzH4Zw5hEWlXQ4yiYL0EjoYo0LZqkcFo4hzIA7r4DOD12IYm0rndWbGXuuhK+e+IwlRJE9iOaT0i6mWXXLZhZDpDdxPYi7Ubd3Mt9u3TgoglqSxDZn2iuU3gCmG5mjxFc1fxNYGpMoxJpJTM1q5rIAYlm6Ox7zewj4BSCbql3uPs/Yx6ZSAu5O7/VrGoiBySakgLAUqDG3d80s45mlufuZbEMTKSlZizTrGoiByqa+RSuAp4B/hCu6gc8H8ugRFpKs6qJNE80P5+uAY4luGgNd18B9IxlUCItNX1pMBLqdScO10ioIgcgmk9LlbtX1y2YWQbhMNoi7VFdj6OB+R01q5rIAYomKbxlZrcSjIH0ReDvwEuxDUuk+V5fsoXFG0v53skqJYgcqGg+MTcDxcBC4NvAK8BPYhmUSHNFIkFbwpDuuZw7ru/+XyAin9LUKKnT3f1k4C53vwn4Y9uFJdI8ry/ZzNJNpdx/0VgyVEoQOWBNdUntY2bHA2eb2VN8MnQ2AO4+N6aRiRygulLC0O65nD1WpQSR5mgqKfyMoOpo76GzIWhoPilWQYk0x6uLNvPx5jIe/Oo4lRJEmqmpobOfAZ4xs5+6+x1tGJPIAaupjXDfG8sY3rMTZ45RKUGkuaIZ5uIOM+sHDGq4vbu/HcvARA7Ec3M3sLp4N7+/dDzpabb/F4hIo/abFMzsbuCrwBKgNlztgJKCtAuVe2r57ZvLGTugK18+vFe8wxFJaNGMfXQecKi7V8U6GJHm+OsH69i4s5JfXzgWM5USRFoimta41UBmrAMRaY5dVTX8z4yVHDvsII4d1j3e4YgkvGhKCuXAfDObDtSXFtz9ey05sJmlA7OBDe5+ppkNAZ4C8oG5wGUNh9cQacwj76xh++5q/uvLI+IdikhSiKak8CJwB/BvYE6DW0tdTzAkd517gAfcfTiwA7iyFY4hSWz77mr++M5qvnx4L8YN6Lr/F4jIfkXT+6jVZ1kzs/7AGcCdwA8tqAg+Cbgk3GQq8HPg4dY+tiSPh2eupLy6hhu+dGi8QxFJGk0NczHN3S8ys4U0Miqqu49pwXF/C9wI5IXLBwEl7l4TLhcSzNvQWFxXA1cDDBw4sAUhSCLbtLOCqe+v5bwj+jO8V97+XyAiUWmqpHB9eH9max7QzM4Eitx9jpmdULe6kU0bHZ7b3acAUwAmTJigIbxT1EPTV+DufP+U4fEORSSpNHVF86bwfm0rH/NYgvGUTgc6AJ0JSg5dzSwjLC30Bza28nElSawu3sW02YVcevRABuR3jHc4IkmlzQeIcfdb3L2/uw8muCjuX+7+dWAG8JVws8nAC20dmySG+15fTlZ6GteepFKCSGtrT6OG3UTQ6LySoI3hkTjHI+3Q7ILtvLxwE1cdN5QeednxDkck6ewzKYTXJWBm98Tq4O4+093PDB+vdveJ7j7M3S/UFdSyt0jEueMfS+jVOZvvHD803uGIJCXNpyAJ44UFG1hQuJPfXDiWjlnRXHcpIgdK8ylIQqioruXe15Yxul8Xzj+i0d7KItIKNJ+CJIQ/vrOaTTsrefCrR5CmobFFYiba+RTOBo4LV81093/ENiyRT2wpreThmas4bVRvJg7Jj3c4Ikltv72PzOwuggvZloS368N1Im3i1/9cRm3Eufk0DXonEmvRtNadAYxz9wiAmU0F5gG3xDIwEYBFG3by7NxCrvrCUAYdlBvvcESSXrTXKTQcgrJLLAIR2Zt70AW1W8csrj1pWLzDEUkJ0ZQU7gLmmdkMgm6px6FSgrSBfy7ewgdrtnPHuaPo3EHzPIm0hWgamp80s5nAUQRJ4SZ33xzrwCS1Ve6p5VevLGV4z0587agB8Q5HJGVEdQVQODjeizGORaTe7/61knXby/nbt44mI709jcYiktz0aZN2Z2VRGX94exXnH9GPYzTvskibUlKQdiUScW59bhEdszK49YyR8Q5HJOU0mRTMLM3MFrVVMCLPzCnkw4Lt3HLaCLp30iioIm2tyaQQXpuwwMw076XE3LZdVfzq1aVMGNSNiyaocVkkHqJpaO4DLDazD4HddSvd/eyYRSUp6VevfMyuyhp+df5ojW8kEifRJIXbYx6FpLz3V23j2bmFfPeEgzmkV168wxFJWdFcp/CWmQ0Chrv7m2bWEUiPfWiSKqpqavnx8wsZkJ/D9zTFpkhcRTMg3lXAM8AfwlX9gOdjGZSklj+8tZrVxbv5xTmjyMnS7w2ReIqmS+o1wLFAKYC7rwB6xjIoSR0rtpTxuxkrOWNMH048VP9WIvEWTVKocvfqugUzyyCYeU2kRaprIlz/1Hw6ZWdw21mHxTscESG6pPCWmd0K5JjZF4G/Ay/FNixJBQ+8uZwlm0q5+/zR9MzrEO9wRIToksLNQDGwEPg28Arwk1gGJclvVsF2fv/WKi6eMIAvHd473uGISCia3keRcGKdDwiqjZa5u6qPpNnKKvfwg6fnM6BbR36qaiORdmW/ScHMzgB+D6wiGDp7iJl9291fjXVwkpxuf2kJG0sq+Pt3Pken7KgG6hWRNhLNJ/I+4ER3XwlgZgcDLwNKCnLAXlu0iWfmFHLticMYPyg/3uGIyF6iaVMoqksIodVAUYzikSRWVFrJLc8tZHS/Llx/ii5SE2mP9llSMLPzw4eLzewVYBpBm8KFwKw2iE2SiLtz47MfUV5dywMXjyNTE+eItEtNVR+d1eDxFuD48HEx0C1mEUlS+uM7q5m5rJjbzz6cYT07xTscEdmHfSYFd7+iLQOR5PXuiq3c/erHnDaqN9/43KB4hyMiTYim99EQ4DpgcMPtNXS2RGP99nKufXIuw3p24jcXjsVMQ2KLtGfR9D56HniE4CrmSGzDkWRSUV3Lt/8yh0jEmXLZBHLV/VSk3YvmU1rp7g+11gHNbADwZ6A3QZKZ4u4Pmlk+8DRBiaQAuMjdd7TWcaVtuTs3P/cRSzeX8ujkoxjcPTfeIYlIFKLpAvKgmd1mZp8zsyPrbi04Zg3wI3cfCUwCrjGzwwiG01Q/Q50AAAyOSURBVJju7sOB6eGyJKhH3l3DC/M38qMvHsKJIzT6qUiiiKakMBq4DDiJT6qPPFw+YO6+CdgUPi4zs6UEczScA5wQbjYVmAnc1JxjSHz9e+VW7nr1Y049vDfXnDgs3uGIyAGIJimcBwxtOHx2azGzwcARBOMq9QoTBu6+ycwa/XlpZlcDVwMMHDiwtUOSFircUc61T85jSPdcfnORGpZFEk001UcLgK6tfWAz6wQ8C3zf3UujfZ27T3H3Ce4+oUePHq0dlrTAtl1VfOPRD9lTG2HKZeM1rpFIAormU9sL+NjMZgFVdStb0iXVzDIJEsJf3f25cPUWM+sTlhL6oKE0EkpZ5R4uf2wWG3ZU8Jcrj2ZoD12gJpKIokkKt7XmAS2oT3gEWOru9zd46kVgMnB3eP9Cax5XYqdyTy1X/Xk2SzeVMuUb45k4RAPdiSSqaOZTeKuVj3ksQcP1QjObH667lSAZTDOzK4F1BGMsSTtXUxvh2r/N44M12/ntxeM4aUSveIckIi0QzRXNZXwyJ3MWkAnsdvfOzTmgu79LMC9DY05uzj4lPiKRYJC7N5du4RfnHM454/rFOyQRaaFoSgp5DZfN7FxgYswikoTg7tzx8hKem7uBH37xEL7xucHxDklEWsEBj1/s7s/TzGsUJDm4Ow+8sZzH3ivgimMHc91JuhZBJFlEU310foPFNGACn1QnSYqJRJzbX1rM1PfXcuH4/vz0jMN0LYJIEomm91HDeRVqCMYlOicm0Ui7Vl0T4Ya/L+DFBRu56gtDuOW0kaSlKSGIJJNo2hQ0r4JQXl3Dd56Yy9vLi7n5tBF8+7ihKiGIJKGmpuP8WROvc3e/IwbxSDu0Y3c1Vzw+i48KS7j3gjFcdNSAeIckIjHSVElhdyPrcoErgYMAJYUUsLGkgm88+iHrtpfz8KXj+fLhveMdkojEUFPTcd5X99jM8oDrgSuAp4D79vU6SR5z1u7gmr/OZXdVDX/+5kQmDT0o3iGJSIw12SXVzPLN7JfARwQJ5Eh3v8ndNS5REnN3HntvDRf/4X2yMtJ4+tufU0IQSRFNtSn8GjgfmAKMdvddbRaVxM2uqhpuevYjXv5oE6eM7MV9F42lS05mvMMSkTbSVJvCjwhGRf0J8OMGPU2MoKG5WcNcSPu1fEsZ33liDgVbd6uHkUiKaqpN4YCvdpbE9fy8Ddzy3EJyszP421WTVF0kkqI0C0qKKyqt5LYXF/Pqos1MHJzP7y45gp6dO8Q7LBGJEyWFFOXuPD1rPXe+spSqmgg3nnooV31hKJnpKiCKpDIlhRS0ungXtzy3kA/WbGfS0HzuOn8MQ7rnxjssEWkHlBRSSHVNhD++s5oHp68gOyONey4YzUUTBqgxWUTqKSmkgNqI88L8DTzw5nLWb6/g9NG9+flZh6vtQEQ+Q0khibk7by4t4jf/XMayLWUc3rczU785muMP6RHv0ESknVJSSFL/Wb2Ne1/7mLnrShjSPZffXXIEp4/qo6GuRaRJSgpJpDbivLFkC4++t4YP12ynV+ds7jp/NF8Z31+9ikQkKkoKSWBnxR6mzVrP1PcLKNxRQb+uOfzkjJFcOmkQHTLT4x2eiCQQJYUEtmJLGX/5z1qemVNIeXUtEwfn85MzRnLKyF5kqGQgIs2gpJBgNu+s5KUFG3l+/gYWbywlKz2Ns8b25YpjBzOqX5d4hyciCU5JIQHsrNjDa4s28fy8jfxnzTbcYWz/LvzszMM4e1xfunfKjneIIpIklBTaIXdnRdEuZnxcxIxlRcwu2EFNxBnSPZfrTx7OOeP66QpkEYkJJYV2Ymf5Hj4s2M7MZUXMXFbMhpIKAEb0zuOq44Zy6uG9GdO/i64+FpGYUlKIA3encEcFs9duZ1bBDmYXbGf5lmAOo9ysdD4/vDvXnTSM4w/tQZ8uOXGOVkRSiZJCjLk7G0oqWLKxlMUbS1myqZSFhTvZXFoJQF52BkcO6sbZY/syYXA+Rw7sRlaGeg6JSHwoKbQSd2fTzkpWF+9m9dZdrC7ezfItZSzZVEpJ+R4AzGBo91wmDslnwuBuTBiUz6G980jXVcYi0k4oKRyAXVU1bNhRwYaScgp3VLBhRwWFOyoo2LabNVt3U15dW79tblY6w3p24rRRfTisb2cO79uZEb3z6Jilt1xE2q+U/4baUxthR3k1O3bvYdvuKrbvrqaotIqisiqKyirDx5VsKa1iZ8WeT702KyON/l1zGJDfkaOHHMTQHrkM7ZHLwT060TMvW43CIpJw2lVSMLNTgQeBdOBP7n53LI7z9Kx1PDxzFdt2V1NWWdPoNpnpRs+8DvTsnM2Q7rkcPeQg+nbNoX+3HPp1C+6752ZrgDkRSSrtJimYWTrwP8AXgUJglpm96O5LWvtYB+VmM6Z/V/Jzs+jWMYv8Tlnkd8yiW24mB+Vm0zMvm64dM/VLX0RSTrtJCsBEYKW7rwYws6eAc4BWTwqnHNaLUw7r1dq7FRFJeO2p72M/YH2D5cJw3aeY2dVmNtvMZhcXF7dZcCIiqaA9JYXG6mr8Myvcp7j7BHef0KOHZhATEWlN7SkpFAIDGiz3BzbGKRYRkZTUnpLCLGC4mQ0xsyzgq8CLcY5JRCSltJuGZnevMbNrgX8SdEl91N0XxzksEZGU0m6SAoC7vwK8Eu84RERSVXuqPhIRkThTUhARkXrm/plenwnDzIqBtc18eXdgayuGkyhS9bwhdc9d551aojnvQe7eaJ/+hE4KLWFms919QrzjaGupet6Quueu804tLT1vVR+JiEg9JQUREamXyklhSrwDiJNUPW9I3XPXeaeWFp13yrYpiIjIZ6VySUFERPaipCAiIvVSMimY2almtszMVprZzfGOJ1bM7FEzKzKzRQ3W5ZvZG2a2IrzvFs8YY8HMBpjZDDNbamaLzez6cH1Sn7uZdTCzD81sQXjet4frh5jZB+F5Px0OOJl0zCzdzOaZ2T/C5aQ/bzMrMLOFZjbfzGaH61r0f55ySaHBtJ+nAYcBXzOzw+IbVcw8Dpy617qbgenuPhyYHi4nmxrgR+4+EpgEXBP+jZP93KuAk9x9LDAOONXMJgH3AA+E570DuDKOMcbS9cDSBsupct4nuvu4BtcmtOj/POWSAg2m/XT3aqBu2s+k4+5vA9v3Wn0OMDV8PBU4t02DagPuvsnd54aPywi+KPqR5OfugV3hYmZ4c+Ak4JlwfdKdN4CZ9QfOAP4ULhspcN770KL/81RMClFN+5nEern7Jgi+PIGecY4npsxsMHAE8AEpcO5hFcp8oAh4A1gFlLh7TbhJsv6//xa4EYiEyweRGuftwOtmNsfMrg7Xtej/vF0Nnd1Gopr2UxKfmXUCngW+7+6lwY/H5ObutcA4M+sK/B8wsrHN2jaq2DKzM4Eid59jZifUrW5k06Q679Cx7r7RzHoCb5jZxy3dYSqWFFJ92s8tZtYHILwvinM8MWFmmQQJ4a/u/ly4OiXOHcDdS4CZBG0qXc2s7gdgMv6/HwucbWYFBNXBJxGUHJL9vHH3jeF9EcGPgIm08P88FZNCqk/7+SIwOXw8GXghjrHERFif/Aiw1N3vb/BUUp+7mfUISwiYWQ5wCkF7ygzgK+FmSXfe7n6Lu/d398EEn+d/ufvXSfLzNrNcM8urewx8CVhEC//PU/KKZjM7neCXRN20n3fGOaSYMLMngRMIhtLdAtwGPA9MAwYC64AL3X3vxuiEZmafB94BFvJJHfOtBO0KSXvuZjaGoGExneAH3zR3/4WZDSX4BZ0PzAMudfeq+EUaO2H10Q3ufmayn3d4fv8XLmYAf3P3O83sIFrwf56SSUFERBqXitVHIiKyD0oKIiJST0lBRETqKSmIiEg9JQUREamnpCASJTPramb/L3zc18ye2d9rRBKNuqSKRCkcR+kf7j4qzqGIxEwqjn0k0lx3AweHA86tAEa6+ygzu5xgJMp0YBRwH5AFXEYwnPXp7r7dzA4mGLa9B1AOXOXuLR6rRqQ1qfpIJHo3A6vcfRzwX3s9Nwq4hGDsmTuBcnc/Angf+Ea4zRTgOncfD9wA/G+bRC1yAFRSEGkdM8K5G8rMbCfwUrh+ITAmHLH1GODvDUZrzW77MEWapqQg0joajqkTabAcIficpRGM7z+urQMTORCqPhKJXhmQ15wXunspsMbMLoRgJFczG9uawYm0BiUFkSi5+zbgPTNbBPy6Gbv4OnClmS0AFpOk08BKYlOXVBERqaeSgoiI1FNSEBGRekoKIiJST0lBRETqKSmIiEg9JQUREamnpCAiIvX+P8aqYyqEPNKQAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(range(0,50),ws)\n",
    "plt.ylabel('Number of infected individuals')\n",
    "plt.xlabel('time')\n",
    "plt.title('Diffusion of infection');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As Scott Page explains in his lecture, what we see here is simply an accelarting process, not a tipping point. This accelaration matches the expectation we developed earlier by studying algebracially the properties of $P_{enc}$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The SIS model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We now turn to the SIS model (Susceptible-Infected-Susceptible). The difference is that now infected individuals can become cured, and reutrn to the pool of susceptible individuals (thus, being sick does not make you immune to subsequent reinfection). The difference consists solely of adding a term to $W_{t+1}$ measuring the numer of cured inviduals. Let $a$ be the rate of being cured."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "a = sm.symbols('a',negative=False, real=True)\n",
    "Wtt_sis = sm.symbols('Wtt_sis',cls=sm.Function)\n",
    "Wtt_sis = Wtt - a*Wt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAARwAAAAsBAMAAABF1xrKAAAAMFBMVEX///8AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAv3aB7AAAAD3RSTlMAEM0ydt27RKvvVGaJIpmJty9HAAAACXBIWXMAAA7EAAAOxAGVKw4bAAAEHElEQVRYCc1YS2gTURQ9E/Np00kaxe/KOCK4a92JIAb8gIgY0IVVS4IIIoLJTihiI+hCVyIogogFPzsxC90I2qkKLgTbhQtRtIOCuBCMn9ZPrfHeN5l03rNpopmBuYs3973z7n3nffJOZgD/7WVrQ9xsrVsrvTZXD+L0xAlcnyj/1V3LA+HCN+D6jxMSGM5NFjuu7E7Fc3u5vWNMQtupxH4BOo04r0RZ+uVMIR4mOknFChnAoAmMFoG11N4PLU0Pbyzx26ajU7p4Rc55n6sDBSJlygDGLZsOtXPQIwX+/6r2HQhNAUcpRTIt5xGLNTJoITYkA1heAnrKiNECcdAyBW6j+hnY9tWef8SU8iSyXLUiFegpCQC60wgPZ8BLykEhouWRnUfi6i+Ex4DjhSM0wqaN2FqtVjN0ii0ewuqaQpIdt3Vm8WKDie1OkOUG2/L7oNP5Wcg5bgGJPdivP3y+/hVVo2UqaJ8mxCpwh7pFzsIct7QMNVAQutJUeGM9uKpN2mf1jVj6M+sQFakjQ/SgfRot8sGSLFSJ9y7P816BghDPSmg7lUIqg898JoF9wMcUO/z7BTq5kfZpPG9yFQsMtpXsRj+F6fwIlhSEeRVu9MQ+PgPOi3lqdP1cEDmfiFLQISQ6bakjJacz6Mya1MxBXtIZvQv03ePxYhXMp0uRbKcoI0P0GKD23yY5kulTQ4hcGqM2CkohlpXQdio78nSFZDiDno730g0UR8Jee3GURwjIFRl2W5iu6hBdD3aQl0d5mOZ4SgyVLC3Ee2hvEbYnq1vUzO6wgN1F1zQplejFQUiabrAt/xhFPxYZEofL0J/cpuwZUact0A6xsnaIqruImbQuRW7hIIR6Z8C5NFn0UvR3JrKZ19+sg4Mvchx6zqXJdjdZf12hTVwhoU36CNgtoY01+Vot1bgl6S+3xjNcNrHoWJMONVhLu/o11mSHjqK/HNsSHS3vGmYOV6bdUJMdOor+cuKW6ODlHBxckLi16nVFk18/3ZqyMYeOor8MtkanPsS/OLImxy6jpxbt0FH0l1Ef6ciavKyMdwodRX/9oVMX2L80eRcP+MAw+gxjNbuK/kI3jFXnDKPEmPemarL2pTaGs1mK/jLqbFY3/fHzxGampWpy/YXAoaPor5vOTBbPPEmT1/xER7r2WujQUfTXZzqyJh/AnZJpT9Who+ivz3RkTV66af3JokxH0V+f6dhjz1I6qzML5BzlWSDfmuY3zqzVLu7GPRxkS7WIeO5iy/2dOL+ehTz9vwoMm9gN+iOq+zXXf86rL6FvHPQ+ERBLIgfw+0QwLInhMqxgcCEWA+jOwgwMnRH6Vme/wgeCkoXEZHB+WLxPhdeBWBgmQZ+bMPghMHT4s2CEP5oEwyJEo6sSDC7A4it7iUqpOZ0/MN50J9LLUkMAAAAASUVORK5CYII=\n",
      "text/latex": [
       "$$- Wt a + Wt + \\frac{Wt c \\tau \\left(N - Wt\\right)}{N}$$"
      ],
      "text/plain": [
       "             Wt⋅c⋅τ⋅(N - Wt)\n",
       "-Wt⋅a + Wt + ───────────────\n",
       "                    N       "
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Wtt_sis"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We are just interested in the change, so let's simplify things a bit:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAOYAAAAyBAMAAAC667qxAAAAMFBMVEX///8AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAv3aB7AAAAD3RSTlMAdt27RKvvEFRmiTIimc3h244VAAAACXBIWXMAAA7EAAAOxAGVKw4bAAAEmUlEQVRYCa1YTYgcRRT+pmdMz/TO7I4/CBJkF8ElRiRD8KKIaUQkeMhuPBgP/oQ1QYygIwl6c9YIXvYw601B3DkIIYrJXgQPwrZKbko2SDxEcPegJ81kkjVEo2Z89Wpmp6vq9Wxn0g+2+r3v/dVWV7+qN4BIQSjCtwr6s0MifDZEZ6nKixYgi8VI4ftkpUJzw+ZjuR0l+c5uiOKuw3VDdUf3YTx75QCOX+FZPa+UpXXDJC68Hxe24P9Q+j01YNxMicJ1IPgTGJsFvAjvcJxHk6L5h5I0Ll5aJazw3r8U31KW/9M5Fd4ktqr0b4VqFCiv4qQk+g8o3jGKv9fy8K/S//438CThvwD+QaUv8qg4i5p1Cxgifqt0FewCTthWN4B7LgMtwnfS36+sV5xED0lgHAtu3467u90uLde9Cq9gaRHO4syhvHAduQiFHzfOAV9whGbED3so0csZSuUdeDz4/MzX35PVI8pyLyam+P8x3KYR0Dt9k7Bih4ZPWOk5U9NwjR/Jg9fCc18hzwY71HiCNmghZDk2XMKCv8FTqawTfJxVpWsxiwG7Eg54kVurK/hL1h1S4yrKG2p70k5pMy0rfk+9ihs8FZoksKIw+PSGBfpJwAzoZZZO8sg5W5TgB8NECWvfAHM8lUZIos6JtlI5dNhBTMCnj53oQR7V2hbmgZkLxJh0+kNg+mOFzcAHntDayUg/jXGMvqnhRAZFlDtspPZQQGvtUcWxaLIGXKoq8BQW+nsIZxctMyXmttq2OA//Z+SmyFZ/AR49Sx0W48NSBDzDwHmfUvcWYtts3KbH5w8KoAEFJz+gb5LCENGKHX3jAWLcUE8T+qmywXcf0fA7s8hPwar+hG+raWW60YvS2cGf0oZBx67+hDdWU0ZhM67xaRx0jaeX8A+s6k/eK3rR3DjvuhAhj4moC+qzjI62y7CqP9k2510HRuSc6sxOQ/25+XQEWdUfWFpOCCHnTHk3KUX9qHTEWdUfmNzUslX5wrlXtLmcsx8q9XODKoVZ/YHTdcP9qWi8k2nO16lSmNWfyqSR038VlVqmOeec6g+8xBn6B0TlGrwqIafa7el2+zelm6DTejTiyNNO9e/nZDUNVCEay1rI6H1STqv622vbaKlLG1NGOe+jbUpvK1b9aRJRLwc/GvN4saiRjHLSHrKqP02it5Q6c75aul9fAZBRTjqerepPkwh1Nj36+7efuUuzGeWkOuRQM3QgDWydU+xbrGi9W4aJriya8qb09iaXyEh9i2VMNd6liVUXS4mIfYvlW/zLApTorQtgOkjsWyzXQLr6VDqWVXpR7ltM/9uk8Ak3bdNTluS+xbQVrz5j0oKbfkmS3LeY1merpqylre7Uko/G5L7FtLeO555yzTS6CWnQtyQ77RZVM6EIpwBbYt9iOkqfJ30so36ghXmxbzFy5uRecOSNm9C3GDnF1oEsjhhW6QWPTIW+xQiQ0NvjZn7DiAVM6ltiJsBFQxoI46O+0EGIJK7Yu7c6+vJrDpQVkPibVO/nhqzyxOO8EBcMvjJriNkJpVpyrP3JqlvS7BviXQmHKEdX+QeU7/+5r07pCVgW3wAAAABJRU5ErkJggg==\n",
      "text/latex": [
       "$$Wt \\left(- a + \\frac{c \\tau \\left(N - Wt\\right)}{N}\\right)$$"
      ],
      "text/plain": [
       "   ⎛     c⋅τ⋅(N - Wt)⎞\n",
       "Wt⋅⎜-a + ────────────⎟\n",
       "   ⎝          N      ⎠"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "change=sm.collect(Wtt_sis -Wt,Wt)\n",
    "change"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now assume $Wt<<N$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAHgAAAAUBAMAAAC0fOTAAAAAMFBMVEX///8AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAv3aB7AAAAD3RSTlMAdt27RKvvEFRmiTIimc3h244VAAAACXBIWXMAAA7EAAAOxAGVKw4bAAAB7ElEQVQ4EY2UP2gUQRTGf5dbnXPcyMYiIBaeTRC1iGAjqDmClSCcFmphIaa1OAu1OzRtBE8rC/+ksAkhiY1gEXAtLASFgNgo6BVaBdYzKuiJOd+bmUCOXfBe8e37vnnfzryZ2YV8vMpLBUolhZ29w5xdO8f0WitUmPGC0gLpAkRdsD9ge13GD2rNjlTx/9GE+K83W6mudNTyoNB3M6faBPNT5voNJ2RwuK0VbnpN+iJvNhdhHXZ9g1kpHVKIqwL5yJv5DJeIZ7qUUjg1cSyB0som55uF0YanG2Y7spvRXq8nlQswhpW+r2jJbYVyS9FHdJmvIQ3meB/H7fzL5dciP0RGZ8x3t2o+auHQDUUf11q8DWkwS2fnn1F24jRMNBLWo5rS/Qpba9DMNFb51GCvaotZNpZlXyQTReKpAs+FvpC+9aAwctzerE+NKcwfnxFmnnJ0zqGYl+5J3/eVRR3kvZuWbbrh6GXQm0XROOBwEvaMS9+JMtuuPOnbMPOLLe1brnBjZrkRFeKO02TDHqVwxpHhuu65XXHEwSHu1mc9Dct+h/lAqeq093BSkseOxEda8oyqjji4OrJ8uuZpMNu5O3ITE6et+pE+HPh6VvtsnhR/GNdzlfJh5KOc5rUipVkkDvozOFpkZrDf0LaUfzfTcv7B3goAAAAAAElFTkSuQmCC\n",
      "text/latex": [
       "$$Wt \\left(- a + c \\tau\\right)$$"
      ],
      "text/plain": [
       "Wt⋅(-a + c⋅τ)"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "change.subs([(N-Wt,N)])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "so the number of infections will increase (i.e., the infection will spread) if $c\\tau-a>0$, or $c\\tau/a>1$. This number $R_0=c\\tau/a$ is called the **Basic reproduction number**. If it is bigger than 1, the disease spreads. This (as opposed to that we saw in the diffusion model) is a true *tipping point*. $R_0$ values for several well-known diseases (via Wikipedia): Measles - 12-18; Mumps - 4-7; HIV - 2-5; Flu - 2-3; Ebola - 1.5-2.5. Ouch!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Vaccines"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Assuming a vaccine is 100% effective, we should check $r_0=R_0(1-v)$, where $v$ is the proportion of the population that is vaccinated (making $1-v$ the proportion of individuals who are not vaccinated) instead of $R_0$. If you are not sure why, review the expression used to derive $R_0$. For vaccination to be effective we need $R_0(1-v)<1$ "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAIYAAAAUBAMAAAC6+VjqAAAAMFBMVEX///8AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAv3aB7AAAAD3RSTlMAZqu7Iu8ydt1UmRDNiUQDe6IrAAAACXBIWXMAAA7EAAAOxAGVKw4bAAABzUlEQVQ4EYWUvUvDQBjGH9vatE1TP3AQcYiCm4Ii4loHByfFxcGl/4HFyc0MujuILkIrCDqKIAo6FAengOJSQcEu7n6hoIK+713OXOJpD3J57ve+91zuK4BWJjRtlOtGSrDzq390hIP28F8piieKSsXfz8CkSzBfjEf0dnaeRqnpRGjGQPoNSLyS6BEwXm1KcHD1QaKuRbdK1JAY1guQoQcrWkIoAw+0sMdhyLt3hBYYyRmgQkuR8gSMVxGPfEOFl26kkh5tLrKzBKyqigNHu8i0y2bEwwly7OtykCw9FvoW30tEkkEnDl7eI+OyACIeLTXB7L2fXOmxvI/efQpl1kScq9TaHrKXshnxyHoCDoap0mMMaOP81gZw6nM5tzGHAufd+v6q759xR5GcfmBp//J4AqZdirCHKoUPWIGOfIf0QCo2l/QjbUuVOmhzgePBMXkUPElja2rR8bor29E1TbjYMHkEa0qhyN7maYqV8hSg9o37OjWbl4iLmovFZyynKBCcMca5i6d2OONlQH2n6DnUIV6hR27gc5uuVElhQJx1iUP4/1kXecdhtlmZ79yJnlzXGyadLJqozgx3Xw+Tbv4Paj4Kmv4Lu2KjyuY3fDF0d0YLbmkAAAAASUVORK5CYII=\n",
      "text/latex": [
       "$$R \\left(- v + 1\\right) < 1$$"
      ],
      "text/plain": [
       "R⋅(-v + 1) < 1"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "v,R = sm.symbols('v,R',negative=False, real=True)\n",
    "R*(1-v)<1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Thus"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAGwAAAAqBAMAAACkSaOPAAAAMFBMVEX///8AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAv3aB7AAAAD3RSTlMAEM0ydpndZolE76siVLs7Lf7LAAAACXBIWXMAAA7EAAAOxAGVKw4bAAABTElEQVRIDWNgoBPgXUiORdLdH8nRxsA+qg093IZzkLCQFd2suz6Fo4cSVfgkpvUrEEtR0/olA0JOgWpDSetlMYR0MWDR5taIQ5dICANbAUQOQxtjmwMOXQwNsxjYNmDXxhgGNQ9TL8+FMAbeBuzaNl3AVC+kBALKjAxLGDhB0k1KSllKStoghdC0zoik7T8UwM3h/MjAAuWg+40HtyMZOBIYOHBoY8ATJNwbGC7j0sbAgDMCOBQYG9C0Iad1nNG9oxyqCxbdqGmd+MQFMwaFrv6/+91zFBEYRxDGwEr/Y2A4BU0XWOWxCzJ/YWDg/oxdDo8oyx8GBjYgJhHwLWBguP+ARE0MDPwbGHhXk6yLwX6791cD0rX5CzAUCZCu7RXQnQ2ka/vNwHCejGj7DgzIAJJtYwHGdL8DI6n6mD4AbXM4S6I2Vv3fBQwcTxwAMCVlvUB5KkIAAAAASUVORK5CYII=\n",
      "text/latex": [
       "$$- v + 1 < \\frac{1}{R}$$"
      ],
      "text/plain": [
       "         1\n",
       "-v + 1 < ─\n",
       "         R"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "1-v<1/R"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAF0AAAAqBAMAAAAjeEqJAAAAMFBMVEX///8AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAv3aB7AAAAD3RSTlMAEHaZ3WaJRO+rIlQyzbupNG4/AAAACXBIWXMAAA7EAAAOxAGVKw4bAAABQElEQVQ4EWNgoBnYIUCK0YzT5pOknoEhflQ93gAebOGjT1p8WfZfDsDrPxIkuRsJK5ZF2CaZ+ZGwesYUBbgiViLUMzCoJMI04FHPqKzA4ApVF2IIZeBRz8bewOAJM3fvLAgLj/rVzA4M22DqGWongpl41AtwKDBUI9Q/RFLPeBcMDgCF+P9DwAcGhnoBBgWwIiBBhHsYUhgYEqDqifEvw0wG7gMQ9USFJ0MugxhYOXJ8MeGJL85USEJASg/MM/6aQp1INsU6/+e8d4gkRtgc/gcMXJ2ElcFVnFdgYJh/AM4lyNBfwMDgD8TEghcMDNz/iFUMVPeTgWHHBOLV8/wMmm9DvHIGYBSyXiRBPV8DA8M3EtTzT2DgAnqBaHA+gIH7M9GqGRj0DzCwfWIQIFrHfWDwf+LaQKx62/9TGRjylBgAd/9RsTO9BHgAAAAASUVORK5CYII=\n",
      "text/latex": [
       "$$v > 1 - \\frac{1}{R}$$"
      ],
      "text/plain": [
       "        1\n",
       "v > 1 - ─\n",
       "        R"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "v>1-1/R"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$v$ has to be above this threshold for vaccination to be effective in stopping the spread othe disease. Below that number there is no population level effect (\"herd immunity\"). That's the key result: the beneifts do not increase linearly, with each vaccinated inidividual improving the protection of the population by the same amount. There's a tipping point, which depends on $R_0$. Thus, for measles, we need:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAALQAAAAPBAMAAAC/7vi3AAAAMFBMVEX///8AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAv3aB7AAAAD3RSTlMAiXZmMs1UEN0i77urRJlR0qN3AAAACXBIWXMAAA7EAAAOxAGVKw4bAAABpUlEQVQ4EbWUPS8EURiFz+zHDGs3VqLTbEgUKpEodOsPsImERiNRKlCJgk6lkOg0JAoFkY0/YBIajfgJGz9gfNsQwnnfhPfNRekm89wz59x7MpkvIBrohYy1nguZJn5BoUXT4NaYqcol0jONlS1OyT52Fjlv8gjRkdIzuNhMVS5hT34PUYM7y3sot4DcEXUADKU0DRY7U2NLRHWliNvc2d1AdAtsz1IH6BxltcHFZqpyiVR0p6g8s269idwLUJXqAFGR1QYXm6nKJVKx3kLlnnUlXvU78lu0AmBMqg0WO1NjS1TtVlGQat6N4gOWweoASY3VBhebqcolWrHbQuGO1TjB3DOqUh0gAqsNLjZTlUu04uuGIB4Zf4zrrA6Aeak2WOxMjS1RJY8xlsfIUW4vgdUBkiqrDS42U5VLtAJdTeTl5eMopgtZ9nodIL7Mbo4Hv1Gz+O9E13x/MvkDTMpXyVf9B0opPYNbY6YqlxxxywZW68kTcjPJFM9w+AtKTZoGt8ZMVS6RnqX+c/1p9GRXPOv7OPuB6PRtGAa3xkxVLtEeucp/GZ+YzwijZ5438gAAAABJRU5ErkJggg==\n",
      "text/latex": [
       "$$0.9444444444444444$$"
      ],
      "text/plain": [
       "0.9444444444444444"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "1-(1.0/18)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "or essentially 95% of the population to be vaccinated to get herd immunity. This has profound policy implications.\n",
    "To be provocative: If you, as a health official, know that this number is unatainable (e.g., because too many people decide not to vaccinate their kids for some reason) -- should you recommend stopping the program to mandate vaccinatons on everyone else?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Plotting the SIS curve"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's compare the spread of the infection in the original model and in the SIS model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEWCAYAAACJ0YulAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3dd3yV9fn/8deVhA2y916CFAURJ9ZZ27pHXV8XqC22X9fX1rpqq/1Z29pWrXZocWKHaNW6Z1GwUsCBUECK7B0SRsIm41y/Pz53QsQQDsk5uZOc9/PxuB/n3Pe5c9/XyYFz5bPN3REREQHIijsAERGpO5QURESknJKCiIiUU1IQEZFySgoiIlJOSUFERMopKUidYGYPm9mPK+x/z8zWmtkWM2tvZqPMbEG0f1YN7vOGmY1OTdT7dN+fmdk6M8ut5LWvmtn8fbjWF343KYzxYjN7O1XXk/rJNE5B0s3MlgKdgRKgFPgMeAoY5+6JSs5vBGwCjnD3WdGxicDL7v5AbcWdKmbWE/gc6O3ueTW81pd+N9W8Th9gCdDI3UtqEpM0LCopSG053d1bAb2BXwI3A4/t4dzOQFNgboVjvXfbr096A+trmhAilf1uRFJGSUFqlbsXuvvLwAXAaDMbCmBmT0ZVLPsDZVUpBWb2rpktAvoBr0RVJk3MbKmZfa3sumZ2p5n9JXre1Mz+YmbrzazAzD4ys87Ra5PM7NvR8ywzu93MlplZnpk9ZWato9f6mJmb2WgzWx5V/fxoT+/LzFpHP58fXe/26PpfA94BukWxP1nJzx5nZisr7C81sxvN7D9mVmhmz0Tv6Uu/m+j8wWb2jpltMLP5ZnZ+hWs1M7N7o5gKzewDM2sGvF/hOlvM7EgzG2NmH1T42aOi311h9HhUhdcmmdldZjbFzDab2dtm1qGKj17qCSUFiYW7fwisBL662/HPga9Eu23c/QR37w8sJ5Q2Wrr7zr1cfjTQGugJtAe+C2yv5Lwx0XY8Iem0BH6/2zlHA4OAE4GfmNkBe7jn76J79gOOBS4DLnf3fwInA6uj2MfsJfYy5wPfBPoCBwFjKvvdmFkLQtL5G9AJ+B/gj2ZWdt5vgEOAo4B2wE1AAjimwnVauvvUijc3s3bAa8CDhN/hfcBru7VhXARcHt23MXBjku9N6jAlBYnTasIXVaoVE77IBrh7qbt/4u6bKjnvYuA+d1/s7luAW4ELzSynwjk/dfftUf39LGDY7hcxs2xCyedWd9/s7kuBe4FLa/AeHnT31e6+AXgFGL6H804Dlrr7E+5e4u4zgOeBc80sC7gCuN7dV0W/i38nkVQBTgUWuPufo+s+DfwXOL3COU+4++fuvh14tooYpR5RUpA4dQc2pOG6fwbeAiaY2Woz+1XUQLu7bsCyCvvLgBxCvX2Zir2FthFKE7vrQPhLefdrda9G7PtyXwjtFYdH1WQFZlZASHZdoriaAouqcf/dfzfw5feUbIxSjygpSCzM7FDCF8wHezt3D7YCzSvsdyl74u7F7v5Tdx9CqDY5jVCds7vVhC/VMr0IPaTW7mMs6wilk92vtWofr1MdK4DJ7t6mwtbS3b8XxbUD6F/Jz+2t2+HuvxuovfckMVJSkFplZvuZ2WnABOAv7j67mpeaSajqaWRmI4FzK9zjeDM7MKrW2UT4wi6t5BpPAzeYWV8zawn8HHhmX7tounspofrkbjNrZWa9ge8Df6nWO9s3rwL7m9ml0e+ikZkdamYHRN19HwfuM7NuZpYdNSg3AfIJbQv99nDd16PrXmRmOWZ2ATAkup80YEoKUlteMbPNhL9sf0RouLy8Btf7MeEv4I3ATwkNrWW6AM8REsI8YDKVf0E/Tqhqep/QZ38HcG0147mWUHpZTCj9/C26flq5+2bg68CFhL/uc4F7gCbRKTcCs4GPCFV19wBZ7r4NuBuYElU7HbHbddcTSlg/ANYTGqhPc/d16X5PEi8NXhMRkXIqKYiISDklBRERKaekICIi5ZQURESkXM7eT6m7OnTo4H369Ik7DBGReuWTTz5Z5+4dK3utXieFPn368PHHH8cdhohIvWJmu49WL6fqIxERKaekICIi5ZQURESknJKCiIiUU1IQEZFyaUsKZvZ4tMThnArH2kXLBi6IHttGx83MHjSzhdEShCPSFZeIiOxZOksKTxKWE6zoFmCiuw8EJkb7EJYrHBhtY4GH0hiXiIjsQdrGKbj7+2bWZ7fDZwLHRc/HA5OAm6PjT3mYsnWambUxs67uviZd8Ulmcnd2FCfYWlTCtp2l7CgpZUdxKTuKExSVJCgqLaWoxCkuTVCSSFBc6pQmnJKEU1qaoNQhkXBK3Um4k0g4CYeEO+7h+v07teTM4TVZdE0kPrU9eK1z2Re9u68xs07R8e6EefbLrIyOfSkpmNlYQmmCXr16pTdaqRcSCWft5h2s2LCdNYXbWVO4g7WbdrBhaxHrtxSxcVsRm3YUU7itmC07S0ikebb4rw/prKQg9VZdGdFslRyr9L+uu48DxgGMHDlSi0FkmM07ipm9qpA5qwqZs2oTn6/dzLL129he/MWF1Vo0zqZ9yya0b9mYzvs1Zf/OrdivaQ6tmjaiRZMcWjbJplnjHJo1yqZpoyya5GTTpFEWjbOzaJSdReMcIzsri5wso1F2FtlZFjYzsrIgO8vIsrINsswwA7PK/imL1B+1nRTWllULmVlXIC86vhLoWeG8HoRVpCTDlSacD5ds4F8L8pmyaD2zVxaU/6XfrXVTBnVpxagBHejToQW92zWnW5umdGndjJZN6srfOyL1S23/z3kZGA38Mnp8qcLxa8xsAnA4UKj2hMzl7ny6ooCXPl3F63Nyyd+8k5wsY3jPNlxz/AAO6dOOod32o33LJnu/mIjsk7QlBTN7mtCo3MHMVgJ3EJLBs2Z2JbAcOC86/XXgFGAhsI2ard0r9VRJaYI35uTy6AdLmLWigCY5WZwwuBOnHdSN4wZ1pIX++hdJu3T2PvqfPbx0YiXnOnB1umKRus3deXNOLj9/Yx4rNmynT/vm3HXmVzh7RA9VA4nUMv2Pk1jNW7OJn74yl2mLNzC4SyseuWwkJw7uRFaWGmxF4qCkILFwdx7512J+9eZ8WjbN4a6zhvI/h/YkJ1szr4jESUlBal3BtiJu/Pss/jkvj5OHduEX5xxIm+aN4w5LRFBSkFq2MG8Lox//kLzNO7jz9CGMPqqP+vaL1CFKClJr5udu5uJHpwHGc989imE928QdkojsRklBasWcVYVc+th0Gudk8bfvHEH/ji3jDklEKqGkIGk3P3czFz0yjVZNG/G37xxO7/Yt4g5JRPZASUHSauPWIr791Ec0bZTNM1cdQY+2zeMOSUSqoKQgaVNSmuCap2ewtnAnE5QQROoFJQVJm5+//l+mLFzPr849iBG92sYdjogkQSOFJC3emL2Gx6cs4fJRfTh/ZM+9/4CI1AlKCpJyBduK+PFLcxjafT9uO+WAuMMRkX2g6iNJubtenUfBtmKeuuJwGmnaCpF6Rf9jJaUmf57P8zNW8t1j+zOk235xhyMi+0hJQVJm684SbnthNv07tuCaEwbEHY6IVIOqjyRlHp68iFUF23nuu0fStFF23OGISDWopCApkb95J499sITTDurKyD7t4g5HRKpJSUFS4o+TFrKzJMH3T9o/7lBEpAb2mhTMbJSZtYieX2Jm95lZ7/SHJvXFyo3b+Ou05Zx3SA/6aaI7kXotmZLCQ8A2MxsG3AQsA55Ka1RSrzw4cQEYXHfiwLhDEZEaSiYplLi7A2cCD7j7A0Cr9IYl9cWi/C0898lKLj2iN93aNIs7HBGpoWR6H202s1uBS4BjzCwbaJTesKS++NPkRTTJyeZ/j+sfdygikgLJlBQuAHYCV7p7LtAd+HVao5J6Yf2Wnbw4czXfOqQ77Vs2iTscEUmBvZYUokRwX4X95ahNQYCnP1xOUUmCMUf1jTsUEUmRPSYFM9sMeGUvAe7umsMggxWVJHhq6jKO2b8jAzqpx5FIQ7HHpODuakyWPXpjzhryNu/knnP7xB2KiKRQ0tNcmFknoGnZflSNJBnq8SlL6dehBccO7Bh3KCKSQskMXjvDzBYAS4DJwFLgjTTHJXXYjOUbmbWigDGj+pCVZXGHIyIplEzvo7uAI4DP3b0vcCIwJa1RSZ3212nLadUkh2+N6BF3KCKSYskkhWJ3Xw9kmVmWu78HDE9zXFJHbSsq4Y05azj1oK60aKJJdkUammT+VxeYWUvgfeCvZpYHlKQ3LKmr3pqby7aiUs5RKUGkQUqmpHAmsB24AXgTWAScXpObmtkNZjbXzOaY2dNm1tTM+prZdDNbYGbPmFnjmtxD0uOFGavo2a4ZI3u3jTsUEUmDvSYFd9/q7qXuXuLu4939wag6qVrMrDtwHTDS3YcC2cCFwD3A/e4+ENgIXFnde0h65Bbu4IOF6zj74B5qYBZpoJLpfbTZzDZF2w4zKzWzTTW8bw7QzMxygObAGuAE4Lno9fHAWTW8h6TYizNX4Q7nHNw97lBEJE2SmebiC4PYzOws4LDq3tDdV5nZb4DlhGqpt4FPgAJ3L2urWEmYY+lLzGwsMBagV69e1Q1D9pG78/wnKzmkd1v6dGgRdzgikib7vPKau79I+Ku+WsysLaGdoi/QDWgBnFzZrfZw/3HuPtLdR3bsqIFTtWXu6k0syNvCOSNUShBpyPZaUjCzcyrsZgEj2cMXdpK+Bixx9/zo+i8ARwFtzCwnKi30AFbX4B6SYi/MWEXj7CxOO7Bb3KGISBol0yW1Yk+jEsKI5jNrcM/lwBFm1pxQfXQi8DHwHnAuMAEYDbxUg3tICrk7b85ZwzH7d6B1cy2lIdKQJdOmcHkqb+ju083sOWAGIcl8CowDXgMmmNnPomOPpfK+Un2zVxWyunAHN5y0f9yhiEiaVTV19u+ooprI3a+r7k3d/Q7gjt0OL6YGDdiSPm/NzSU7y/jaAZ3jDkVE0qyqhuaPCb2CmgIjgAXRNhwoTX9oUle8OSeXw/u2o20LjScUaeiqWk9hPICZjQGOd/fiaP9hQjdSyQAL8zazKH8ro4/qE3coIlILkumS2g2oOFahZXRMMsBbc9cC8PUhXWKORERqQzK9j34JfGpm70X7xwJ3pi0iqVPenJPL8J5t6NK66d5PFpF6L5neR0+Y2RvA4dGhW9w9N71hSV2wqmA7s1cVcsvJg+MORURqyR6rj8xscPQ4glBdtCLaukXHpIF7a07I/d/4iqqORDJFVSWF7xPmGLq3ktecGkx1IfXDO5+tZWCnlvTVXEciGaOq3kdjo8fjay8cqSu27Czh42UbuOLovnGHIiK1KJmps2eZ2a1m1r82ApK6YcrCdRSXOsft3ynuUESkFiXTJfUMwmC1Z83sIzO70cw0Z3UDN2l+Pi2b5DCyj1ZYE8kkyay8tszdf+XuhwAXAQcBS9IemcTG3Zk8P49RA9rTKHufZ1cXkXosqf/xZtbHzG4izGA6GLgprVFJrBbkbWF14Q6OG6SqI5FMk8x6CtOBRsDfgfPcfXHao5JYTZqfB8Bxg7SIkUimSWZE82h3/2/aI5E6Y9L8fAZ1bkXX1s3iDkVEallVU2df4u5/AU4xs1N2f93d70trZBKLLTtL+GjpBq4Ypa6oIpmoqpJC2YilVlWcIw3Mv6OuqMeq6kgkI1U1eO1P0eNPay8cidukz/Np0Tibkb3bxR2KiMSgquqjB6v6wZqsvCZ1178W5HPUgA40zlFXVJFMVNX//E/QymsZZcWGbazYsJ2jB3SIOxQRiYlWXpNyUxauA2DUgPYxRyIicdHKa1Lug4Xr6NSqCf07tow7FBGJiVZeEwASCWfqovUcu39HzCzucEQkJlp5TQCYv3Yz67cWcZTaE0QyWrJdTLKBfGAjsL+ZHZO+kCQOak8QEUhu7qN7gAuAuUAiOuzA+2mMS2rZlIXr6Nexhaa2EMlwybQpnAUMcved6Q5G4lFcmmD6kg18a0SPuEMRkZglU320mDBLqjRQM1cUsK2oVFVHIpJUSWEbMNPMJgLlpQWNaG44pixchxkc2U+NzCKZLpmk8HK0SQP174XrObB7a1o3V4FQJNMl0yV1fG0EIvHYXlTKpys2aqpsEQGqnhDvWXc/38xmE3obfYG7H5TWyKRWzFi+keJS54j+ak8QkapLCtdHj6el+qZm1gZ4FBhKSDhXAPOBZ4A+wFLgfHffmOp7yxdNW7ye7CxjZO+2cYciInXAHnsfufua6HFZZVsN7/sA8Ka7DwaGAfOAW4CJ7j4QmBjtS5pNW7yeod1b06qp2hNEJLmG5pQys/2AY4AxAO5eBBSZ2ZnAcdFp44FJwM21HV8m2V5UyswVBVx5dL+4Q5F0KymBHTt2bTt3hseioi9vxcVhKynZ9Vi2lZbueizbEokvP3cPzytuZccqPpZtu++XbfDFxz09L7On4xXtfv7ezknmeLJq+vMVjR0LX/966q4XqfWkAPQjTJnxhJkNI6zZcD3QuULpZI2Zdarsh81sLDAWoFevXrUTcQNV3p7QT6us1VnuUFAAeXmwfj1s2BAeCwp2bZs27dq2bIGtW8Pjtm27tuLi2onXLGzZ2ZCVtet52fGsrF3HK26VHSvbyq5b8fqVPd89ht2P7x5nZc/3dE4yx5OVqgknN6andj2OpJBDWLTnWnefbmYPsA9VRe4+DhgHMHLkyBSm3cxT3p7QR0khFokErF4NixfD8uWwbBmsWBGOrVkTHvPzq/5Cb9kSWreG/faDVq3C1qEDtGgRXmveHJo1C49Nm4atWTNo0mTX1rhx2Bo12vU8Jyfs5+R8ccvO3vVYccvK2vXlL/VaVb2PKu11VKYGvY9WAivdfXq0/xwhKaw1s65RKaErkFfN60uSpi4K4xNaNonjb4MMsnMnfPYZzJkD8+aFbf58WLIkVOFU1L49dO8O3brBgQdC587QqRN07Bi+7Nu3h3btoG3bkAhy9NlJalX1L6qs19HV0eOfo8eLCaOcq8Xdc81shZkNcvf5wInAZ9E2mrB+w2jgpereQ/ZuW1EJs1aqPSHlSkrCl//06TBtGsyYERJCSUl4PScHBg6EQYPg1FNhwADo1w9694YePcJf9CIxqmo5zmUAZjbK3UdVeOkWM5sC/L8a3Pda4K9m1pgwt9LlhJ5Qz5rZlcBy4LwaXF/2YsayArUnpEIiATNnwsSJMGkS/OtfsHlzeK1jRxg5Mnz5Dx8e/vIfMCBUy4jUUcmUPVuY2dHu/gGAmR0FtKjJTd19JjCykpdOrMl1JXlqT6iB7dvhzTfh1Vfh9dchN1pzatAguPhiOPpoOPJI6NtXdexS7ySTFK4EHjez1oQ2hkLCYDOpx6YtXs9BPdSekLTi4pAAnnkGXn459PBp3Rq+8Y1QEjjpJOjaNe4oRWosmbmPPgGGReMLzN0L0x+WpNP2olK1JyRr3jx47DH4859Dt9D27UNp4Lzz4NhjVRUkDU4yK691Bn4OdHP3k81sCHCkuz+W9ugkLcrGJxyu9oTKucNbb8H998Pbb4fG4dNPhyuvDCUD9fiRBiyZRXaeBN4CukX7nwP/l66AJP2mL15PlqH5jnZXWgp/+xsMHQonnwyzZ8Pdd8PKlfDCC6GaSAlBGrhkkkIHd3+WaH1mdy8BStMalaTVtCUbNN9RRYlEaCs48MBQNZSVFaqLli6F224LYwVEMkQySWGrmbUnGshmZkcQGpulHtpRHOY7Oryvqo4AmDo19BS68MLQU+jZZ2HWLLjkkjCyVyTDJFMW/j5h5bX+0fiEjmgMQb01c0UBRSUJDu+b4esnrF4NP/xhqC7q2hWefDIkguzsuCMTiVUySWEucCwwCDDCugfJlDCkDpq+eANmcGimlhTcQ2+iG28MU0zcfjvcfHOYJ0hEkkoKU919BCE5AGBmMwiT2kk9M33Jeg7osh+tm2Vge8LSpaEH0bvvhu6kjzwSppwQkXJVTYjXBegONDOzgwmlBID9AE3QUg8VlSSYsXwjFx3WO+5Qat8zz8BVV4VG5T/9Cb797dCgLCJfUFVJ4RuEhXB6APdVOL4JuC2NMUma/GdlATuKE5k1PmHrVrj++lBldMQRoQ2hb9+4oxKps6qaEG88MN7MvuXuz9diTJIm05dsAOCwTJnvaMkSOOusMN7gttvgzjs1AllkL5IpPx9iZm3KdsysrZn9LI0xSZpMW7yewV1a0bZFBnS1fPddOPTQsHjN66+HQWhKCCJ7lUxSONndC8p23H0jcEr6QpJ0KC5N8MmyjZkxPuHhh8PatZ06wYcfwje/GXdEIvVGMkkh28yalO2YWTOgSRXnSx00e1Uh24pKOaJfAx6f4B66mH7veyERTJum3kUi+yiZLql/ASaa2ROEUc1XAOPTGpWk3LTF6wE4rKGWFIqL4bvfhccfDz2LHnpI8xSJVEMyU2f/ysz+A3yN0C31Lnd/K+2RSUpNW7yBQZ1b0b5lAyzk7dwJF1wAL70EP/lJaFDW4jYi1ZLsn1LzgBJ3/6eZNTezVu6+OZ2BSeoUlyb4eOkGzjukR9yhpN6OHfCtb4XG5N/9Dq65Ju6IROq1vbYpmNl3gOeAP0WHugMvpjMoSa0G256wfTucfXZICA8/rIQgkgLJNDRfDYwiDFrD3RcAndIZlKRWg2xPKCqCc84Ji+E8+mgYrSwiNZZMUtjp7kVlO2aWQzSNttQP0xZvYP/OLRtOe0JpaZjR9M03Ydy4MJ+RiKREMklhspndRpgD6STg78Ar6Q1LUqWsPaHBVB25h15Gf/87/OY3oaeRiKRMMknhFiAfmA1cBbwO3J7OoCR1Glx7wq23huqiH/0IfvCDuKMRaXCqmiV1orufCPzC3W8GHqm9sCRVGlR7wsMPwz33hJLCXXfFHY1Ig1RVl9SuZnYscIaZTWDX1NkAuPuMtEYmKVHWntChvrcnvP46XH01nHpq6HqqcQgiaVFVUvgJoepo96mzITQ0n5CuoCQ1ytoTzq3v4xNmzgyD04YNgwkTNFJZJI2qmjr7OeA5M/uxu6usXg/9Z2UB24pKObI+tyfk5sJpp0HbtvDqq1o2UyTNkpnm4i4z6w70rni+u7+fzsCk5qYsXI8Z9beRuagIzj0XNmyAf/8bunWLOyKRBm+vScHMfglcCHwGlEaHHVBSqOOmLFzHkK771c/1E9zDCOUpU8JSmsOHxx2RSEZIpnL2bGCQu+9MdzCSOtuLSvl0eQFjRvWJO5TqefhheOSR0AX1/PPjjkYkYyQzTmExoCWr6pmPlm6gqDTBUf3rYdXRtGlhXeVTTlHXU5FalkxJYRsw08wmAuWlBXe/riY3NrNs4GNglbufZmZ9gQlAO2AGcGnF6TVk30xZtI5G2Vb/xiesWxdKBj16wF/+AtnZcUckklGSSQovR1uqXU+Yknu/aP8e4H53n2BmDwNXAg+l4b4ZYeqi9Rzcsy3NG9ej7puJRJjTKC8vNCy3bRt3RCIZJ5neRylfZc3MegCnAncD3zczI4x7uCg6ZTxwJ0oK1VK4rZjZqwq57oR6thTl3XeHWU//9CcYMSLuaEQyUlXTXDzr7ueb2WwqmRXV3Q+qwX1/C9wEtIr22wMF7l4S7a8krNtQWVxjgbEAvXr1qkEIDdfUxetxh1EDOsQdSvImTYI77gglhe98J+5oRDJWVSWF66PH01J5QzM7Dchz90/M7Liyw5WcWun03O4+DhgHMHLkSE3hXYl/L1pHs0bZDO/ZJu5QkrNuHVx8Mey/f1hbWVNYiMSmqhHNa6LHZSm+5yjCfEqnAE0JbQq/BdqYWU5UWugBrE7xfTPGlIXrOKxvOxrnJNO5LGbuMGZMSAyvvaYRyyIxq/VvDXe/1d17uHsfwqC4d939YuA94NzotNHAS7UdW0OQW7iDRflb609X1AceCMngN7/RADWROqAu/Sl5M6HReSGhjeGxmOOpl95fkA/AMft3jDmSJMycCTffDGecofWVReqIPSaFaFwCZnZPum7u7pPc/bTo+WJ3P8zdB7j7eRpBXT2TP8+nU6smDO7Sau8nx2nbNrjoImjfHh57TO0IInWE1lNoQEoTzgcL1nHSkM5YXf+SvekmmDcP3n4bOtSjXlIiDZzWU2hAZq0soHB7McfW9aqj11+HP/wBbrgBTjop7mhEpAKtp9CAvP95PmZwdF0en5CXB5dfDgceCD//edzRiMhukl1P4QzgmOjQJHd/Nb1hSXVM/jyfYT3a1N2pst1h7FgoKIB//hOaNo07IhHZzV57H5nZLwgD2T6LtuujY1KHFGwrYtaKgrpddfTkk/DSS6GEcOCBcUcjIpVIZra0U4Hh7p4AMLPxwKfArekMTPbNBwvXkfA63BV1yRK47jo47rjQliAidVKy4xQqzpfQOh2BSM1Mnp9P62aNGNajDn48paUwejRkZYXSQlZdGh4jIhUlU1L4BfCpmb1H6JZ6DCol1CnuzvsL8jl6YAdysuvgF+5998G//gXjx0Pv3nFHIyJVSKah+WkzmwQcSkgKN7t7broDk+T9N3czazft5NiBdbDqaM4cuP12OPtsuPTSuKMRkb1IagWWaHK8dCy0Iynwz8/WYgbHDa5jSaGoCC67DFq3Dmsk1PUBdSKSXFKQuu2deWsZ3rMNnVrVsS6eP/sZfPop/OMf0LGOJSwRqVQdrICWfZFbuIP/rCzkpCGd4w7liz78MHQ9vewyOOusuKMRkSRVmRTMLMvM5tRWMLLv3pm3FoCTDqhDSWH79pAMunYNU2OLSL1RZfWRuyfMbJaZ9XL35bUVlCTvnc/W0qd9cwZ0qkOL09x6K8yfH0Ytt6knq7+JCJBcm0JXYK6ZfQhsLTvo7mekLSpJyuYdxUxdtI4xR/WpO7OivvtuKB1cey2ceGLc0YjIPkomKfw07VFItbz/+TqKS52ThnSJO5SgsDBMdrf//vDLX8YdjYhUQzLjFCabWW9goLv/08yaA9npD0325p3PcmnbvBEjetWRKprrr4eVK2HKFGjePO5oRKQakpkQ7zvAc8CfokPdgRfTGZTsXXFpgnf/m8cJgzvXjVHML7wQRiz/6EdwxBFxRyMi1ZTMt8nVwChgE4C7LwA6pTMo2bsPl2xg09Ywq3IAABE5SURBVI6SutEVNTc3TIl9yCHw4x/HHY2I1EAySWGnuxeV7ZhZDmHlNYnRK7NW06JxdvxTZbvDt78NW7fCn/8MjRrFG4+I1EgySWGymd0GNDOzk4C/A6+kNyypSlFJgjfm5HLSkM40axxz884jj8Brr8E998ABB8Qbi4jUWDJJ4RYgH5gNXAW8DtyezqCkah8szKdwezFnDO8WbyDz54e1Eb72NbjmmnhjEZGUSKb3USJaWGc6odpovrur+ihGL89cTetmjTh6QIxVR0VFcPHFYUnN8eO1RoJIA7HXpGBmpwIPA4sIU2f3NbOr3P2NdAcnX7a9qJR3PlvL6cO60Tgnxi/iO++ETz4JvY66xVxiEZGUSWbw2r3A8e6+EMDM+gOvAUoKMXhvfh5bi0o5Y1iMX8STJ4fBaVdeGdZJEJEGI5k/NfPKEkJkMZCXpnhkL16euZqOrZpweL/28QSwfn2oNurfH37723hiEJG02WNJwczOiZ7ONbPXgWcJbQrnAR/VQmyym807inl3fh4XHdaL7KwY5jpyhyuugLw8mDYNWtahSfhEJCWqqj46vcLztcCx0fN8oG3aIpI9enNOLkUlCU6Pq+ro97+Hl1+G+++HESPiiUFE0mqPScHdL6/NQGTvJny0gv4dW8Qz19HMmXDjjXDqqWGOIxFpkJLpfdQXuBboU/F8TZ1duz5fu5lPlm3kR6ccUPvTZBcWwnnnQYcO8MQTWmtZpAFLpvfRi8BjhFHMifSGI3vy9IfLaZRtnDOie+3e2D30MlqyBN57T2stizRwySSFHe7+YKpuaGY9gaeALoQkM87dHzCzdsAzhBLJUuB8d9+YqvvWZzuKS/nHp6v4+le60L5lk9q9+YMPwvPPw69/DV/9au3eW0RqXTJdUh8wszvM7EgzG1G21eCeJcAP3P0A4AjgajMbQphOY6K7DwQmRvsCvDU3l4JtxVx0WK/avfHUqaEd4cwz4Qc/qN17i0gskikpHAhcCpzAruojj/b3mbuvAdZEzzeb2TzCGg1nAsdFp40HJgE3V+ceDc3THy6nV7vmHFmbYxNyc+Hcc6FnT3jySbUjiGSIZJLC2UC/itNnp4qZ9QEOJsyr1DlKGLj7GjOrdM0GMxsLjAXo1auW/3KOweL8LUxbvIEffmMQWbU1NmHnTvjWt6CgAP79b2hTR1Z2E5G0S6b6aBaQ8m8FM2sJPA/8n7tvSvbn3H2cu49095EdM6DR86mpy8jJMs47pEft3fS660IyeOIJGDas9u4rIrFLpqTQGfivmX0E7Cw7WJMuqWbWiJAQ/uruL0SH15pZ16iU0BVNpcHGrUU889EKzhzenU77Na2dmz78MIwbB7fcAuefXzv3FJE6I5mkcEcqb2ihk/1jwDx3v6/CSy8Do4FfRo8vpfK+9dFfpi1je3EpY4/pVzs3fOedsC7CySfDz35WO/cUkTolmfUUJqf4nqMIDdezzWxmdOw2QjJ41syuBJYT5ljKWDuKSxk/dSnHD+rIoC6t0n/DuXNDw/IBB8CECZAd84puIhKLZEY0b2bXmsyNgUbAVnffrzo3dPcPCOsyVObE6lyzIXp+xkrWbSniqmP7p/9ma9eG6SuaNYNXX4X9qvXRikgDkExJ4Qt/pprZWcBhaYtIKE04j/5rCcN6tObwvu3Se7MtW+CMM8LMp5MnQ+/e6b2fiNRp+7x0l7u/SDXHKEhy3p6by5J1W7nq2P7pneeoqCh0Pf34Y/jb3+DQQ9N3LxGpF5KpPjqnwm4WMJJd1UmSYiWlCe5753P6dWzBN77SJX03SiRg9Gh4+2149FE466z03UtE6o1keh9VXFehhDAv0ZlpiUZ4fsZKFuRt4eFLRqRvIR33MP31hAm7ltUUESG5NgWtq1BLtheVcv87Czi4V5v0lRLcwzxGv/99mNfoppvScx8RqZeqWo7zJ1X8nLv7XWmIJ6M9+e+l5G7awQMXDk9PW4J7SAL33x9GLf/qV5rTSES+oKqSwtZKjrUArgTaA0oKKbRxaxF/nLSQEwd34vB0THznDrfeCr/5DVx9Nfz2t0oIIvIlVS3HeW/ZczNrBVwPXA5MAO7d089J9TwwcQFbdpZw0zcHp/7iiUQoGfzhD/Dd78LvfqeEICKVqrJNIVr45vvAxYTprEdo4ZvUm7F8I+OnLuWSw3unfvRycTGMGRO6nP7wh3DPPUoIIrJHVbUp/Bo4BxgHHOjuW2otqgxSVJLgluf/Q5f9mnLTNwel9uJbt8IFF8Brr8EvfhEmuRMRqUJVJYUfEGZFvR34UYWGTyM0NGsuhBR4aNIiPl+7hcfHjKRV00apu/CqVXD66TBrVpj59KqrUndtEWmwqmpT2OfRzrJvFqzdzO/fW8AZw7pxwuDOqbvwp5/CaafBpk3wyitwyimpu7aINGj64o/JjuJS/u+ZmbRsksMdpw9J3YWffhqOPjrMcjplihKCiOwTJYWY3PHSXOau3sS95w+jfcsmNb9gURFcey1cdBEcfDBMnw4HHVTz64pIRlFSiMGzH63gmY9XcM3xA1JTbbR0KRx7bBilfMMN8N570LVrza8rIhknmbmPJIXmrCrkxy/NYdSA9txw0v41u5g7PPVUKCEAPPOMltAUkRpRSaEWrdiwjW+P/5i2zRvz4IUH12zCu7y8kADGjIHhw+E//1FCEJEaU1KoJXmbd3DpY9PZVlTCE5cfWv12BHd44omwbOZLL4VZTt97D/r0SWm8IpKZlBRqQeG2Yi577EPWbtrJE5cfxgFdqznEY+5cOOEEuOIKGDIEZs6Em2/WesoikjJKCmm2bstOLnt8OovztzLuskM4pHfbfb9IXh5873uhN9HMmTBuXFg6c0gKu7KKiKCG5rRanL+FMU98RN7mHTx0yQi+OrDjvl1g0yZ44IEws+nWrWF20zvugPZpmEVVRAQlhbT5eOkGvv3Ux2Sb8fR3juDgXvtQQti0Kcxkeu+9sHEjnHlmaDsYnIYZVEVEKlBSSLFEwvnT+4u59+359GzXnCcvP5Te7Vsk98PLl4eSwSOPwObNYe6iO++EESPSGrOISBklhRRaU7id7z8zi6mL13PqgV35+dkH0rr5Xia5SyRg4sTQTvCPf4RjF1wQlsxUMhCRWqakkALFpQmemrqM3/7zc0oTzq/OPYjzDulR9ZKaCxeGNQ6efBKWLIEOHeD73w8D0Xr2rLXYRUQqUlKoAXfn/QXruOvVz1iYt4Vj9u/IT8/4Cn077KG6aPFiePHFMPL4ww/DYjfHHQc//zmcfTY0ScEcSCIiNaCkUA2JhPPOvLX8cdIiZq0ooFe75jx62UhOPKDTF0sHRUUwbRq89VaYwnr27HD84IPh17+GCy+EHj3ieRMiIpVQUtgH67fs5B+frmLCRytYmLeFnu2a8bOzhnLuIT1o2ig7LH05Ywa8/34YRzB5MmzZEgaXjRoF998fehL17Rv3WxERqZSSwl4Ubivm3flreWN2Lu/+N4+ShDOsZxseOO9ATm26hZzZU+HvH4XqoE8+ge3bww8OGgSXXALf+AYcfzy0bh3vGxERSYKSwm6KShLMWlnA9MXrmbJwPTMW5dFtYy7Dd+TzcKMCDt2WS+v3Pocfzt6VAJo0CT2FrroqlAi++lXonMKV1EREaklGJ4WtO4pZsjSXFXMWsXbeIjYtWEJi+XI6b1zLsMI8ztmaT5cNuWQlSnf9UOfOYTK6sWPD7KTDhsHQodAohesri4jEpE4lBTP7JvAAkA086u6/TMd9Prz9V3R76AHabV7P0OKdDN3t9Z3tO5Ldry85/b8KAwaEbeDAMKK4Xbt0hCQiUifUmaRgZtnAH4CTgJXAR2b2srt/lup7NenWlbWDD2Jd1y407d2DNv1703Fwf7J79YTu3WnSrFmqbykiUi/UmaQAHAYsdPfFAGY2ATgTSHlSGPa/l8L/Xprqy4qI1Ht1aers7sCKCvsro2NfYGZjzexjM/s4Pz+/1oITEckEdSkpVDYnhH/pgPs4dx/p7iM7dtzHqahFRKRKdSkprAQqTvrTA1gdUywiIhmpLiWFj4CBZtbXzBoDFwIvxxyTiEhGqTMNze5eYmbXAG8RuqQ+7u5zYw5LRCSj1JmkAODurwOvxx2HiEimqkvVRyIiEjMlBRERKWfuX+r1WW+YWT6wrJo/3gFYl8Jw6otMfN+Z+J4hM993Jr5n2Pf33dvdK+3TX6+TQk2Y2cfuPjLuOGpbJr7vTHzPkJnvOxPfM6T2fav6SEREyikpiIhIuUxOCuPiDiAmmfi+M/E9Q2a+70x8z5DC952xbQoiIvJlmVxSEBGR3SgpiIhIuYxMCmb2TTObb2YLzeyWuONJBzPraWbvmdk8M5trZtdHx9uZ2TtmtiB6bBt3rKlmZtlm9qmZvRrt9zWz6dF7fiaacLFBMbM2Zvacmf03+syPzJDP+obo3/ccM3vazJo2tM/bzB43szwzm1PhWKWfrQUPRt9t/zGzEft6v4xLChWW/TwZGAL8j5kNiTeqtCgBfuDuBwBHAFdH7/MWYKK7DwQmRvsNzfXAvAr79wD3R+95I3BlLFGl1wPAm+4+GBhGeP8N+rM2s+7AdcBIdx9KmEjzQhre5/0k8M3dju3psz0ZGBhtY4GH9vVmGZcUqLDsp7sXAWXLfjYo7r7G3WdEzzcTviS6E97r+Oi08cBZ8USYHmbWAzgVeDTaN+AE4LnolIb4nvcDjgEeA3D3IncvoIF/1pEcoJmZ5QDNgTU0sM/b3d8HNux2eE+f7ZnAUx5MA9qYWdd9uV8mJoWklv1sSMysD3AwMB3o7O5rICQOoFN8kaXFb4GbgES03x4ocPeSaL8hft79gHzgiaja7FEza0ED/6zdfRXwG2A5IRkUAp/Q8D9v2PNnW+Pvt0xMCkkt+9lQmFlL4Hng/9x9U9zxpJOZnQbkufsnFQ9XcmpD+7xzgBHAQ+5+MLCVBlZVVJmoHv1MoC/QDWhBqD7ZXUP7vKtS43/vmZgUMmbZTzNrREgIf3X3F6LDa8uKk9FjXlzxpcEo4AwzW0qoFjyBUHJoE1UvQMP8vFcCK919erT/HCFJNOTPGuBrwBJ3z3f3YuAF4Cga/ucNe/5sa/z9lolJISOW/Yzq0h8D5rn7fRVeehkYHT0fDbxU27Gli7vf6u493L0P4XN9190vBt4Dzo1Oa1DvGcDdc4EVZjYoOnQi8BkN+LOOLAeOMLPm0b/3svfdoD/vyJ4+25eBy6JeSEcAhWXVTMnKyBHNZnYK4S/IsmU/7445pJQzs6OBfwGz2VW/fhuhXeFZoBfhP9V57r57I1a9Z2bHATe6+2lm1o9QcmgHfApc4u4744wv1cxsOKFxvTGwGLic8Edfg/6szeynwAWE3nafAt8m1KE3mM/bzJ4GjiNMj70WuAN4kUo+2yg5/p7QW2kbcLm7f7xP98vEpCAiIpXLxOojERHZAyUFEREpp6QgIiLllBRERKSckoKIiJRTUhBJUjQT6f9Gz7uZ2XN7+xmR+kZdUkWSFM0h9Wo0I6dIg5Sz91NEJPJLoL+ZzQQWAAe4+1AzG0OYpTIbGArcSxhEdimwEzglGljUnzBte0fCwKLvuPt/a/9tiOyZqo9EkncLsMjdhwM/3O21ocBFhKnZ7wa2RZPTTQUui84ZB1zr7ocANwJ/rJWoRfaBSgoiqfFetG7FZjMrBF6Jjs8GDopmqz0K+HuYiQCAJrUfpkjVlBREUqPi3DqJCvsJwv+zLMI8/8NrOzCRfaHqI5HkbQZaVecHo7UslpjZeVC+lu6wVAYnkgpKCiJJcvf1wJRoAfVfV+MSFwNXmtksYC4NcBlYqf/UJVVERMqppCAiIuWUFEREpJySgoiIlFNSEBGRckoKIiJSTklBRETKKSmIiEi5/w9me6GG1nQGTAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "Wupdate_sis=Wtt_sis.subs([(N,100),(tau,0.5),(c,0.5),(a,0.1)])\n",
    "ws_sis=[1]\n",
    "for t in range(0,99):\n",
    "    ws_sis.append(Wupdate_sis.subs([(Wt,ws_sis[-1:][0])]))\n",
    "plt.plot(range(0,50),ws)\n",
    "plt.plot(range(0,100),ws_sis,color='r')\n",
    "plt.ylabel('Number of infected individuals')\n",
    "plt.xlabel('time')\n",
    "plt.title('Diffusion of infection');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# SIR\n",
    "\n",
    "We now change the previous model, so that people who recover become immune to the disease and do not get it again. This version is called the SIR model (S for the number of susceptible, I for the number of infectious, and R for the number recovered individuals). The derivation is essentially the same as before."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [],
   "source": [
    "Rt = sm.symbols('Rt')\n",
    "Rtt, Penc1 = sm.symbols('Rtt,Penc1',cls=sm.Function)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Penc = (Wt/N)*((N - Wt)/N)\n",
    "Penc1= (Wt/N)*((N - Rt - Wt)/N)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAK0AAAAsBAMAAADlbSUNAAAAMFBMVEX///8AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAv3aB7AAAAD3RSTlMAdt27RKvvEFRmiTIimc3h244VAAAACXBIWXMAAA7EAAAOxAGVKw4bAAADjUlEQVRIDe1XQWgbVxB9WktZaSXZS0sPgRTrEpOWgBQaUgwJXkIOIiFYFNIQyMGY5tBTVJyzbdxLwAfrUmghxHso6JQ4l5y9SaD4EqwQcumh9aHQQkusWHWDEzvqm1nZWq1ssFn5EMiA/s68mf92/vy//38B3ZKpdGMHQJKeBMfyDSdx83g5mf9su+8ElY+bDpL5G+VtzH/+0LzyxXAnJNZHzdP4eu0qZtY0o+saMeoCCw7wkJbhscE/0owUgP4QLdAAfikBJyWiLfE3gPUvkKaLFJPqWa75vC6tOUFSBBD/fpOxYgYl/QpIvEWyHgSBzJbPK/GksGxxT5WA1QriDvXfBNCkrVuMLYoZlBhfZmwiuxIEAXOdY9sALhAmhTkm7oEVxOZtP7fPBXgsTRZ5oCpaUPrGmFABhhsEqb8DjnIoAgvFH/zhSA5Liy6OcezPG88IfCpoFvMVSEE6ZKCE5A1cHjmnI227xpGZfYOY16K4Jx5jDO5yzZRIv2464UUM5PT97c7UFq99svWgNQ1BxxAs1vi7bYq74uuvJx9MFaTkrbqdELXKCdaSi7Ejg2XcLbemYQeksopZs6FpaOlnxNf3OsYaS8lbdftG1BoyDX0XJl+K/C0gznA6XL+GMBV/yfS5KMs23mkaWvpFwbKbNmvsijrlSKu8LoOfitUhnPflEkwu1U75/REwrmkohfJaG9MwvvUkcBQmW6lDfJrWr1Q6JM21NFdDvA4WIygLPwJDdwRRivOixfgN9XORUO5jlq3Mm8WORldasbf8hGzTWknq6KWLL4MF1tgWXSl03lJc64mc+l/oqpBVYvCXqisYaPqJzNlL2ZLMfFDmPeArBZRCBxp3mZ6j4JOf5MFRTNyUTaikYLvJrq5XYJ2yM8OVNqjaRbY/q6YU/iSHYvzvOAwexDZzu0brvrOrZ5+gv+90B5/thg6ETO4RLft6FNkrr4jnUMqLktSHvvurQPNQpGsL2F8y70/UvXPlw0jW8hIrh8FreGlu172T2AgPj5n/rhpOpusUifSWPp5W+JO/5OtIPOHOxRHPvycZdtgVya6O8hyeJoXekfjskdSMuh7SkQ+TUD611AbvgsBFXAp5IpnxaazJVTb91+0vIxGFOvOesuBc4KW22ezplsISLBfc0Mt6YLIEfZu1HhCFKIq8Cm65IbAHZpUceacHRCGKHO35EBbdNM/IH8nEQYj+B9n7MR/tbUynAAAAAElFTkSuQmCC\n",
      "text/latex": [
       "$$\\frac{Wt \\left(N - Rt - Wt\\right)}{N^{2}}$$"
      ],
      "text/plain": [
       "Wt⋅(N - Rt - Wt)\n",
       "────────────────\n",
       "        2       \n",
       "       N        "
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Penc1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [],
   "source": [
    "Wtt_sir = sm.symbols('Wtt_sir',cls=sm.Function)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAU0AAAAsBAMAAAATAfe8AAAAMFBMVEX///8AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAv3aB7AAAAD3RSTlMAEM0ydt27RKvvVGaJIpmJty9HAAAACXBIWXMAAA7EAAAOxAGVKw4bAAAE8UlEQVRYCdVYXWgcVRT+Zrs/ye7M7hqKPm5MER+zj0WQDLQqIpIFX1pr2cUHpYjuvhQxFLMVFVQkUagU+tCFEMEnF+mTopko6EPB7INIsWiWtkgfAt2EJKak6Xruudm5c2fXLdtJNut5uHPP+c4998y9Z+58M8D+y/VgU8RrwcZ7Rz/TfB0fbZzH/EbVa+a+UQDCxX+A+bvndTAxs/Xha+0DwvnN0tDlE+l4/hT7f6qPCqLF7gEmpXKoTFFO65FCYj2im9Q8rgNAqgDrIhDxrdi0AyyVgKPkT9iIf9jD69Z9madJIeINPc4PQp0qUi6ODgATNjCTxdc++3Jd5umQnTDTpuveiLEFhLaBtylcclSPycu7OF1HrKIDwCQln6/hhs+eKQPjVcRoSQVmzPrwAOoa8MK6XLGIo8WxckKtRxow0xpAyhlafSqIV3z21CjCCzbE5jD2pA8PoF6ANXcP4RrwbvEtmuL4MTzfbDZteoTqImw9sY2k6GhCu/BVOfb35i3NiuEc/phw8CIgsV91OIh2EibV6GERgirKehmvmj9de/pPUqNVamjLN+T6CI+WHNq6MnOzvaARmYWzXDds8uNi/6Y1IPh1HHPGpnxQqKJo6z9+ClEOG6nQhbZ8qSSKV5PwKhJ0BvgLGqFGPJsp8LYzNq+NCqQU0zbWuPBFRd3hSjzKEYdLdKEtXy44rI+MCTki+tFZgA40uisSr301jNQo3xZjEwLfG7nzG3CBV8CgY/QLDnqVW86TkOhO3T9VqgyLSjQj7sQryR2batQRJsb2MM+lb4GT34vIsQYeoUUieYnbSIUuU2S/71BHk4kq4nTwTsPQzDC3K4hcqgkjY2d1OIg2WaAjzxYRzNF4lk7SOKyGUOVztEidvH/ZMJmFuY70j5hjT7cJ01kVolOOhLE9fI4W6O4/4MjJ8mHchnET4RzrZp0uorvAqrcZp9tZtyq3DdtrBRI7wJAczdiKDgfRztHgXziA9UYV5tUr9OzYrMdyMM4IejLEqqe50fwE+Pk4fv/OYxTdmEPbUmIjY0e4K5tujIc9HujgCaZ1T2vaQyiGXFo50st4/PyFPbwOfkrUdXLmIV09HgCatsfBy3jojaIkbsu+16FFiVqY8u7Qi9Y6GHsxjXidvYxH4y+tXLwOrTdIC/MGausbhTZTb4Y3Nfc1xXjojaLEzcXjwG8JcnEx5d6hd72DrQdTQt8Pl/EwR7FWbn0mY7m5uA6SEgnUxXqYNqirYjyi+t6phagV4uaiHFwC7mLStS+tYjxUfcYlJAtyWjcX5eAScBfb/wxdmqIYD1VfchURmyY3x8ae+HxsrCzyUA6SgHsxgfdJFOMhjjJcQCYrJ3bXTDkISsTSwlJE3fdVdmtQTKoYD3GUjNNeg8qBKFFay5OVPjWK8RBHyVTw3u5/iNaaQTkISnRgeSrGQxwlaidOmHKB3DyVg6BEB5anYjzEUYz3j117zpeLchCU6MDylDO3t+56tkPqbO2A9dtkyEem47TdsI4D2Phss4R4/mKXuP89tq9IsUDfFYOfZuxL+q7YfUr7uj49TmY+Rh+Z9AU66JJEHhBfoAMuSSxUUR/wJCm9KaRy8lfUYOe6SL/I5a+owc6zDmvzf/C4iy0v/jXYSymyo3+2mF4Z/DzFb/pIi2kPcLoRyi3RGOAEZWqPXj5FnXKPef4Lulu3NdGiGX0AAAAASUVORK5CYII=\n",
      "text/latex": [
       "$$- Wt a + Wt + \\frac{Wt c \\tau \\left(N - Rt - Wt\\right)}{N}$$"
      ],
      "text/plain": [
       "             Wt⋅c⋅τ⋅(N - Rt - Wt)\n",
       "-Wt⋅a + Wt + ────────────────────\n",
       "                      N          "
      ]
     },
     "execution_count": 30,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Wtt_sir=N*c*Penc1*tau + Wt - a*Wt\n",
    "Wtt_sir"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAFkAAAAQBAMAAACPcs44AAAAMFBMVEX///8AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAv3aB7AAAAD3RSTlMAZqu7Iu8ydt1UmRDNiUQDe6IrAAAACXBIWXMAAA7EAAAOxAGVKw4bAAABl0lEQVQoFXWRvy8sURTHP1i7O3PN2icaicSQyJNQrLyI5BVMiIgKr5JQ7H9gvEJUJEIpUW10Xig0xHbKVSo22f0DBI1GgmXFE7/PnWEmE3GKO/d+z+ee+z1ngKa3jl+9fImFtwE6H7rZeEhHclUYsWElIpJ6AesV4pJKOkEu/h9q7zEqvmK4/jf26NOWHDcDmMSdXL7DPPSlT1rdQsMNXIm8G9J1E3I5Q7IcpXmGOXlWy2sh3Whj/OFi8Ldv4bM2U8RyLyQcUtvVHWJ7O9P6zmz738ds6C2gV7HE+4EQuqVzp0FWWMrTmg+9BXSJnKp6RqQlNYmZ0XQfNJY/vFnF4tFMsWhrfTDv8pw6k520ZD6RdLUqvY/bKBmtF0Ht/VGY0vOjcEZ9hkJWtnGZ0uY/UhXEjkRAn7bAarOWrlGFsj/0xD2cuso6NLI6E9JL4rPkvX5MrrBOj+HIH6hIbXfMtHX3EkHtNge6PGlLuXVuzbK4Mku3aax+N9ab9jIhfSnnE08b/oH6OT+06AORNagdUb87qPx3GXgHPm5q0+Jkt/MAAAAASUVORK5CYII=\n",
      "text/latex": [
       "$$Rt + Wt a$$"
      ],
      "text/plain": [
       "Rt + Wt⋅a"
      ]
     },
     "execution_count": 31,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Rtt = Rt + a * Wt\n",
    "Rtt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAEWCAYAAABhffzLAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3dd3hUZdrH8e+dTkIJoXdCB+mEagMrCoiIYheVXXcta69rede14Vp21bWuKNgFREXFyoKIIEjovQYSWkIJBEL6/f4xh91shGQSMnOm3J/rOtfMHGbO+WUyuTnznOc8j6gqxhhjwkeE2wGMMcb4lxV+Y4wJM1b4jTEmzFjhN8aYMGOF3xhjwowVfmOMCTNW+I1fichrIvJwqcc3ishuETkkIvVE5GQR2eA8vvAE9vO1iIytntSV2u/jIrJHRHYd499OFZF1ldjW/7w31ZjxShH5rrq2Z4KPWD9+U11EJA1oBBQBxcBq4B3gDVUtOcbzo4GDwABVXeasmwlMV9UX/JW7uohIC2A90EpVM09wW795b6q4ndbAFiBaVYtOJJMJHXbEb6rbCFWtBbQCxgP3AROO89xGQBywqtS6VmUeB5NWwN4TLfqOY703xlQLK/zGJ1T1gKpOBy4FxopIVwARmeg0h3QAjjZ7ZIvIv0VkE9AG+MJp3ogVkTQROevodkXkLyLynnM/TkTeE5G9IpItIr+KSCPn32aLyO+c+xEi8pCIbBWRTBF5R0TqOP/WWkRURMaKyDanmebB4/1cIlLHeX2Ws72HnO2fBXwPNHWyTzzGaweLSEapx2kicreILBeRAyLysfMz/ea9cZ7fSUS+F5F9IrJORMaU2lYNEXnOyXRAROaKSA1gTqntHBKRgSJyrYjMLfXaQc57d8C5HVTq32aLyGMi8rOI5IjIdyJSv5xfvQkCVviNT6nqQiADOLXM+vXASc7DRFU9Q1XbAtvwfGuoqar5FWx+LFAHaAHUA/4IHDnG8651liF4/mOpCfyzzHNOAToCZwKPiEjn4+zzJWefbYDTgWuA61T1B+A8YIeT/doKsh81BhgKJAPdgWuP9d6ISAKe/1g+ABoClwOviMjR5z0L9AEGAUnAvUAJcFqp7dRU1fmldy4iScBXwIt43sPnga/KnFO4ArjO2W8McLeXP5sJUFb4jT/swFOMqlshnmLVTlWLVTVVVQ8e43lXAs+r6mZVPQQ8AFwmIlGlnvOoqh5x2tOXAT3KbkREIvF8g3lAVXNUNQ14Drj6BH6GF1V1h6ruA74Aeh7necOBNFV9W1WLVHUx8AlwsYhEANcDt6nqdue9mOfFf5wAw4ANqvqus90PgbXAiFLPeVtV16vqEWByORlNkLDCb/yhGbDPB9t9F/gW+EhEdojI35yTomU1BbaWerwViMLTjn5U6V44uXi+FZRVH88Rb9ltNatC9srsFzznD/o7TVrZIpKN5z+0xk6uOGBTFfZf9r2B3/5M3mY0QcIKv/EpEemLp4jMrei5x3EYiC/1uPHRO6paqKqPqmoXPE0cw/E0vZS1A0/hPKolnp5HuyuZZQ+ebxllt7W9ktupinTgR1VNLLXUVNUbnVx5QNtjvK6ibntl3xvw389kXGKF3/iEiNQWkeHAR8B7qrqiiptaiqdZJlpEUoCLS+1jiIh0c5pgDuIpysXH2MaHwB0ikiwiNYEngY8r271RVYvxNHU8ISK1RKQVcCfwXpV+ssr5EuggIlc770W0iPQVkc5OV9m3gOdFpKmIRDoncWOBLDxt/W2Os90ZznavEJEoEbkU6OLsz4QoK/ymun0hIjl4jlAfxHOy8LoT2N7DeI5k9wOP4jm5eVRjYCqeor8G+JFjF+G38DQLzcHTpz0P+FMV8/wJz7eQzXi+xXzgbN+nVDUHOAe4DM9R+i7gaSDWecrdwArgVzzNak8DEaqaCzwB/Ow0EQ0os929eL4p3QXsxXNSeLiq7vH1z2TcYxdwGWNMmLEjfmOMCTNW+I0xJsxY4TfGmDBjhd8YY8JMVMVPcV/9+vW1devWbscwxpigkpqaukdVG5RdHxSFv3Xr1ixatMjtGMYYE1REpOxV2YA19RhjTNixwm+MMWHGCr8xxoQZK/zGGBNmrPAbY0yYscJvjDFhxgq/McaEmaDox29MVakqGfuPsGDLPnYdOEJSQiz1asZQv2YM7RrUok78sSbsMia0WeE3ISl16z7emb+VhVv2sfNA3jGfExMZwbldG3N53xYMaFOPiAjxc0pj3GGF34SUnQeOMP7rtXy+dAd146MZ1K4+/ZOT6J9cj1b14tmfW8DeQwVkHcrnx3VZTFucwRfLdtCqXjw3D2nHJX2aI2L/AZjQFhQTsaSkpKgN2WDKU1hcwus/buLlWZsoVuWGU9tw4+C2JMSWf2yTV1jM1yt3MmneVpamZzOkYwPGj+5Oo9pxfkpujO+ISKqqpvxmvRV+E+yOFBRz0/upzFqXxbknNeLB87vQsl58xS8spaREeWd+GuO/WUtMZASPjjyJC3s2s6N/E9SOV/itV48JageOFHLNWwuYvT6LJ0Z15fWrUypd9AEiIoRrT07m69tOo32jWtzx8TL+/OlKiksC/8DImMqywm+CVmZOHpe+Pp+l6dn88/LeXNm/1QlvM7l+ApP/MJAbB7flw4XbuPXDJRQUlVRDWmMCh53cNUEpKyefS16bT+bBfCaM7ctpHX4z5HiVRUYI9w3tRN34aJ6csZac/CJeu6o38TH252JCgx3xm6BTWFzCze8vZvfBPN77Xf9qLfql3XBaW54e3Y25G7K4esJCDuYV+mQ/xvibFX4TdJ74ag0L0/bx9Oju9GlV16f7urRvS16+ojfL0rO56b3FFBZbs48Jflb4TVCZtjiDifPSGHdKMiN7NvPLPs/r1oTxo7szd+Me/jxtBcHQE86Y8lijpQkaK7cf4IFpKxjQJokHzuvk131f3Kc56ftyeWHmBlomxfOnM9v7df/GVCcr/CYoHMov4o/vpVIvIYaXr+hNVKT/v6zeflZ70vfl8tz362meVINRvZr7PYMx1cEKvwkKz367ju3ZR5j6x0HUqxnrSgYRYfzo7uw8kMe9U5fTul4CvVr69hyDMb5gbfwm4C3Ztp9J89O4ZkArn5/MrUhMVASvXdWHRrXjuOWDJRzItZ4+JvhY4TcBrbC4hAemraBx7TjuGerfdv3jqRMfzUuX92L3wTzumbrMTvaaoGOF3wS0N+ZsZu2uHP46sis1KxhwzZ96tazL/ed14rvVu5k4L83tOMZUihV+E7C27DnMCzM3cH63xpzdpZHbcX5j3CnJnNW5IU/OWMPyjGy34xjjNSv8JiCpKg9/tpLYqAj+MuIkt+Mck4jw7CU9aFAzlls+WMKh/CK3IxnjFSv8JiDN2bCHuRv3cOfZHWgYwGPjJ8bH8MLlvUjfn8v4r9e4HccYr1jhNwGnpEQZ//VaWibFV8uIm77Wt3US405O5r1ftvHzxj1uxzGmQlb4TcD5YvkO1uw8yF3ndCAmKjg+onef25E29RO4d+pya/IxAS84/qpM2CgoKuHZ79bRpUltRnRv6nYcr8VFR/LMJd3ZceAIT86wJh8T2Kzwm4DywYKtpO87wn3ndSIiIrimPezTKonfnZLMBwu2MXeDNfmYwGWF3wSMQ/lFvPTvjQxsU4/T2td3O06V3HWOp8nnvk+Wk1tgTT4mMFVY+EXkZBFJcO5fJSLPi0jgn3EzQWfCT1vYe7iA+87rFLSTnMdFRzJ+dHe2Zx/hxZkb3Y5jzDF5c8T/KpArIj2Ae4GtwDs+TWXCTk5eIRPmbuacLo3o2SLR7TgnpF9yEpf0ac6bP21m3a4ct+MY8xveFP4i9QxGMhJ4QVVfAGr5NpYJNx8s2MbBvCJuOaOd21GqxQPnd6ZWXBQPfbaCkhIby8cEFm8Kf46IPABcBXwlIpFAtG9jmXCSV1jMm3O3cEq7+nRvHtxH+0clJcTwwPmd+TVtP1NTM9yOY8z/8KbwXwrkA+NUdRfQDHimoheJSAsRmSUia0RklYjc5qxPEpHvRWSDc2sDmoe5TxZnkJWTz02D27odpVpd3Ls5/Von8eTXa9h3uMDtOMb8R4WFX1V3qerzqvqT83ibqnrTxl8E3KWqnYEBwM0i0gW4H5ipqu2Bmc5jE6aKikt4/cfN9Gheh4Ft67kdp1pFRAiPj+rKobwinrK+/SaAHLfwi0iOiBw8xpIjIgcr2rCq7lTVxc79HGANnm8LI4FJztMmARee+I9hgtVXK3aybV8uNw5uF7Q9ecrToVEtxp2azJTUDJam2wieJjAct/Crai1VrX2MpZaq1q7MTkSkNdALWAA0UtWdzj52Ag2P85obRGSRiCzKysqqzO5MkFBVXp29iXYNa3JOAA67XF3+dEZ7GtSK5S/TV9mJXhMQvL6AS0QaikjLo0slXlcT+AS4XVUr/KZwlKq+oaopqprSoEEDb19mgsjsdVms3ZXDH09vG3RX6VZGzdgo7hvaiaXp2Uxbst3tOMZ4dQHXBSKyAdgC/AikAV97s3ERicZT9N9X1WnO6t0i0sT59yZAZhVymxDwr58206ROHCN7Bs+YPFV1Ua9m9GyRyNPfrLVB3IzrvDnifwzPydn1qpoMnAn8XNGLxNNgOwFYo6rPl/qn6cBY5/5Y4PNKJTYhYd2uHOZt2ss1A1sTHRn6I4dERAh/ueAksnLyeenfG9yOY8KcN39xhaq6F4gQkQhVnQX09OJ1JwNXA2eIyFJnOR8YD5ztfIs423lswszEeWnERkVwWd8Wbkfxm54tErm4T3PemruFzVmH3I5jwpg3s1dnO+30c4D3RSQTT1fNcqnqXOB4Dbdneh/RhJrs3AI+XZLBqF7NqJsQ43Ycv7p3aEe+WbmLJ2es4c2xfd2OY8KUN0f8I4EjwB3AN8AmYIQvQ5nQ9vGv6eQVljB2UGu3o/hdw1px3DykHT+syWSezdZlXOLNBVyHVbVYVYtUdZKqvug0/RhTaUXFJbwzfysD2iTRuUmlegWHjOtObk2zxBo8/tUaiq17p3GBN716Sl/IlScixd5cwGXMsfywJpPt2Ue4dlCy21FcExcdyX3ndWL1zoN8stjG8TH+580Rf+kLueKA0cA/fR/NhKKJ87bQLLEGZ3U+5nV7YWNE9yb0apnIs9+u47B17zR+Vul+dKr6GXCGD7KYELdm50F+2byPawa2IioMunCWR0R4aFgXMnPyeWPOZrfjmDBTYa8eEbmo1MMIIAWwhklTae/9spXYqAguDaMunOXp06ouw7s34fU5m7isXwua1KnhdiQTJrw57BpRajkXyMHT08cYrx3OL+LzpTsY1r0JifHh1YWzPPcN7URJCTz/3Xq3o5gwUuERv6pe548gJrR9sWwHh/KLuLK/18M8hYUWSfGMHdSKN+du4fpTksO2p5Pxr+MWfhF5iXKadFT1Vp8kMiHpg4Xb6NCoJr1b2rw7Zd0ypD2TF2Uw/uu1TLq+n9txTBgor6lnEZAKxAG9gQ3O0hMo9n00EypWbj/A8owDXNGvZUiOuX+i6sRHc8uQdvy4Pou5G+yiLuN75Y3HP0lVJwHtgSGq+pKqvoRnuAVvxuoxBoD3F2wjLjqCUb2bux0lYF0zqBXN69bgyRlrbMx+43PenNxtCtQq9bims86YCh3KL2L60u0M796UOjWi3Y4TsGKjIrnn3I6s3nmQz5fZmP3Gt7wp/OOBJSIyUUQmAouBJ32ayoSM6Ut3cLigmCvspG6FRnRvSrdmdXj22/XkFVprqvEdb67cfRvoD3zqLAOdJiBjKvTBwq10alyLXi0S3Y4S8CIihAfO78T27CNMmpfmdhwTwsqbbL2Tc9sbT9NOurM0ddYZU64VGQdYuf0gV/S3k7reGtS2PoM7NuCV2Zs4kFvodhwToso74r/TuX3uGMuzPs5lQsDHi7YRGxXByJ7N3I4SVO4b2omDeYW8Mnuj21FMiDpuP35VvcG5HeK/OCZU5BUW8/nSHZzfrYmd1K2kzk1qc1Gv5rw9L41rBnmGcDamOnkzLPMyEXlARNr6I5AJDd+s3EVOXhGXpFgXzqq485wOgA3lYHzDm149F+C5YGuyiPwqIneLiHXRMOWavCidFkk1GJBcz+0oQalZYg2uHdSaaUsyWLvLpr8w1cubXj1bVfVvqtoHuALoDmzxeTITtNL35TJv017G9GlBRISd1K2qmwa3pVZsFE9/vdbtKCbEeDUouoi0FpF7gY+ATsC9Pk1lgtqURemIwOg+1sxzIhLjY7h5SDtmrcti/iab7dRUH2/a+BcA04BI4BJV7aeqz/k8mQlKxSXK1NQMTmvfgKZ2UvKEjR3UmqZ14hj/9RpUbSgHUz28OeIfq6q9VfUpVbWpgky5ft64hx0H8hiTYpOtVIe46EjuOLsDyzIOMGPFLrfjmBBR3rDMV6nqe8D5InJ+2X9X1ed9mswEpY8XpVM3PpqzuoT3nLrV6aLezXnzpy088+1azjmpEdFhPm2lOXHlfYISnNtax1mM+R/ZuQV8v2o3F/ZqRmxUpNtxQkZkhHDfeR1J25vLRwu3uR3HhIDyLuB63bl91H9xTDCbvmwHBcUlXGwndavdkI4N6Z+cxAszNzCqd3NqxlY4eZ4xx1VeU8+L5b3QZuAyZU1ZlEGXJrU5qWkdt6OEHBHhgfM7c+HLP/PmT5u5/awObkcyQay8pp5UbAYu46W1uw6yYvsBO9r3oZ4tEjm/W2PemLOZzJw8t+OYIGYzcJlq8UlqBlERwsieNkePL91zbicKikp4ceYGt6OYIGYzcJkTVlhcwqdLdnBm54bUqxnrdpyQllw/gSv6t+TDhelsyjrkdhwTpGwGLnPC5qzPYs+hfC7uY333/eHWM9sTFxXBM9+sczuKCVI2A5c5YVNTM6iXEMPgjg3cjhIW6teM5Y+nt+WbVbtI3brP7TgmCHl7JUgkkAXsBzqIyGm+i2SCyb7DBfywxtN33y4s8p9xpybTsFYsT85Ya0M5mEqrsDOwiDwNXAqsAkqc1QrM8WEuEySmL91OYbFabx4/i4+J4o6zO/DAtBV8t3o3557U2O1IJoh4cxXIhUBHVc33dRgTfKYuzuCkprXp3KS221HCziV9mjNh7hae/notZ3RqaN+4jNe8+aRsBio9d56IvCUimSKystS6v4jIdhFZ6iy/GQPIBI81Ow+ycvtBRve2o303REVGcP/QTmzec9iGcjCV4s0Rfy6wVERmAv856vfiyt2JwD+Bd8qs/7uq2mTtIeCT1AyiI4ULe9lk6m45s3NDBrRJ4u8/bGBkr2bUjrP5jU3FvDninw48Bszjv1fzplb0IlWdA1iXgxBVWFzCZ0t3cEanhiQlxLgdJ2yJCA+e34V9hwt4bfYmt+OYIFHhEb8Pum7eIiLXAIuAu1R1/7GeJCI3ADcAtGxpU/wGmqN9962Zx33dmtdhVK9mTJi7hasGtLIJcEyFjnvELyKTndsVIrK87FLF/b0KtMUz5MNO4LgzeanqG6qaoqopDRpY//BA88liT9/9IZ1s3P1AcPe5HVHg2e/soi5TsfKO+G9zbodX185UdffR+yLyL+DL6tq28Z/s3AJ+WJ3JVQNaWU+SANEssQbXn5zM63M2cf3JyXRtZiOkmuMrb5C2nc7t1mMtVdmZiDQp9XAUsPJ4zzWB6wtn3P3RfeykbiC5aUhb6sbH8MRXNj+vKZ/PDtdE5ENgPtBRRDJEZBzwt6NNR8AQ4A5f7d/4ztTUDDrbuPsBp3ZcNLef1Z75m/fy/erdFb/AhC2fTeOjqpcfY/UEX+3P+MeG3TksyzjAw8O7uB3FHMMV/VryzvytPPX1WgZ3bEhMlDXFmd+yT4WplKmLbdz9QBYVGcGDwzqzZc9h3v2lSi2yJgyUN/XiCjxj8hyTqnb3SSITsIqKS/h08XYGd2xAfRt3P2AN7tCAU9vX54Uf1nNRr2bUtessTBnlHfEPB0YA3zjLlc4yA5jq+2gm0Py0YQ+ZOTbufqATER4a1oVD+UW8YDN1mWMor1fP0d47J6vqvaq6wlnuB871X0QTKKakppOUEMMZ1nc/4HVsXIvL+7Xk3V+2sjHTZuoy/8ubNv4EETnl6AMRGQQk+C6SCUT7D3v67o/s2dROGAaJO87uQHx0JE98tdrtKCbAePMXPA54WUTSRGQL8ApwvW9jmUAz3em7f4k18wSN+jVjufXM9sxal8WstZluxzEBxJupF1NVtQfQHeipqj1VdbHvo5lAMiU1nZOa1qZLUxt3P5iMHdSaNvUTeOzL1RQUlVT8AhMWKiz8ItJIRCYAH6vqARHp4lyMZcLE0XH3bZat4BMTFcHDw7uwec9hJs1LczuOCRDeNPVMBL4FjnbcXg/c7qtAJvBMdcbdH9nThmgIRkM6NWRIxwa8OHMDWTk2kZ7xrvDXV9XJOPPtqmoRUOzTVCZgFBaX8NmS7ZzZqZGNux/EHh7ehSOFxTz7rY3eabwr/IdFpB7OxVwiMgA44NNUJmDMWpvJ3sMFXJJizTzBrE2Dmlx3cmsmp6azPCPb7TjGZd4U/jvxzMLVVkR+xjOVYkXTLpoQMSU1g/o1Yzm9g82JEOxuPbM99RJieeTzVZSU2Oid4cybwr8KOB0YBPwBOAlY68tQJjBkHszj32szGd2nGVE27n7QqxUXzZ/P78TS9GympKa7Hce4yJu/5vmqWqSqq1R1paoW4hlu2YS4TxZvp7hEGZNiffdDxahezejbui5Pf7OO7NwCt+MYl5Q39WJjEekD1BCRXiLS21kGA/F+S2hcoapMWZRO39Z1adugpttxTDURER69oCvZuQU2TWMYK288/nOBa4HmwPOl1h8E/uzDTCYA/Jq2n817DnPTkHZuRzHVrEvT2lwzsDWT5qdxWd+WNk1jGCpvkLZJqjoEuFZVh5RaRqrqND9mNC74+Nd0asZGcX63xm5HMT5wx9kdqJcQw8Ofr7QTvWHImzb+PiKSePSBiNQVkcd9mMm47GBeIV+t2MGIHk2Jj/HZJG3GRXVqRHP/eZ1Zsi2byYvsRG+48abwn6eq/+n4q6r7gfN9F8m47YtlO8grLOHSvnZSN5SN7t2MfslJPPX1WvYcsit6w4k3hT9SRP4z3ZKI1ABs+qUQNvnXdDo1rkWP5tb2G8pEhCdHdSW3oIgnv1rjdhzjR94U/veAmSIyTkSuB74HJvk2lnHLmp0HWZZxgDEpLRARt+MYH2vXsBZ/OK0t05ZsZ97GPW7HMX7izbDMfwMeBzrjuXjrMWedCUEf/5pOTGQEo3rZgGzh4pYz2tGqXjwPfbaSvEIbhisceHs55hrgG1W9C/hJRGr5MJNxyZGCYj5ZnMHQro1tgu4wEhcdyWMju7J5z2Fenb3J7TjGD7wZj//3eCZXf91Z1Qz4zJehjDu+XL6DnLwirujf0u0oxs9O69CAET2a8ursTTZHbxjw5oj/ZuBkPBduoaobAJttOwR9sHAbbRsk0D85ye0oxgWPDO9CjZhI7v9kufXtD3HeFP58Vf3PoB4iEoUzRLMJHWt2HmTJtmwu79fSTuqGqQa1Ynl4eBcWbd3Pewu2uh3H+JA3hf9HEfkznjF7zgamAF/4Npbxtw8WbCMmKsKmVwxzo3s349T29Xn667Vszz7idhzjI94U/vuBLGAFnmGZZwAP+TKU8a/cgiI+W7KdYd2akBhvJ3XDmadvfzcUePDTFajal/tQVN7onDOdu0+p6r9U9RJVvdi5b5+GEPLFsh3k5NtJXePRIimeu8/pyOx1WXy+dIfbcYwPlHfE30RETgcuKDMsc28R6e2vgMb3PliwjfYNa5LSqq7bUUyAGDuoNb1aJvLoF6tsgvYQVF7hfwRPM8/RYZmfK7U86/toxh9Wbj/AsowDXNHfTuqa/4qMEJ65uDuHC4p56DNr8gk15Q3LPFVVzwP+VmZY5iGqeoYfMxofeu+XrcRFR3BRLzupa/5Xu4a1uPucDny7arc1+YSYCsfcVdXHRKQZ0Kr081V1ji+DGd/Lzi3gs6XbGdWrGXXio92OYwLQuFPa8O2q3Tzy+UoGtq1Ho9pxbkcy1cCbK3fHAz/j6clzj7Pc7eNcxg8mL0onr7CEawa2djuKCVBHm3wKikt4YJo1+YQKb2bZGAV0VFU7wxNCikuUd3/ZSr/kJDo3qe12HBPA2jSoyb3nduKvX65mSmoGY1JsnoZg500//s1ApdsBROQtEckUkZWl1iWJyPcissG5tW4kLpm9LpP0fUcYa0f7xgvXDmpN/+Qk/vrFatL35bodx5wgbwp/LrBURF4XkRePLl68biIwtMy6+4GZqtoemOk8Ni6YOC+NxrXjOOekRm5HMUEgIkJ4bkwPBLjj46UU21g+Qc2bwj8deAyYB6SWWsrlnPzdV2b1SP47icsk4EKvk5pqsynrED9t2MOV/VsSHentyNwm3DWvG89jF3Zl0db9vDp7o9txzAnwpldPdc621UhVdzrb3Skixx3lU0RuAG4AaNnSriitTu/O30pMZASX25W6ppJG9mzKzLWZ/OOHDZzavgE9WiS6HclUQXlDNkx2bleIyPKyi6+DqeobqpqiqikNGjTw9e7CxqH8IqamZjCsexPq17Spk03liAiPX9iVhrViuf3jpRzOL3I7kqmC8r7n3+bcDgdGHGOpit0i0gTAuc2s4nZMFU3+NZ1D+UVcM7CV21FMkKpTI5rnxvQkbe9hHvtytdtxTBWUd+Xu0SaZrcdaqri/6cBY5/5Y4PMqbsdUQVFxCW/9vIWUVnXp1dI6VJmqG9i2Hjee3paPfk1n+jK7qjfY+OzMnoh8CMwHOopIhoiMA8YDZ4vIBuBs57Hxk29X7SZj/xF+d2obt6OYEHDH2R3o06ouf562grQ9h92OYyrBZ4VfVS9X1SaqGq2qzVV1gqruVdUzVbW9c1u214/xEVXlXz9tplW9eM7uYl04zYmLjozgxct7ERkh3PLhYvKLit2OZLxU4Xj8IvK0/+IYX0ndup+l6dmMOyWZyAgbhdNUj2aJNXjm4u6s3H6Qp2asdTuO8VJ53TlLj8f/EfA/1UJVF/s0malWb8zZTGJ8tE2taKrdOSc15tpBrZk4L40BbeoxtGtjtyOZCt1IfW4AABaDSURBVJRX+MuOx1+aAjY0c5DYsucw36/Zzc2D2xEf483wTMZUzgPnd2Lxtv3cM2UZHRvXIrl+gtuRTDlsPP4w8NbcLURHRHDNIOvCaXwjNiqSV67sTWSkcON7qeQWWP/+QFbhyV1nPP4LRORZZxnuj2Cmeuw7XMCU1HRG9mxKw1o2lrrxneZ143nhsl6s253Dg5+utCGcA5g34/E/hedirtXOcpuzzgSBt+ZuIb+ohBtOsy6cxvdO79CAO87qwKdLtvPuL1W93Mf4mjcNvsOAnqpaAiAik4AlwAO+DGZO3IEjhUyal8bQkxrTvlEtt+OYMHHLkHYsTc/msS9Xc1LT2vRpleR2JFOGt/34S4/EVMcXQUz1e3d+Gjn5Rdw8pJ3bUUwYiYgQ/j6mJ00Ta/CHdxezI/uI25FMGd4U/qeAJSIy0TnaTwWe9G0sc6JyC4qYMHcLQzo2oGsz+7/a+Fed+GjevCaFvMJibnh3EUcK7OKuQOLNyd0PgQHANGcZqKof+TqYOTEfLNjG/txCbjmjvdtRTJhq36gWL1zWk1U7DnLP1GV2sjeAeNXUo6o7VXW6qn6uqrt8HcqcmLzCYt6Ys5lBbevRp5UNxmbcc2bnRtx7bie+XL6TV2ZvcjuOcdj0SyFoSmoGmTn53GJt+yYA/PH0NlzYsynPfLuOr1fsdDuOwQp/yCkoKuG12Zvo3TKRgW3ruR3HGESE8aO707tlIrd/vJTF2/a7HSnslVv4RSRCRFb6K4w5cR/9uo3t2Ue47awOiNhgbCYwxEVH8q9rUmhcJ47fT1rE1r02jLObyi38Tt/9ZSJik7MGgSMFxbz07430S07itPb13Y5jzP+oVzOWt6/tS7Eq1739K9m5BW5HClveNPU0AVaJyEwRmX508XUwU3mT5qeRlZPPPed2tKN9E5DaNKjJG1enkLH/CDe8k0peoXXzdIM3V+4+6vMU5oQdzCvk1dmbGNKxAX1b25WSJnD1S07i2TE9uPXDJdz64RJeubI3UZF2utGfvOnH/yOQBkQ7938FbCz+APPmT1s4cKSQu87p6HYUYyp0QY+m/N+ILny3ercN6OaCCo/4ReT3wA1AEtAWaAa8Bpzp22jGW3sP5TPhp80M69bErtI1QeO6k5PZe6iAf87aSFLNGO4b2sntSGHDm6aem4F+wAIAVd0gIg19mspUyquzN3GksJg7zu7gdhRjKuWuczqwL7eAV2dvol5CDL871UaR9QdvCn++qhYcPVkoIlF4ZuAyAWDb3lzemb+V0b2b065hTbfjGFMpIsJjI7uSnVvA41+toUZMJFf2twmDfM2bwv+jiPwZqCEiZwM3AV/4Npbx1lNfryEqUrj7XGvbN8EpMkL4x6W9yCtM5cFPVxIdGcGYlBZuxwpp3pxKvx/IAlYAfwBmAA/5MpTxzi+b9/L1yl3cNLgtjWrb7FomeMVERfDKlb05tX197vtkOZ8v3e52pJBW4RG/qpY4wzEvwNPEs07tFLzrikuUx75cTbPEGtYuakJCXHQkb1ydwnUTF3Ln5GVERUQwrHsTt2OFJG+mXhwGbAJeBP4JbBSR83wdzJTvk9QMVu04yH3ndSIuOtLtOMZUixoxkUwY25deLRK59aMlTF+2w+1IIcmbpp7ngCGqOlhVTweGAH/3bSxTnkP5Rfzt23X0bpnICDsiMiEmITaKidf3I6VVXW7/aAlTUzPcjhRyvCn8maq6sdTjzUCmj/IYL7w8ayN7DuXzyIiTbGgGE5JqxkYx8bp+DGpbn3umLuODBdvcjhRSjtvGLyIXOXdXicgMYDKeNv5L8Fy9a1ywfncO/5qzmdG9m9OzRWLFLzAmSNWIieTNsSnc+F4qf/50BXmFxVx/SrLbsUJCeSd3R5S6vxs43bmfBdi0Ti4oKVH+PG0FteKieHBYZ7fjGONzcdGRvHZ1H279cAl//XI1+3MLuPNsG3L8RB238Kvqdf4MYir28aJ0Fm3dzzMXdycpIcbtOMb4RWxUJC9f0ZuHPlvJS//eyJ5DBTx+YVciI6z4V5U3Y/UkA38CWpd+vqpe4LtYpqysnHyemrGG/slJXNynudtxjPGrqMgInrqoG/VqxvDyrE3sP1zAPy7raT3aqsibK3c/AybguVq3xLdxzPE8/tVq8gpLeGJUN/uaa8KSiHDPuZ1ISojlsS9Xc/WEBbx+dYp9+60Cbwp/nqq+6PMk5rh+XJ/F50t3cNuZ7W08HhP2xp2STMNasdw1ZRkXvfIzb1/Xj+T6CW7HCiredOd8QUT+T0QGikjvo4vPkxkADhwp5P5PltO2QQI3Dm7rdhxjAsKIHk358Pf9OZhXxKhXfmbhln1uRwoq3hT+bsDvgfF4LuZ6DnjWl6HMf/1l+ioyc/J5foy1ZxpTWp9WSXx60yCSEmK46s0FdqFXJXjT1DMKaKOq1TYzsoikATlAMVCkqinVte1QMmPFTj5dsp3bz2pPD+uzb8xvtKqXwLQbB3HT+4u5e8oyVm4/wIPDOhNtUzmWy5t3Zxngi6ozRFV7WtE/tsyDeTz46Qq6N6/DzUPauR3HmICVGB/DO9f34/qTk5k4L41rJixk3+FqO04NSd4U/kbAWhH5VkSmH118HSycqSr3fbKc3IJinh/T045ejKlAVGQEj4zownOX9CB1235GvDSXFRkH3I4VsLxp6vk/H+xXge9ERIHXVfWNsk8QkRvwzPVLy5YtfRAhcL37y1Zmrcvi/0Z0sV48xlTC6D6emehufC+V0a/O4+HhnblqQCvrAl2GuDG0vog0VdUdzty93wN/UtU5x3t+SkqKLlq0yH8BXbQ0PZtLXpvHKe3qM2FsXyLs6kRjKm3/4QLunLyUWeuyGN69CeNHd6dmrDfHuaFFRFKP1ZzuzXj8OSJy0FnyRKRYRA6eSBhV3eHcZgKf4pnMPeztP1zAze8vpmGtOP5+aU8r+sZUUd2EGCaM7cu9QzsyY8VORrw0l+UZ2W7HChgVFn5VraWqtZ0lDhiNZ0KWKhGRBBGpdfQ+cA6wsqrbCxUlJcqdk5eSmZPHK1f2JjHerkY05kRERAg3DW7Hh78fwJGCYi56ZR6vzt5EcYlNIFjps4aq+hlwxgnssxEwV0SWAQuBr1T1mxPYXkh49cdNzFqXxSPDu1jXTWOqUf829fjm9lM556RGPP3NWq588xd2ZB9xO5arvBmk7aJSDyOAFDwnZ6tEVTcDPar6+lA0e10mz323jgt6NOWqAa3cjmNMyEmMj+HlK3ozNTWDv0xfxbn/mMMjw7twcZ/mYXni15uzHaXH5S8C0oCRPkkThtbuOsgtHyyhU+PaPHWRDcBmjK+ICJektKBfchL3TFnOPVOXM2PFTp66qDuN68S5Hc+vXOnVU1mh2qsnMyePC//5M8WqfHbzyTSpU8PtSMaEhZISZdL8NJ7+Zi3RkRE8PKwLl6SE3tH/8Xr1lDf14iPlbE9V9bFqSRamjhQU8/tJi9ifW8iUPw60om+MH0VECNednMyQjg25d+py7v1kOVMXZ/DkqK60a1jL7Xg+V97J3cPHWADGAff5OFdIKy5R7vh4Kcu3H+CFy3rStVkdtyMZE5Za10/goxsGMP6ibqzblcN5L/zEc9+tI6+w2O1oPuVVU4/T/fI2PEV/MvCc0wffL0KpqaekRLl/2nImL8rgoWGd+d2pbdyOZIwB9hzK54mv1vDpku20SKrBQ8O6cE6XRkHd/FOlC7hEJElEHgeW42kW6q2q9/mz6IcSVeWvX65m8qIM/nRGOyv6xgSQ+jVj+fulPXn/d/2pER3JH95N5eoJC9mwO8ftaNXuuIVfRJ4BfsUzfHI3Vf2Lqu73W7IQ9My365g4L41xpyRz59kd3I5jjDmGk9vVZ8atp/KXEV1YnpHN0Bd+4pHPV7LnUL7b0arNcZt6RKQEyMfThbP0kwTPyd3avo/nEQpNPS/N3MBz36/niv4teeLCrkH99dGYcLH3UD5//2E9Hy5MJy4qgj+e3pZxpyYTHxMc4/4cr6nHunP6mKryt2/X8ersTVzUqxnPXtLDxuAxJshsyjrE375Zy7erdtOwVix/OqMdY/q2IDYqsGfFs8LvguIS5ZHPV/L+gm1c0b8lj43sSqQVfWOC1qK0fTz9zVp+TdtPs8Qa3HJGOy7u0zxg58ywwu9nBUUl3DVlGV8s28GNg9ty77kdrXnHmBCgqvy0YQ/Pf7+epenZtEyK58bBbbmod7OA+wZghd+PDhwp5JYPFvPThj3cf14n/nh6W7cjGWOqmaoya10m//hhA8szDtCodiy/P7UNl/drSUKAjP1vhd9Ptuw5zLhJv7Jtby5PjurGmL4t3I5kjPEhVeXnjXt5edZG5m/eS2J8NFf2b8k1A1vTqLa7YwBZ4feDnzZkcfP7i4mKjODVK3vTv009tyMZY/wodet+3pizie9W7yYqQhjRoynXn5zs2tX5Vvh9SFWZMHcLT329lnYNavLm2BRaJMW7HcsY45Ktew/z9s9pTF6UTm5BMX1a1eXqAa04r1tjv54HsMLvI3sP5XP3lGXMWpfFOV0a8fylPcNybk9jzG8dyC1kSmo67/2ylbS9udRLiGFM3xaMSWlBcv0En+/fCr8PzNu0h9s/Wkp2biEPDuvMNQNbWc8dY8xvlJQoczfu4d1ftvLvtZkUlyj9kpO4rG8LzuvahBoxvvkWYIW/Gh0pKOb579fx5twtJNdP4KXLe3FSUxth0xhTscyDeUxdnMHkX9NJ25tLQkwkQ7s2YVSvZgxsW69ar/Wxwl9N5m3cw/3TVrBtXy6X92vJw8M7B83l28aYwKGqLNyyj2mLtzNjxU5y8otoVDuW4d2bMqx7E3q1SDzhFgQr/Cdo/+ECnv5mLR/9mk7revE8dVF3Bra1XjvGmBOXV1jMD2t289mS7cxZv4eC4hKaJdZgWPcmXNW/FS3rVa2zSKVn4DIeBUUlvPvLVl74YT2HC4r5w+ltuOOsDsRFB9YVesaY4BUXHcnw7k0Z3r0pB44U8v3q3Xy1fAdvzd3CmZ0aVrnwH48V/uNQVWauyeTJGWvYvOcwp7avz0PDutCxcehPy2aMcU+dGtFc3Kc5F/dpTnZuAbXioqt9H1b4y1BV5mzYwz9+WM+Sbdm0aZDA29f2ZXDHBtZjxxjjV4nxMT7ZrhV+h6oye30WL87cwJJt2TRLrMETo7oyJqVFwI68Z4wxVRH2hT+vsJhpi7fz9s9b2JB5iGaJNXhyVDcu7tOcmCgr+MaY0BO2hX/LnsNMXpTORwu3sT+3kJOa1ub5MT0Y3r2pFXxjTEgLq8KfW1DEjBW7mLwonYVb9hEhcFbnRow7JZl+yUnWhm+MCQshX/iPFBQze10mXy7fycy1u8krLCG5fgL3Du3I6N7NXR821Rhj/C2kC/+LMzfw2o+byC0opn7NGC7p04ILejYlpVVdO7o3xoStkC78jevEcWGvZgzv1oT+bap3DAxjjAlWIV34x6R4hj81xhjzX9Z9xRhjwowVfmOMCTNW+I0xJsxY4TfGmDDjSuEXkaEisk5ENorI/W5kMMaYcOX3wi8ikcDLwHlAF+ByEeni7xzGGBOu3Dji7wdsVNXNqloAfASMdCGHMcaEJTcKfzMgvdTjDGfd/xCRG0RkkYgsysrK8ls4Y4wJdW5cwHWsy2d/M/Gvqr4BvAEgIlkisrWK+6sP7Knia30tULMFai4I3GyBmgsCN1ug5oLAzVbZXK2OtdKNwp8BlL6ctjmwo7wXqGqDqu5MRBYda7LhQBCo2QI1FwRutkDNBYGbLVBzQeBmq65cbjT1/Aq0F5FkEYkBLgOmu5DDGGPCkt+P+FW1SERuAb4FIoG3VHWVv3MYY0y4cmWQNlWdAczw0+7e8NN+qiJQswVqLgjcbIGaCwI3W6DmgsDNVi25RPU351WNMcaEMBuywRhjwowVfmOMCTMhXfgDaUwgEXlLRDJFZGWpdUki8r2IbHBu67qQq4WIzBKRNSKySkRuC4RsIhInIgtFZJmT61FnfbKILHByfez0DPM7EYkUkSUi8mWA5UoTkRUislREFjnrXP+cOTkSRWSqiKx1Pm8D3c4mIh2d9+roclBEbnc7V6l8dzif/5Ui8qHzd3HCn7WQLfwBOCbQRGBomXX3AzNVtT0w03nsb0XAXaraGRgA3Oy8T25nywfOUNUeQE9gqIgMAJ4G/u7k2g+M83Ouo24D1pR6HCi5AIaoas9S/b3d/l0e9QLwjap2Anrgef9czaaq65z3qifQB8gFPnU7F4CINANuBVJUtSueXpCXUR2fNVUNyQUYCHxb6vEDwAMuZ2oNrCz1eB3QxLnfBFgXAO/b58DZgZQNiAcWA/3xXLUYdazfsR/zNMdTDM4AvsRzNbrruZx9pwH1y6xz/XcJ1Aa24HQoCaRspbKcA/wcKLn47/A2SXh6YH4JnFsdn7WQPeLHyzGBXNZIVXcCOLcN3QwjIq2BXsACAiCb05yyFMgEvgc2AdmqWuQ8xa3f6T+Ae4ES53G9AMkFnuFPvhORVBG5wVnn+u8SaANkAW87TWRvikhCgGQ76jLgQ+e+67lUdTvwLLAN2AkcAFKphs9aKBd+r8YEMh4iUhP4BLhdVQ+6nQdAVYvV8xW8OZ5RXTsf62n+zCQiw4FMVU0tvfoYT3Xrs3ayqvbG08R5s4ic5lKOsqKA3sCrqtoLOIx7TU6/4bSTXwBMcTvLUc55hZFAMtAUSMDzey2r0p+1UC78lR4TyAW7RaQJgHOb6UYIEYnGU/TfV9VpgZQNQFWzgdl4zkEkisjRCw/d+J2eDFwgIml4hhQ/A883ALdzAaCqO5zbTDxt1f0IjN9lBpChqgucx1Px/EcQCNnAU1AXq+pu53Eg5DoL2KKqWapaCEwDBlENn7VQLvzBMCbQdGCsc38snvZ1vxIRASYAa1T1+UDJJiINRCTRuV8Dzx/BGmAWcLFbuVT1AVVtrqqt8Xym/q2qV7qdC0BEEkSk1tH7eNqsVxIAnzNV3QWki0hHZ9WZwOpAyOa4nP8280Bg5NoGDBCReOfv9Oh7duKfNbdOpPjp5Mj5wHo8bcMPupzlQzztdIV4jn7G4WkbnglscG6TXMh1Cp6visuBpc5yvtvZgO7AEifXSuARZ30bYCGwEc/X8lgXf6eDgS8DJZeTYZmzrDr6mXf7d1kqX09gkfM7/QyoGwjZ8HQe2AvUKbXO9VxOjkeBtc7fwLtAbHV81mzIBmOMCTOh3NRjjDHmGKzwG2NMmLHCb4wxYcYKvzHGhBkr/MYYE2as8BtThjOK5E3O/aYiMtXtTMZUJ+vOaUwZzphFX6pnRERjQo4rc+4aE+DGA22dAeI2AJ1VtauIXAtciGd43K7Ac0AMcDWeYaTPV9V9ItIWz5DgDfAM8/t7VV3r/x/DmGOzph5jfut+YJN6Boi7p8y/dQWuwDMGzhNArnoGHZsPXOM85w3gT6raB7gbeMUvqY3xkh3xG1M5s1Q1B8gRkQPAF876FUB3Z5TTQcAUz/AqgOcye2MChhV+Yyonv9T9klKPS/D8PUXgGS+9p7+DGeMta+ox5rdygFpVeaF65jLYIiKXgGf0UxHpUZ3hjDlRVviNKUNV9wI/i8hK4JkqbOJKYJyIHB0lc2R15jPmRFl3TmOMCTN2xG+MMWHGCr8xxoQZK/zGGBNmrPAbY0yYscJvjDFhxgq/McaEGSv8xhgTZv4fAivD169uwGkAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "Wupdate_sir=Wtt_sir.subs([(N,100),(tau,0.5),(c,0.5),(a,0.1)])\n",
    "Rupdate= Rtt.subs([(a,0.1)])\n",
    "ws_sir=[1]\n",
    "Rt_v = 0\n",
    "for t in range(0,79):\n",
    "    ws_sir.append(Wupdate_sir.subs([(Wt,ws_sir[-1:][0]),(Rt,Rt_v)]))\n",
    "    Rt_v = Rupdate.subs([(Wt,ws_sir[-1:][0]),(Rt,Rt_v)])\n",
    "plt.plot(range(0,80),ws_sir)\n",
    "plt.ylabel('Number of infected individuals')\n",
    "plt.xlabel('time')\n",
    "plt.title('Diffusion of infection');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following image uses the qualitative conclusions from this model to argue why protective measures during an outbreak (like washing hands, quarantine, canceling events) is heplful. It delays the spread of the disease and reduces the height of the peak. This reduces the risk of overrunning the capacity of the health care system (which, itself, can increase spread) and giving it more time to prepare. \n",
    "\n",
    "![The value of protective measures](https://pbs.twimg.com/media/ESas8tHVAAAhr_I?format=jpg&name=4096x4096)\n",
    "\n",
    "See the [following thread](https://twitter.com/CT_Bergstrom/status/1235865328074153986) for discussion."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAEWCAYAAABhffzLAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nOzdd3hUxfrA8e+bQkIJHULoHQKk0JvA0kEBKYIiKFjAhh2Va0G9WLj+QAWVi2IBQSyg0kSQIE2KECCR3gMEkpAQEhLSk/n9cTZ7A6Rskt1Nm8/znCc5u+fMvLuE2bNzZt4RpRSapmla2eFU1AFomqZpjqUbfk3TtDJGN/yapmlljG74NU3Tyhjd8GuappUxuuHXNE0rY3TDrzmUiCwUkTey7D8hIhEiEi8iNUSkp4icMu+PLEQ9v4vIJNtEna963xGRKBEJz+a5XiJyIh9l3fTe2DDGCSLyh63K00oe0eP4NVsRkRDAE0gD0oGjwLfAF0qpjGyOdwWuA92UUsHmxzYDa5RS8xwVt62ISAPgJNBIKXWlkGXd9t4UsJzGwDnAVSmVVpiYtNJDX/FrtjZcKeUBNAJmA68AX+VwrCfgDhzJ8lijW/ZLkkbA1cI2+mbZvTeaZhO64dfsQikVq5RaA9wLTBKRdgAistjcHdISyOz2iBGRP0XkDNAUWGvu3nATkRARGZBZroi8JSLLzL+7i8gyEbkqIjEisk9EPM3PbRWRR82/O4nI6yJyXkSuiMi3IlLF/FxjEVEiMklELpi7aV7L6XWJSBXz+ZHm8l43lz8A2ATUNce+OJtzTSISmmU/RESmi8g/IhIrIj+aX9Nt7435+NYisklEokXkhIiMy1JWeRGZa44pVkT+EpHywPYs5cSLSHcRmSwif2U5t4f5vYs1/+yR5bmtIjJLRHaKSJyI/CEiNXP5p9dKAN3wa3allNoLhAK9bnn8JNDWvFtVKdVPKdUMuIDxraGSUio5j+InAVWABkAN4HEgMZvjJpu3vhgfLJWAT2855g6gFdAfmCki3jnU+Ym5zqZAH+BB4CGlVAAwFLhsjn1yHrFnGgcMAZoAvsDk7N4bEamI8cGyHKgNjAcWiEjmcXOAjkAPoDrwMpAB9M5STiWl1O6slYtIdeA3YD7Ge/gh8Nst9xTuBx4y11sOmG7la9OKKd3wa45wGaMxsrVUjMaquVIqXSm1Xyl1PZvjJgAfKqXOKqXigX8B94mIS5Zj3lZKJZr704MBv1sLERFnjG8w/1JKxSmlQoC5wAOFeA3zlVKXlVLRwFrAP4fjhgEhSqlvlFJpSqkDwM/APSLiBDwMPKuUumR+L3ZZ8cEJcBdwSim11Fzu98BxYHiWY75RSp1USiUCP+USo1ZC6IZfc4R6QLQdyl0KbAR+EJHLIvKB+aboreoC57PsnwdcMPrRM2UdhZOA8a3gVjUxrnhvLateAWLPT71g3D/oau7SihGRGIwPtDrmuNyBMwWo/9b3Bm5/TdbGqJUQuuHX7EpEOmM0In/ldWwObgAVsuzXyfxFKZWqlHpbKdUGo4tjGEbXy60uYzScmRpijDyKyGcsURjfMm4t61I+yymIi8A2pVTVLFslpdQT5riSgGbZnJfXsL1b3xtw3GvSiohu+DW7EJHKIjIM+AFYppQ6VMCigjC6ZVxFpBNwT5Y6+oqIj7kL5jpGo5yeTRnfA8+LSBMRqQS8B/yY3+GNSql0jK6Od0XEQ0QaAS8Aywr0yvJnHdBSRB4wvxeuItJZRLzNQ2W/Bj4Ukboi4my+iesGRGL09TfNodz15nLvFxEXEbkXaGOuTyuldMOv2dpaEYnDuEJ9DeNm4UOFKO8NjCvZa8DbGDc3M9UBVmI0+seAbWTfCH+N0S20HWNMexLwdAHjeRrjW8hZjG8xy83l25VSKg4YBNyHcZUeDvwHcDMfMh04BOzD6Fb7D+CklEoA3gV2mruIut1S7lWMb0ovAlcxbgoPU0pF2fs1aUVHT+DSNE0rY/QVv6ZpWhmjG35N07QyRjf8mqZpZYxu+DVN08oYl7wPKXo1a9ZUjRs3LuowNE3TSpT9+/dHKaVq3fp4iWj4GzduTGBgYFGHoWmaVqKIyK2zsgHd1aNpmlbm6IZf0zStjNENv6ZpWhlTIvr4Na20S01NJTQ0lKSkpKIORSuB3N3dqV+/Pq6u2SWnvZ1u+DWtGAgNDcXDw4PGjRsjIkUdjlaCKKW4evUqoaGhNGnSxKpzdFePphUDSUlJ1KhRQzf6Wr6JCDVq1MjXt0Xd8GtaMaEbfa2g8vu3o7t6ypjw2CTW/XOZyu6uNKtdCf8GVXF20g2OppUl+oq/jLielMrzPwZxx3/+5J3fjvHyz/8w5r+7GP/FHkKvJRR1eFoRu3jxIn379sXb25u2bdsyb948y3PR0dEMHDiQFi1aMHDgQK5duwbAzz//TNu2benVqxdXr14F4MyZM9x33312i/Py5cvcc889uR4TEhJCu3bt8ixr8uTJrFy5EoBHH32Uo0eP5njs4sWLuXz5cv6CtUJe9dqLbvjLgISUNB7+Zh9rgy8zqUdjtk43se0lE++MbMfRsOsMnbeDXaf1uhtlmYuLC3PnzuXYsWPs2bOHzz77zNIgzZ49m/79+3Pq1Cn69+/P7NmzAZg7dy579uzhwQcfZPlyY32c119/nVmzZtktzrp161oaa1v68ssvadOmTY7P26vhz6tee9ENfymXkpbBY0v3c+DCNeaPb88bw9rQuGZFGtWoyMRujVj/TC/qVHbnie8OcC7qRlGHqxURLy8vOnToAICHhwfe3t5cumQsu7t69WomTZoEwKRJk1i1ahUATk5OJCcnk5CQgKurKzt27MDLy4sWLVpkW0d6ejqTJ0+mXbt2+Pj48NFHHwEQFBREt27d8PX1ZdSoUZZvFKdPn2bAgAH4+fnRoUMHzpw5c9PVfEhICL169aJDhw506NCBXbt25foalVJMmzaNNm3acNddd3HlyhXLcyaTicDAwGxjXLlyJYGBgUyYMAF/f38SExP597//TefOnWnXrh1Tp04lc0Erk8nEK6+8QpcuXWjZsiU7duywvPbp06fj4+ODr68vn3zyyU31AlSqVInXXnsNPz8/unXrRkSEsST0mTNn6NatG507d2bmzJlUqlT4te51H38pt2jHWXaciuKDe3y508frtucb1qjA15M7M+LTv5jybSC/PtkDD3frxgJrdvL7DAgv6BLFOajjA0NnW3VoSEgIBw8epGvXrgBERETg5WX87Xh5eVkazDfffJPBgwdTt25dli1bxrhx4/jhhx9yLDcoKIhLly5x+PBhAGJiYgB48MEH+eSTT+jTpw8zZ87k7bff5uOPP2bChAnMmDGDUaNGkZSUREZGxk2Nde3atdm0aRPu7u6cOnWK8ePH55rT69dff+XEiRMcOnSIiIgI2rRpw8MPP5xnjFWrVuXTTz9lzpw5dOrUCYBp06Yxc+ZMAB544AHWrVvH8OHDAUhLS2Pv3r2sX7+et99+m4CAAL744gvOnTvHwYMHcXFxITo6+rb4bty4Qbdu3Xj33Xd5+eWXWbRoEa+//jrPPvsszz77LOPHj2fhwoU5vr780Ff8pdilmEQ+/fM0g9t6Mq5TgxyPa1C9AgsmdCQk6gaz1jm+v1ErPuLj4xkzZgwff/wxlStXzvXYgQMHsn//ftauXcuqVau48847OXHiBPfccw9TpkwhIeHme0dNmzbl7NmzPP3002zYsIHKlSsTGxtLTEwMffr0AYxvFNu3bycuLo5Lly4xatQowJigVKFChZvKS01NZcqUKfj4+DB27Ng8+8q3b9/O+PHjcXZ2pm7duvTr1++2Y7KLMTtbtmyha9eu+Pj48Oeff3LkyBHLc6NHjwagY8eOhISEABAQEMDjjz+Oi4txrV29evXbyixXrhzDhg277dzdu3czduxYAO6///5cX6O19BV/KfbOuqMoFG8My7sPsXuzGjxyRxM+336W+7o0pEPDag6IUMuWlVfmtpaamsqYMWOYMGGCpfEC8PT0JCwsDC8vL8LCwqhdu/ZN5yUkJLBkyRI2btzIoEGDWL16NcuXL+e7775jypQpluOqVatGcHAwGzdu5LPPPuOnn36ydPfcypq1wD/66CM8PT0JDg4mIyMDd3f3PM/Ja9hjdjF+/fXXNx2TlJTEk08+SWBgIA0aNOCtt966aQy9m5sbAM7OzqSlpVleT151u7q6Wo7Jeq496Cv+UmrvuWh+PxzOtL7NqV+tQt4nAE/3b4FnZTfeXH2E9Iy8/+NppYdSikceeQRvb29eeOGFm54bMWIES5YsAWDJkiXcfffdNz3/wQcf8Oyzz+Lq6kpiYiIigpOT021X/FFRUWRkZDBmzBhmzZrFgQMHqFKlCtWqVbP0hS9dupQ+ffpQuXJl6tevb7mfkHkvIavY2Fi8vLxwcnJi6dKlpKen5/oae/fuzQ8//EB6ejphYWFs2bLltmOyixGM+x5xcXEAlka+Zs2axMfHW3WzedCgQSxcuNDSmGfX1ZOTbt268fPPPwPk2pWWH7rhL6W+3HGWahVceeSOplafU8nNhVfv9ObQpVhWBF60Y3RacbNz506WLl3Kn3/+ib+/P/7+/qxfvx6AGTNmsGnTJlq0aMGmTZuYMWOG5bzLly8TGBho+TB48cUX6datG0uWLLmtW+LSpUuYTCb8/f2ZPHky77//PmB8mLz00kv4+voSFBRk6TtfunQp8+fPx9fXlx49ehAeHn5TeU8++SRLliyhW7dunDx5kooVK+b6GkeNGkWLFi3w8fHhiSeesHQvWRPj5MmTefzxx/H398fNzc3SxTRy5Eg6d+6c5/v76KOP0rBhQ3x9ffHz87OMgrLGxx9/zIcffkiXLl0ICwujSpUqVp+bE7HmK1VR69Spk9ILsVjv/NUbmOZs5SlTc6YPbpWvc5VSjFqwi8i4ZLa+ZMLVWV8bOMKxY8fw9vYu6jC0YighIYHy5csjIvzwww98//33rF69+rbjsvsbEpH9SqlOtx6r/1eXQt/sDMHFSXige6N8nysiPN2vOZdiElkdZPtxy5qm5c/+/fvx9/fH19eXBQsWMHfu3EKXqW/uljLXk1JZEXiRYb518ayc982u7PRrXRtvr8os2HqaUe3r6ZQOmlaEevXqRXBwsE3L1Ff8pcyaoMvcSElnco/GBS5DRHiqbzPORt5gw+HwvE/QbKIkdLtqxVN+/3Z0w1/KrAm6TPPalfCtX7gbQEPbedG4RgW+3nnORpFpuXF3d+fq1au68dfyLTMfvzXDWTPprp5S5FJMIntDonlxYMtCp/h1dhImdmvEO78d4+jl67Spm/tkHq1w6tevT2hoKJGRkUUdilYCZa7AZS3d8Jcia4ONm7Ej/OvapLyxHRsw548TLPv7PO+N8rFJmVr2XF1drV49SdMKS3f1lCKrgy7j36AqjWrkPp7ZWlUquDLcty6rDl7ielKqTcrUNK3o6Ya/lDgVEcexsOvcbaOr/UwPdG9EQko6vx64ZNNyNU0rOrrhLyU2HjFG39yVTQbOwvCtXxWfelX4Sc/k1bRSQzf8pcTm41fwq1+F2gUcu5+bezrW58jl6xwPv27zsjVNczzd8JcCUfHJBF2MoV9rT7uUP9yvLq7Ows/7Q+1SvqZpjqUb/lJg64lIlIL+3rXzPrgAqlcsR99WtVkVdJm09Ay71KFpmuPohr8U+PN4BJ6V3Whrx7H2ozvUJzIumR16bV5NK/F0w1/CpaRlsP1kFP1a1y70pK3c9Gtdm2oVXPlFj+7RtBJPN/wl3L6QaOKT0+zWv5+pnIsTQ328CDgaQUKK/VYG0jTN/vJs+EWkp4hUNP8+UUQ+FJH85/vV7GL7qUhcnYWezWvYva7hvnVJTE3nz+NX8j5Y07Riy5or/v8CCSLiB7wMnAe+tWtUmtV2nb5K+wbVqFDO/tk3ujSpTi0PN0tqCE3TSiZrGv40ZaQMvBuYp5SaB3jkdZKINBCRLSJyTESOiMiz5seri8gmETll/qlX9S6gmIQUDl+OpYcDrvbBSNx2l48XW05EEqdTOGhaiWVNwx8nIv8CJgK/iYgz4GrFeWnAi0opb6Ab8JSItAFmAJuVUi2AzeZ9rQD2nL2KUtCzeU2H1Tncz4uUtAw2HY1wWJ2aptmWNQ3/vUAy8IhSKhyoB/xfXicppcKUUgfMv8cBx8zn3g0sMR+2BBhZgLg1YOfpq1Qo54xf/aoOq7N9g2rUq1qedf+EOaxOTdNsK8+GXykVrpT6UCm1w7x/QSmVrz5+EWkMtAf+BjyVUmHmssIA+8w6KgN2nomiS5PqlHNx3OAsJydhcNs6/HUqivhkPbpH00qiHFsMEYkTkevZbHEiYnXSFhGpBPwMPKeUys95U0UkUEQC9eIUtwuPTeJs5A16NnNcN0+mwW09SUnPYOsJPbpH00qiHBt+pZSHUqpyNpuHUsqqKaIi4orR6H+nlPrF/HCEiHiZn/cCsm09lFJfKKU6KaU61apVK3+vqgzYdcaYQdu9mWNu7GbVqXF1alQsp9fj1bQSyuo+AhGpLSINMzcrjhfgK+CYUurDLE+tASaZf58ErM5PwJph77loKru70MbL8UsiOjsJA7w92XoikuS0dIfXr2la4VgzgWuEiJwCzgHbgBDgdyvK7gk8APQTkSDzdicwGxhoLnOgeV/Lp70h0XRqXB0nJ/ulacjNkHZ1iE9OY9fpq0VSv6ZpBWfNrJ9ZGMMxA5RS7UWkLzA+r5OUUn8BObVK/a0PUbtVVHwyZyNvMLZjgyKLoUfzGlRyc2HjkXD6ttb35zWtJLGmqydVKXUVcBIRJ6XUFsDfznFpuQgMiQagS5Oim/vm5uKMqVUtNh2NID1DFVkcmqblnzUNf4x5ZM524DsRmYcxOUsrInvPXcPNxQmfeo4bv5+dwW3rcPVGCvvPXyvSODRNyx9rGv67gUTgeWADcAYYbs+gtNztC4nGv0FVh47fz46pVS3KOTtZ1vvVNK1ksGYC1w2lVLpSKk0ptUQpNd/c9aMVgfjkNI5cjqVLk+pFHQoe7q70bF6DjUfCMdI5aZpWElgzqifrRK4kEUnPzwQuzbYOnL9GhoLOjYu+4Qejuyf0WiJHw/SfhKaVFNZc8WedyOUOjAE+tX9oWnYCz1/DSaBDo+KR1HRAG0+cBDYe0UnbNK2kyHcnsVJqFdDPDrFoVjh44Rqt6lSmkpv98+9bo2YlNzo1qs4fup9f00qMPFsPERmdZdcJ6AToDt0ikJGhCLoQwwj/ukUdyk0GtfXknd+OcTE6gQbVKxR1OJqm5cGaK/7hWbbBQBzGSB/NwU5diScuOY0ODYtHN0+m/t7Ger+bj+nuHk0rCfK84ldKPeSIQLS8HbhgjJcvLv37mZrUrEizWhUJOHaFyT2bFHU4mqblIceGX0Q+IZcuHaXUM3aJSMvRgfPXqF6xHI1rFL/ulAFtPPlqxzmuJ6VS2d2aBdo0TSsquXX1BAL7AXegA3DKvPkDOiVjEThw4RrtG1TFSHxavAz09iQtQ7HthF47QdOKu9zy8S9RSi0BWgB9lVKfKKU+wUiwpnP1OFhMQgpnIm8Uu26eTO0bVqN6xXIE6H5+TSv2rLm5WxfwyLJfyfxYiRASEkK7du1sUtbixYuZNm0aAKtWreLo0aOW50wmE4GBgTapJztBF2MAaN/QPvl53nvvvUKd7+wk9Gtdmy3Hr5CanmGjqApu5syZBAQEAPDxxx+TkJBQxBFpWvFhTcM/GzgoIotFZDFwAChcK1EK3Nrw20Na2v9y4R24EIOTYLeF1Qvb8AMM8PbkelIagSFFn7Tt3//+NwMGDAB0w69pt7Jm5u43QFfgV/PW3dwFVGKkp6czZcoU2rZty6BBg0hMTATgzJkzDBkyhI4dO9KrVy+OHz8OwNq1a+natSvt27dnwIABRETc3H2xa9cu1qxZw0svvYS/vz9nzpwBYMWKFXTp0oWWLVuyY8cOS93Tp0/Hx8cHX19fPvnkE8BomDp37ky7du2YOnWqJdeNyWTi1VdfpU+fPsybN4/IyEjGjBnDe1NHcvW7FwkK/Pu213fkyBG6dOmCv78/vr6+nDp1ijfeeIN58+ZZjnnttdeYP38+YWFh9O7dG39/f9q1a8eOHTuYMWMGiYmJ+Pv7M2HCBACWLVtmKfOxxx4jPd24rVOpUiVeeeUVOnbsyIABA9i7dy8mk4mmTZty/cRuyjk73dbd88EHH+Dj44Ofnx8zZswAYNGiRXTu3Bk/Pz/GjBljaZgnT57M448/Tq9evWjZsiXr1q0DjG9uvXr1okOHDnTo0IFdu3blWv7kyZNZuXIl8+fP5/Lly/Tt25e+ffvy1Vdf8fzzz1vOXbRoES+88II1f0aaVnoopbLdgNbmnx2y23I6zx5bx44dVUGdO3dOOTs7q4MHDyqllBo7dqxaunSpUkqpfv36qZMnTyqllNqzZ4/q27evUkqp6OholZGRoZRSatGiReqFF15QSin1zTffqKeeekoppdSkSZPUihUrLPX06dPHctxvv/2m+vfvr5RSasGCBWr06NEqNTVVKaXU1atXb/qplFITJ05Ua9assZTzxBNPWJ4bP3682r59u/J/e6N6fMEG1bp169te47Rp09SyZcuUUkolJyerhIQEde7cOdW+fXullFLp6emqadOmKioqSs2ZM0e98847Siml0tLS1PXr15VSSlWsWNFS3tGjR9WwYcNUSkqKUkqpJ554Qi1ZskQppRSg1q9fr5RSauTIkWrgwIEqJSVFBQUFKT8/PzXp679V7w/+tLx/69evV927d1c3bty46XVHRUVZ6nvttdfU/PnzLe/r4MGDVXp6ujp58qSqV6+eSkxMVDdu3FCJiYlKKaVOnjypMv8mcio/679Po0aNVGRkpFJKqfj4eNW0aVPLa+vevbv6559/bntPNa00AAJVNm1qbuP4XwCmAnOz+7ygBKVtaNKkCf7+xv3ojh07EhISQnx8PLt27WLs2LGW45KTkwEIDQ3l3nvvJSwsjJSUFJo0sW5s+ujRo2+qAyAgIIDHH38cFxfjra5e3UiutmXLFj744AMSEhKIjo6mbdu2DB9uZLu+9957LWUGBAQQfOgwp67EE13FnbQb14mLi8PD43+3Xbp37867775LaGgoo0ePpkWLFjRu3JgaNWpw8OBBIiIiaN++PTVq1KBz5848/PDDpKamMnLkSMv7ktXmzZvZv38/nTt3BiAxMZHatY1VtsqVK8eQIUMA8PHxwc3NDVdXV3x8fAgJCWG6tyevrzrM6SvxtPD0ICAggIceeogKFSrc9PoPHz7M66+/TkxMDPHx8QwePNhS/7hx43BycqJFixY0bdqU48eP06RJE6ZNm0ZQUBDOzs6cPHnS8v5kV35OKlasSL9+/Vi3bh3e3t6kpqbi4+OT+z+sppUyOTb8Sqmp5p99HReOfbi5uVl+d3Z2JjExkYyMDKpWrUpQUNBtxz/99NO88MILjBgxgq1bt/LWW2/lqx5nZ2dL/7xS6rbhl0lJSTz55JMEBgbSoEED3nrrLZKSkizPV6xY0fJ7RkYGb3/5K9N/Pc66p++gXb0qt9V7//3307VrV3777TcGDx7Ml19+Sb9+/Xj00UdZvHgx4eHhPPzwwwD07t2b7du389tvv/HAAw/w0ksv8eCDD95UnlKKSZMm8f77799Wl6urq+X1ODk5WV6zk5MTaWlp9PeuzeurYNOxCFp4emT7+sHoilm1ahV+fn4sXryYrVu3Wp679XgR4aOPPsLT05Pg4GAyMjJwd3fP8f3Ny6OPPsp7771H69ateeghPT9RK3usScscLCL/EpFmjgjIUSpXrkyTJk1YsWIFYDQgwcHBAMTGxlKvXj0AlizJ/naGh4cHcXFxedYzaNAgFi5caPkgiI6OtjTyNWvWJD4+npUrV+Z6/pef/xc3Fyda1fHI9oPq7NmzNG3alGeeeYYRI0bwzz//ADBq1Cg2bNjAvn37LFfU58+fp3bt2kyZMoVHHnmEAwcOAEaDnpqaCkD//v1ZuXIlV65cscR8/vz5PF8rgFeV8vjUq8LmY1cs8X/99deWPvzoaGPZyLi4OLy8vEhNTeW77767qYwVK1aQkZHBmTNnOHv2LK1atSI2NhYvLy+cnJxYunSp5Z5DTuVndeu/VdeuXbl48SLLly9n/Pg8l4/WtFLHmlE9IzAmbP0kIvtEZLqINLRzXA7x3Xff8dVXX+Hn50fbtm1ZvXo1AG+99RZjx46lV69e1KxZM9tz77vvPv7v//6P9u3bW27uZufRRx+lYcOG+Pr64ufnx/Lly6latSpTpkzBx8eHkSNHWrpUsjN//nyO/BNE+DdP4+fTjoULF952zI8//ki7du3w9/fn+PHjliv4cuXK0bdvX8aNG4ezszMAW7duxd/fn/bt2/Pzzz/z7LPPAjB16lR8fX2ZMGECbdq04Z133mHQoEH4+voycOBAwsLCrHtTMUb3HLhwjaj4ZIYMGcKIESPo1KkT/v7+zJkzB4BZs2bRtWtXBg4cSOvWrW86v1WrVvTp04ehQ4eycOFC3N3defLJJ1myZAndunXj5MmTlm9FOZWf1dSpUxk6dCh9+/7vy+u4cePo2bMn1aoVz3kRmmZX2XX857RhTOb6FkjPz3mF3Qpzc7ekS01LV61f/129ufpwvs9NT09Xfn5+lhvYjnIoNEY1emWd+nHfhXyfe+tNc3u56667VEBAgN3r0bSiRA43d63Kxy8ijUXkZeAHoDXwsr0+iLSbnY6MJzE1Hf8G+Ru/f/ToUZo3b07//v1p0aKFnaLLXtu6lfGq4k7A0eI3izcmJoaWLVtSvnx5+vfvX9ThaFqRsCYf/9+AK7ACGKuUOmv3qDSLYPOMXd/6t9/UzU2bNm04e7Zo/qlEhAHenqzcH0pSajrurs5Wn7t48WL7BQZUrVrVMiJI08oqa674JymlOiil3teNvuMFh8ZS2d2FxjUq5n1wMTKgjSeJqensOhNV1KFomnaL3NIyT1RKLQPuFJE7b31eKfWhXSPTAOOK37d+VZycil9Gztx0a1qdiuWcCTh2hX6tPYs6HE3Tssjtij/zEtMjh02zs6TUdE6Ex7mmkJAAACAASURBVOHXIH/dPMWBm4szfVrVYvOxCDIy9Eqdmlac5JaW+XPzz7ez2xwXYuGYTCZLv3Fqaiomk4lly5YBkJCQgMlk4scffwSM8fsmk4lffvkFgKioKEwmE2vXrgUgPDwck8nEhg0bALh48SImk8mSBfLs2bOYTCa2bdsGwIkTJzCZTJa8MocPH8ZkMrFv3z4AgoKCMJlMlrH5+/btw2QycfjwYQCWr91E6LJXqJ5qdJds27YNk8lk6bsPCAjAZDJx8eJFADZs2IDJZCI83Fj4fO3atZhMJqKijPN/+eUXTCYTsbGxgDEM1GQyWcbAL1u2DJPJZBnPv3jxYkwmk+W9XLRokSXxGcCCBQsYOnSoZX/evHmMGDHCsn9j368cXvImhy8b9c2ePZv77rvP8vysWbOYOHGiZX/mzJk3Taj617/+xdSpUy3706dP56mnnrLsP/fcczz33HOW/aeeeorp06db9qdOncq//vUvy/5DDz3EzJkzLfsTJ05k1qxZlv377ruP2bNnW/bHjBlz0/DQESNG3JT/aOjQoSxYsMCyP2DAABYtWmTZL8l/e7t27cJkMnHixAmg5P3tzZkzhzFjxlj2Hf23V9zl1tUzP7cTlV6By+7OXIkHoLVX5SKOpGCa1a4EQMDRCHztlFVU07T8E6Wy/xouIpPMv/YE2gA/mvfHAvuVUs9ne6IddOrUSdkz131x9dwPB9lzNpo9r5bcYYfjFu4mLjmN35/tVdShaFqZIyL7lVKdbn08t1w9S8wnTsZYgSvVvL8Q+MNOcWpZ/BMam+9hnA6VmgjntsPpAIg8DgnXQASqNQbPduA9nAHetXjv9xOEXkugfrXit1awppVFeY7j538rcGUmQSlRK3CVVLGJqZyNusGYjvWLOpTbJcfDvkWw6xNIuAquFYyGvkp9yEiFK0fh2FrY+h6Taniz32kwfx5rw4M9rMtyqmmafVnT8GeuwLXFvN8HeMtuEWkAHAo1boLZa8WtAju3HX59HK5fgmb9ofuT0OgOcHW/+bi4CDi+Frc9C/m83Mec3roBWnwLtVoVTdyaplmUiRW4SqLgUGPGrk9x6erJyIDNs2DJCHAtDw9vhAd+geYDbm/0ATw8ofOj8OQe1jV5nerJl1Cf94a/v4Ac7itpmuYYVuXqAZyBSOAa0FJEeud1goh8LSJXRORwlsfeEpFLIhJk3m6bGKYZgi/G0LRmRaqUdy3qUCA9FVY/BTvmQPsJ8Nh2aNjNunOdXfDs8wiDk//DlZpd4feXYM3TRpmaphUJa3L1/Ae4FzgCZJgfVsD2PE5dDHyKkc0zq4+UUrfnztVuEhwaQ49m2aeEdqj0VPhpEpz4Dfq+Br1fMm7g5kOHhtVIq1CL/1R9iw9broft/wexF+G+5VCuZKWi0LTSwJo+/pFAK6VUcn4KVkptF5HGBQmqrAuPTSLienLRj+hRCtY+azT6Qz+Aro8VqBhnJ6Fv69psPnaFtLGv4lKtsXHVv/xeuP9H3fhrmoNZ09VzFiM7p61ME5F/zF1BOa6CISJTRSRQRAIjIyNtWH3xl9m/75fPVMw2t+VdCPoO+swocKOfaaC3J7GJqew/fw3aT4RRX8D5nUbjn5qUdwGaptmMNQ1/AhAkIp+LyPzMrYD1/RdoBvgDYWS/kDsASqkvlFKdlFKdatWqVcDqSqZ/QmNwcRLaFOWM3aNrjC6Z9g+AaUahi+vVshblnJ0IOGbO0e87FkYuhJAdsOpx4+axpmkOYU1XzxrzVmhKKcvKHCKyCFhni3JLm+CLsbT28shXHnubunrGuJlbtwPcNTffffrZqeTmQvdmNdh0NIJX7/Q2Fkj3uxfiI2DTG8YcgEHv2CB4TdPykmfDb8uhmyLipZTKXLx1FHA4t+PLoowMRXBoDMP9imiOXFoKrJgE4gTjloCLm82KHtDGkzdWHeZM5A2am/P40ONpiLlgTAar5W2MGtI0za5yS9L2k1JqnIgcwhjFcxOllG9uBYvI94AJqCkiocCbgElE/M3lhQCF6zguhUKu3iAuKQ3/opq4tf0DCD9kjLip2tCmRfdvXZs3gIBjEf9r+EVgyGyIOgnrngfPNlC3vU3r1TTtZrld8T9r/jmsIAUrpcZn8/BXBSmrLCnSG7uXDsCOD8Hvfmh9l82Lr1u1PG3rVibgaASP92n2vyecXeCeb+CLPvDjA8Y8gQrVbV6/pmmG3PLxh5l/ns9uc1yIZUvwxVgqlHP+3xWxo6SlwKonoZInDHnfbtUM8PbkwIVrXI2/ZXRwxRpw71KICzeGkOrZvZpmN9bO3NUcJDg0hnZ1q+Ds6KUW93wGkcdg+MdQ3n7fNga28SRDwZYT2QzRrdse+r8Bx9bAwWV2i0HTyjrd8BcjKWkZHLl8Hf+GDu7mibkI2z6A1sOg5WC7VtW2bmXqVHYn4GhE9gd0fxoa94LfXzFGF2maZnO64S9GToTHkZKW4fiMnBtmGF0rduziySQiDGhTm+2nIklKTb/9ACcnGLXQ6Pf/ZarO6aNpdpBjwy8ih8wzbLPdHBlkWRF08RqAY6/4z26D4+ug93Sbj+LJyaA2dUhISWfHqajsD6hSH4Z9DJcCjW8imqbZVG6jejJH82SuMLzU/HMCxmxezcaCLsZSs5Ibdatkk+bYHjIy4I/XoEpD6D7NMXUC3ZvVoEp5V34/HMbANp7ZH9RuNJz6w8gI2mIQNOjssPg0rbTLbVRP5uidnkqpl5VSh8zbDMC+HcFlVHBoDP4NqhizWh3hnx+MMfsD3sw+p76duDo7MbCNJ5uORpCSlkuqhqEfgIeXkdAtLcVh8WlaaWdNH39FEbkjc0dEegA6naKNXU9K5UxkvOP691MTjYVV6nWEdmMcU2cWd/rUIS4pjV1ncujuAXCvDHd9aIw22jnPccFpWilnTcP/CPCZiISIyDlgAfCwfcMqew6HxqKUAydu7V0EcZdh4Cyb5OLJr57Na+Lh5sLvh8JzP7DVEGg72phRHHnSMcFpWilnzdKL+5VSfoAv4K+U8ldKHbB/aGXLwYvmGbuOuOJPug5/fWSsmdu4p/3ry4abizP9vWvzx9Fw0tLzyMw59D/Ggu5rn9FZPDXNBvJs+EXEU0S+An5USsWKSBsRecQBsZUplqUWKzhgqcU9CyAxGvq9bv+6cjGknRfXElL5+1x07gdWqg2D34ULu2H/N44JTtNKMWu6ehYDG4HMdJEngefsFVBZFRwa45hunoRo2PUpeA+Heh3sX18uTK1qUaGcM+sPheV9sP8EaNIbNr0J1y/bPzhNK8WsafhrKqV+wrzerlIqDchm5o1WUJlLLfo5YqnFvz+HlDgwvWr/uvLg7upM31a12XgkgvSMPHLziMDweZCRChtfc0yAmlZKWdPw3xCRGphTM4tINyDWrlGVMUEXHZSRM+k6/P1fIzWDZxv71mWloT51iIpPJjAkj+4egOpN4Y7n4cgvcG67/YPTtFLKmob/BYwVuJqJyE7gW+AZu0ZVxgRdjMHVWfC291KLgV9BUiz0etG+9eRD31a1cXNx4vfDeYzuydTzWajaCNa/pNM5aFoBWdPwHwH6AD0wFk5pCxy3Z1BlTfDFGLy9Ktt3qcWUBNj9mTGSp4j79rOq6OZCn5a12HA4nIy8unsAXMsbo3wij8PfC+0foKaVQtY0/LuVUmlKqSNKqcNKqVRgt70DKyvSMxSHLsXib+9ungPfwo1IIydPMXOnjxfh15M4cOGadSe0GgotBsPW2XDdihvDmqbdJLckbXVEpCNQXkTai0gH82YCKjgswlLubGQ88clp9h2/n5YMu+ZDwx7QqIf96imgAW08cXd1Yk1wPkbrDJ1tdPVsesN+gWlaKZXbFf9gYA5QH/gQmGvengeKfkhIKeGQG7vB38P1S8Xyah+gkpsLA7w9+e2fMFLzmsyVqXpTo7//0AoI+cu+AWpaKZNbkrYlSqm+wGSlVN8s291KqV8cGGOpFnQxBg83F5rWtFP6o/Q0Y5Zu3fbQrJ996rCBEX51uXojhZ2nc8ndc6s7njcyi+obvZqWL9b08XcUEcvlqIhUE5F37BhTmXLwQgy+DargZK+lFo+ugmsh0Gt6keTksVafVrWo7O7CmqB8dPeUq2B0+Vw5auQe0jTNKtY0/EOVUjGZO0qpa8Cd9gup7IhPTuN4+HU6NqpunwqUgt2fQo3m0Kp4/5O5uThzp48XG4+Ek5iSj/mBre6E5gNh6/sQl8Nyjpqm3cSaht9ZRNwyd0SkPOCWy/GalYIuxJChoGOjavap4PwuuHwQuj1pLGlYzI3wq8uNlHT+PH7F+pNEjOGdaUkQ8Kb9gtO0UsSa1mAZsFlEHhGRh4FNwBL7hlU27D9/DRFob6+lFnd/CuWrg994+5RvY12b1qC2hxurgy7l78QazYwVxIK/h4t77ROcppUi1qRl/gB4B/DGmLw1y/yYVkiB56Np5elBZXc7ZOSMOg0nfofOjxp94SWAs5Mw3K8uW09EEpuQz5u1vV4Ej7rGjd4MnUpK03Jj7ff/Y8AGpdSLwA4R8bBjTGVCeoYi6EKM/bp59nwGzuWgyxT7lG8nd/vXJSU9gw1H8jkxy60SDJoFYUFwcJl9gtO0UsKafPxTgJXA5+aH6gGr7BlUWXAyIo645DT7NPw3rkLQcvAdZ+SyL0F86lWhcY0KrM7P6J5M7cYYk9Q2vw2JMXkfr2lllDVX/E8BPYHrAEqpU0DJak2Kof3njfQEnewxoifwK+NmZ/dpti/bzkSEEf712H32KhHXk/J7snGjN/GaMcpH07RsWdPwJyulUjJ3RMQFc4pmreD2n79GzUpuNKhe3rYFpybB3i+MIY61W9u2bAcZ6V8XpeDXg/m8yQvg5QsdHzLG9UcctX1wmlYKWNPwbxORVzFy9gwEVgBr7RtW6bf//DU6NqqK2HpS1aGfjGRsPUre1X6mprUq0alRNX4KvIhSBbjG6Pc6uFeG31825jJomnYTaxr+GUAkcAgjLfN6oGgXay3hrsQlcSE6wfbdPEoZqZc9faBJH9uW7WDjOjXgbOQN6zN2ZlWhutH4h+yAo6ttH5ymlXC5ZefcbP71faXUIqXUWKXUPebf9WVUIRww9+93sPWN3dMBRp76HtOKdXoGa9zp60WFcs6sCAwtWAEdHzI+AP943ViLQNM0i9yu+L1EpA8w4pa0zB1EpPis5FEC7T9/jXIuTrSrZ+MVt3Z9Ah5e0Ha0bcstApXcXLjTx4u1wZdJSEnLfwFOzsaN3tiLsPNj2weoaSVYbg3/TIxunlvTMs/FSNecKxH5WkSuiMjhLI9VF5FNInLK/NNOg9iLt8Dz1/CtVwU3FxuuuBX2D5zbBl0fA5dytiu3CI3r1IAbKemsP2Tlsoy3atwT2t0Df31sJKrTNA3IPS3zSqXUUOCDW9Iy91VKWZPfdzEw5JbHZgCblVItgM3m/TIlKTWdw5di6djYxp95uz8D14rQcbJtyy1CnRtXo3GNCqwIvFjwQgb+27j6/0PfltK0TNakbJglIvVEpIeI9M7crDhvOxB9y8N38788P0uAkfmOuIQ7dCmW1HRFx4Y2bPivX4bDK6HDA1C+9HyJEhHGdmrA3+eiCYm6UbBCqtQz0jkcWwtnttg2QE0roayZuTsb2Ikxkucl81bQpZw8lVJhAOafOU4EE5GpIhIoIoGRkZEFrK74yZy4ZdMZu39/DioDuj1huzKLiTEd6uMksHJ/AW/ygjGRrVoT+P0VvWCLpmHdcM5RQCul1J1KqeHmbYS9A1NKfaGU6qSU6lSrVi17V+cwgSHXaFKzIjUq2SizdXI87P8GvIdDtca2KbMYqVPFnd4ta7FyfyjpGQUcTObqDkPeh6gTxuQ2TSvjrGn4zwK2Sh8ZISJeAOaf+Ui8XvJlZCgCz0fb9mr/4DJIioXuT9uuzGJmXKcGhF9PYvvJQnzzazkEmg+ArbMhvkz92Wnabaxp+BOAIBH5XETmZ24FrG8NMMn8+ySgTM2uORERR0xCKt2a1rBNgelpsGcBNOgGDTrbpsxiaIC3JzUrubF0z/mCFyICQ2ZDaiIEvGWz2DStJLKm4V8DzAJ2AfuzbLkSke+B3UArEQkVkUeA2cBAETkFDDTvlxl7zl4FoGsTG83YPb4WYs5Dj9J7tQ9QzsWJ+7s0YMuJK1y4WojJWDVbQPenIOg7CNlpuwA1rYRxyesApVSBVttSSuW07FP/gpRXGvx9Npr61crToLoNFkZRCnbOh+rNoNXQwpdXzN3ftRGfbT3Dsr/P8+qd3gUvqM/LcPgX+O0FeGxHqZnzoGn5kVvKhp/MPw+JyD+3bo4LsXTIyFD8fe6q7bp5LuyGyweMK1gnG04EK6bqVHFncFtPftx3MX+Lsd+qXEW4a46R2mL3p7YLUNNKkNy6ep41/xwGDM9m0/Lh5JU4rtmyf3/XJ1ChRolZT9cWHuzemNjEVNYGF2CRlqxaDjZGQW37QM/o1cqk3GbuZo63P5/d5rgQS4c9Z2zYvx91Ck6sh85TSsx6urbQtUl1Wnl6sGR3SMHSNWc15D/GN6X1L+nUzVqZY+2au1oh/X0umnpVbdS/v/tTcHE3FlIvQ0SEB7o34sjl6xy4UMilFavUg76vwqk/4Nga2wSoaSWEbvgdICNDseesjfr34yMh6Huji6dS6ZnYZq1R7evh4ebCt7tDCl9Yl8egjo8xozc5rvDlaVoJkWc+fhH5j+PCKZ2Ohl3nWkIqd7SwQcO/bxGkpxg3dcugim4u3NOpPusPhREZl1y4wpxdYNjHEBcOW96zTYCaVgLofPwOsONUFAA9m9csXEEpCcZasq3uNMakl1EPdGtEaroq3ISuTPU7QaeH4e+FEBZc+PI0rQSwWz5+7X92no6ilacHtT3cC1dQ8HJIjC71E7by0rRWJQa28eTb3SEFW6TlVv1nGiOk1j5rzIbWtFLOnvn4NYz8+3tDogt/tZ+RbuTcr9cJGnazTXAl2ON9mhGTkMoPewuRqz9T+arGal2XD8KezwpfnqYVc9bm4x8hInPM2zBHBFZaBIZcIyUtg14tCtnwH/8Nos8aV/slfD1dW+jYqBpdmlTnyx1nSU3PKHyBbUdD62FGX3/U6cKXp2nFmDX5+N/HmMx11Lw9a35Ms8Jfp6NwdRa6FGb8vlLw14dG2mVvPXcu0xN9mnE5Nok1QYWc0AXGh+ldc8HFDdZMgwwbfJhoWjFlzXDOu4CBSqmvlVJfYyyneJd9wyo9/jodSfuG1ajolmdapJyd+dPohrjj+TKRnsFapla1aF3Hg8+3nyGjoLn6s/KoA4PfN9Jh7Puy8OVpWjFl7Tj+qll+r2KPQEqjyLhkDl+6Tq/C9u/vmAsedctUegZriAiP92nGyYh4/jxuoxz7/vdDs/5G6madzkErpaxp+N8HDorIYhFZgpGSWQ96tsI288IhfVvnuMJk3s7vhvM7oeczRjeEdpNhvl7Ur1aehdvO2KZAERg+z/i55hnd5aOVStbc3P0e6Ab8Yt66K6V+sHdgpcGWE1eo7eFG27qVC17IjjlQoSZ0mJT3sWWQi7MTU3o1JfD8Nfaei7ZNoVUbwKBZcG6bMWFO00oZq7p6lFJhSqk1SqnVSqlwewdVGqSlZ7D9ZCSmVrWQgo7CuXwQTgdA9yfLVDK2/BrXqQE1K5Vj3uaTtiu040PQYhBsmgmRJ2xXrqYVAzpXj50cuBBDXFIafVsVoptnx1xwr2Jk4dRyVL6cM0+amrPz9FV2no6yTaEiMOJTcK0Av0yBtBTblKtpxYBu+O1ky4kruDgJPQs6fv/KcTi21kgk5l6IrqIy4v6uDalbxZ3/23ii8CmbM3l4Gv39YcGw/QPblKlpxUCuDb+IOInIYUcFU5psOX6FTo2rUdndtWAF7JgDrhWh2xO2DayUcnd15pn+LQi6GMPmYzYa4QPQZgT4TzC+fV3ca7tyNa0I5drwK6UygGARaeigeEqFi9EJHA+Po19BR/NEHIVDK6HrVKhgo4XZy4AxHevTpGZF5vxxwjbj+jMNmQ1V6htdPknXbVeuphURa7p6vIAjIrJZRNZkbvYOrCT742gEAIPb1ilYAVveBTcP6PGMDaMq/VydnXh+YEuOh8ex9h8bzObN5F4ZRi+CmItGIje9YpdWwlkznfRtu0dRymw8Ek7rOh40qlEx/ydfOgDH14HpVX21XwDDfLxYsOU0H206yZ0+Xrg62+g2VsNu0O812PxvaNIbOj1km3I1rQhYM45/GxACuJp/3wccsHNcJdbV+GQCQ6IZ1MazYAX8+Q6Ur6779gvIyUmYPqgVIVcTWBEYatvCez4PzfrBhhkQrm99aSWXNUnapgArgc/ND9UDVtkzqJJs87ErZCgYVJBunvO74MxmuOM5PZKnEPp716ZTo2rM/eMEsYmptivYyQlGfWEMsV0xGZLjbVe2pjmQNd+DnwJ6AtcBlFKngEIMTi/dNh4Jp17V8vmfrauUcbVfyVOP2y8kEeGtEW2JTkhhXsAp2xZeqRaM+RKiz8BvL+r+fq1EsqbhT1ZKWWaviIgLoP/asxGXlMqO01EMblsn/7N1z24xcvL0fknP0rWBdvWqcF/nhizZHcKpCBsvpN6kN/R5Bf75AQK/sm3ZmuYA1jT820TkVaC8iAwEVgBr7RtWybTpaAQpaRnc5euVvxMzMoxskFUaQIcH7RJbWfTS4FZULOfMW2uP2G5SV6beL0OLwfD7K0YXnaaVINY0/DOASOAQ8BiwHnjdnkGVVGuCL1Ovank6NKya98FZBX9vzA7t/6bOwGlD1SuW48VBrdh5+iobj9g4xZSTE4z+Aqo2gp8ehNhLti1f0+zImlE9GcASYBbG0M4lyuaXTyVf9I0U/joVxXC/uvnr5kmOh81vQ/3O4HOP/QIsoyZ0bUjrOh7MWneMxJR02xZeviqM/x5Sk+DHicZPTSsBrBnVcxdwBpgPfAqcFpGh9g6spPn9cBhpGYrhfvns5tn5McRHGCs/6bV0bc7F2Yk3h7flUkwi/7VVzv6sarWC0Z/D5QPw2wv6Zq9WIljT1TMX6KuUMiml+gB9gY/sG1bJszb4Ms1qVaSNVz5G88RchF2fQLt7oEFn+wVXxnVvVoO7/evy362nORZmh5QLre+CPjMg6DvYOc/25WuajVnT8F9RSp3Osn8WsGEWrJLvckwif5+Lzn83T8Bbxs8Bb9khKi2rN4e3pUp5V15aGUxquh1W1erzCrQbAwFvwuFfbF++ptlQjg2/iIwWkdEYeXrWi8hkEZmEMaJnn8MiLAF+3h+KUjCmQ33rT7q4Fw6vhB5PGys+aXZVvWI53hnZjsOXrvO5Pbp8nJzg7gXQsDv8+jhc2GP7OjTNRnK74h9u3tyBCKAPYMIY4VOtMJWKSIiIHBKRIBEJLExZRS0jQ7Fifyg9mtWgQXUrx99npBvDACvVgZ7P2TdAzWJIOy+G+Xoxb/MpToTbeGw/gKs73Lfc+CD//j6IOp33OZpWBHJM0qaUsncWqr5KKRstl1R0/j4XzYXoBF4Y2NL6k/YuMm4GjvkK3CrZLzjtNm+PaMvuM1eZviKYX5/sgYutkrhlqlAdJqyALwfCd2PgkQBjtq+mFSPWjOppIiIfisgvOi3z7VYEXsTD3YUh7azMzRNz0cjw2Hyg0SesOVSNSm7MGtmOQ5diWWiPLh+A6k1h/A8QfwWWjoLEa/apR9MKyJrLnVUY2Tk/wRjhk7kVhgL+EJH9IjI1uwNEZKqIBIpIYGRkZCGrs4/YxFTWHw5jhF9d3F2d8z5BKVg/HVBw11w9fLOI3OljdPl8FHCKwJBo+1TSoDPc9x1EnYDvxkKyHbqWNK2ArGn4k5RS85VSW5RS2zK3QtbbUynVARgKPCUivW89QCn1hVKqk1KqU61axfOr8sr9oSSlZjC+i5ULlB1dBSc3QL/XoVoj+wan5eq90T7Ur1aeacsPEn3DTgupN+sH93xtrLHw/XhITbRPPZqWT9Y0/PNE5E0R6S4iHTK3wlSqlLps/nkF+BXoUpjyikJGhmLp7hA6NqpGu3pV8j4h8Rqsfxm8/I0F1LUiVdndlc/u70D0jRSe/zHItks1ZuU9HEb+F0L+gp8mQZqdPmQ0LR+safh9gCnAbP7XzTOnoBWKSEUR8cj8HRgElLhVLbadiiTkagIPdrfyyn3Tm5BwFUbMB2drFj7T7K1dvSq8MbwN205G2mdWbya/e42uvVMbjbw+OrWDVsSsaYFGAU2zpmYuJE/gV/NEJxdguVJqg43Kdphvd4VQy8ONoe2sSNFw4nc4sAR6PgtefvYPTrPaxK4N+fvsVeb+cYJOjarRtWkN+1TU+RFQGcY9nh/Gw73f6fTbWpGx5oo/GMhnusmcKaXOKqX8zFtbpdS7tirbUc5ExrP1ZCTjuzSknEseb2FcOKx+Cur4Qt/XHBOgZjUR4f3RPjSqUZGnlh8k9FqC/SrrMgVGfApntugbvlqRsqbh9wSOi8hGPZzT8Pm2M5Rzdsq7mycjw5jFmZJgjNnXKZeLJQ93VxY92JHktHQeXryP60k2XK7xVh0eMFbwurDbPNQzxn51aVoOrOnqedPuUZQgYbGJ/HrwEuO7NKRmpTwa8j2fGStrDZ8HtfIxwUtzuOa1Pfh8Ykce/HovT313gK8nd8bV1pO7MvncY1wErHgIvh4CE1dClXyk+9C0QrImH/+27DZHBFccLdp+jgwFU3o1zf3AsGAIeBtaD4MOkxwTnFYoPZrX5L1RPuw4FcXM1Ydtv2pXVt7DYeLPcP0SfDkAwv6xX12adgtrZu7Gich185YkIukiYofctsVf9I0Uvt97gbv96uaelyfpOqx8BCrWghGf6IlaJci4zg14qm8zvt97kc+3n7VvZU37wMMbQZzgm6FwerN969M0M2uu+D2UUpXNmzswBmNBljJnwZbTJKel82Tf5jkflJEBq56A6LMwZpGRrZd1uQAAE/lJREFUu0UrUV4c2IrhfnWZ/ftxlu05b9/KPNvAowFQrTEsHweBX9u3Pk3Dupu7N1FKrQL62SGWYu1yTCLf7jnPmA71aV47l8Rqf30Ix9fBoHeg8R2OC1CzGScnYe5YPwZ41+b1VYf5Ye8F+1ZYuS489Ds07Qvrnoc1z0Basn3r1Mq0PG/umnPyZ3ICOmHk2ilT5gWcAgXP5ZaF89g6+PMd8BkL3Z5wXHCazZVzceKzCR14bOl+/vXrIVycnbinox1vwLpXhvt/hC3vwo65EHEE7l1qfChomo1Zc8U/PMs2GIgD7rZnUMXNqYg4Vuy/yMRujahXtXz2B10+CL9MgXoddb9+KeHm4szCiR25o3lNXloZzK8HQ+1boZMz9J8J45ZC5HH4vI+R6kHTbCzPK34H5OUv1pRSvLnmCB7urkzrl0Pf/rXzsPw+qFATxn8Prjl8OGgljrurM1880ImHF+/jhZ+CiUtK48Huje1baZsRULMl/DgBlgyHXi8aSzs6u9q3Xq3MyLHhF5GZuZynlFKz7BBPsbP+UDi7zlxl1t1tqV6x3O0HxEcaE3HSEuHBVVCptuOD1OyqfDlnvp7cmae/P8jM1UcIj03ipcGt8re+cn7Vbg1TtxkrtW3/Pzi7FUYvgupN7FenVmbk1tVzI5sN4BHgFTvHVSzcSE7jnd+O0sarMvd3zWaWblIsfHePMRb7/p+gtrfjg9Qconw5ZxZO7MD4Lg1ZsPUML66w06LtWbn9f3v3Hh5VeSdw/PvLTCbJkAuEBAmXEEREkSIiyE1qV+uKtPXSYqt1tXTZ2m1ra31qrdS9dNvtPrVWW7deutRatWu1FaGl1mp5KK5SEYmWqxAlgNwCCSQkIbeZyfz2j/cAQ0gQJJkzyfw+z3Oec03mx0vO7z3nPWfeNxeuech17VzzDvxsJqz5tRvXwZjTcKKhF48MtuL1pnkb8HngGU5/IJZe4YcvbqaqvpWf3nABgYwOV3etDfC/n4J9G1yHW6VT/QnSJE0wkMF/XTuOkoJs7l/6DjWNbTx4w0QKwj3cBDPuUzBsMiy6xb0qvOE5+PhP3Ni+xnwAJ3y4KyKFIvKfwDpcJTFRVb/l9aPfp62sPMATK99j7vQyJpV1eBe/tcF1srX7LbjucRgzy5cYTfKJCF+7bDT3fOpDrKw8wCceXMHGPfU9/8H9S2HuC3DlD+G9lfDwVDd2c7yH7zpMn9Rl4heRe4HVuLd4PqSq31HVtBg8tKktxp3PraVsYJg7Z43psPMAPHkV7C6HOb9wX703aeczk0v5zRen0RZr55MPv8az5Tt7/kMzMmDKF+HLK90dwAt3wGNXuDfKjDkFJ7ri/wYwBPgXYE9Ctw2NfbnLBlXl7sXr2VXXwr3XnU84lNAadnAHPD4bqje55p3zrvUvUOO7C0cM4I9fm8nE0gF8c+E65i9aT2u0vec/eMAIuGmxG9mrbhss+DtY8lX3ooExJ6HLxK+qGaqa06HLhvzD68kMMpmeWb2T363Zw+0fPZvJiU08O9+An18KDVVw40Jr3jEAFOVm8at5F/HPl4zi6Td2MPuBV3lrRxJujEVgwmfhq2/CtK+4h74/vRBee9BG+DLvq4f6ne2d1u+q59+XbGTm6CJuTeyPZ+1v4PGPQSjX9asycqZ/QZqUEwxkcNeV5/DUP02hNdrOnEde454XN9MWS8LVf3YBXPF9+NJKGDYJ/ny3qwDeehLaYz3/+aZXssTv2XOwhXlPrKY4N4sff2YCGRkC8XZY9l1YfAsMnwJf+Iv1q2+6NOOsIl68/cPMuXAYj7xcyVU//Svl22uT8+HFZ8NNi+DmJZB3hmv6eXiKewMonoQKyPQq0qN9jneTSZMmaXl5eY/9/kNtMeY88hq761pY+KXpjBmcB/W73OhZ2191/enP/hEEO/kClzGd+MvmfXx70Qb2NrRy7QVDmX/lOQzKz07Oh6tCxQuw7HtQswkKR7nxns+/3kaBSzMi8qaqTjpue7on/uZIjLmPrebNHXU8Nncyl5xdDBsXwx9uc7fKs38IE260vnfMKWtqi/Hwy1v4+SvbyAy410DnzigjKxhITgDxdtdT7Kv3Q9UayB0M074ME2+GnAHJicH4yhJ/J1qj7cx7YjUrKw/wwPUX8Imzw/DSt2HNUzB0EnxyAQwc1e2fa9LL9v1NfO/5t1m2uZohBdnceulorps0rOeGduxI1Q0B+ur97g42mAPjP+0Gfx/8oeTEYHxhib+D+pYoX3iynNXba7lvzng+GVgBS/8Vmg/AzDvgkjutUyzTrVa8u5/7llbwtx0HGV6Yw22Xnc01E4YQTFYFAG6Ix9U/h3XPuv6lhk91A8CPvcZ1EWH6FEv8CarqW/j8L1dTWXOIR6/I5pJ374Gdr7svxcz+EQyZ0G2fZUwiVeXlihruW1rBht0NDO2fw9zpZXzmouHkZyfxQqO51t3Zlv8Saishsx+cd417RbR0uvuymOn1LPF7Xt96gFt//RZF0SqeGPV/nLFtkWvvvPy7cP5n7Q/eJIWqsmxTNY+u2MrrW2vpFwrw6cnDuWnqCM4sTuKVtyrsXOUqgQ2LIdLongWc+wkYezWMmO7GCTC9Uton/va4suCVrTz95xV8K7yE2fGXEQnA5HmuWccedhmfbNhdzy9WbOMPa/cQiyuTywbw6UnD+dj4kmO/Od7TIk1Q8Sd4+/fw7lLXFNSv2FUC514FI2bYm229TFon/h0Hmnnkqd8yqeZZrgmsJCMQQC78PFx8O+SXdGOkxnxw1Q2tPPfWbp4t38nW/U30CwWYNa6Ej40fzMVnFRMKJvFuNNLkkv/bv4d3XoJok/sCY9lMGHUpnHUZFJ5pb7uluLRM/K0tzbyyeAGDK55kvFQSDYQJXvgPyIyvQ8HQHojUmNOnqpS/V8dvV+/kxY17aWyNkZ8d5PKxg7nivDOYcVYR/bKSeCcQbYHK5VC5DLYsc/0DAfQf4SqBsouhdJqdUykoLRN/+f1zmNSwlL2Zwwlf/CXyp9zkBrU2ppeIxOKs2FLDH9ft5c9vu0ogFMhgypmFfGTMID48uoizBuX27GhgHdVudRVA5XLY9op7LgBQUAqlU9zYFKXToPgcez7gs7RM/NvXv0ZDXTXjZ15tt6Sm14vE4pRvr2V5RTXLK2rYUn0IgMJ+ISaXDeCikQOZMrKQc0vyjx84qKe0x2DfetixCnashB2vw6G9bl9mP/c9gZLzoWS8mxefY69JJ1FaJn5j+rKdtc2srDzAqm21rN5ey47aZgDysoJMKO3PuKEFjBtSwHlD8iktDLv+p3qaKhx8z1UAe9ZA1VrYuw4irpIikAXFY1wFcGR+Dgwog0ASm6/ShCV+Y/q4qvoW3thWy6pttazdeZB39jUSbXfnd15WkLFD8hk7JJ+zBuUyqthNRbmhnm8misdd81DVGjdVb4KaCqhPGLwmEHJ9Cg0oc1PhyKPL/UshM6dnY+yjLPEbk2baYu28u+8QG3bXs3FPAxv21LO5qpGWhMFi8rODjBqUy8iifgwbEGZY/xyGDshhaP8cSvpn92y/Qm2NbhD5ms2uM7kDW6Fuu5uiTccemzfEVQD5JW45vwTySiB/iJvnlUBmkjrB60Us8RtjiMeVvQ2tVNYcorL6EFtqDlFZ3cS2/U3sa2wlMR2IQHFu1pGKYFBeNkV5IYpysyjOy6I4N4ui3CwG5oa6t98hVWja794eOlwR1G13I+A17IHGKog2H/9zOYXQrwjCRRBOXB7oLSfMswvc66l9/NlfSiV+EZkFPAAEgEdV9QcnOt4SvzE9LxKLs7e+lV0Hm9ld18Lugy1H5nsOtlDT2EZTpPO+/fuHMxkQDpGfk0l+dpCCnEwKcjLJ9+YFOZnkZ7t5bnaQcChATmaAcChAOBQkOzPj5JucVKG13lUAhyuChio3b97vuqNo2u8tHwDtYkB6yYCsPMgqcG/7ZeV3PQ/luuamUBgyw2654zyYk3Lf/O8q8Sf9aYqIBICHgMuBXcBqEVmiqm8nOxZjzFGhYAalA8OUDgx3eUxzJMb+xgg1h9rY7001jW5e3xKjviVKfUuU3XUtR5Zj8fe/uBThmIogHAqQEwocqSBCwQxCgQw3D2YQCnjbgsPJCo5w+wozCA069rigKFmxRrKjdWRF6siKHCSrrZZgtIHM6CGCsUMEIo0Eo40EIo1k1O0iI9JIRlsjEmlA4qc4ilkwp0OlkO0eaAdC7m2m4OFlbwqGjl0/ZluW+5nRfw/9h59aHO8XZrf+tpNzEbBFVbcCiMgzwNWAJX5jUlw4FKR0YPCElUMiVaUl2k59S5QGr2I41BalOdJOc6SdliPzGM2RdpoSlluibl9dU5RIe5xIzJs6LJ+aPG86qejJJkIezYSljRwi5NBGtkQI00qORAhLhDARcsTty4m3EY56y0TIlgiZxAjRQibthIh669Ej60FvPUiMLKLHRVFx+eOMmdH7E/9QIOFxPruAKR0PEpFbgFsASktLkxOZMaZbiYh3BR+kpKD7f7+qdlkptMXixFWJxZX2uBJr14T1OO1xaI/Hj+xvj+sxy0fX3bGx9jhxBUWJK6BuHldFgSZVmrz1uLoWqbjXlO62qbcNQInHj/5sXBUU4vE4aDsBjZARjxKMR7lp+PhuLzc/En9nDXnH3Quq6gJgAbg2/p4OyhjT+4gIWcFA8kY16yP8eBKxC0i8bxkG7PEhDmOMSUt+JP7VwGgRGSkiIeB6YIkPcRhjTFpKelOPqsZE5FbgJdzrnI+p6sZkx2GMMenKl84xVPUF4AU/PtsYY9Jdan3bwBhjTI+zxG+MMWnGEr8xxqQZS/zGGJNmekXvnCJSA7z3AX+8CNjfjeH0BIuxe1iMpy/V4wOL8VSMUNXijht7ReI/HSJS3lnvdKnEYuweFuPpS/X4wGLsDtbUY4wxacYSvzHGpJl0SPwL/A7gJFiM3cNiPH2pHh9YjKetz7fxG2OMOVY6XPEbY4xJYInfGGPSTJ9O/CIyS0QqRGSLiNyVAvEMF5HlIrJJRDaKyG3e9kIRWSoi73rzASkQa0BE/iYiz3vrI0VklRfjb7wutf2Mr7+ILBSRzV55Tku1chSR273/5w0i8rSIZPtdjiLymIhUi8iGhG2dlps4/+2dP+tEZKKPMd7r/V+vE5HFItI/Yd98L8YKEbnCrxgT9t0hIioiRd66L+V4In028ScM6n4lMBa4QUTG+hsVMeAbqnouMBX4ihfTXcAyVR0NLPPW/XYbsClh/R7gx16MdcA8X6I66gHgRVU9BzgfF2vKlKOIDAW+BkxS1XG4Lsivx/9yfByY1WFbV+V2JTDam24BHvExxqXAOFUdD7wDzAfwzp/rgfO8n3nYO/f9iBERGQ5cDuxI2OxXOXZNVfvkBEwDXkpYnw/M9zuuDjH+HvdHUgGUeNtKgAqf4xqGSwCXAs/jhsvcDwQ7K1sf4ssHtuG9nJCwPWXKkaNjSxfiuj9/HrgiFcoRKAM2vF+5Af8D3NDZccmOscO+a4GnvOVjzmvcOB/T/IoRWIi7ENkOFPldjl1NffaKn84HdR/qUyzHEZEy4AJgFXCGqlYBePNB/kUGwE+AO4G4tz4QOKiqMW/d77I8E6gBfuk1Rz0qIv1IoXJU1d3Aj3BXflVAPfAmqVWOh3VVbql6Dv0j8CdvOWViFJGrgN2qurbDrpSJ8bC+nPhPalB3P4hILvAc8HVVbfA7nkQi8nGgWlXfTNzcyaF+lmUQmAg8oqoXAE2kRvPYEV47+dXASGAI0A93y99RSvxNdiHV/t8RkbtxTaZPHd7UyWFJj1FEwsDdwL91truTbb6WY19O/Ck5qLuIZOKS/lOqusjbvE9ESrz9JUC1X/EBM4CrRGQ78AyuuecnQH8ROTxim99luQvYpaqrvPWFuIoglcrxo8A2Va1R1SiwCJhOapXjYV2VW0qdQyLyOeDjwI3qtZmQOjGOwlXya71zZxjwlogMJnViPKIvJ/6UG9RdRAT4BbBJVe9P2LUE+Jy3/Dlc278vVHW+qg5T1TJcmf1FVW8ElgNzvMP8jnEvsFNExnibLgPeJoXKEdfEM1VEwt7/++EYU6YcE3RVbkuAm723UqYC9YebhJJNRGYB3wKuUtXmhF1LgOtFJEtERuIeoL6R7PhUdb2qDlLVMu/c2QVM9P5WU6Ycj/DzAUMSHr7Mxr0BUAncnQLxXIy7xVsHrPGm2bg29GXAu9680O9YvXg/AjzvLZ+JO6G2AM8CWT7HNgEo98ryd8CAVCtH4D+AzcAG4FdAlt/lCDyNe+YQxSWneV2VG66J4iHv/FmPe0PJrxi34NrJD583P0s4/m4vxgrgSr9i7LB/O0cf7vpSjiearMsGY4xJM325qccYY0wnLPEbY0yascRvjDFpxhK/McakGUv8xhiTZizxG9OB1/Pnl73lISKy0O+YjOlO9jqnMR14/Sg9r65XTWP6nOD7H2JM2vkBMEpE1uC+1HSuqo4TkbnANbgulscB9wEh4CagDZitqrUiMgr3hZ1ioBn4gqpuTv4/w5jOWVOPMce7C6hU1QnANzvsGwd8FrgI+D7QrK6juJXAzd4xC4CvquqFwB3Aw0mJ2piTZFf8xpya5araCDSKSD3wB2/7emC81/PqdOBZ10UP4LpqMCZlWOI35tS0JSzHE9bjuPMpA9fn/oRkB2bMybKmHmOO1wjkfZAfVDe+wjYRuQ6OjLd6fncGZ8zpssRvTAeqegD4qzeQ9r0f4FfcCMwTkbXARtyALMakDHud0xhj0oxd8RtjTJqxxG+MMWnGEr8xxqQZS/zGGJNmLPEbY0yascRvjDFpxhK/Mcakmf8HwBCeRgdHZG8AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "Wtt_sir,Penc1 = sm.symbols('Wtt_sir,Penc1',cls=sm.Function)\n",
    "Penc1= (Wt/N)*((N - Rt - Wt)/N)\n",
    "Wtt_sir=N*c*Penc1*tau + Wt - a*Wt\n",
    "Wupdate_sir=Wtt_sir.subs([(N,100),(tau,0.5),(c,0.5),(a,0.1)])\n",
    "Rupdate= Rtt.subs([(a,0.1)])\n",
    "ws_sir=[1]\n",
    "Rt_v = 0\n",
    "for t in range(0,149):\n",
    "    ws_sir.append(Wupdate_sir.subs([(Wt,ws_sir[-1:][0]),(Rt,Rt_v)]))\n",
    "    Rt_v = Rupdate.subs([(Wt,ws_sir[-1:][0]),(Rt,Rt_v)])\n",
    "#plt.plot(range(0,80),ws_sir)\n",
    "#plt.ylabel('Number of infected individuals')\n",
    "#plt.xlabel('time')\n",
    "#plt.title('Diffusion of infection')\n",
    "\n",
    "Wtt_sir2,Penc2 = sm.symbols('Wtt_sir2,Penc2',cls=sm.Function)\n",
    "Penc2= (Wt/N)*((N - Rt - Wt)/N) * 0.80\n",
    "Wtt_sir2=N*c*Penc2*tau + Wt - a*Wt\n",
    "Wupdate_sir2=Wtt_sir2.subs([(N,100),(tau,0.5),(c,0.5),(a,0.1)])\n",
    "\n",
    "ws_sir2=[1]\n",
    "Rt_v = 0\n",
    "for t in range(0,149):\n",
    "    ws_sir2.append(Wupdate_sir2.subs([(Wt,ws_sir2[-1:][0]),(Rt,Rt_v)]))\n",
    "    Rt_v = Rupdate.subs([(Wt,ws_sir2[-1:][0]),(Rt,Rt_v)])\n",
    "plt.plot(range(0,150),ws_sir)\n",
    "plt.plot(range(0,150),ws_sir2)\n",
    "plt.gca().get_lines()[1].set_label('20% social distancing')\n",
    "plt.plot((0,140),(16,16),color='k',linestyle=':')\n",
    "plt.text(0,17,\"healthcare system capacity\")\n",
    "plt.ylabel('Number of infected individuals')\n",
    "plt.xlabel('time')\n",
    "plt.title('Diffusion of infection')\n",
    "plt.legend();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Our next two steps are for better communicating the results of the model. We first show how to make the graph interactive, so the user can change the degree of social isolation and see the effects. Next we create a gif animation showing how the spread of the epidemic changes as social isolation measures are increased. This gif can then be posted to social media, for example."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "423633d66eb6434fb47023d5be891126",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "interactive(children=(FloatSlider(value=0.8, description='percent', max=1.0), Output()), _dom_classes=('widget…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import ipywidgets as widgets\n",
    "from ipywidgets import interact, interact_manual\n",
    "\n",
    "@interact(percent=(0.,1.))\n",
    "def compare(percent=0.8):\n",
    "    Wtt_sir,Penc1 = sm.symbols('Wtt_sir,Penc1',cls=sm.Function)\n",
    "    Penc1= (Wt/N)*((N - Rt - Wt)/N)\n",
    "    Wtt_sir=N*c*Penc1*tau + Wt - a*Wt\n",
    "    Wupdate_sir=Wtt_sir.subs([(N,100),(tau,0.5),(c,0.5),(a,0.1)])\n",
    "    Rupdate= Rtt.subs([(a,0.1)])\n",
    "    ws_sir=[1]\n",
    "    Rt_v = 0\n",
    "    for t in range(0,149):\n",
    "        ws_sir.append(Wupdate_sir.subs([(Wt,ws_sir[-1:][0]),(Rt,Rt_v)]))\n",
    "        Rt_v = Rupdate.subs([(Wt,ws_sir[-1:][0]),(Rt,Rt_v)])\n",
    "    #plt.plot(range(0,80),ws_sir)\n",
    "    #plt.ylabel('Number of infected individuals')\n",
    "    #plt.xlabel('time')\n",
    "    #plt.title('Diffusion of infection')\n",
    "\n",
    "    Wtt_sir2,Penc2 = sm.symbols('Wtt_sir2,Penc2',cls=sm.Function)\n",
    "    Penc2= (Wt/N)*((N - Rt - Wt)/N) * percent\n",
    "    Wtt_sir2=N*c*Penc2*tau + Wt - a*Wt\n",
    "    Wupdate_sir2=Wtt_sir2.subs([(N,100),(tau,0.5),(c,0.5),(a,0.1)])\n",
    "\n",
    "    ws_sir2=[1]\n",
    "    Rt_v = 0\n",
    "    for t in range(0,149):\n",
    "        ws_sir2.append(Wupdate_sir2.subs([(Wt,ws_sir2[-1:][0]),(Rt,Rt_v)]))\n",
    "        Rt_v = Rupdate.subs([(Wt,ws_sir2[-1:][0]),(Rt,Rt_v)])\n",
    "    plt.plot(range(0,150),ws_sir)\n",
    "    plt.plot(range(0,150),ws_sir2)\n",
    "    plt.gca().get_lines()[1].set_label('{:.2f} social distancing'.format(1-percent))\n",
    "    plt.plot((0,140),(16,16),color='k',linestyle=':')\n",
    "    plt.text(0,17,\"healthcare system capacity\")\n",
    "    plt.ylabel('Number of infected individuals')\n",
    "    plt.xlabel('time')\n",
    "    plt.title('Diffusion of infection')\n",
    "    plt.legend()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {},
   "outputs": [],
   "source": [
    "import gif"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {},
   "outputs": [],
   "source": [
    "@gif.frame \n",
    "def compare(percent=0.8):\n",
    "    Wtt_sir,Penc1 = sm.symbols('Wtt_sir,Penc1',cls=sm.Function)\n",
    "    Penc1= (Wt/N)*((N - Rt - Wt)/N)\n",
    "    Wtt_sir=N*c*Penc1*tau + Wt - a*Wt\n",
    "    Wupdate_sir=Wtt_sir.subs([(N,100),(tau,0.5),(c,0.5),(a,0.1)])\n",
    "    Rupdate= Rtt.subs([(a,0.1)])\n",
    "    ws_sir=[1]\n",
    "    Rt_v = 0\n",
    "    for t in range(0,149):\n",
    "        ws_sir.append(Wupdate_sir.subs([(Wt,ws_sir[-1:][0]),(Rt,Rt_v)]))\n",
    "        Rt_v = Rupdate.subs([(Wt,ws_sir[-1:][0]),(Rt,Rt_v)])\n",
    "    #plt.plot(range(0,80),ws_sir)\n",
    "    #plt.ylabel('Number of infected individuals')\n",
    "    #plt.xlabel('time')\n",
    "    #plt.title('Diffusion of infection')\n",
    "\n",
    "    Wtt_sir2,Penc2 = sm.symbols('Wtt_sir2,Penc2',cls=sm.Function)\n",
    "    Penc2= (Wt/N)*((N - Rt - Wt)/N) * percent\n",
    "    Wtt_sir2=N*c*Penc2*tau + Wt - a*Wt\n",
    "    Wupdate_sir2=Wtt_sir2.subs([(N,100),(tau,0.5),(c,0.5),(a,0.1)])\n",
    "\n",
    "    ws_sir2=[1]\n",
    "    Rt_v = 0\n",
    "    for t in range(0,149):\n",
    "        ws_sir2.append(Wupdate_sir2.subs([(Wt,ws_sir2[-1:][0]),(Rt,Rt_v)]))\n",
    "        Rt_v = Rupdate.subs([(Wt,ws_sir2[-1:][0]),(Rt,Rt_v)])\n",
    "    plt.plot(range(0,150),ws_sir)\n",
    "    plt.plot(range(0,150),ws_sir2)\n",
    "    plt.gca().get_lines()[1].set_label('{:.2f} social distancing'.format(1-percent))\n",
    "    plt.plot((0,140),(16,16),color='k',linestyle=':')\n",
    "    plt.text(0,17,\"healthcare system capacity\")\n",
    "    plt.ylabel('Number of infected individuals')\n",
    "    plt.xlabel('time')\n",
    "    plt.title('Diffusion of infection')\n",
    "    plt.legend()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "metadata": {},
   "outputs": [],
   "source": [
    "frames = []\n",
    "for p in range(11):\n",
    "    frame=compare(1 - 0.1*p)\n",
    "    frames.append(frame)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 56,
   "metadata": {},
   "outputs": [],
   "source": [
    "gif.save(frames, \"contagion.gif\", duration=350)"
   ]
  },
  {
   "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.6.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}