Doenças infecciosas e equações diferenciais não lineares

cupom com desconto - o melhor site de cupom de desconto cupomcomdesconto.com.br


[Esteartigofoipublicadopelaprimeiravezem[Thisarticlewasfirstpublishedon Fabian Dablander, e gentilmente contribuiu para os R-blogueiros]. (Você pode relatar um problema sobre o conteúdo desta página aqui)


Deseja compartilhar seu conteúdo com R-blogueiros? clique aqui se você tiver um blog ou aqui se não tiver.

No verão passado, escrevi sobre casos de amor e equações diferenciais lineares. Enquanto o tópico é alegre, as equações diferenciais lineares são severamente limitadas nos tipos de comportamento que eles podem modelar. Nesta postagem do blog, que passei escrevendo em auto-quarentena para impedir a disseminação adicional do SARS-CoV-2 – entenda isso, alegria -, introduzo equações diferenciais não lineares como um meio para modelar doenças infecciosas. Em particular, discutiremos os modelos simples de SIR e SIRS, os alicerces de muitos dos modelos mais complicados usados ​​em epidemiologia.

Antes de fazê-lo, no entanto, discuto algumas das ferramentas básicas da dinâmica não linear aplicada à equação logística como modelo para o crescimento populacional. Se você já está familiarizado com isso, pode avançar. Se você não teve experiência anterior com equações diferenciais, sugiro que você verifique primeiro meu post anterior sobre o assunto.

Eu devo preceder isso dizendo que não sou um epidemiologista e que nenhuma análise que apresento aqui está especificamente relacionada à atual pandemia de SARS-CoV-2, nem que algo que eu diga seja interpretado como dando conselhos ou fazendo previsões. Estou apenas interessado em equações diferenciais e, como nos casos de amor, as doenças infecciosas são um bom exemplo. Portanto, sem mais delongas, vamos mergulhar!

Modelando o Crescimento Populacional

Antes de começarmos a modelar doenças infecciosas, vale a pena estudar os conceitos necessários para estudar equações diferenciais não lineares em um exemplo simples: modelar o crescimento populacional. Vamos $ N> 0 $ denotar o tamanho de uma população e assumir que seu crescimento depende de si mesmo:

frac {dN} {dt} = ponto {N} = r N espaço em espaço.

Como mostrado em uma postagem anterior, isso leva a um crescimento exponencial de $ r> 0 $:

N

onde $ N_0 = N (0) $ é o tamanho inicial da população no momento $ t = 0 $. A figura abaixo visualiza a equação diferencial (painel esquerdo) e sua solução (painel direito) para $ r = 1 $ e uma população inicial de $ N_0 = 2 $.

trama do pedaço sem nome-pedaço-1

Claramente, este não é um modelo realista, pois o crescimento de uma população depende de recursos finitos. Para modelar recursos finitos, escrevemos:

ponto {N} = rN esquerda (1 – frac {N} {K} direita) espaço em branco,

onde $ r> 0 $ e $ K $ são os chamados capacidade de carga, ou seja, a população de tamanho máximo que pode ser sustentada pelos recursos disponíveis. Observe que, à medida que $ N $ cresce e se $ K> N $, $ (1 – N / K) $ diminui, diminuindo a taxa de crescimento $ dot {N} $. Se, por outro lado, $ N> K $, a população precisa de mais recursos do que os disponíveis, e a taxa de crescimento se torna negativa, resultando em diminuição da população.

Para simplificar, deixe $ K = 1 $ e interprete $ N in [0, 1]$ como a proporção em relação à capacidade de carga; isto é, $ N = 1 $ implica que estamos com capacidade de carga. A figura abaixo visualiza a equação diferencial e sua solução para $ r = 1 $ e uma condição inicial $ N_0 = 0,10 $.

lote de pedaço sem nome-pedaço-2

Em contraste com o crescimento exponencial, a equação logística leva ao crescimento sigmoidal que se aproxima da capacidade de carga. Esse é um comportamento muito mais interessante do que a equação diferencial linear acima permite. Em particular, a equação logística tem dois pontos fixos – pontos nos quais a população não aumenta nem diminui, mas permanece fixa, ou seja, onde $ dot {N} = 0 $. Isso ocorre em $ N = 0 $ e em $ N = 1 $, como pode ser deduzido no painel esquerdo na figura acima.

Analisando a estabilidade de pontos fixos

Qual é a estabilidade desses pontos fixos? Intuitivamente, $ N = 0 $ deve ser instável; se há indivíduos, eles procriam e a população aumenta. Da mesma forma, $ N = 1 $ deve ser estável: se $ N 0 $ e a população crescer em direção a $ N = 1 $, e se $ N> 1 $, então $ dot {N} <0 $ e os indivíduos morrerão até $ N = 1 $.

Para tornar esse argumento mais rigoroso e obter uma avaliação mais quantitativa da rapidez com que as perturbações se afastam ou se aproximam de um ponto fixo, derivamos uma equação diferencial para essas pequenas perturbações próximas ao ponto fixo (ver também Strogatz, 2015, p. 24) Deixe $ N ^ { star} $ denotar um ponto fixo e defina $ eta

frac {d eta} {dt} = frac {d} {dt} left (N

já que $ N ^ { star} $ é uma constante. Isso implica que a dinâmica da perturbação é igual à dinâmica da população. Se $ f (N) $ denota a equação diferencial para $ N $, observe que $ N = N ^ { star} + eta $ é tal que $ dot {N} = dot { eta} = f (N ) = f (N ^ { estrela} + eta) $. Lembre-se de que $ f $ é uma função não linear e funções não-lineares são difíceis de lidar. Assim, simplesmente fingimos que a função é linear próxima ao ponto fixo. Mais precisamente, aproximamos $ f $ em torno do ponto fixo usando uma série de Taylor (veja este excelente vídeo para obter detalhes) escrevendo:

f (N ^ { estrela} + eta) = f (N ^ { estrela}) + eta f ‘(N ^ { estrela}) + mathcal {O} ( eta ^ 2) enspace,

onde ignoramos termos de ordem superior. Observe que, por definição, não há alterações no ponto fixo, ou seja, $ f (N ^ { star}) = 0 $. Supondo que $ f ‘(N ^ { star}) neq 0 $ – caso contrário, os termos de ordem superior importam, pois não haveria mais nada – temos isso perto de um ponto fixo

ponto { eta} approx eta f ‘(N ^ { star}) espaço,

que é uma equação diferencial linear com solução:

eta

Usando esse truque, podemos avaliar a estabilidade de $ N ^ { star} $ da seguinte maneira. Se $ f ‘(N ^ { star}) 0 $, então a pequena perturbação $ eta

%

Conectando nossos dois pontos fixos $ N ^ { star} = 0 $ e $ N ^ { star} = 1 $, descobrimos que $ f ‘(0) = r $ e $ f’ (1) = -r $. Como $ r> 0 $, isso confirma nossa suspeita de que $ N ^ { star} = 0 $ seja instável e $ N ^ { star} = 1 $ seja estável. Além disso, essa análise nos diz com que rapidez as perturbações crescem ou decaem; para a equação logística, isso é dado por $ r $.

Em suma, linearizamos um sistema não linear próximo a pontos fixos para avaliar a estabilidade desses pontos fixos e a rapidez com que as perturbações próximas a esses pontos fixos crescem ou decaem. Essa técnica é chamada análise de estabilidade linear. Nas próximas duas seções, discutiremos duas maneiras de resolver equações diferenciais usando a equação logística como exemplo.

Solução Analítica

Ao contrário das equações diferenciais lineares, que foram o tópico de uma postagem anterior, as equações diferenciais não lineares geralmente não podem ser resolvidas analiticamente; isto é, geralmente não podemos obter uma expressão que, dada uma condição inicial, nos diga o estado do sistema a qualquer momento $ t $. A equação logística pode, no entanto, ser resolvida analiticamente e pode ser instrutivo ver como. Nós escrevemos:

%

Olhando para isso um pouco, percebemos que podemos usar frações parciais para dividir a integral. Nós escrevemos:

%

Os expoentes e os logs se cancelam muito bem. Nós escrevemos:

%

Um último truque é multiplicar por $ e ^ {- rt + Z} $, o que gera:

N = frac { left (e ^ {- rt + Z} right) left (e ^ {rt – Z} right)} { left (e ^ {- rt + Z} right) + { left (e ^ {- rt + Z} right) left (e ^ {- rt + Z} right)}} = frac {1} {1 + e ^ {- rt + Z}} enspace ,

onde $ Z $ é a constante de integração. Para resolver isso, precisamos da condição inicial. Suponha que $ N (0) = N_0 $, que, usando a terceira linha na derivação acima e o fato de que $ t = 0 $, leve a:

%

Conectar isso à nossa solução a partir dos resultados acima:

N

Embora isso tenha sido um aborrecimento, outras equações diferenciais não lineares são muito, muito mais difíceis de resolver, e a maioria não admite uma solução em forma fechada – ou pelo menos se o fizer, a expressão resultante geralmente não é muito intuitiva. Felizmente, podemos calcular a evolução no tempo do sistema usando métodos numéricos, conforme ilustrado na próxima seção.

Solução Numérica

Uma equação diferencial codifica implicitamente como o sistema que modelamos muda com o tempo. Especificamente, dado um estado específico (potencialmente de alta dimensão) do sistema no momento $ t $, $ mathbf {x} _t $, sabemos em que direção e com que rapidez o sistema mudará, porque é exatamente isso que está codificado na equação diferencial $ f = frac { mathrm {d} mathbf {x}} { mathrm {d} t} $. Isso sugere a seguinte aproximação numérica: Suponha que conheçamos o estado do sistema em um ponto de tempo (discreto) $ n $, denotado $ x_n $, e que a alteração no sistema seja constante em um pequeno intervalo $ Delta_t $. Então, a posição do sistema no momento $ n + 1 $ é dada por:

Leia Também  As operações de fusões e aquisições da Fintech atingem a maioridade com um valor de transação de US $ 40 bilhões no primeiro semestre de 2018

mathbf {x} _ {n + 1} = mathbf {x} _n + Delta t cdot f ( mathbf {x} _n) enspace.

$ Delta t $ é um parâmetro importante, codificando em que período de tempo assumimos que a alteração $ f $ seja constante. Podemos codificar isso em R para a equação logística:

solve_logistic <- function(N0, r = 1, delta_t = 0.01, times = 1000) {
  N <- rep(N0, times)
  dN <- function(N) r * N * (1 - N)
  
  for (i in seq(2, times)) {
    # Euler
    N[i] <- N[i-1] + delta_t * dN(N[i-1])
    
    # Improved Euler
    # k <- N[i-1] + delta_t * dN(N[i-1])
    # N[i] <- N[i-1] + 1 /2 * delta_t * (dN(N[i-1]) + dN(k))
    
    # Runge-Kutta 4th order
    # k1 <- dN(N[i-1]) * delta_t
    # k2 <- dN(N[i-1] + k1/2) * delta_t
    # k3 <- dN(N[i-1] + k2/2) * delta_t
    # k4 <- dN(N[i-1] + k3) * delta_t
    #
    # N[i] <- N[i-1] + 1/6 * (k1 + 2*k2 + 2*k3 + k4)
  }
  
  N
}

Claramente, a precisão dessa aproximação é uma função de $ Delta t $. Para ver como, o painel esquerdo mostra a aproximação para vários valores de $ Delta t $, enquanto o painel direito mostra o erro absoluto (log) como uma função de (log) $ Delta t $. O erro é definido como:

E = | N (10) – hat {N} (10) | enspace,

onde $ hat {N} $ é a aproximação de Euler.

lote de pedaço sem nome-pedaço-4

O painel direito mostra aproximadamente o relacionamento:

%

Portanto, o erro diminui linearmente com $ Delta t $. Outros métodos, como o método Euler aprimorado ou os solucionadores Runge-Kutta (consulte o código comentado acima), são melhores. No entanto, não é aconselhável escolher $ Delta t $ extremamente pequeno, porque isso leva a um aumento no tempo de computação e pode levar a erros de precisão que são exacerbados ao longo do tempo.

Em resumo, vimos que equações diferenciais não lineares podem modelar comportamentos interessantes, como vários pontos fixos; como classificar a estabilidade desses pontos fixos usando análise de estabilidade linear; e como resolver numericamente equações diferenciais não lineares. No restante deste post, estudamos equações diferenciais não lineares acopladas – os modelos SIR e SIRS – como uma maneira de modelar a propagação de doenças infecciosas.

Modelando Doenças Infecciosas

Muitos modelos foram propostos como ferramentas para entender epidemias. Nas seções a seguir, enfocarei as duas mais simples: o modelo SIR e SIRS (ver também Hirsch, Smale, Devaney, 2013, cap. 11).

O modelo SIR

Usamos o modelo SIR para entender a propagação de doenças infecciosas. O modelo SIR é o mais básico compartimental modelo, o que significa que agrupa a população geral em subpopulações distintas: uma população suscetível $ S $, uma população infectada $ I $ e uma população recuperada $ R $. Fazemos várias outras suposições simplificadoras. Primeiro, supomos que a população geral seja $ 1 = S + I + R $, de modo que $ S $, $ I $ e $ R $ sejam proporções. Assumimos ainda que a população geral não muda, ou seja,

frac {d} {dt} esquerda (S + I + R direita) = 0 espaço em branco.

Segundo, o modelo SIR pressupõe que, uma vez que uma pessoa tenha sido infectada e se recupere, ela não poderá mais ser infectada novamente – iremos relaxar mais tarde essa suposição. Terceiro, o modelo assume que a taxa de transmissão da doença é proporcional ao número de encontros entre pessoas suscetíveis e infectadas. Modelamos isso definindo

frac {dS} {dt} = – beta IS enspace,

onde $ beta> 0 $ é a taxa de infecção. Quarto, o modelo assume que o crescimento da população recuperada é proporcional à proporção de pessoas infectadas, ou seja,

frac {dR} {dt} = gama I espaço,

onde $ gama> 0 $ é a taxa de recuperação. Como a população geral é constante, essas duas equações naturalmente levam à seguinte equação para os infectados:

begin {alinhado}
frac {d} {dt} left (S + I + R right) = 0 \[0.50em]
frac {dI} {dt} = – frac {dS} {dt} – frac {dR} {dt} \[0.50em]
frac {dI} {dt} = beta IS – gamma I enspace.
end {alinhado}

onde $ beta I S $ fornece a proporção de indivíduos recém-infectados e $ gama I $ fornece a proporção de indivíduos recém-recuperados. Observe que, uma vez que assumimos que a população geral não muda, precisamos apenas focar em dois desses subgrupos, pois $ R

%

Antes de analisarmos este modelo matematicamente, vamos implementar o método de Euler e visualizar algumas trajetórias.

solve_SIR <- function(S0, I0, beta = 1, gamma = 1, delta_t = 0.01, times = 8000) {
  res <- matrix(NA, nrow = times, ncol = 3, dimnames = list(NULL, c('S', 'I', 'R')))
  res[1, ] <- c(S0, I0, 1 - S0 - I0)
  
  dS <- function(S, I) -beta * I * S
  dI <- function(S, I)  beta * I * S - gamma * I
  
  for (i in seq(2, times)) {
    S <- res[i-1, 1]
    I <- res[i-1, 2]
    
    res[i, 1] <- res[i-1, 1] + delta_t * dS(S, I)
    res[i, 2] <- res[i-1, 2] + delta_t * dI(S, I)
  }
  
  res[, 3] <- 1 - res[, 1] - res[, 2]
  res
}
 
 
plot_SIR <- function(res, main = '') {
  cols <- brewer.pal(3, 'Set1')
  matplot(
    res, type = 'l', col = cols, axes = FALSE, lty = 1, lwd = 2,
    ylab = 'Subpopulations
    ylim = c(, 1), main = main, cex.main = 1.75, cex.lab = 1.5,
    font.main = 1, xaxs = 'i', yaxs = 'i'
  )
  
  axis(1, cex.axis = 1.25)
  axis(2, las = 2, cex.axis = 1.25)
  legend(
    3000, 0.65, col = cols, legend = c('S', 'I', 'R'),
    lty = 1, lwd = 2, bty = 'n', cex = 1.5
  )
}

A figura abaixo mostra trajetórias para uma taxa de recuperação fixa de $ gama = 1/8 $ e uma taxa crescente de infecção $ beta $ para a condição inicial $ S_0 = 0,95 $, $ I_0 = 0,05 $ e $ R_0 = 0 $.

lote de pedaço sem nome-pedaço-6

Para $ beta = 1/8 $, não ocorre nenhum surto (painel esquerdo). Em vez disso, a proporção de pessoas suscetíveis e infectadas diminui monotonicamente, enquanto a proporção de pessoas recuperadas aumenta monotonicamente. O painel do meio, por outro lado, mostra um pequeno surto. A proporção de pessoas infectadas aumenta, mas depois cai novamente. Da mesma forma, o painel direito também mostra um surto, mas mais grave, já que a proporção de pessoas infectadas aumenta mais acentuadamente antes que, eventualmente, diminua novamente.

Como as coisas mudam quando alteramos a taxa de recuperação $ gamma $? A figura abaixo mostra novamente três casos de trajetórias para a mesma condição inicial, mas para uma taxa de recuperação menor $ gama = 1/12 $.

lote de pedaço sem nome-pedaço-7

Novamente, não observamos surtos no painel esquerdo e surtos de gravidade crescente no painel do meio e do lado direito. Em contraste com os resultados de $ gama = 1/8 $, o surto é mais grave, como seria de esperar, já que a taxa de recuperação com $ gama = 1/12 $ agora é menor. De fato, se um surto ocorre ou não, e quão grave será, não depende apenas de $ beta $ e $ gamma $, mas da proporção. Essa proporção é conhecida como $ R_0 = beta / gamma $, pronuncia-se “R-zero”. (Observe a escolha infeliz de terminologia bem estabelecida neste contexto, pois $ R_0 $ também denota a proporção inicial de pessoas recuperadas; no entanto, deve ficar claro a partir do contexto a que se destina.) Podemos pensar em $ R_0 $ como o número médio de pessoas que uma pessoa infectada infectará antes de melhorar. Se $ R_0> 1 $, ocorre um surto. Na próxima seção, procuramos os pontos fixos desse sistema e avaliamos sua estabilidade.

Analisando pontos fixos

Uma olhada nas figuras acima sugere que o modelo SIR permite vários estados estáveis. Os painéis da esquerda, por exemplo, mostram que, se não houver surto, a proporção de pessoas suscetíveis permanece acima da proporção de pessoas recuperadas. No entanto, se houver um surto, ele sempre desaparece e a proporção de pessoas recuperadas será maior que a proporção de pessoas suscetíveis; quanto maior depende da gravidade do surto.

Embora pudéssemos brincar um pouco mais com visualizações, vale a pena fazer uma análise formal. Observe que, ao contrário da equação logística, que modelou apenas uma única variável – tamanho da população -, uma análise do modelo SIR exige que lidemos com duas variáveis, $ S $ e $ I $; o terceiro, $ R $, decorre da suposição de um tamanho populacional constante. Nos pontos fixos, nada muda, ou seja, temos:

%

Isso só pode acontecer quando $ I = 0 $, independentemente do valor de $ S $. Em outras palavras, todos os $ (I ^ { star}, S ^ { star}) = (0, S) $ são pontos fixos; se ninguém está infectado, a doença não pode se espalhar – e todo mundo fica suscetível ou recuperado. Para avaliar a estabilidade desses pontos fixos, derivamos novamente uma equação diferencial para as perturbações próximas ao ponto fixo. No entanto, observe que, ao contrário do caso unidimensional estudado acima, as perturbações agora podem ser em relação a $ I $ ou $ S $. Sejam $ u = S – S ^ { star} $ e $ v = I – I ^ { star} $ as respectivas perturbações, e sejam $ dot {S} = f (S, I) $ e $ ponto {I} = g (S, I) $. Primeiro derivamos uma equação diferencial para $ u $, escrevendo:

ponto {u} = frac {d} {dt} esquerda (S – S ^ { estrela} direita) = ponto {S} espaço,

já que $ S ^ { star} $ é uma constante. Isso implica que $ u $ se comporta como $ S $. Em contraste com o caso unidimensional acima, temos dois acoplado equações diferenciais e, portanto, temos que levar em consideração como $ u $ muda em função de $ S $ e $ I $. Nós Taylor expandimos no ponto fixo $ (S ^ { star}, I ^ { star}) $:

%

desde $ f (S ^ { estrela}, I ^ { estrela}) = 0 $ e eliminamos termos de ordem superior. Observe que usar a derivada parcial de $ f $ em relação a $ S $ (ou $ I $) gera uma função e os subscritos $ (S ^ { estrela}, I ^ { estrela}) $ significam que avaliamos esta função no ponto fixo $ (S ^ { star}, I ^ { star}) $. Da mesma forma, podemos derivar uma equação diferencial para $ v $:

Leia Também  nnlib2Rcpp: um pacote R (nother) para redes neurais

ponto {v} aprox u frac { parcial g} { parcial S} _ {(S ^ { estrela}, I ^ { estrela})} + v frac { parcial g} { parcial Eu} _ {(S ^ { estrela}, eu ^ { estrela})} espaço.

Podemos escrever tudo isso de forma concisa usando álgebra matricial:

%

Onde

%

é chamado de Matriz jacobiana no ponto fixo $ (S ^ { star}, I ^ { star}) $. O jacobiano dá a dinâmica linearizada perto de um ponto fixo e, portanto, nos diz como as perturbações evoluirão perto de um ponto fixo.

Ao contrário dos sistemas unidimensionais, nos quais simplesmente verificamos se a inclinação é positiva ou negativa, ou seja, se $ f ‘(x ^ star) 0 $, o teste para determinar se um ponto fixo é estável é um pouco mais complicado em configurações multidimensionais . De fato, e não surpreendentemente, já que temos linearizado Nesta equação diferencial não linear, a verificação é igual à dos sistemas lineares: calculamos os autovalores $ lambda_1 $ e $ lambda_2 $ de $ J $, observando que autovalores negativos significam decaimento exponencial e autovalores positivos significam crescimento exponencial ao longo das direções de os respectivos autovetores. (Observe que isso não funciona para todos os tipos de pontos fixos, consulte Strogatz (2015, p. 152).)

O que isso significa para o nosso modelo SIR? Primeiro, vamos derivar o jacobiano:

%

Avaliando isso no ponto fixo $ (S ^ { estrela}, I ^ { estrela}) = (S, 0) $ resulta em:

%

Como essa matriz é triangular superior – todas as entradas abaixo da diagonal são zero – os valores próprios são dados pela diagonal, ou seja, $ lambda_1 = 0 $ e $ lambda_2 = beta S – gamma $. $ lambda_1 = 0 $ implica uma solução constante, enquanto $ lambda_2> 0 $ implica crescimento exponencial e $ lambda_2 0 $ para $ S> gamma / beta $, o que resulta em um ponto fixo instável. Por outro lado, temos esse $ lambda_2 <0 $ para $ S < gamma / beta $, o que resulta em um ponto fixo estável. Na próxima seção, usaremos campos vetoriais para obter mais intuição para a dinâmica do sistema.

Nullclines e campo de vetor

Um campo vetorial mostra para qualquer posição $ (S, I) $ em que direção o sistema se move, o que indicamos pela ponta de uma seta e com que rapidez, o que indicamos pelo comprimento de uma seta. Usamos o código R abaixo para visualizar um campo vetorial e trajetórias selecionadas nele.

cupom com desconto - o melhor site de cupom de desconto cupomcomdesconto.com.br
library('fields')
 
plot_vectorfield_SIR <- function(beta, gamma, main = '', ...) {
  S <- seq(, 1, 0.05)
  I <- seq(, 1, 0.05)
  
  dS <- function(S, I) -beta * I * S
  dI <- function(S, I)  beta * I * S - gamma * I
  
  SI <- as.matrix(expand.grid(S, I))
  SI <- SI[apply(SI, 1, function(x) sum(x) <= 1), ] # S + I <= 1 must hold
  dSI <- cbind(dS(SI[, 1], SI[, 2]), dI(SI[, 1], SI[, 2]))
  
  draw_vectorfield(SI, dSI, main, ...)
}
 
draw_vectorfield <- function(SI, dSI, main, ...) {
  S <- seq(, 1, 0.05)
  I <- seq(, 1, 0.05)
  
  plot(
    S, I, type = 'n', axes = FALSE, xlab = '', ylab = '', main = '',
    cex.main = 1.5, xlim = c(-0.2, 1), ylim = c(-0.2, 1.2), ...
  )
  
  lines(c(-0.1, 1), c(, ), lwd = 1)
  lines(c(, ), c(-0.1, 1), lwd = 1)
  
  arrow.plot(
    SI, dSI,
    arrow.ex = .075, length = .05, lwd = 1.5, col = 'gray82', xpd = TRUE
  )
  
  cx <- 1.5
  cn <- 2
  
  text(0.5, 1.05, main, cex = 1.5)
  text(0.5, -.075, 'S', cex = cn, font = 1)
  text(-.05, 0.5, 'I', cex = cn, font = 1)
  text(-.03, -.04, , cex = cx, font = 1)
  text(-.03, .975, 1, cex = cx, font = 1)
  text(0.995, -0.04, 1, cex = cx, font = 1)
}

Para $ beta = 1/8 $ e $ gama = 1/8 $, sabemos de cima que nenhum surto ocorre. O campo vetorial mostrado no painel esquerdo abaixo ilustra ainda mais que, como $ S leq gamma / beta = 1 $, todos os pontos fixos $ (S ^ { estrela}, I ^ { estrela}) = (S, 0) $ são estáveis. Por outro lado, sabemos que $ beta = 3/8 $ e $ gamma = 1/8 $ resultam em um surto. O campo vetorial mostrado no painel direito abaixo indica que os pontos fixos com $ S> gama / beta = 1/3 $ são instáveis, enquanto os pontos fixos com $ S <1/3 $ são estáveis; a linha pontilhada é $ S = 1/3 $.

lote de pedaço sem nome-pedaço-9

Podemos encontrar alguma estrutura em tais campos vetoriais? Uma maneira de “organizá-los” é desenhando os chamados nullclines. No nosso caso, o $ I $ -nullcline fornece o conjunto de pontos pelos quais $ dot {I} = 0 $ e o $ S $ -nullcline fornece o conjunto de pontos pelos quais $ dot {S} = 0 $ . Encontramos esses pontos de maneira semelhante à localização de pontos fixos, mas, em vez de definir $ dot {S} $ e $ dot {I} $ como zero, os enfrentamos um de cada vez.

Os nullclines $ S $ são dados pelos eixos $ S $ – e $ I $, porque $ dot {S} = 0 $ quando $ S = 0 $ ou quando $ I = 0 $. No eixo do eixo $ I $, temos $ ponto {I} = – gama I $ desde $ S = 0 $, resultando em deterioração exponencial da população infectada; isso é indicado pelas setas cinza ao longo do eixo $ I $, que são de comprimento progressivamente menor à medida que se aproximam da origem.

Os nullclines $ I $ são dados por $ I = 0 $ e por $ S = gamma / beta $. Para $ I = 0 $, temos $ dot {S} = 0 $ e, portanto, isso gera pontos fixos. Para $ S = gama / beta $, temos $ ponto {S} = – gama I $, resultando em deterioração exponencial da população suscetível, mas desde que $ dot {I} = 0 $, a proporção de infectados pessoas não mudam; isso é indicado no campo vetorial esquerdo acima, onde temos setas horizontais na linha tracejada fornecidas por $ S = gamma / beta $. No entanto, isso é válido apenas por um breve momento, uma vez que $ S $ diminui e por $ S < gamma / beta$ we again have $dot{I} gamma / beta$, which results in $dot{I} > 0 $ e, portanto, a proporção de pessoas infectadas cresce.

Em resumo, vimos como o modelo SIR permite surtos sempre que a taxa de infecção é maior que a taxa de recuperação, $ R_0> beta / gamma $. Se isso ocorrer, temos uma proporção crescente de pessoas infectadas enquanto $ S> gamma / beta $. Como ilustrado pelo campo vetorial, a proporção de pessoas suscetíveis $ S $ diminui com o tempo. Em algum momento, portanto, temos esse $ S < gamma / beta $, resultando em uma diminuição na proporção de pessoas infectadas até finalmente $ I = 0 $. Observe que, no modelo SIR, as infecções sempre desaparecem. Na próxima seção, estendemos o modelo SIR para permitir o estabelecimento de doenças na população.










O modelo SIRS

O modelo SIR pressupõe que as pessoas infectadas ficam imunes à doença para sempre e, portanto, qualquer doença ocorre apenas uma vez e nunca mais volta. Dinâmicas mais interessantes ocorrem quando permitimos a reinfecção de pessoas recuperadas; podemos então perguntar, por exemplo, em que circunstâncias a doença se estabelece na população. O modelo SIRS estende o modelo SIR, permitindo que a população recuperada se torne suscetível novamente (daí o ‘S’ extra). Pressupõe que a população suscetível aumente proporcionalmente à população recuperada, de modo que:

%

onde, como adicionamos $ mu R $ à mudança na proporção de pessoas suscetíveis, tivemos que subtrair $ mu R $ da mudança na proporção de pessoas recuperadas. Novamente assumimos que a população em geral não muda e, portanto, basta estudar o seguinte sistema:

%

desde $ R

solve_SIRS <- function(
  S0, I0, beta = 1, gamma = 1, mu = 1, delta_t = 0.01, times = 1000
) {
  res <- matrix(NA, nrow = times, ncol = 3, dimnames = list(NULL, c('S', 'I', 'R')))
  res[1, ] <- c(S0, I0, 1 - S0 - I0)
  
  dS <- function(S, I, R) -beta * I * S + mu * R
  dI <- function(S, I, R)  beta * I * S - gamma * I
  
  for (i in seq(2, times)) {
    S <- res[i-1, 1]
    I <- res[i-1, 2]
    R <- res[i-1, 3]
    
    res[i, 1] <- res[i-1, 1] + delta_t * dS(S, I, R)
    res[i, 2] <- res[i-1, 2] + delta_t * dI(S, I, R)
    res[i, 3] <- 1 - res[i, 1] - res[i, 2]
  }
  
  res
}
 
 
plot_SIRS <- function(res, main = '') {
  cols <- brewer.pal(3, 'Set1')
  matplot(
    res, type = 'l', col = cols, axes = FALSE, lty = 1, lwd = 2,
    ylab = 'Subpopulations
    main = main, cex.main = 1.75, cex.lab = 1.25, font.main = 1,
    xlim = c(, 4000), font.main = 1, xaxs = 'i', yaxs = 'i'
  )
  
  axis(1, cex.axis = 1.5)
  axis(2, las = 2, cex.axis = 1.5)
  legend(
    3000, 0.95, col = cols, legend = c('S', 'I', 'R'),
    lty = 1, lwd = 2, bty = 'n', cex = 1.5
  )
}

A figura abaixo mostra trajetórias para uma taxa de recuperação fixa de $ gama = 1/8 $, uma taxa de reinfecção fixa de $ mu = 1/8 $ e uma taxa crescente de infecção $ beta $ para a condição inicial $ S_0 = 0,95 $, $ I_0 = 0,05 $ e $ R_0 = 0 $.

lote de pedaço sem nome-pedaço-11

Quanto ao modelo SIR, descobrimos novamente que nenhum surto ocorre por $ R_0 = beta / gamma <1 $, que é o caso do painel esquerdo. O mais interessante, porém, é que descobrimos que a proporção de pessoas infectadas não, em contraste com o modelo SIR, diminua para zero nos outros painéis. Em vez disso, a doença se estabelece na população quando $ R_0> 1 $, e o painel do meio e o direito mostram diferentes pontos fixos.

Leia Também  Obtenha seu ingresso (gratuito) para a conferência virtual e-Rum2020

Como as coisas mudam quando variamos a taxa de reinfecção $ mu $? A figura abaixo mostra novamente três casos de trajetórias para a mesma condição inicial, mas para uma taxa de reinfecção menor $ mu $.

lote de pedaço sem nome-pedaço-12

Mais uma vez, não encontramos surtos no painel esquerdo e surtos de gravidade crescente no painel médio e direito. Ambos os surtos são menos graves em comparação com os surtos nas figuras anteriores, como seria de esperar, devido à diminuição da taxa de reinfecção. Da mesma forma, o sistema parece se estabilizar em diferentes pontos fixos. Na próxima seção, forneceremos uma análise mais formal dos pontos fixos e de sua estabilidade.

Analisando pontos fixos

Para encontrar os pontos fixos do modelo SIRS, novamente buscamos soluções para as quais:

%

onde substituímos $ R = 1 – S – I $ e daí resulta também também $ dot {R} = 0 $, pois assumimos que a população geral não muda. Vimos imediatamente que, em contraste com o modelo SIR, $ I = 0 $ não pode ser um ponto fixo para qualquer $ S $ devido ao termo adicionado que depende de $ mu $. Em vez disso, é um ponto fixo apenas para $ S = 1 $. Para obter o outro ponto fixo, observe que a última equação fornece $ S = gamma / beta $, que conectado à primeira equação produz:

%

Portanto, os pontos fixos são:

%

Observe que o segundo ponto fixo não existe quando $ gamma / beta> 1 $, pois a proporção de pessoas infectadas não pode ser negativa. Outra perspectiva mais intuitiva é escrever $ gamma / beta> 1 $ como $ R_0 = beta / gamma <1 $. Isso nos permite ver que o segundo ponto fixo, que teria uma proporção diferente de zero de pessoas infectadas na população, não existe quando $ R_0 <1 $, pois não ocorre surto. Voltaremos a isso daqui a pouco.

Para avaliar a estabilidade dos pontos fixos, derivamos a matriz jacobiana para o modelo SIRS:

%

Para o ponto fixo $ (S ^ { star}, I ^ { star}) = (1, 0) $, temos:

%

que é novamente triangular superior e, portanto, possui valores próprios $ lambda_1 = – mu $ e $ lambda_2 = beta – gamma $. Isso significa que é instável sempre que $ beta> gama $ desde então $ lambda_2> 0 $, e qualquer indivíduo infectado espalha a doença. O jacobiano no segundo ponto fixo é:

%

o que é mais assustador. No entanto, sabemos desde o post anterior que, para classificar a estabilidade do ponto fixo, basta examinar o traço $ tau $ e o determinante $ Delta $ do jacobiano, dados por

%

O rastreamento pode ser escrito como $ tau = lambda_1 + lambda_2 $ e o determinante pode ser escrito como $ Delta = lambda_1 lambda_2 $, conforme mostrado em uma postagem anterior do blog. Aqui, temos esse $ tau 0 $ porque os dois termos acima são positivos. Isso restringe $ lambda_1 $ e $ lambda_2 $ a serem negativos e, portanto, o ponto fixo é estável.

Campos vetoriais e nulos

Como feito anteriormente para o modelo SIR, podemos visualizar novamente as direções nas quais o sistema muda a qualquer momento usando um campo vetorial.

plot_vectorfield_SIRS <- function(beta, gamma, mu, main = '', ...) {
  S <- seq(, 1, 0.05)
  I <- seq(, 1, 0.05)
  
  dS <- function(S, I) -beta * I * S + mu * (1 - S - I)
  dI <- function(S, I)  beta * I * S - gamma * I
  
  SI <- as.matrix(expand.grid(S, I))
  SI <- SI[apply(SI, 1, function(x) sum(x) <= 1), ] # S + I <= 1 must hold
  dSI <- cbind(dS(SI[, 1], SI[, 2]), dI(SI[, 1], SI[, 2]))
  
  draw_vectorfield(SI, dSI, main, ...)
}

A figura abaixo mostra o campo vetorial para o modelo SIRS, várias trajetórias e os nulos para $ gama = 1/8 $ e $ mu = 1/8 $ para $ beta = 1/8 $ (painel esquerdo) e $ beta = 3/8 $ (painel direito). O painel esquerdo mostra que existe apenas um ponto fixo estável em $ (S ^ { star}, I ^ { star}) = (1, 0) $ para o qual todas as trajetórias convergem.

lote de pedaço sem nome-pedaço-14

O painel direito, por outro lado, mostra dois pontos fixos: um ponto fixo instável em $ (S ^ { star}, I ^ { star}) = (1, 0) $, que só atingimos quando $ I_0 = 0 $ e estável em

(S ^ { estrela}, I ^ { estrela}) = left ( frac {1/8} {3/8}, frac {1/8 left (1 – frac {3/8} {1/8} right)} {1/8 + 1/8} right) = (1/3, 1/3) enspace.

Em contraste com o modelo SIR, portanto, existe um ponto fixo estável que constitui uma população que inclui pessoas infectadas e, portanto, a doença não é erradicada, mas permanece na população.

As linhas tracejadas fornecem os nullclines. O nullcline $ I $ fornece o conjunto de pontos em que $ dot {I} = 0 $, que são – como no modelo SIR acima – dados por $ I = 0 $ e $ S = gamma / beta $. A cláusula $ S $ é dada por:

%

que é uma função não linear em $ S $. Os nullclines nos ajudam novamente a “organizar” o campo vetorial. Isso pode ser visto melhor no painel direito acima. Em particular, e semelhante ao modelo SIR, novamente teremos uma diminuição na proporção de pessoas infectadas à esquerda da linha, dada por $ S = gamma / beta $, ou seja, quando $ S gamma / beta $. Da mesma forma, a proporção de pessoas suscetíveis aumenta quando o sistema está “abaixo” da nulina $ S $, enquanto aumenta quando o sistema está “acima” da nulina $ S $.

Bifurcações

Nos campos vetoriais acima, vimos que o sistema pode deixar de ter apenas um ponto fixo e ter dois pontos fixos. Sempre que um ponto fixo é destruído ou criado ou altera sua estabilidade, como um parâmetro interno é variado – aqui a proporção de $ gamma / beta $ – falamos de um bifurcação.

Como apontado acima, o segundo ponto de equilíbrio existe apenas para $ gamma / beta leq 1 $. Contanto que $ gamma / beta <1 $, tenhamos dois pontos fixos distintos. Em $ gamma / beta = 1 $, o segundo ponto fixo passa a ser:

%

que é igual ao primeiro ponto fixo. Assim, em $ gamma / beta = 1 $, os dois pontos fixos se fundem em um; este é o ponto de bifurcação. Isso faz sentido: se $ gamma / beta 1 $, ocorre um surto, que estabelece a doença na população, pois permitimos reinfecções.

Podemos visualizar essa mudança em pontos fixos no chamado diagrama de bifurcação. Um diagrama de bifurcação mostra como os pontos fixos e sua estabilidade mudam conforme variamos um parâmetro interno. Como lidamos com pontos fixos bidimensionais, dividimos o diagrama de bifurcação em dois: o painel esquerdo mostra como a parte $ I ^ { star} $ do ponto fixo muda à medida que variamos $ gamma / beta $ e o painel direito mostra como a parte $ S ^ { star} $ do ponto fixo muda à medida que variamos $ gamma / beta $.

lote de pedaço sem nome-pedaço-15

O painel esquerdo mostra que, enquanto $ gamma / beta 1 $, temos dois pontos fixos em que o ponto fixo estável é aquele com uma proporção diferente de zero de pessoas infectadas – a doença se estabelece. Esses pontos fixos estão na linha diagonal, indicam como pontos pretos. Curiosamente, isso mostra que a proporção de pessoas infectadas nunca pode ser estável em um valor superior a US $ 1/2 $. Também existem pontos fixos instáveis ​​para os quais $ I ^ { star} = 0 $. Esses pontos fixos são instáveis ​​porque, se existir apenas uma pessoa infectada, ela espalhará a doença, resultando em mais pessoas infectadas. No ponto em que $ beta = gamma $, os dois pontos fixos se fundem: a doença não pode mais ser estabelecida na população e a proporção de pessoas infectadas sempre chega a zero.

Da mesma forma, o painel direito mostra como os pontos fixos $ S ^ { star} $ mudam em função de $ gamma / beta $. Como a infecção se espalha por $ beta> gamma $, o ponto fixo $ S ^ { star} = 1 $ é instável, pois a proporção de pessoas suscetíveis deve diminuir desde que sejam infectadas. Para surtos que se tornam cada vez mais leves como $ gamma / beta rightarrow 1 $, a proporção estável de pessoas suscetíveis aumenta, atingindo $ S ^ { star} = 1 $ quando, finalmente, $ gamma = beta $.

Em resumo, vimos como o SIRS estende o modelo SIR, permitindo reinfecções. This resulted in possibility of more interesting fixed points, which included a non-zero proportion of infected people. In the SIRS model, then, a disease can become established in the population. In contrast to the SIR model, we have also seen that the SIRS model allows for bifuractions, going from two fixed points in times of outbreaks ($beta > gamma$) to one fixed point in times of no outbreaks ($beta < gamma$).

Conclusão

In this blog post, we have seen that nonlinear differential equations are a powerful tool to model real-world phenomena. They allow us to model vastly more complicated behaviour than is possible with linear differential equations, yet they rarely provide closed-form solution. Luckily, the time-evolution of a system can be straightforwardly computed with basic numerical techniques such as Euler’s method. Using the simple logistic equation, we have seen how to analyze the stability of fixed points — simply pretend the system is linear close to a fixed point.

The logistic equation has only one state variable — the size of the population. More interesting dynamics occur when variables interact, and we have seen how the simple SIR model can help us understand the spread of infectious disease. Consisting only of two parameters, we have seen that an outbreak occurs only when $R_0 = beta / gamma > 1$. Moreover, the stable fixed points always included $I = 0$, implying that the disease always gets eradicated. This is not true for all diseases because recovered people might become reinfected. The SIRS model amends this by introducing a parameter $mu$ that quantifies how quickly recovered people can become susceptible again. As expected, this led to stable states in which the disease becomes established in the population.

On our journey to understand these systems, we have seen how to quantify the stability of a fixed point using linear stability analysis, how to visualize the dynamics of a system using vector fields, how nullclines give structure to such vector fields, and how bifurcations can drastically change the dynamics of a system.

The SIR and the SIRS models discussed here are without a doubt crude approximations of the real dynamics of the spread of infectious diseases. There exist several ways to extend them. One way to do so, for example, is to add an exposed population which are infected but are not yet infectious; see here for a visualization of an elaborated version of this model in the context of SARS-CoV-2. These basic compartment models assume homogeneity of spatial-structure, which is a substantial simplification. There are various ways to include spatial structure (e.g., Watts, 2005; Riley, 2007), but that is for another blog post.


I would like to thank Adam Finnemann, Anton Pichlere Oísin Ryan for very helpful comments on this blog post.


Referências

  • Strogatz, S. H. (2015). Nonlinear Dynamics and Chaos: With applications to Physics, Biology, Chemistry, and Engineering. Colorado, US: Westview Press.
  • Hirsch, M. W., Smale, S., & Devaney, R. L. (2013). Differential equations, dynamical systems, and an introduction to chaos. Boston, US: Academic Press.
  • Riley, S. (2007). Large-scale spatial-transmission models of infectious disease. Science, 316(5829), 1298-1301.
  • Watts, D. J., Muhamad, R., Medina, D. C., & Dodds, P. S. (2005). Multiscale, resurgent epidemics in a hierarchical metapopulation model. Proceedings of the National Academy of Sciences, 102(32), 11157-11162.

var vglnk = {key: ‘949efb41171ac6ec1bf7f206d57e90b8’};

(função (d, t) {
var s = d.createElement
s.src = ‘//cdn.viglink.com/api/vglnk.js’;
var r = d.getElementsByTagName
} (documento, ‘script’));

Para Deixe um comentário para o autor, siga o link e comente no blog: Fabian Dablander.

R-bloggers.com oferece atualizações diárias por email sobre notícias e tutoriais do R sobre o aprendizado do R e muitos outros tópicos. Clique aqui se você deseja publicar ou encontrar um emprego em ciência da dados / R.


Deseja compartilhar seu conteúdo com R-blogueiros? clique aqui se você tiver um blog ou aqui se não tiver.



cupom com desconto - o melhor site de cupom de desconto cupomcomdesconto.com.br