{ "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", "