MM4 - Aula 7 - Espalhamento de Epidemias.sagews

%auto
typeset_mode(True, display=False)

Modelando Doenças Infecciosas

Entendendo Epidemias como uma Reação Química

A modelagem do espalhamento de doenças infecciosas ou epidemias começou clássicamente como extensão dos princípios empregados na modelagem de reações químicas. Assumindo uma população "bem misturada", ou seja, em que a probabilidade de encontro entre qualquer par de indivíduos é igual. 

Tais populações são então divididas em classes imunológicas, e a partir da interação entre indivíduos pertencentes a estas classes, a dinâmica se origina.

O modelo SI

No modelo SI temos apenas duas classes de indivíduos, os saudáveis, mas suscetíveis a contrair uma doença, e os que já foram infectados pela doença e pode portanto, espalhá-la por meio do contato direto. Note que nem todas as doenças se transmitem desta forma, mas neste modelo inicial vamos considerar que sim. O indivíduos transitam entre estes dois estados por meio de eventos de infecção ($S\rightarrow I$) e recuperação ($I\rightarrow S$)

$$S\leftrightharpoons_{\beta}^{\mu}I$$

em termos de reações temos:

$$S+I\,\, \rightarrow \,\, 2I$$

$$I \rightarrow S$$

Escrevendo as Equações

Novamente lançaremos mão da lei de ação de massas e do fato de que o sistema é fechado, ou seja o número de indivíduos não varia.

$$\frac{dS}{dt}= \mu I -\beta S I$$

$$\frac{dI}{dt}= \beta SI - \mu I$$

Utilizamos ps parâmetros $\beta$ e $\mu$ para representar as taxas de infecção e recuperação, respectivamente.

Exercícios:

1. Mostre que a população total é conservada no modelo SI.

2. O Modelo SI apresentado acima poderia ser simplificado para uma única equação diferencial? Se sim, escreva esta equação.

Aplicando Análise Dimensional ao Modelo

Vamos adotar as unidades de tempo em dias. Logo, o lado esquerdo das equações terá unidades de [número de pessoas]/[dias]. Logo, concluímos que $\mu$ possui unidades dias$^{-1}$ enquanto $\beta$ tem unidades de pessoas$^{-1} \times$ dias$^{-1}$, ou seja, é uma taxa de infecção per capita. 

Exercício 3:

Dada a análise dimensional acima, que parâmetros (ou combinação de parâmetros) carrega a informação sobre a estala de tempo? e qual seria uma escala de tempo adequada para descrever os eventos descritos pelo modelo?

Adimensionalizando o modelo

Seja $x^*(t)$ a fração da população no estado infeccioso e  $y^*(t)$ a fração da população total suscetível no tempo $t$. Vamos assumir que estas variáveis adimensionais somam $1$, ou seja, $x^* + y^* = 1$. temos então:

$$y^* = \frac{S}{N},$$

$$x^* = \frac{I}{N},$$

$$t^* = \frac{t}{1/\mu} = \mu t.$$

se $S(t) +I(t) = N,$ então:

$$\frac{S}{N}+\frac{I}{N}=y^* + x^* = \frac{N}{N} = 1.$$

Agora podemos substituir as novas variáveis adimensionais no modelo e obter:

$$\frac{d(y^* N)}{d(t^*/\mu)} = \mu x^* N - \beta (y^*N)(x^* N)$$

$$\frac{d(x^* N)}{d(t^*/\mu)} = \beta(y^* N)(x^* N) - -\mu x^* N$$

Cancelando os fatores comuns $N$ e $\mu$ em ambos os lados das duas equações, chegamos a:

$$\frac{dy^*}{dt^*} =  x^*  - \left(\frac{\beta N}{\mu}\right)x^*y^*$$

$$\frac{dx^*}{dt^*} =  \left(\frac{\beta N}{\mu}\right) x^*y^* -  x^*$$

