{ "cells": [ { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [], "source": [ "import matplotlib.pyplot as plt\n", "plt.rcParams[\"animation.html\"] = \"jshtml\"\n", "import matplotlib.animation\n", "import numpy as np\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Compartmental models\n", "\n", "Comparmental models are a category of infectious diseases in which a population is divided into compartments, with the assumption that all elements in the same compartment present the same characteristics" ] }, { "attachments": { "image.png": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlgAAAA1CAYAAACdk5VxAAAABHNCSVQICAgIfAhkiAAAIABJREFUeF7tXQd8VNXSn00gIUACIXRCCL036V0pIgj6IaCAKL1XRREfYgGeIg9UiiiItCeggFIEUZ4gHWlKEykhhBJKaIaeut9/7u5N7m7uhuTuJmw2M7/fJHfPPefcM//b5szMmUuUTKWxORd8ARwPNgsLBgavgYdo9xd4AjgvWMgxAl7Y9Sp4B/guWO47wcCZa+AGrqE14GaOLznZY0WgFP5/AT4Plnee3HfO3Hf8zjsOfg/sr95hJuvG0/j/PThvzpxkLl3a2+Ttre5Sq8p/QSBtCMTGmik8PMFsNhNfRKfAbcD8EBOyRcAHP1eAn+fifOA8gpAg4AQC19E21tKeXxbvgD90ojtPbtoSwrEi6k9eXmYqXNiE/54sr8iWkQjEQz+PiuJ7jt95YWB+50Xwj1DwUZOJ8owdG2zq3DmIcuSQCy0jz0V26Ds6OoFmzrxAq1ffYnEPghuAE7KD7OmQ8TPUHVUSf57HnVgwHQ2lqiDgCIEz2LHSTGZMqfn53gXMk2ehZASCsckWdn/TM8+YctSrR+TtLfgIAk4hYH7wgBI2bybzgQPcz2FwXb6qPgI3fuutYNNLLxWCEi+WK6dQlsYKArlyeVGLFvnp1KnbFBERVxxFx8BsQhWyIMCYLMuNl+BAE5kCBBVBwEUIFEA/5XBN4THPM+pa4Nku6tpTupkEQZqb2rY15WjUiMRy5Smn9fHKYYL7z6tCBTJfu0Z07VpRjOYkm6ra+/iQmS1XQoKAqxHo3buE2mV7V/edxftrjfHnaAwFyzeLCyLDdz8E+Ole3mLBqojNMu43wsc6ovZkMpkVy5WQIOBiBLyaNFF7bMcKVgmOuRK3oItRlu4UBEqV8lORSNK0BBoFAQWPIDEYy+WQQQhoXM7sEhNKRiCYChY0IRZGMBEEXI6AV8GkOy+YFSxvCWh3OcbSoRUBTWiDPM1srwoFD9Gv5FbJKAQ0kbRy79mC7CVuwYy66qRfzbUl0exyOQgCgoAgIAgIAoKAIOBqBGS5oKsRlf4EAUFAEBAEBAFBINsjIApWtr8EBABBQBAQBAQBQUAQcDUComC5GlHpz50Q4NiTsuDC7jQoGYsgIAgIAoKA5yMgCpbnn+PsKmFzCM5Z5PnzT/vAnO9NS/x5Gi6PBvPnRZaAq1kr8KqrKG1l2RYEBAGXIiD3mEvhlM7cEQFRsNzxrMiYnEWALVffgnuBOd9UJ3A3TaecaPAlcH8wW7eeAHMS1NetdZB9kPZYt+WfICAIuB4Bucdcj6n06GYIyPJdNzshMhyXINAUvdwE8weUmWqA+QsiKg3HRjkwW66YzoGnWLf5nzz8NWDIpiBgAAHOQMJWY57oHNK074vts2C5xwyAKk2yFgJiwXLT8/XEE38Ssz05KrevZ//baDv7frLIb/7uIbv/SoP7gKeCtR+9ZbfgWDCSXesSP/x36+6RQkFAEEgLAvwh8/zgLzWVq2J7MpgVLrnH0oKi1MnSCDhlwUpISKTvv79Jv/xylU6ciCV865By4+Nq5cvnpLp1A6lNmwJUoUJSJu8sDZQMPkshwC6/neD/gPmr5jxjPq2RoAO2/wXeBc4J/h/4AzB/AJa/XFMdrHyxU0gQEAQMIRCDVu+A2TrM1iz+LuIs8ATwfbDcYwBByLMRMKxg3b+fSMOGHafDh+NsELqPW4fLDh+Ooq+/jqI//qjt2QgakE61TAk2BsBLWxO+6GaC+YHO1/hiMMdd9bY25w9P97Buc7Atuy3WgtltyMrZCTC/BIR0EGg4ajTVHTKEAkuXJi984JTpAzf7SPx7ifw+d79x6cDpyUXXIRym3VQK3BCcB7zAui33mCef+SwgW+x77ymj9PmA59YZQ4YVrLlzLyuKlL+/iUaMCKamTf0pKCgnxcSY6ezZB7R//z3FsiVkDAFRvozhhlb+4DLgw9Ye4vGf468cfc38IvbNBo+01vco14WrFY16Q4ZS208/tUIl/wSBRyLAihR/VZnd9C+AWfP1qHuMEVBf1lo07sfFUcQ//9AvZ87QtF276Nrdu9rdsp0NEDCsYG3adE2BZ8qUstSoEb/TLMQT2urV8yrct2+RbAChiOhmCLD1KhHMs2Z2+RUH9wKrFqvl2N4K3gSOBIeAOcD9ezATP/x/sG7LPzsEnujfXynZNGYM7Z85g+ITEjIdI1crjZkuQPY6ICtYM8DrwKrbPVvcY7nxMqxSqJDCPapXpybz59N5KFxC2QcBwwrW1asWE3zNmmz1TR89ykXmaP+ePbfpm28u0pkzsXTtmpkKFDBRrVp5qEuXotSgQbKSp44mOjqOVqy4Tjt2XKewsHh6+JCoSBET1a6dlzp3Lkp16uS1GXhkZCwtWXKF9uy5RZGRicTKYs2aPtStWzA99VS+FEKq4zx4sDatXn0D8WiX6NSpeOVbj7Vr+1KfPiVTjEsbuG4fxK61WjnCQB2EGfCn9ZgpBm5XkF65H9XfY97PCtYl8O9gjv8IALN7kOOtmNhtOBQ8HswzAK67Aqzaifnh/yZYSAeBgpUqKaUHZs98LMqVzpB0i9zNZak7yOxRyApWJzDfbyp57D2mdTcVQEByo5IlaXrbtlQmMJAmtmxJvX+QuVv2uOwtUhpWsFhRYSWLlZ5WrXixSMbSqlXX6cMPL9gc5MYNM23efBccliLW68iRe/Taa6fp1i2LIqg25DH//PMdhbUKzZ49d2jUqDCKZ4eSlWJjCa7OWHA49e9fhIYOZWNISpoy5TytXKmu+CfiSf2+fTHgMHr77ZLUtWvBlI2cLHHVMZ2R20kRMqo5x1DNAfOsmbVi+4ShvDowtRWCHJMl5ACBHH6WRStxcZobxUFdKRYEgAAH6bF1OPkBSZQt7rGbCEjecPIkXcf/HX37UusyHLkglJ0QMKxgdehQSAliHzv2LHXqFEgtWgRStWq5KTDQEvTqahAXLuRQGaLnnw+kfv2KwhLlQ/fvm+nIkbv07bfs6UmmGzfioAydVvZXq5aDBgwIoRo18lKePCbifQcP3oO16XJSg6tX42j0aIty1a6dP73ySnEKCclFd+7E0aZN/9CMGZdo/vyr1LBhAFIn2Fq9uBNWrl54IZB69SpKRYv60JUrsbRo0RVYmG7BhXpBsZSVKZNLOR4rdY+yTmlEcbiZnmM66sRZuR31+5jLWcH6BsyrmOyVK6ND41VQ74JXgjkhaZYnrZut5iuvUv0RI6hQlSqUiBnGBcSLbP7Xv+jK4eT0RWp9VXDtb3trUf7Q0tTkzTep7NNPU0BwMMXDdHz5jz9o78yZdGLtGl3s/Arkp3pDR1CFDh0oqEIF8smbl+5euULnd+ygA3Pn0rnt25R22uPaj0k7jtTciCUbNaZGcHGGNGlCfgUK0IObN+n8zp20e9o0uriXDZ+2lFpf2jHZ41D26bbU6PXXqXDVqpSncGG6f/06Xdi9mw58+SWFb/7V/jCe+pvvxy1OCNcKbVkhYyszB8xnOTp01RKLnC+X5R1gL0BI/vw0Btfi02XLUnBAAD3EPfjH5cs0a+9eWn+CDYApKRATnUH16lEH3CsVgoIorw/eO4jx2nH+PM07cIB2nWPjfTLVhyXt9UaNqElICBVA25tY8r8Tdafz9XjR8m5Va+8bNIhqFS1KL65cSWuO6z/uXsA1/W2XLvQn7tEGuD+1lF55tMHm/erUoX5PPEGVChYkP7iQ/CZOTOo6vf1yw6eg1L7VtCnVLV6cvEwmOoF7cC7wWfxnyhRINkK46IdhBWvgwGKwYMXQ+vXR9MMPtxRmYstWnToBxApYw4Yp3XZGx626JEeNKkH581uUuHywTzRrlk9hLS1ZEqUoV7Vq+dC8eZUpR47kdF9FivhS+/bMBZKaLFsWRYhHxJjz0cSJybOM3Ll9oWwVgavQi6ZOvUjffXcZylH5FCK0aJGb3nknNKm8ZMlcNGFCKJS5GNq+/T4tXx5F48dzqI/ryBXHdFZu10nj0p44seF+l/ZoCcxtiT7HgNmluBD8NZhXSWVpajv9E2r42ms2MpRr145CmjWjeXjY3TjNXxtKO5Vt8zS9+P33ioKkkrevL4U+9ZTC2ydPpt/e5ZX6yRTcsBF1W72a8hSxjdlk5axa9+4K2ysvaR+Rbc26g4dQ+9mzycR+fCvxcSt37kyVOnWi9YMH0x/zvzLafVK7uoMG07NffGHTT95ixZTjMLtKHqcHmvEdbMYhON7RKHGw7wbwv8H8VuTUK5xA2NY1YbT3TGhXB+ed6eLt2ymO1hJK1aoXX1QUJJV8vb3pqdBQhSdv304Tf/vNpl093Bc/dOtGRfLYhuewcta9WjWFta7KAXXr0qz27RUFQyVu27lyZeoEl/+Q9etpISZAKi2A8jETz4C+tWs7VLD6YB/TAk07/m1EHvW4szC5GoRnjkraE2ykX5Z79rPPKjlCVGJFq+5zz1Ft6znR7MqQzeSnTDq7Z6WDlZEVKyohmL0wlJmcUGRIcRv+9FM0LEhh9Oabp6G4cLyx81SjhuUCHDjwb1q8+CodPXqXYmP1+96502KNHjaslI1y5WgUu3db6nfrVlS3Sps2Fhfo4cP3dPf36qVv8e7d21K+f79F+dRtbLDQFcd0Vm6DQ8/oZvwgzohI0iHol2fQFcEc08WpHtjc0Rmc/HTEj6xEbLna8eGHNKt8efrQNycthGJ149QpRUFq+vbbSaKwQqBVCtTf2rKAEsHUFbNebnt02TJFQfsolw99ilnzpjfeoHjMmpu/8w6VatY8qd+8UG56/PijolxF7ttHyzt2pKmB+WiSt0lpt/qVVxQrlkqOxpAWhaVIjZrUDlY0Vq4OLVxIszD7n+ztpfw/vHixUv7s55/D4lTN6VPYdNw4pY8/FyygmeXKKcdhuVi+M5uc0TecHlpmd/AlDhjuxEGPoO1OcAlwB/CP4PPgmeCyTvSb4U0DEYPVEQrMov/7P+VYS4+wKMlUDArRiq5dFeVq2dGjVHfePPL/6CMKxSrdN3CNPIAl653mzalxqVJJjQrh3lrXo4eiXO2LjKTnli+nwlOnku+kSUq7VzFRYSuWSlVxX82AssTK1cJDh6jSrFmUC5Mc/r/48GGlnJWQSrCwqrQc44xBnAtb1IpjjPbEZW1gGWJL27cYt0pG5NH2PQCWq+l79lBlTIB8YbnytaZPMNJvRSws+Axys3I19+BBqoD7njGqALnnQykcCutfZpBhC5Y6uHLl/Gj4cL72SxAHXp869YB+/fUWXGRXlfior7664jB2KT0Cvv12GRznJILVExSXHRMr5E2b5qKePUtSvXrJM+Zz5ywrm6pWRdbTNFBEhKV+z54nU63NgfV6pLr/7PeVLm0xCV+8qK8I2tdPz29XHNNZudMzXg+oy7byaDA/idiEyuYW5kXWctYCPgG72nqGLjOOduKB/tt77P200PldO+mnYcPolf/9j0rD4pQeajh6NPni4Xt4yRJa07tXUtPYixdozyfTKQFBjazg1B8+nM7t2K7sb/zGm+QHFwe76Ba3aEbW9FXKvttod2TpNwq7ghqMHElemAWewox9bb++SV3eDDtNa/r0VpS8cs88Qw1GjaIfBw5w6pDsEmT6ddxYuAYtE7gH0bfp1Ib1CgulCwFWphqC2WTDb3zm4eAXwWxFZh/VEjDfn4+V9NI18ICWHztG/4Y1SksjGzakAFh3l0DR6b8m2XV+CffJTCgacVByWDkaVr8+7ba6/F5v3JiC4OJj914rTArMicnvlkuwkH0L5YhZpRENGlAOTBzWY9I0aO3apPJwuMUH4JisqD2DCcBI1BuKiQ5TNFz6q//+m7rBEvZqrVo0xW7cvWC9YsVsNdyXXFclI/IkNcbGF3Ddva0z+TDSL8udE3LzGEfgflcpAnKznEWgqHbExCqjyWkFSztAVngqVvRTuEqV3DRmzFnasCEqXQoWK2l6VL68H61bV4O2bo1G8Hg0HTt2m06fTsAKwYfg0/Txx6HIHB+o1/SRZZprNNW6jsbmqJHGIuuoisvL03PMjJLboFDPoh1r6u5MyU+T5FGyZs/cDcwzbLZy7QWz+YetXG5NhxYtSjG+yN1sMIBQ6TSjs3LCtA+zRD06vnKFomCVxEtCpfJwXTBtgWVLq1zptXe2LLRFC6WLXR9/rNsVK5ssQ+iTT+ruT0/hJcyaQxD70WvLVjry3/8ihmw7Xd6/lxIyWsj0DNJSl1futAN7p79pprXgGBD7VRVsnFAnOTyxmQjmYKdFYHYjZn7+EBxUj85HR9OkbdtslCGux4oN02xYbvVoFeKfWMFqjPgpldrD0sw0YcuWFP3p9dECbkamqYir1KMpmNjwOJ601lPrsJuQFSx2E9orWH2gdDHZuweNyKMdE8dG6ZGRflsiCTITx5jp0VTIneUULK0gjRpZTItXrqTUmFgJYGXl3r1EBJ572ch/6ZLeO8xSxdfXi9q2DVSY6fbtBFq48ApchlE0d+6FJAWrVClvJDtNoL/+up8iFYMe2CEhXnTuXCKC0qtQqVL8pZT0UXj4QyVdhD1xOVNwsK2M9vWM/HbFMZ2V28i4HbQJRbknTOtVZYsVLVYYmTc6kNktiv85m9J7E/PAct16a+JC0jLY/NaHWv/frYHiVm3fpGr91v+qdYf7zG99sF+yKnVpOY7ROv4lLPr7tWO2rhq1v6gjlqD+AGs9o8fhdmwFfHnjRiqMl1Rrq0LHwf5nN2+m3dOnU8TW35zp3pVtJ6OzPq7s8DH0xYYCjuNg5g9MdwfXfAzjSIp98kYcVXVYRD9u04Y4lmo93HrV5szBKvVkPbE0gtuZ9lhzy7HGyKTeL+rvwppYq1Brm32XLF4caxOH/0r4W+KgT1y7plvneJRlHVAJO1fg1vBwOgfFkI/XAvf1trNnlfZPwjXIZZxAVS1TOzYij3ZQYbAu6ZGRfjkejemkA7lPItg9M8ilFiztgI8etcQrFS6sXibJezkQnhWvv/66R/Xr2wbCr12rD7IeGAEB3oj/KqooWKq7i+s1bx4EBSuKvvjiHBSvyuSN+IfUqFmzglCwotDPZXr33dDUquruW7z4IhQsDs2xpUWLLioF/F1GLaWmYNr34eh3eo+p14+zcuv1abAsAu16gvWD2Qx2mgHNXkWfVVLp9z728erFw+C3wPpT01Q6yOxdKac/xkegBo6b8HJJjdTP66RWJ6vvu3r0CM0sFUKVunSl0q1aUQm4eVjZKo94F7barUKQ8l+w6LkBsYLFCXkz7F3gAhl54jISbHlrpuyQrVX8wuGXx3fgZJ93yrqZUpIA994hKEEdEIt4YOBAqoxVca9hFd9/NPGEatC59yPcDuzqehy0EFas92HN7YvYKFWZYosWE8dz2ZOz8miVT23fzvZrP07+bU6vO0qvkzSUGb6punQ5guSbQVhV548UBL5KegYec1RUrOLG40/pMLVuXSjFMOrVC6Aff4xGCoOz9N57ZRWX4vXrscTK1ddfW5a02jfq2fMoPftsUSXWKjjYV0nmeelSHGK9IpWqJUokX4Q9exZG6oRrSIkQi5QOx5HDKgQJQ/PgQ9QmHIfTNNzHqkdOvVBVafvqq4Xx+xqtWXOLbt6MoZdeKk7lyuVCIlNvfMDaTOfPx6AN586KoqVLq9sPjbZtu0+TJ0coaRqKFfOhy5djoaxdUVYQ8r3TvXtyACE3LlnSC30mQt5rSHFRiPz80n8DpfeYKQbtArn1+nSibKkTbTOrqd5Mn12CvIqB7dvsmmBbvCv1lsySzenjRCMuhFMszK5YMc2rD/+JiKCCWM1UvHHTpFQMTg/EQQd3EBQciMDdglWrI11CSpdJoeoWo8dt1NNSIpYYs1Lo4+tDsTFIjqehwNJlHBwNPi28ZI99963CTLkC/KnZ+AnUGCksWuA7aG6iYLEJc7pDIdxjB/ueR+kMhYPbmBeA51u3dao9vqI4XDuv//wzbezZk0Yj5uoTxFYlWK1Y7DrkFAuVENQdfsMSp/eokbLliJW1+lgNZ5+KQa9t5J07VBZJTiuizd4LF1JUqYxgcKZInRWOi6FEvQu3+gu4P0cixQRb1jphOxEv+iU6CpYReVIMSKfASL+8YrMcUrBwsPs+HbkrWeXWOZxLiwwrWOHhCRQebvmgs6MR1aiRkwYPLppid58+JWjjxmjF6tSnj+0y8B49CtGyZSnNmcePx9Px4xaLUIoOUTB8eKmkYv4m4pw55ZE49DTyZMXRyJFn9JoklRUsmJO+/LIcEpOGKUrR9u1hqda339m1a5CSC0tNVaHdP25cSSpb1hLsrpa3bl2QFiyIomnTLimsUnq+P5jeY9qPmX87K7denx5cxhGRFpu+5VM8/ETkJ9Zn4JVgx75tDwZFKxoHj3PepyZjx9K6Af3TJDW3YQXrKawaWvJkC+UbR4+ihJgY4tQPegpPam0jEAfDClbTt96i5c8/l6KquvIvYutWm313kJMoH1Y0lmjSjM5u4awDyVSrj57OnaJrpeDh7Tu0Y9IHioLFH8oWSjMCrFyprg7OdcDMqRtmgP9Ocy+PqeJmfIvwAKxZnCKgd82a9DXi85g48JxzU72JHFhD1q1L0+i4DStYH2ABShssJtEGuet1sA0TGFawxiIesDNWHNoT54hi2op69hQJBXAzXIO8YrB7jRrKijxOIcHfVuR99mREHvs+9H4b6XcLxs0K1hjEe770HRs1bYnxyAxKv+nEOqpVqyrTsGHF8R1CPyX3FVtqmIOCTNS4sR+9/34IlIgqiLFK6S4IDcXKiSWVqEkTP+VzNGyNqlTJG9asEHrjDX0v0bJlFenllwtThQrecPkhIhPMMUQdO2Lp8/JKiL9S332WAdaokQcxVVWRZLSI0jcfh8dXtKgJObD8Yb0qb4NvtWp5kXy0uiJTtWo5Sc0Jx/85WemgQUWIZdajceNClIztfBweF6erqFcvF5S8crpZ3DmHWJ8+hRVL1iOsw3qHU8rSe0xHHTkjt6M+PbR8DuTiQDvWvt8H1wBz0pb/grO9cgUMkKjzP/QAM/Ha/fpR97XriHNi+SNQnh8yvnn8qHidulDAxtDA/cnBrLy6kNuUwnL0Prt2U/l27SlX3jxKm4DgklS9x8vUe+s27j6J2FLG9MTgobiv0z5H5ESnnEi1AlIldJw7jwqULYfVUKT8f27+14rrjq1VXE9LEdY8RJw/i5OU5szhTYFlylLLSZNJVcpsGuDHgH37qcGIkXALVlfqe+NAQRUqUtuZ/F1xaAh2VjL79vI7CQF2MzwJZrMi57foAubZ9GCw2ytXGKNC06zB1rwKUI2x+gRlN5C6pB/cbt8j1xvneiqCmCl2tefFBKIWFLKRUMD2wMWo0mewgHGb5kjdsBXK/dMIevdn6xLacPqEl/DNw197906qz8lK47GaiQO6P8d1Hwqlg+vy/y+RD4qD5uOwn+vpkRrIzq5BdhUysetQj4zIo9ePfZmRfnnhAMvFeb5mwi0fCiVTlZtxyIwAd5aDlVJzlSo58I2/lK4ve0Hlty0CrsjI7umY3ruXgESwSlAxR/W2zKLy8n3CwbPs60kZfGBMqAloNrEbeq5krH2aWjnKQu6oXO3U0X5H5Wq7EvUbKElDH7UCUZu3Skk0iiXj2uB3e+G09Vv9+0ObHF1qXW0dR+OsN2SokmhUb2bD1oANQ4bQwa/m2RyeFaOBWOGkTZ6qVtg7Y4aS1oFJ7/j2cvBvM1yHq/BCPb6KDZ8ZR/+Do9rqCG2Fo2zJuCNlaM+8wpFn3Wx+ueuiI8Xwp0By4ly7irTZyPX6ZKXqb+Sc428SvgCLipqhvQ4WVHDS0GKaxLx67bWJQznR6Gq00Qa/27fR1h+InE+caJQfYvbE7r6hGzbQAqtVzX4/B+tH4qsHnP2diZW7YCzS4BgzPTIiz6Ow4+MY6Zflnm1dpWw/1jn79yflwtJiZV/P0G9M4uKQawy007AFy9CBpZEgkDUR4JiqcWBXKVdZE4U0jDpy3176vHJF2jJ+PEViVvwQMSOsUPB/TiS6DckD5+AzG1q6+PsemlOtCm1HIsArmB3HIm6ELUlsqTq6dCktsqZXUNtsnTCedk6ZQjfDwpTcWumh/V/MoQVIpnoCSuA9rKBiixb//xsf4V0At4G9csV93zh1khbCwhaGVYExiO3ghKk8znWw1P382mjdw3OS1d+R+PEqchzF4Vt0zDdPn1YSms5DhumMVq50B5U1C3kVLqfWd5Vy9VhQ4KDqT2F9YnpDk6bkICyZNZDcdjzSLuzF9j9YaZqAuvyfE4l+ALc2rz7U0n582qY6yiYh9Qd/quYO7gG21nCs0lIk/nxy0SKb+vOgTDRHwlvOCRV1755i0eL/PyDXVTOUO1KuuBNWpL7R5NXiZKmOlCuub0Qem8E6+GGkX5b7mW++IXYX3gVG9/FM4U8Qceb60T/95OBIri0WC5YTeIoF69HgeYgF69GCpr9Gpliw0j8saeEpCHiIBSsjTofLLVgZMUjpM4siIBasLHriZNiCgCAgCAgCgoAgkCUQEBdhljhNMkhBQBAQBAQBQUAQyEoIpH0JTlaSKpPGmp60Cpk0JDmMICAICAKCgCAgCLgBAmLBcoOTIEMQBAQBQUAQEAQEAc9CQBQszzqfIo0gIAgIAoKAICAIuAEComC5wUmQIQgCgoAgIAgIAoKAZyEgCpZnnU+RRhAQBAQBQUAQEATcAAFRsNzgJMgQBAFBQBAQBAQBQcCzEBAFy7POp0gjCAgCgoAgIAgIAm6AgChYbnASZAiCgCAgCAgCgoAg4FkIiILlWedTpBEEBAFBQBAQBAQBN0BAFCw3OAkyBEFAEBAEBAFBQBDwLAREwfKs8ynSCAKCgCAgCAgCgoAbIMAKVkJiotkNhiJD8EQEEhOTpErwRPmckCme28qd5wSC0jRVBJJvPVKuNaEkBBLJLHeeXA8ZhEDySy+BFaxL4eEJ5vh4ze2YQceVbrMfAufOPVCFjsx+0qd9KCtFAAAByElEQVQq8SXee0Oe86mCJDuNI3Ajuance7YwRtK1a2aKF73T+NUlLR0hkHgj6c67yArWxthYMq1efdNRfSkXBAwjsHjxZbXtRsOdeGbDXyFWwu8wYsV4pnwi1WNE4CqOfcpiIMU/OvMYh+KOh94IC5Yp/sABdxybjCmLI5C4a5cqwc/e2PoL3H/Xrts5g4JymipWzEVeXqYsLqIM/3EjcOdOAk2bdpE2bIjmofwJfg0s9prkE3MHmwVjiRqcx0YIbrncj/ukyfE9AoGzkGKpmcywz/CDfAj4uEcI5johjqGr/hQW5kO5c5u8ihUDUvLOcx282bMnc0wMJfzyC5kPHWIA+BoboV5V7fBjJTiPjw+Zy5bNIVdb9rxGXCJ1fLwZz64EM8Ic+DoKA7cBR7ikc8/qxBfi8H3XkcUKAnOBkCBgFIFbaGh1yvNk5n3wRKN9eXi71pBvNTgveXubqUgReed5+AnPUPE47urqVUxrzHwdhYP5ncf/k6gMtuaD2afDN6ewYGD0GojD9XMS/D7YHyzkGAF20/cG7wE/BBvFXNoJdnwN3AavB7cAC6WOQGnsngfmeEi5fwQDZ64BDuhjdzxPaALUy+7/AWUS+CmJqoVxAAAAAElFTkSuQmCC" } }, "cell_type": "markdown", "metadata": {}, "source": [ "## SIR model\n", "\n", "The SIR model is one of the most basic examples of compartmental models. The model consists of three compartments: \n", "- Susceptible (S)\n", "- Infectious (I)\n", "- Recovered (R) (or better not infectious, includes deceased)\n", "\n", "The assumptions underlying this model require the disease to be able to infect through human to human transmission, and to have the recovered individuals acquire lasting immunity.\n", "\n", "![image.png](attachment:image.png)\n", "\n", "A full modeling requires the specification of the transition rates. The S->I infection rate depends on the average number of contacts per person per time, the probability of transmission when in contact with an infected individual, and the number of infected instances (or the probability that any individual is infected).\n", "The I->R rate depends on the rate of recovery and mortality, and the total number of infected individuals." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### ODE formulation of SIR model\n", "\n", "Compartmental models can be formulated in a deterministic manner through differential equations. This makes them less realistic than stochastic models, but this formalism is easier to analyse.\n", "\n", "The ODE SIR model for diseases whose dynamics are faster than the birth and death dynamics of a population (i.e. without vital dynamics) is modeled by the following equations.\n", "\n", "$$\\begin{aligned}\n", "&\\frac{d S}{d t}=-\\frac{\\beta I S}{N}\\\\\n", "&\\frac{d I}{d t}=\\frac{\\beta I S}{N}-\\gamma I\\\\\n", "&\\frac{d R}{d t}=\\gamma I\n", "\\end{aligned}$$\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "If we discretize these equations in order to simulate it we obtain:\n", "$$\\begin{aligned}\n", "&S_{k+1}= S_K (1 -\\frac{\\beta I_k}{N})\\\\\n", "&I_{k+1}=I_k (1 +\\frac{\\beta S_k}{N}-\\gamma)\\\\\n", "&R_{k+1}= R_k + \\gamma I_k\n", "\\end{aligned}$$\n", "\n", "Where $\\beta$ is a product of number of encounters and infection probability." ] }, { "cell_type": "code", "execution_count": 26, "metadata": {}, "outputs": [], "source": [ "def SRI_step(s, i, r, beta=0.5, gamma=0.5):\n", " N = int(s+i+r)\n", " s_new = s*(1 - beta*i/N)\n", " r_new = r + gamma * i\n", " #i_new = i*(1 - gamma + beta*s/N)\n", " i_new = N - s_new - r_new\n", " return s_new, i_new, r_new" ] }, { "cell_type": "code", "execution_count": 29, "metadata": { "scrolled": false }, "outputs": [ { "data": { "text/html": [ "\n", "\n", "\n", "\n", "\n", "\n", "
\n", " \n", "
\n", " \n", "
\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
\n", "
\n", " \n", " \n", " \n", " \n", " \n", " \n", "
\n", "
\n", "
\n", "\n", "\n", "\n" ], "text/plain": [ "" ] }, "execution_count": 29, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAmIAAAHSCAYAAABPdKcOAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3df7CdVWHv/89qQMJPISTttxLaxLkBIXByAocUisiP0CQoDU4nFCpfCEoHFb3gl9sK6FB6aZmxylwByy11JAo2t+GXF1K/WMMPBQs0eCKBImCTKGKQSiTCJaj8kHX/OA/pIRyI5OSwkvB6zWTOftZ+nn3WWcPsefPsZ+9daq0BAOCN9xutJwAA8GYlxAAAGhFiAACNCDEAgEaEGABAI0IMAKCRrVpPYEONHTu2TpgwofU0AADWa8mSJT+ttY5bd3yzDbEJEyakv7+/9TQAANarlPLDoca9NAkA0IgQAwBoRIgBADSy2V4jBgC8tueffz4rV67ML3/5y9ZTedMYPXp0xo8fn6233vrX2l+IAcAWauXKldlxxx0zYcKElFJaT2eLV2vNE088kZUrV2bixIm/1jFemgSALdQvf/nL7LrrriLsDVJKya677vq6zkAKMQDYgomwN9brXW8hBgCMmAsuuCCTJ09OT09Pent7s3jx4ibzWLp0aW688ca12wsXLsynPvWpJMnJJ5+ca6+99hXHfPOb38zRRx89ovNyjRgAMCLuuuuufPWrX813vvOdbLPNNvnpT3+a5557rslcli5dmv7+/rz73e9OksyePTuzZ89uMpfBnBEDAEbEY489lrFjx2abbbZJkowdOzZve9vbMmHChPz0pz9NkvT39+ewww5Lktx2223p7e1Nb29vpk6dmqeffjpJ8ulPfzr77rtvpkyZkrPPPjtJsmLFisyaNSv7779/DjnkkDz00ENJBs5ufehDH8ohhxySPfbYI1/96lfz3HPP5S/+4i9y1VVXpbe3N1dddVW+9KUv5aMf/ejaud58880vO2ZdzzzzTD7wgQ/kgAMOyNSpU3PDDTdslDVyRgwA3gQ+9rFk6dKN+5i9vclFF736/TNmzMj555+fPfbYI0ceeWSOO+64HHrooa+6/4UXXphLL700Bx98cNasWZPRo0fna1/7Wq6//vosXrw42223XVavXp0kOfXUU3PZZZdl0qRJWbx4cU477bTceuutSZKHH344t912W1asWJHDDz88y5cvz/nnn5/+/v787d/+bZLkS1/60st+91DHDHbBBRfkiCOOyLx58/Lkk09m2rRpOfLII7P99ttvwMr9JyEGAIyIHXbYIUuWLMm3vvWtfOMb38hxxx239rqsoRx88ME588wzc8IJJ+SP/uiPMn78+Nx88815//vfn+222y5JMmbMmKxZsyZ33nlnjj322LXHPvvss2tv//Ef/3F+4zd+I5MmTcrb3/72tWfLXsv6jlm0aFEWLlyYCy+8MMnAO1IfeeSR7LXXXq9rTdYlxADgTeC1zlyNpFGjRuWwww7LYYcdln333TdXXHFFttpqq7z44otJ8rKPejj77LPznve8JzfeeGMOPPDA3Hzzzam1vuKdiC+++GJ23nnnLH2VU3zr7v/rvJNxfcfUWnPddddlzz33XO9jvR6uEQMARsT3vve9LFu2bO320qVL87u/+7uZMGFClixZkiS57rrr1t6/YsWK7LvvvjnrrLPS19eXhx56KDNmzMi8efPy85//PEmyevXq7LTTTpk4cWKuueaaJAORdO+99659nGuuuSYvvvhiVqxYke9///vZc889s+OOO6695mwoQx0z2MyZM/O5z30utdYkyT333DPM1Rmw3hArpcwrpTxeSrl/0NiYUspNpZRl3c9duvFSSrmklLK8lHJfKWW/QcfM7fZfVkqZO2h8/1LKv3XHXFJ84AkAbBHWrFmTuXPnZu+9905PT08eeOCB/OVf/mXOO++8nHHGGTnkkEMyatSotftfdNFF2WeffTJlypRsu+22OeqoozJr1qzMnj07fX196e3tXfvS4Pz583P55ZdnypQpmTx58ssunt9zzz1z6KGH5qijjspll12W0aNH5/DDD88DDzyw9mL9dQ11zGDnnntunn/++fT09GSfffbJueeeu1HWqLxUdq+6QynvSrImyZW11n26sU8nWV1r/VQp5ewku9RazyqlvDvJf03y7iS/l+TiWuvvlVLGJOlP0pekJlmSZP9a689KKXcnOSPJvya5MckltdavrW/ifX19tb+/f8P+agB4E3jwwQeHfQ3T5ubkk0/O0UcfnTlz5jSbw1DrXkpZUmvtW3ff9Z4Rq7XenmT1OsPHJLmiu31FkvcOGr+yDvjXJDuXUn47ycwkN9VaV9daf5bkpiSzuvt2qrXeVQeK8MpBjwUAsEXb0Iv1f6vW+liS1FofK6X8Zje+W5IfDdpvZTf2WuMrhxgfUinl1CSnJsnv/M7vbODUAYAt1bofS7Gp29gX6w91fVfdgPEh1Vo/X2vtq7X2jRs3bgOnCACwadjQEPtJ97Jiup+Pd+Mrk+w+aL/xSX68nvHxQ4wDAGzxNjTEFiZ56Z2Pc5PcMGj8pO7dkwcmeap7CfPrSWaUUnbp3mE5I8nXu/ueLqUc2L1b8qRBjwUAsEVb7zVipZR/THJYkrGllJVJzkvyqSRXl1JOSfJIkpc+2vbGDLxjcnmSnyd5f5LUWleXUv4qybe7/c6vtb70BoAPJ/lSkm2TfK37BwCwxVtviNVa/+RV7po+xL41yUde5XHmJZk3xHh/kn3WNw8AYPOzww47ZM2aNa+5z7e+9a186EMfytZbb5277ror22677a/9+Ndff3322GOP7L333ht9Xm8En6wPADQ1f/78/Nmf/VmWLl36uiIsGQixBx54YIRmNvKEGAAw4r75zW/msMMOy5w5c/KOd7wjJ5xwQmqt+cIXvpCrr746559/fk444YQkyWc+85kccMAB6enpyXnnnbf2Ma688sr09PRkypQpOfHEE3PnnXdm4cKF+fM///P09vZmxYoVWbFiRWbNmpX9998/hxxyyNov7/7BD36Qgw46KAcccMBG+1T8jcGXfgPAm8HHPpa8ypdkb7De3tf1beL33HNPvvvd7+Ztb3tbDj744Nxxxx350z/90/zLv/zL2k/DX7RoUZYtW5a77747tdbMnj07t99+e3bddddccMEFueOOOzJ27NisXr06Y8aMyezZs1/2SfrTp0/PZZddlkmTJmXx4sU57bTTcuutt+aMM87Ihz/84Zx00km59NJLN+46DIMQAwDeENOmTcv48QOfWtXb25uHH34473znO1+2z6JFi7Jo0aJMnTo1ycD3VS5btiz33ntv5syZk7FjxyZJxowZ84rHX7NmTe68884ce+yxa8eeffbZJMkdd9yx9gvGTzzxxJx11lkb/w/cAEIMAN4MXseZq5GyzTbbrL09atSovPDCC6/Yp9aac845Jx/84AdfNn7JJZdk4JOuXt2LL76YnXfeOUtf5czf+o5vwTViAMAmY+bMmZk3b97adzQ++uijefzxxzN9+vRcffXVeeKJJ5Ikq1cPfArWjjvumKeffjpJstNOO2XixIm55pprkgxE3b333pskOfjgg7NgwYIkA28O2FQIMQBgkzFjxoy8733vy0EHHZR99903c+bMydNPP53Jkyfnk5/8ZA499NBMmTIlZ555ZpLk+OOPz2c+85lMnTo1K1asyPz583P55ZdnypQpmTx5cm64YeBz4i+++OJceumlOeCAA/LUU0+1/BNfpgx89Nfmp6+vr/b397eeBgBssh588MHstdderafxpjPUupdSltRa+9bd1xkxAIBGhBgAQCNCDACgESEGANCIEAMAaESIAQA0IsQAgBEzatSo9Pb2Zp999skf/uEf5sknn2w9pQ1y2GGHZSQ+NkuIAQAjZtttt83SpUtz//33Z8yYMZvUF24P9RVLbzQhBgC8IQ466KA8+uija7c/85nP5IADDkhPT0/OO++8teNXXnllenp6MmXKlJx44olJkh/+8IeZPn16enp6Mn369DzyyCN56qmnMmHChLz44otJkp///OfZfffd8/zzz2fFihWZNWtW9t9//xxyyCF56KGHkiQnn3xyzjzzzBx++OE566yz8swzz+QDH/hADjjggEydOnXtJ/H/4he/yPHHH5+enp4cd9xx+cUvfjEia+JLvwHgTeBj//yxLP2Pob8Me0P1/j+9uWjWr/dl4r/61a9yyy235JRTTkmSLFq0KMuWLcvdd9+dWmtmz56d22+/PbvuumsuuOCC3HHHHRk7duza75T86Ec/mpNOOilz587NvHnzcvrpp+f666/PlClTctttt+Xwww/PP/3TP2XmzJnZeuutc+qpp+ayyy7LpEmTsnjx4px22mm59dZbkyT//u//nptvvjmjRo3KJz7xiRxxxBGZN29ennzyyUybNi1HHnlk/v7v/z7bbbdd7rvvvtx3333Zb7/9NuravUSIAQAj5he/+EV6e3vz8MMPZ//9988f/MEfJBkIsUWLFmXq1KlJkjVr1mTZsmW59957M2fOnIwdOzZJMmbMmCTJXXfdla985StJkhNPPDEf//jHkyTHHXdcrrrqqhx++OFZsGBBTjvttKxZsyZ33nlnjj322LXzePbZZ9fePvbYYzNq1Ki181i4cGEuvPDCJMkvf/nLPPLII7n99ttz+umnJ0l6enrS09MzIusjxADgTeDXPXO1sb10jdhTTz2Vo48+OpdeemlOP/301Fpzzjnn5IMf/ODL9r/kkktSSlnv4760z+zZs3POOedk9erVWbJkSY444og888wz2XnnnbN06dBnALfffvu1t2utue6667Lnnnu+6u8YSa4RAwBG3Fvf+tZccsklufDCC/P8889n5syZmTdvXtasWZMkefTRR/P4449n+vTpufrqq/PEE08kydqXJn//938/CxYsSJLMnz8/73znO5MkO+ywQ6ZNm5YzzjgjRx99dEaNGpWddtopEydOzDXXXJNkILbuvffeIec1c+bMfO5zn0utNUlyzz33JEne9a53Zf78+UmS+++/P/fdd99ILIsQAwDeGFOnTs2UKVOyYMGCzJgxI+973/ty0EEHZd99982cOXPy9NNPZ/LkyfnkJz+ZQw89NFOmTMmZZ56ZZOBM2Re/+MX09PTky1/+ci6++OK1j3vcccflH/7hH3LcccetHZs/f34uv/zyTJkyJZMnT157Ef66zj333Dz//PPp6enJPvvsk3PPPTdJ8uEPfzhr1qxJT09PPv3pT2fatGkjsiblpQLc3PT19dWR+DwPANhSPPjgg9lrr71aT+NNZ6h1L6UsqbX2rbuvM2IAAI0IMQCARoQYAEAjQgwAtmCb67Xgm6vXu95CDAC2UKNHj84TTzwhxt4gtdY88cQTGT169K99jA90BYAt1Pjx47Ny5cqsWrWq9VTeNEaPHp3x48f/2vsLMQDYQm299daZOHFi62nwGrw0CQDQiBADAGhEiAEANCLEAAAaEWIAAI0IMQCARoQYAEAjQgwAoBEhBgDQiBADAGhEiAEANCLEAAAaEWIAAI0IMQCARoQYAEAjQgwAoBEhBgDQiBADAGhEiAEANCLEAAAaEWIAAI0IMQCARoQYAEAjQgwAoBEhBgDQiBADAGhEiAEANCLEAAAaEWIAAI0IMQCARoQYAEAjQgwAoBEhBgDQiBADAGhEiAEANCLEAAAaEWIAAI0IMQCARoQYAEAjQgwAoJFhhVgp5f8rpXy3lHJ/KeUfSymjSykTSymLSynLSilXlVLe0u27Tbe9vLt/wqDHOacb/14pZebw/iQAgM3DBodYKWW3JKcn6au17pNkVJLjk/xNks/WWicl+VmSU7pDTknys1rrf0ny2W6/lFL27o6bnGRWkv9ZShm1ofMCANhcDPelya2SbFtK2SrJdkkeS3JEkmu7+69I8t7u9jHddrr7p5dSSje+oNb6bK31B0mWJ5k2zHkBAGzyNjjEaq2PJrkwySMZCLCnkixJ8mSt9YVut5VJdutu75bkR92xL3T77zp4fIhjAAC2WMN5aXKXDJzNmpjkbUm2T3LUELvWlw55lftebXyo33lqKaW/lNK/atWq1z9pAIBNyHBemjwyyQ9qratqrc8n+UqS30+yc/dSZZKMT/Lj7vbKJLsnSXf/W5OsHjw+xDEvU2v9fK21r9baN27cuGFMHQCgveGE2CNJDiylbNdd6zU9yQNJvpFkTrfP3CQ3dLcXdtvp7r+11lq78eO7d1VOTDIpyd3DmBcAwGZhq/XvMrRa6+JSyrVJvpPkhST3JPl8kv8/yYJSyl93Y5d3h1ye5MullOUZOBN2fPc43y2lXJ2BiHshyUdqrb/a0HkBAGwuysBJqc1PX19f7e/vbz0NAID1KqUsqbX2rTvuk/UBABoRYgAAjQgxAIBGhBgAQCNCDACgESEGANCIEAMAaESIAQA0IsQAABoRYgAAjQgxAIBGhBgAQCNCDACgESEGANCIEAMAaESIAQA0IsQAABoRYgAAjQgxAIBGhBgAQCNCDACgESEGANCIEAMAaESIAQA0IsQAABoRYgAAjQgxAIBGhBgAQCNCDACgESEGANCIEAMAaESIAQA0IsQAABoRYgAAjQgxAIBGhBgAQCNCDACgESEGANCIEAMAaESIAQA0IsQAABoRYgAAjQgxAIBGhBgAQCNCDACgESEGANCIEAMAaESIAQA0IsQAABoRYgAAjQgxAIBGhBgAQCNCDACgESEGANCIEAMAaESIAQA0IsQAABoRYgAAjQgxAIBGhBgAQCNCDACgESEGANCIEAMAaESIAQA0IsQAABoRYgAAjQgxAIBGhBgAQCNCDACgESEGANDIsEKslLJzKeXaUspDpZQHSykHlVLGlFJuKqUs637u0u1bSimXlFKWl1LuK6XsN+hx5nb7LyulzB3uHwUAsDkY7hmxi5P8c631HUmmJHkwydlJbqm1TkpyS7edJEclmdT9OzXJ3yVJKWVMkvOS/F6SaUnOeyneAAC2ZBscYqWUnZK8K8nlSVJrfa7W+mSSY5Jc0e12RZL3drePSXJlHfCvSXYupfx2kplJbqq1rq61/izJTUlmbei8AAA2F8M5I/b2JKuSfLGUck8p5QullO2T/Fat9bEk6X7+Zrf/bkl+NOj4ld3Yq40DAGzRhhNiWyXZL8nf1VqnJnkm//ky5FDKEGP1NcZf+QClnFpK6S+l9K9ater1zhcAYJMynBBbmWRlrXVxt31tBsLsJ91Ljul+Pj5o/90HHT8+yY9fY/wVaq2fr7X21Vr7xo0bN4ypAwC0t8EhVmv9jyQ/KqXs2Q1NT/JAkoVJXnrn49wkN3S3FyY5qXv35IFJnupeuvx6khmllF26i/RndGMAAFu0rYZ5/H9NMr+U8pYk30/y/gzE3dWllFOSPJLk2G7fG5O8O8nyJD/v9k2tdXUp5a+SfLvb7/xa6+phzgsAYJNXah3ycqxNXl9fX+3v7289DQCA9SqlLKm19q077pP1AQAaEWIAAI0IMQCARoQYAEAjQgwAoBEhBgDQiBADAGhEiAEANCLEAAAaEWIAAI0IMQCARoQYAEAjQgwAoBEhBgDQiBADAGhEiAEANCLEAAAaEWIAAI0IMQCARoQYAEAjQgwAoBEhBgDQiBADAGhEiAEANCLEAAAaEWIAAI0IMQCARoQYAEAjQgwAoBEhBgDQiBADAGhEiAEANCLEAAAaEWIAAI0IMQCARoQYAEAjQgwAoBEhBgDQiBADAGhEiAEANCLEAAAaEWIAAI0IMQCARoQYAEAjQgwAoBEhBgDQiBADAGhEiAEANCLEAAAaEWIAAI0IMQCARoQYAEAjQgwAoBEhBgDQiBADAGhEiAEANCLEAAAaEWIAAI0IMQCARoQYAEAjQgwAoBEhBgDQiBADAGhEiAEANCLEAAAaEWIAAI0IMQCARoQYAEAjQgwAoBEhBgDQyLBDrJQyqpRyTynlq932xFLK4lLKslLKVaWUt3Tj23Tby7v7Jwx6jHO68e+VUmYOd04AAJuDjXFG7IwkDw7a/pskn621TkrysySndOOnJPlZrfW/JPlst19KKXsnOT7J5CSzkvzPUsqojTAvAIBN2rBCrJQyPsl7knyh2y5JjkhybbfLFUne290+pttOd//0bv9jkiyotT5ba/1BkuVJpg1nXgAAm4PhnhG7KMnHk7zYbe+a5Mla6wvd9soku3W3d0vyoyTp7n+q23/t+BDHAABssTY4xEopRyd5vNa6ZPDwELvW9dz3Wses+ztPLaX0l1L6V61a9brmCwCwqRnOGbGDk8wupTycZEEGXpK8KMnOpZStun3GJ/lxd3tlkt2TpLv/rUlWDx4f4piXqbV+vtbaV2vtGzdu3DCmDgDQ3gaHWK31nFrr+FrrhAxcbH9rrfWEJN9IMqfbbW6SG7rbC7vtdPffWmut3fjx3bsqJyaZlOTuDZ0XAMDmYqv17/K6nZVkQSnlr5Pck+TybvzyJF8upSzPwJmw45Ok1vrdUsrVSR5I8kKSj9RafzUC8wIA2KSUgZNSm5++vr7a39/fehoAAOtVSllSa+1bd9wn6wMANCLEAAAaEWIAAI0IMQCARoQYAEAjQgwAoBEhBgDQiBADAGhEiAEANCLEAAAaEWIAAI0IMQCARoQYAEAjQgwAoBEhBgDQiBADAGhEiAEANCLEAAAaEWIAAI0IMQCARoQYAEAjQgwAoBEhBgDQiBADAGhEiAEANCLEAAAaEWIAAI0IMQCARoQYAEAjQgwAoBEhBgDQiBADAGhEiAEANCLEAAAaEWIAAI0IMQCARoQYAEAjQgwAoBEhBgDQiBADAGhEiAEANCLEAAAaEWIAAI0IMQCARoQYAEAjQgwAoBEhBgDQiBADAGhEiAEANCLEAAAaEWIAAI0IMQCARoQYAEAjQgwAoBEhBgDQiBADAGhEiAEANCLEAAAaEWIAAI0IMQCARoQYAEAjQgwAoBEhBgDQiBADAGhEiAEANCLEAAAaEWIAAI0IMQCARoQYAEAjQgwAoJENDrFSyu6llG+UUh4spXy3lHJGNz6mlHJTKWVZ93OXbryUUi4ppSwvpdxXStlv0GPN7fZfVkqZO/w/CwBg0zecM2IvJPlvtda9khyY5COllL2TnJ3kllrrpCS3dNtJclSSSd2/U5P8XTIQbknOS/J7SaYlOe+leAMA2JJtcIjVWh+rtX6nu/10kgeT7JbkmCRXdLtdkeS93e1jklxZB/xrkp1LKb+dZGaSm2qtq2utP0tyU5JZGzovAIDNxUa5RqyUMiHJ1CSLk/xWrfWxZCDWkvxmt9tuSX406LCV3dirjQMAbNGGHWKllB2SXJfkY7XW//Nauw4xVl9jfKjfdWoppb+U0r9q1arXP1kAgE3IsEKslLJ1BiJsfq31K93wT7qXHNP9fLwbX5lk90GHj0/y49cYf4Va6+drrX211r5x48YNZ+oAAM0N512TJcnlSR6stf6PQXctTPLSOx/nJrlh0PhJ3bsnD0zyVPfS5deTzCil7NJdpD+jGwMA2KJtNYxjD05yYpJ/K6Us7cY+keRTSa4upZyS5JEkx3b33Zjk3UmWJ/l5kvcnSa11dSnlr5J8u9vv/Frr6mHMCwBgs1BqHfJyrE1eX19f7e/vbz0NAID1KqUsqbX2rTvuk/UBABoRYgAAjQgxAIBGhBgAQCNCDACgESEGANCIEAMAaESIAQA0IsQAABoRYgAAjQgxAIBGhBgAQCNCDACgESEGANCIEAMAaESIAQA0IsQAABoRYgAAjQgxAIBGhBgAQCNCDACgESEGANCIEAMAaESIAQA0IsQAABoRYgAAjQgxAIBGhBgAQCNCDACgESEGANCIEAMAaESIAQA0IsQAABoRYgAAjQgxAIBGhBgAQCNCDACgESEGANCIEAMAaESIAQA0IsQAABoRYgAAjQgxAIBGhBgAQCNCDACgESEGANCIEAMAaESIAQA0IsQAABoRYgAAjQgxAIBGhBgAQCNCDACgESEGANCIEAMAaESIAQA0IsQAABoRYgAAjQgxAIBGhBgAQCNCDACgESEGANCIEAMAaESIAQA0IsQAABoRYgAAjQgxAIBGhBgAQCNCDACgESEGANDIJhNipZRZpZTvlVKWl1LObj0fAICRtkmEWCllVJJLkxyVZO8kf1JK2bvtrAAARtYmEWJJpiVZXmv9fq31uSQLkhzTeE4AACNqUwmx3ZL8aND2ym4MAGCLtamEWBlirL5ip1JOLaX0l1L6V61a9QZMCwBg5GwqIbYyye6Dtscn+fG6O9VaP19r7au19o0bN+4NmxwAwEjYVELs20kmlVImllLekuT4JAsbzwkAYERt1XoCSVJrfaGU8tEkX08yKsm8Wut3G08LAGBEbRIhliS11huT3Nh6HgAAb5RN5aVJAIA3HSEGANCIEAMAaESIAQA0IsQAABoRYgAAjQgxAIBGhBgAQCNCDACgESEGANCIEAMAaESIAQA0IsQAABoRYgAAjQgxAIBGSq219Rw2SCllVZIfjvCvGZvkpyP8O96MrOvIsK4jw7qODOs6MqzryNgY6/q7tdZx6w5utiH2Riil9Nda+1rPY0tjXUeGdR0Z1nVkWNeRYV1Hxkiuq5cmAQAaEWIAAI0Isdf2+dYT2EJZ15FhXUeGdR0Z1nVkWNeRMWLr6hoxAIBGnBEDAGhEiA2hlDKrlPK9UsryUsrZreezOSulzCulPF5KuX/Q2JhSyk2llGXdz11aznFzU0rZvZTyjVLKg6WU75ZSzujGreswlFJGl1LuLqXc263rf+/GJ5ZSFnfrelUp5S2t57o5KqWMKqXcU0r5ardtXYeplPJwKeXfSilLSyn93ZjngWEqpexcSrm2lPJQ9zx70EiuqxBbRyllVJJLkxyVZO8kf1JK2bvtrDZrX0oya52xs5PcUmudlOSWbptf3wtJ/lutda8kByb5SPffqHUdnmeTHFFrnZKkN8msUsqBSf4myWe7df1ZklMaznFzdkaSBwdtW9eN4/Baa++gj1bwPDB8Fyf551rrO5JMycB/tyO2rkLslaYlWV5r/ZyNuj8AAAL/SURBVH6t9bkkC5Ic03hOm61a6+1JVq8zfEySK7rbVyR57xs6qc1crfWxWut3uttPZ+BJYrdY12GpA9Z0m1t3/2qSI5Jc241b1w1QShmf5D1JvtBtl1jXkeJ5YBhKKTsleVeSy5Ok1vpcrfXJjOC6CrFX2i3JjwZtr+zG2Hh+q9b6WDIQFUl+s/F8NlullAlJpiZZHOs6bN3LZ0uTPJ7kpiQrkjxZa32h28XzwYa5KMnHk7zYbe8a67ox1CSLSilLSimndmOeB4bn7UlWJfli91L6F0op22cE11WIvVIZYsxbS9nklFJ2SHJdko/VWv9P6/lsCWqtv6q19iYZn4Gz43sNtdsbO6vNWynl6CSP11qXDB4eYlfr+vodXGvdLwOX0nyklPKu1hPaAmyVZL8kf1drnZrkmYzwy7tC7JVWJtl90Pb4JD9uNJct1U9KKb+dJN3PxxvPZ7NTStk6AxE2v9b6lW7Yum4k3UsR38zANXg7l1K26u7yfPD6HZxkdinl4Qxc6nFEBs6QWddhqrX+uPv5eJL/nYH/efA8MDwrk6ystS7utq/NQJiN2LoKsVf6dpJJ3Tt63pLk+CQLG89pS7Mwydzu9twkNzScy2anu77m8iQP1lr/x6C7rOswlFLGlVJ27m5vm+TIDFx/940kc7rdrOvrVGs9p9Y6vtY6IQPPp7fWWk+IdR2WUsr2pZQdX7qdZEaS++N5YFhqrf+R5EellD27oelJHsgIrqsPdB1CKeXdGfg/tlFJ5tVaL2g8pc1WKeUfkxyWgW+u/0mS85Jcn+TqJL+T5JEkx9Za172gn1dRSnlnkm8l+bf85zU3n8jAdWLWdQOVUnoycBHuqAz8T+rVtdbzSylvz8CZnDFJ7kny/9Zan203081XKeWwJH9Waz3aug5Pt37/u9vcKsn/qrVeUErZNZ4HhqWU0puBN5a8Jcn3k7w/3XNCRmBdhRgAQCNemgQAaESIAQA0IsQAABoRYgAAjQgxAIBGhBgAQCNCDACgESEGANDI/wUJhCdOsKV6kgAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "%matplotlib inline\n", "\n", "\n", "initial_infected = 5\n", "encounters_per_person = 2\n", "infection_probability = 0.3\n", "beta = encounters_per_person * infection_probability\n", "\n", "D = 5 # Duration of the disease\n", "gamma = 1/D\n", "\n", "N = 10000\n", "\n", "timeline = np.arange(60)\n", "susceptible = [N-initial_infected]\n", "infected = [initial_infected]\n", "recovered = [0]\n", "\n", "\n", "for step in timeline[1:]:\n", " s, i, r = SRI_step(susceptible[-1], infected[-1], recovered[-1], beta=beta, gamma=gamma)\n", " susceptible.append(s)\n", " infected.append(i)\n", " recovered.append(r)\n", "\n", "fig, ax = plt.subplots(figsize=(10, 8))\n", "l1, = ax.plot(timeline, susceptible, color=\"blue\", label=\"Susceptible\")\n", "l2, = ax.plot(timeline, infected, color=\"red\", label=\"Infected\")\n", "l3, = ax.plot(timeline, recovered, color=\"green\", label=\"Recovered\")\n", "ax.legend()\n", "\n", "\n", "def animate_step(i):\n", " l1.set_data(timeline[:i], susceptible[:i])\n", " l2.set_data(timeline[:i], infected[:i])\n", " l3.set_data(timeline[:i], recovered[:i])\n", "\n", "ani = matplotlib.animation.FuncAnimation(fig, animate_step, frames=len(timeline))\n", "ani" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.7.6" }, "toc": { "base_numbering": 1, "nav_menu": {}, "number_sections": true, "sideBar": true, "skip_h1_title": false, "title_cell": "Table of Contents", "title_sidebar": "Contents", "toc_cell": false, "toc_position": {}, "toc_section_display": true, "toc_window_display": false }, "varInspector": { "cols": { "lenName": 16, "lenType": 16, "lenVar": 40 }, "kernels_config": { "python": { "delete_cmd_postfix": "", "delete_cmd_prefix": "del ", "library": "var_list.py", "varRefreshCmd": "print(var_dic_list())" }, "r": { "delete_cmd_postfix": ") ", "delete_cmd_prefix": "rm(", "library": "var_list.r", "varRefreshCmd": "cat(var_dic_list()) " } }, "types_to_exclude": [ "module", "function", "builtin_function_or_method", "instance", "_Feature" ], "window_display": false } }, "nbformat": 4, "nbformat_minor": 4 }