{ "cells": [ { "cell_type": "code", "execution_count": 1, "metadata": { "ExecuteTime": { "end_time": "2021-08-31T16:28:35.316574Z", "start_time": "2021-08-31T16:28:35.300597Z" } }, "outputs": [], "source": [ "%display typeset" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Usando Modelos Populacionais Logísticos para Epidemias\n", "Nesta aula vamos conhecer o modelo logístico de Richards, e estuda sua aplicabilidade no estudo de epidemias. Este notebook se baseia [neste artigo](https://userweb.ucs.louisiana.edu/~xxw6637/papers/JTB2012.pdf)\n", "\n", "O modelo de Richards se baseia em um modelo proposto por Verhulst em 1838.\n", "\n", "$$\\frac{dN(t)}{dt}=r N(t)\\left[1-\\frac{N(t)}{K}\\right]$$\n", "\n", "Richards generalizou este modelo em 1959, adicionando o parâmetro $a$ permitindo o desvio da dependência estrita da densidade." ] }, { "cell_type": "code", "execution_count": 6, "metadata": { "ExecuteTime": { "end_time": "2021-08-30T12:09:51.582016Z", "start_time": "2021-08-30T12:09:51.571761Z" } }, "outputs": [ { "data": { "text/html": [ "" ], "text/latex": [ "$$\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\frac{\\partial}{\\partial t}C\\left(t\\right) = r_{1} C\\left(-\\left(\\frac{C\\left(t\\right)}{K}\\right)^{a} + 1\\right)$$" ], "text/plain": [ "diff(C(t), t) == r1*C(-(C(t)/K)^a + 1)" ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" } ], "source": [ "var('C r1 K a t')\n", "C = function('C')(t)\n", "dcdt = diff(C,t)==r1*C(1-(C/K)^a)\n", "dcdt" ] }, { "cell_type": "code", "execution_count": 12, "metadata": { "ExecuteTime": { "end_time": "2021-08-30T12:12:03.274132Z", "start_time": "2021-08-30T12:12:03.097801Z" } }, "outputs": [ { "data": { "text/html": [ "" ], "text/latex": [ "$$\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\frac{\\int \\frac{\\frac{\\partial}{\\partial t}C\\left(t\\right)}{C\\left(-\\left(\\frac{C\\left(t\\right)}{K}\\right)^{a} + 1\\right)}\\,{d t}}{r_{1}} = C + t$$" ], "text/plain": [ "integrate(1/C(-(C(t)/K)^a + 1), C(t))/r1 == _C + t" ] }, "execution_count": 12, "metadata": {}, "output_type": "execute_result" } ], "source": [ "desolve(dcdt, C, ivar=t)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "O modelo de Richards admite uma solução explícita:\n", "\n", "$$C(t)=K\\left[1+a e^{-a r1(t-t_c)}\\right]^{-1/a},$$\n", "\n", "onde $t_c$ é o ponto no tempo em que a segunda derivada de $C(t)$ se torna 0, ou seja, quando $C(t)=K(1+a)^{-1/a}$. Seja $r=a r1$, então a equação acima se torna:\n", "$$C(t)=K\\left[1+a e^{-r(t-t_c)}\\right]^{-1/a}.$$\n", "\n", "## Relação com o modelo SIR\n", "Seja o modelo SIR,\n", "$$\n", "\\frac{dS}{dt} = -\\beta S I/N \\label{dsdt}\n", "$$\n", "$$\\frac{dI}{dt} = \\beta S I/N -\\gamma I$$\n", "$$\\frac{dR}{dt} = \\gamma I$$\n", "\n", "Onde podemos considerar que $N=S+I$. A partir das equações acima temos que\n", "\n", "$$\\frac{d(S+I)}{ds}=\\frac{\\gamma (S+I)}{\\beta S}$$\n", "\n", "cuja solução é:\n", "\n", "$$S(t)+I(t) = cS(t)^{\\gamma/\\beta} \\label{eq10}$$\n", "\n", "onde\n", "\n", "$$c:=[S(0)+I(0)]S(0)^{\\gamma/\\beta}$$\n", "\n", "combinando (\\ref{dsdt}) e (\\ref{eq10}) obtemos:\n", "\n", "$$\\frac{dS}{dt}=-\\beta S[1-(S/L)^\\alpha] \\label{eq12}$$\n", "\n", "onde $\\alpha:= 1-\\gamma/\\beta$ e $L:=c^{1/(1-\\gamma/\\beta)}$.\n", "\n", "Como a equação (\\ref{eq12}) tem a mesma forma que a equação (\\ref{dsdt}), pode ser resolvida como:\n", "\n", "$$S(t)=L[1+\\alpha e^{b(t-t_j)}]^{-1/\\alpha}$$\n", "\n", "onde $b=\\alpha\\beta$ e $t_j$ o tempo finito em que a segunda derivada de $S(t)$ é $0$.\n", "\n", "Como estamos interessados no numero acumulado de de casos, podemos definir:\n", "\n", "\\begin{align}\n", "J(t) :=& I(t)+R(t)\\nonumber\\\\\n", "=& N-S(t)\\nonumber\\\\\n", "=& N-L[1+\\alpha e^{b(t-t_j)}]^{-1/\\alpha}\n", "\\end{align}\n", "\n", "Sendo $N=S(0)+I(0)$ então $N \\approx L$. Como $I(0)/S(0)\\rightarrow 0$,\n", "\n", "$$J(t)\\approx L-L[1+\\alpha e^{b(t-t_j)}]^{-1/\\alpha}.$$\n", "\n", "Esta última equação nos dá uma forma eficiente de mapear uma curva de casos acumulados ao modelo SIR sem utilizar equações diferenciais.\n", "\n", "Na prática podemos igualar $L$ ao tamanho final da epidemia, e $t_j$ ao ponto de inflexão da curva. Temos ainda que $\\alpha=1-1/{\\cal R}_0$, e $b=\\beta -\\gamma$ a taxa de geração de infeções.\n", "\n", "### Instante do pico, $t_i$\n", "O pico da epidemia ocorre quando $\\dot{I}(t_i)=0$ o que implica que \n", "$$\\frac{\\beta S(t_i)}{S(t_i)+I(t_i)}=\\gamma$$\n", "\n", "Considerando as equações (\\ref{eq10}) e $\\alpha:= 1-\\gamma/\\beta$ e $L:=c^{1/(1-\\gamma/\\beta)}$, podemos obter\n", "\n", "$$S(t_i)=L (\\beta/\\gamma)^{-1/\\alpha}$$\n", "\n", "combinando as duas equações acima temos\n", "\n", "$$I(t_i)=L{\\cal R}_0^{-{\\cal R}_0/({\\cal R}_0-1)}({\\cal R}_0-1)$$\n", "\n", "com \n", "\n", "$${\\cal R}_0=e^{b(t_i-t_j)}$$." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Ajustando o modelo de Richards a Dados\n", "Agora vamos ver como o modelo de Richards se presta para representar uma epidemia com múltiplas ondas. Vamos utilizar os dados da Epidemia da COVID-19 em 2020 e 21" ] }, { "cell_type": "code", "execution_count": 18, "metadata": { "ExecuteTime": { "end_time": "2021-08-31T22:16:09.589023Z", "start_time": "2021-08-31T22:16:09.575167Z" } }, "outputs": [], "source": [ "import pandas as pd\n", "import numpy as np\n", "from matplotlib import pyplot as plt\n", "%matplotlib inline" ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "ExecuteTime": { "end_time": "2021-08-31T16:39:21.134371Z", "start_time": "2021-08-31T16:39:21.100928Z" } }, "outputs": [ { "data": { "text/html": [ "
\n", " | total_cases | \n", "new_cases | \n", "location | \n", "total_deaths | \n", "new_deaths | \n", "total_vaccinations | \n", "new_vaccinations | \n", "stringency_index | \n", "population | \n", "
---|---|---|---|---|---|---|---|---|---|
date | \n", "\n", " | \n", " | \n", " | \n", " | \n", " | \n", " | \n", " | \n", " | \n", " |
2020-01-31 | \n", "2 | \n", "2 | \n", "United Kingdom | \n", "NaN | \n", "NaN | \n", "NaN | \n", "NaN | \n", "8.33 | \n", "68207114 | \n", "
2020-02-01 | \n", "2 | \n", "0 | \n", "United Kingdom | \n", "NaN | \n", "NaN | \n", "NaN | \n", "NaN | \n", "8.33 | \n", "68207114 | \n", "
2020-02-02 | \n", "2 | \n", "0 | \n", "United Kingdom | \n", "NaN | \n", "NaN | \n", "NaN | \n", "NaN | \n", "11.11 | \n", "68207114 | \n", "
2020-02-03 | \n", "8 | \n", "6 | \n", "United Kingdom | \n", "NaN | \n", "NaN | \n", "NaN | \n", "NaN | \n", "11.11 | \n", "68207114 | \n", "
2020-02-04 | \n", "8 | \n", "0 | \n", "United Kingdom | \n", "NaN | \n", "NaN | \n", "NaN | \n", "NaN | \n", "11.11 | \n", "68207114 | \n", "
... | \n", "... | \n", "... | \n", "... | \n", "... | \n", "... | \n", "... | \n", "... | \n", "... | \n", "... | \n", "
2021-08-24 | \n", "6586181 | \n", "30762 | \n", "United Kingdom | \n", "132174.0 | \n", "174.0 | \n", "89865264.0 | \n", "186086.0 | \n", "43.98 | \n", "68207114 | \n", "
2021-08-25 | \n", "6621799 | \n", "35618 | \n", "United Kingdom | \n", "132323.0 | \n", "149.0 | \n", "90095045.0 | \n", "229781.0 | \n", "43.98 | \n", "68207114 | \n", "
2021-08-26 | \n", "6659916 | \n", "38117 | \n", "United Kingdom | \n", "132465.0 | \n", "142.0 | \n", "90295121.0 | \n", "200076.0 | \n", "43.98 | \n", "68207114 | \n", "
2021-08-27 | \n", "6697770 | \n", "37854 | \n", "United Kingdom | \n", "132566.0 | \n", "101.0 | \n", "90466529.0 | \n", "171408.0 | \n", "43.98 | \n", "68207114 | \n", "
2021-08-28 | \n", "6729912 | \n", "32142 | \n", "United Kingdom | \n", "132699.0 | \n", "133.0 | \n", "NaN | \n", "NaN | \n", "NaN | \n", "68207114 | \n", "
576 rows × 9 columns
\n", "