Agora podemos notar que restou uma razão de parâmetros que vamos denotar por $R_0 = \frac{\beta N}{\mu}$. Este novo parâmetro é muito importante e controla completamente o comportamento qualitativo do modelo, como veremos mais adiante. Por ora, vamos re-escrever o modelo em trmos do $R_0$ e deixar de lado as $*$:

$$\frac{dy}{dt} =  x  - R_0 xy$$

$$\frac{dx}{dt} = R_0 xy -  x$$

Se repetirmos o processo de adimensionalização para o modelo reduzido a uma única equação, derivado no exercício 2, acima, chegamos à seguinte equação:

$$\frac{dx}{dt} =  \left(\frac{\beta N}{\mu}\right) (1-x)x-x$$

ou, em termos do R_0:

$$\frac{dx}{dt} = R_0 (1-x)x-x$$

Exercício 4:

Considere que $1/\mu$ é o tempo típico de recuperação, o que significa o tempo em que a pessoa está doente e pode transmitir a doença para outros. Suponha que inicialmente tenhamos uma única pessoa infectada em uma população de $N$ indivíduos suscetíveis (para $N$ grande $N-1 \approx N$). Explique o significado de R_0 para a dinâmica da epidemia.

Analisando o Modelo

Agora que temos o modelo construído e devidamente simplificado, pessemos à sua análise aplicando as ferramentas que já conhecemos.

Encontrando os Equilíbrios

Vamos partir do modelo reduzido a uma equação, expresso em tremo do $R_0$. Lembre-se que os equilíbrios devem satisfazer $dx/dt=0$.

Exercício 5

Encontre os equilíbrios, se estes existirem, e interprete os resultados.

var('R_0')
f(x) = R_0*(1-x)*x-x
solve(f,x)

$\renewcommand{\Bold}[1]{\mathbf{#1}}\left[x = \frac{R_{0} - 1}{R_{0}}, x = 0\right]$

O equilíbrio em que $x$, nosso número de infectados é 0, é chamado de equilíbrio livre de doença, e neste caso $y=1$.

Exercício 6:

No caso em que $x=(R_0-1)/R_0$, qual o valor de $y$?

var('y')
x=(R_0-1)/R_0
html('se $x$ é')
show(expand(x))
y=1-x
show(html('então $y$ no "equilíbrio endêmico" torna-se:'))
show(expand(y))

<html><font color='black'>se <script type="math/tex">x</script> é</font></html>

$$\renewcommand{\Bold}[1]{\mathbf{#1}}-\frac{1}{R_{0}} + 1$$

<html><font color='black'>então <script type="math/tex">y</script> no "equilíbrio endêmico" torna-se:</font></html>

$$\renewcommand{\Bold}[1]{\mathbf{#1}}$$

$$\renewcommand{\Bold}[1]{\mathbf{#1}}\frac{1}{R_{0}}$$

É importante notar também que o "equilíbrio endêmico" só faz sentido, biologicamente falando, se $x>0$, logo $R_0$ também precisa ser maior que zero. O que nos leva ao seguinte teorema:

Teorema:

Em um modelo SI, uma doença só pode tornar-se endêmica, se $R_0>1$, onde $R_0$ é o numero reprodutivo básico da doença, $R_0=\beta N/\mu$.

Comportamento Qualitativo

Para aplicar nosso inspeção gráfica dos equilíbrios do sistema, precisamos plotar $f(x) = dx/dt$

pl = plot(f(R_0=1.5),(x,-0.05,0.4),ymin=-0.01)
e1 = point((1-1/1.5,0),size=50,color='red') # 1-1/R_0
show(pl+e1)

Simulações

def fun(t,x):
    R_0=1.5
    x=x[0]
    return [R_0*(1-x)*x-x]

T=ode_solver()
T.function = fun
t_span = [0,10]
ic = [1]
T.ode_solve(t_span,ic,num_points=100)

sol = T.plot_solution()
ee = plot(1-1/1.5,color='green')
show (sol+ee)

Error in lines 3-3 Traceback (most recent call last): File "/cocalc/lib/python2.7/site-packages/smc_sagews/sage_server.py", line 995, in execute exec compile(block+'\n', '', 'single') in namespace, locals File "", line 1, in <module> File "/ext/sage/sage-8.0/local/lib/python2.7/site-packages/sage/plot/graphics.py", line 1031, in __radd__ raise TypeError TypeError

Bifurcações

Com vimos acima, o valor de $R_0=1$ parece ser um ponto de bifurcação para este modelo dada a alteração na natureza dos seus equilíbrios. Vamos construir um diagrama de bifurcações em função de $R_0$:

import numpy as np
def drawbif(func,l,u):
    pts = []
    for v in np.linspace(l,u,100):
        g = func(R_0=v)
        xvals = solve(g,x)
        pts.extend([[v,n(i.rhs().real_part())] for i in xvals if n(i.rhs().real_part())>=0])

    show(points(pts),axes_labels=['$R_0$','$x$'],gridlines=True, xmin=0)
var('R_0')
f(x) = R_0*(1-x)*x - x
drawbif(f,0,2)

R_0

Exercício 7:

Interprete a bifurcação acima e diga que tipo de bifurcação dentre as estudadas nós observamos acima

O Modelo SIR

Se extendermos o cenário epidemico que deu origem ao modelo SI, com a introdução da possibilidade dos indíviduos acometidos pela doença tornarem-se imunes à doença, obtemos outro modelo epidemiológico clássico, o modelo SIR., que pode ser descrito por meio do seguinte sistema de EDOs:

$$\frac{dS}{dt} = - \beta I S$$

$$\frac{dI}{dt} = \beta I S - \mu I$$

$$\frac{dR}{dt} = \mu I$$

Exercício 8:

Interprete as equações acima e esboçe um diagrama de blocos para o modelo SIR

Exercício 9:

Explique porque o modelo acima pode ser estudado como um sistema de apenas duas variáveis. Qual variável não encontra-se acoplada às outras duas?

Exercício 10:

Divida $dI/dt$ por $dS/dt$ e simplifique a equação para obter uma EDO para I em função de S. Resolva a equação, obtendo a seguinte solução: $I(S)=-S+\frac{\mu}{\beta}ln S + K$,  plote a solução e interprete.

var('beta I S mu t')
S = function('S')(t)
I = function('I')(t)
dsdt = diff(S,t) == -beta*I*S
didt = diff(I,t) == beta*I*S - mu*I
dids = didt/dsdt
show(dids)
html("Agora vamos simplificar esta equação:")
dids2 = simplify(expand(dids))
show(dids2)

(beta, I, S, mu, t)

$\displaystyle \frac{\frac{\partial}{\partial t}I\left(t\right)}{\frac{\partial}{\partial t}S\left(t\right)} = -\frac{\beta I\left(t\right) S\left(t\right) - \mu I\left(t\right)}{\beta I\left(t\right) S\left(t\right)}$

Agora vamos simplificar esta equação:

$\displaystyle \frac{\frac{\partial}{\partial t}I\left(t\right)}{\frac{\partial}{\partial t}S\left(t\right)} = \frac{\mu}{\beta S\left(t\right)} - 1$

Agora podemos resolver facilmente a equação simplificada:

sol = integrate(dids2,S)
show(sol)

$\displaystyle I\left(t\right) = c_{2} + \frac{\mu \log\left(S\left(t\right)\right)}{\beta} - S\left(t\right)$

E, após atribuir valores a $\mu$ e $\beta$ podemos plotar a solução e interpretá-la

mu=0.9
beta = .1
c=10
plot(mu*log(x)/beta -x +c, (1,100), axes_labels=["S","I"])

Simulando o Modelo SIR

Podemos também explorar o modelo por meio de simulações noméricas

def fun (t,y, pars):
    S,I = y
    beta, mu = pars
    return [-beta*I*S,
            beta*I*S - mu*I
            ]

T = ode_solver()
T.function = fun
inits = 500,1
tspan = [0,10]
T.ode_solve(tspan, inits, num_points=500, params=[.008,.09])

T.plot_solution(0,interpolate=True)

T.plot_solution(1,interpolate=True)

Analisando o final da Epidemia

Será que os suscetíveis são necessáriamente consumidos completamente durante uma epidemia?

Para tentar responder a esta pergunta, podemos dividir a 1ª equação do modelo SIR pela 3ª.

A partir de agora vamos interpretar os nossas variáveis $S$, $I$, $R$ como Frações de N, ou seja, $S/N$, $I/N$, $R/N$, Isto faz com que a nossa expressão para $R_0$ se reduza a $\beta/\mu$, assumindo $S \approx N = 1$.

var('beta mu S I R t')
S = function('S',t)
R = function('R',t)
f1 = diff(S,t) == -beta*S*I
f2 = diff(R,t) == mu*I
f1/f2

(beta, mu, S, I, R, t) diff(S(t), t)/diff(R(t), t) == -beta*S(t)/mu

Podemos então reescrever a equação obtida acima como

$$\frac{dS}{dR}==R_0 S$$

f3 = f1/f2


html("Resolvendo a equação, obtemos")

var ('beta X mu Y R_0 S_0')
X = function('X')(Y)
solution = desolve(diff(X,Y)== -R_0*X, X, ivar=Y)
solution

Resolvendo a equação, obtemos

(beta, X, mu, Y, R_0, S_0) _C*e^(-R_0*Y)

Onde $c$ é $S_0$

solution.subs(c=S_0)

_C*e^(-R_0*Y)

Sabendo que $R_0 = \beta S/\mu$, podemos reescrever a solução acima como

$$S = S_0 e^{\frac{R_0 R}{S}}$$

mu=0.9
beta = .5
s0 = 100
R_0 = 1
plot(s0*exp(-R_0*x), (x, 0,100),axes_labels=["R","S"])

f(S,I) = -beta*I*S
g(I,S) = beta*I*S - mu*I
IN = plot3d(g(beta=0.008,mu=0.1),(S,0,10),(I,0,500),alpha=.5, color="red")
SN = plot3d(f(beta=0.008,mu=0.1),(S,0,10),(I,0,500))
show(IN+SN)

3D rendering not yet implemented

var('beta I S mu t')
jack = jacobian([-beta*I*S,beta*I*S - mu*I],[S,I])
show(jack)

($\beta$, $I$, $S$, $\mu$, $t$)

$\displaystyle \left(\begin{array}{rr} -I \beta & -S \beta \\ I \beta & S \beta - \mu \end{array}\right)$

jack.eigenvalues()

[$-\frac{1}{2} \, {\left(I - S\right)} \beta - \frac{1}{2} \, \mu - \frac{1}{2} \, \sqrt{{\left(I^{2} - 2 \, I S + S^{2}\right)} \beta^{2} - 2 \, {\left(I + S\right)} \beta \mu + \mu^{2}}$, $-\frac{1}{2} \, {\left(I - S\right)} \beta - \frac{1}{2} \, \mu + \frac{1}{2} \, \sqrt{{\left(I^{2} - 2 \, I S + S^{2}\right)} \beta^{2} - 2 \, {\left(I + S\right)} \beta \mu + \mu^{2}}$]

matrix(SR,2,2,[[diff(-beta*I*S,S),diff(-beta*I*S,I)],[diff(beta*I*S-mu*I,S),diff(beta*I*S-mu*I,I)]])

$\left(\begin{array}{rr} -I \beta & -S \beta \\ I \beta & S \beta - \mu \end{array}\right)$