💀 🏏 👩🏽‍🚒 Pandemia de COVID-19 através dos olhos de um matemático, ou por que o modelo SEIRD clássico não funciona 💄 💻 🥪

Resumo, ou sobre o lazer de jovens cientistas

Nas últimas semanas, meus colegas e eu terminamos o dia útil competindo na precisão da previsão para o desenvolvimento da epidemia de COVID-19 na Rússia, usando vários métodos de regressão não linear. E se a previsão para amanhã inevitavelmente for boa, a previsão para um período de mais de uma semana reflete a realidade apenas em termos gerais. Parece que tudo está claro: existem modelos epidemiológicos , existem métodos de otimização , existem dados suficientemente detalhados , basta combinar isso e obter uma previsão precisa de um mês ou até seis meses de antecedência. Neste artigo, compartilharei meus pensamentos sobre o que há de errado com o modelo SEIRD clássico e como corrigi-lo. E, claro, vou abrir o véu de segredo que envolve nosso futuro com você.

Sente-se, um matan furioso nos espera para quem sabe quais são as equações diferenciais (para o resto, belas imagens são anexadas).

A figura acima mostra o número total de casos confirmados de COVID-19 em uma escala logarítmica para a Rússia e os três países europeus no top 5 em termos de número de infectados. A explicação está mais adiante no texto.

Minuto de Cuidados com OVNI

A pandemia de COVID-19, uma infecção respiratória aguda potencialmente grave causada pelo coronavírus SARS-CoV-2 (2019-nCoV), foi anunciada oficialmente no mundo. Há muitas informações sobre Habré sobre esse tópico - lembre-se sempre de que pode ser confiável / útil e vice-versa.

Pedimos que você seja crítico com qualquer informação publicada.

Fontes oficiais
C
C
()

, .

Lave as mãos, cuide de seus entes queridos, fique em casa sempre que possível e trabalhe remotamente.

Leia publicações sobre: coronavírus | trabalho remoto

Modelo SEIRD

O modelo epidêmico SEIRD pertence à classe dos chamados modelos compartimentais , cuja essência é dividir a população em vários grupos ( compartimentos ingleses ), no nosso caso:

$inline$ ( Suscetível em inglês ) - suscetível,

$inline$ ( Eng. Exposto) - aqueles que têm a doença no período de incubação,

$inline$ ( Inglês infeccioso) - doente,

$inline$ ( Inglês recuperado) - recuperado,

$inline$ ( Inglês morto) - os mortos. Em seguida, o tamanho de cada grupo é comparado com uma variável no sistema de equações diferenciais, resolvendo quais, você pode prever a dinâmica da epidemia. Existem muitas modificações no modelo SEIRD, por exemplo, o SEIR é um modelo simplificado que não considera separadamente os casos de recuperação e morte. Para se familiarizar com outros modelos, posso recomendar um bom artigo sobre este tópico.

Pouco de teoria

Pela primeira vez, um modelo epidêmico na forma de um sistema de três equações diferenciais para variáveis

$inline$ apareceu no trabalho de W. Kermak e A. McKendrick em 1927.
Essas equações diferenciais têm a forma:

$\begin{align} \frac{dS}{dt }&=-\beta \frac{SI}{N},\\ \frac{dI}{dt}&= \beta \frac{SI}{N}-\gamma I,\\ \frac{dR}{dt}&= \gamma I, \end{align}$

onde, além das variáveis familiares, aparecem as seguintes constantes:

$inline$ - tamanho total da população,

$\beta$ - taxa de transmissão de infecções,

$\gamma$ - velocidade de recuperação.

O significado da equação de Kermak e McKendrick é o seguinte: o número de suscetíveis diminui proporcionalmente ao seu número vezes a proporção média de infectados na população

$inline$ , o número de infectados cresce no mesmo ritmo, ajustado pelo fato de que alguns deles

$\gamma I$ se recuperando e o número de convalescentes aumenta devido a uma diminuição no número de infectados. Vale ressaltar que o modelo SIR contém não linearidade

$inline$ , pelo qual a solução analítica do sistema de equações se torna geralmente impossível, mas, felizmente, os métodos de diferenciação numérica podem lidar facilmente com essa tarefa.

Adicionando outra variável aqui

$inline$ (o número de pessoas com a doença no período de incubação), obtemos o modelo SEIR:

$\begin{align} \frac{dS}{dt }&=-\beta \frac{SI}{N},\\ \frac{dE}{dt}&= \beta \frac{SI}{N}-\kappa E,\\ \frac{dI}{dt}&= \kappa E-\gamma I,\\ \frac{dR}{dt}&= \gamma I, \end{align}$

onde outra constante aparece

$\kappa$ - a taxa de transição da doença do estágio de incubação para o aberto. Figura tirada do artigo .

É imediatamente claro que o modelo SEIR não é muito adequado para descrever o COVID-19, apenas porque neste modelo existem portadores ocultos de infecção

$inline$ não é contagioso. Essa deficiência pode ser corrigida através da introdução, seguindo Pengpeng et al. parâmetro opcional

$\theta$ , caracterizando o grau de infectividade dos portadores latentes de infecção em comparação com os doentes. O modelo SEIR modificado, que tentaremos aplicar à epidemia atual, será semelhante a:

$\begin{align} \frac{dS}{dt }&=-\beta \frac{S(I + \theta E)}{N},\\ \frac{dE}{dt}&= \beta \frac{S(I + \theta E)}{N}-\kappa E,\\ \frac{dI}{dt}&= \kappa E-\gamma I,\\ \frac{dR}{dt}&= \gamma I. \end{align}$

À primeira vista, o modelo resultante promete ser completamente crível.

Experiência numérica com o modelo SEIR

Para modelagem, tentaremos seguir os seguintes parâmetros, focando em dados abertos . Assumindo que a doença dura em média 14 dias (pelo menos quanto tempo dura a forma branda, responsável por até 80% dos casos), encontramos o valor

$\gamma=1/14=0,0714$ . Aceitará

$\beta=3/14=0,2143$ . Magnitude

$\theta=0,6$ emprestado de Pengpeng et al . Dado o período médio de incubação de 3 dias, tome

$\kappa=1/3 = 0,33$ . Aceitamos a população da Rússia igual

$N = 144,5\cdot10^6$ pessoa.

Como condições iniciais, usamos os dados para a Rússia em 2 de abril, quando as medidas adotadas no final de março para limitar a propagação da infecção teriam efeito, a saber:

$\begin{align} &S_0=3548,\\ &I_0= 3283,\\ &E_0=0,5 I_0. \end{align}$

Avaliação

$inline$ a tomamos de forma relativamente arbitrária e isso não importa, porque, como você sabe, algo vai dar errado.

Como resultado da modelagem pelo método de Euler, com uma etapa de 1 dia de 2 a 24 de abril, inclusive, obtemos os gráficos mostrados abaixo: à esquerda em escala linear, à direita em logarítmica.

Marcadores redondos marcam dados reais sobre o número total de casos na Rússia, marcadores quadrados sobre o número de pacientes. À primeira vista, os resultados parecem bons, exceto por um: com os parâmetros do modelo, claramente não adivinhávamos. E aqui os métodos de otimização vêm em nosso auxílio.

Otimize-o

Os métodos de otimização, se o leitor não estiver familiarizado com eles, são algoritmos que permitem encontrar o mínimo de alguma função objetiva. No nosso caso, temos diante de nós o problema da regressão não linear: como escolher o vetor de parâmetros da equação diferencial

$\mathbf{x} = (\beta, \gamma, \kappa, E_0)^\top$ de modo que o conjunto de pontos de solução da equação diferencial

$inline$ estava o mais próximo possível do conjunto de pontos de observação

$\mathcal{F}$ .

Usamos o desvio padrão como uma medida do erro do modelo. A função objetivo assumirá a forma

$f(\mathbf{x}) = \frac{1}{M}\sqrt{\sum_{i=1}^{M}(F_{i} - \mathcal{F}_{i})^2 + \sum_{i=1}^{M}(G_{i} - \mathcal{G}_{i})^2},$

Onde

$inline$ - o número de pontos

$inline$ - o número total de casos de infecção que o modelo apresenta,

$\mathcal{F}$ - o número total real de casos,

$inline$ - o número de pacientes no momento, que fornece o modelo,

$\mathcal{G}$ - o total real dos casos ativos atuais.

Usando a caixa de ferramentas de otimização no MATLAB, ajustaremos os parâmetros do modelo para os dados de observação. Como resultado, obtemos a solução mostrada na figura abaixo.

À primeira vista, está tudo bem. A discrepância acabou sendo igual

$f(\mathbf{x}) = 131,98$ , e "no olho do mar convexo" o ajuste da solução correu bem. Vamos dar uma olhada nos parâmetros recebidos:

$\begin{align} &\beta = 0,374,\\ &\gamma = 0,0117,\\ &E_0 = 7,84\cdot 10^6,\\ &\kappa = 4,81\cdot 10^{-5}. \end{align}$

O valor de quase 8 milhões de pacientes latentes, com cerca de 60 mil casos registrados
em 24 de abril, é algo duvidoso. Também descobrimos que o tempo médio de transição para a fase ativa da doença é

$1/\kappa = 2079$ dias.

Por quê isso aconteceu? Tudo ficará claro se analisarmos o formato da curva em uma longa escala de tempo. Para fazer isso, pegue nosso modelo SEIR com parâmetros "plausíveis" e simule-o por um longo período de tempo (neste experimento, assumi um novo significado

$\beta=0,186$ ):

A curva correspondente ao número total de casos tem um formato S característico em uma escala linear. O programa de otimização tentou fazer uma curva neste formulário. Além disso, a previsão em si com parâmetros "plausíveis" é horrível - segundo ela, quase 90% da população do país ficará doente em setembro - é obviamente irrealista se você olhar para os resultados de outros países (a mesma imagem do começo do artigo, apenas em uma escala linear ):

Aqui eu comparo os três países europeus no top 5 no número de casos e a Rússia. Pode-se ver que estamos cerca de um mês atrás do ritmo da epidemia e que, durante um mês, o crescimento no número total de casos nos três países é quase linear (e até mais lento que o linear), em contraste com os resultados obtidos no modelo SEIR. Isso levanta três questões:

?
SEIR, ?
, , ?

Começarei respondendo à terceira pergunta. Quando prevemos algo, enfrentamos uma tarefa bastante desagradável: os dados nos quais construímos o modelo não são ideais - eles contêm erros, ruído e o modelo construído com base em eles também contém algum erro. Quando continuamos a série temporal com nossos pontos de modelo, o erro se acumula - e rapidamente, se prevermos uma função que aumenta no tempo. E este é o caso no nosso caso. Além disso, o modelo é um modelo que reflete a situação real é muito limitado. O súbito desenvolvimento da epidemia em uma nova cidade grande, o uso de um método de tratamento mais eficaz, uma mudança no método de coleta de informações - tudo isso pode cometer tantos erros em dados reais que uma previsão a longo prazo estará completamente longe da realidade.

Muletas e bicicletas: modificamos o modelo SEIR

Vamos tentar responder à pergunta por que o crescimento da epidemia está diminuindo para um linear. Com o número de pessoas infectadas que temos agora, um papel significativo é desempenhado pelas economias de escala associadas à velocidade limitada da comunicação entre as pessoas.

Mais precisamente, lembramos: o número de casos no modelo SEIR é diretamente proporcional ao número médio de casos na população

$inline$ . Essa regra funciona bem em pequenas populações, onde todos podem se comunicar com todos e os pacientes são distribuídos igualmente. Na realidade, especialmente na escala de dezenas e centenas de milhares de pessoas, se você leva duas pessoas aleatoriamente doentes, acontece que elas não apenas nunca se comunicavam e não se viam, nem sequer andavam no mesmo vagão do metrô. E de fato eles vivem em cidades diferentes. Tudo o que os une é uma cadeia de laços sociais, que levou ao fato de que o vírus lhes foi transmitido.

Como exemplo, construí um modelo epidêmico na forma de um autômato celular, em que cada célula interage com apenas quatro vizinhos. Isso é equivalente ao fato de que cada indivíduo da população possui 4 contatos sociais - esse é um número muito pequeno para a população humana, mas quanto mais rápido o efeito de limitar as conexões sociais se manifestará. A cada iteração com uma probabilidade de 0,1, cada um dos 4 vizinhos da célula infectada pode ser infectado. A doença dura em média 14 dias. Os resultados da simulação para um conjunto de células de 200x200 são apresentados na figura abaixo, onde

$inline$ É o número da iteração.

A cor azul indica suscetível, amarelo - doente, verde - recuperado. O mais interessante é como são os gráficos do número de casos. E eles parecem ter sido planejados: após uma curta fase de crescimento subexponencial, assim como no modelo SEIR, há uma fase prolongada de crescimento linear - exatamente como na realidade.

Meu objetivo não era obter uma imagem que se assemelhasse à realidade quantitativamente. Se você quiser mais credibilidade, posso recomendar o projeto de Sergey Potekhin, publicado recentemente em Habré . Para leitores meticulosos, a prova mais rígida de crescimento linear é apresentada abaixo.

Prova do teorema da epidemia de crescimento linear em uma grande população

$inline$ - . . , 20 ( )

$d \approx 4$ .

$inline$ - .

$n^{\frac{1}{d}}$ , ,

$inline$ ,

$inline$ . , ,

$n_{k+1}=\left(n_k^{\frac{1}{d}} + 2P \right)^d,$

$inline$ :

$n_k = n(t_k); n_{k+1} = n(t_k + 1)$ . , :

$n'_{k+1}=n'_{k}\left(1 + \frac{2P}{n_k^{\frac{1}{d}}} \right)^{d-1}.$

$inline$

$n'_{k+1}=n'_{k}$ , , .

4% 16- .

No que diz respeito às estatísticas globais da epidemia, veremos o mesmo: há um mês, o crescimento tem sido linear, apesar de os modelos clássicos nos prometerem uma continuação do aumento exponencial do número de casos.

Tarefas para os curiosos

COVID-19, offline .
.1 ?
.1 2, .

Agora, sobre a modificação do modelo SEIR. A coisa mais simples que podemos fazer é multiplicar o componente não linear por alguma função, dependendo do número de pacientes. Com um pequeno número de casos, essa função deve estar próxima de 1, e com um grande número, deve assintoticamente tender a zero. O candidato adequado mais simples é

$\mathcal{\varphi}(I,E) = e^{-\alpha(I +\theta E)^{K_0}}.$

Seleção de parâmetros

$\alpha$ e

$inline$ pode compensar o crescimento exponencial no modelo original.

Adicione ao modelo, para obter mais informações, e o componente

$inline$ - o número de mortes. Obtenha uma versão modificada do modelo SEIRD:

$\begin{align} \frac{dS}{dt }&=-\beta \frac{S(I + \theta E) \mathcal{\varphi}(I,E)}{N},\\ \frac{dE}{dt}&= \beta \frac{S(I + \theta E) \mathcal{\varphi}(I,E)}{N}-\kappa E,\\ \frac{dI}{dt}&= \kappa E-\gamma I-\mu D,\\ \frac{dR}{dt}&= \gamma I,\\ \frac{dD}{dt}&= \mu I. \end{align}$

Os resultados da simulação são mostrados na figura abaixo.

O erro padrão em comparação com o modelo original é quase inalterado. Os valores dos parâmetros já são realistas. Por conveniência, designei o número inicial de infectados na fase ativa como

$inline$ .

$\ begin {align} & \ beta = 0,219, \\ & \ gamma = 0,0102, \\ & E_0 = 0,13 \ cdot I_0, \\ & \ kappa = 1/3, \\ & \ mu = 1 , 13 \ cdot 10 ^ {- 3}. \ end {align}$

O modelo interpola muito bem os derivativos - os valores do aumento diário no número de casos e no número de mortes.

Vamos tentar fazer uma previsão. Temos um horizonte de previsão de 2 meses e continuamos a modelar soluções com os parâmetros encontrados pelo programa de otimização.

À primeira vista, não é ruim, mas você não pode desejar essa previsão para sua amada terra natal: o número de novos casos continuará a diminuir, mas o número total continuará a crescer. Nesse caso, uma epidemia pode ser interrompida apenas com a ajuda de uma vacina ou aguardando até que quase toda a população esteja doente. O número de novas mortes é fixado em aproximadamente 200 por dia. Esta é uma ilustração clara do que acontecerá se não fortalecermos medidas para combater a epidemia. É isso que nos espera? E por causa desse futuro não muito brilhante, muitos de nós sentamos diligentemente em casa, comprando trigo sarraceno e papel higiênico?

Abaixo, considerarei dois cenários e, olhando para a nebulosa distância dos próximos meses de 28 de abril de 2020, não posso dizer com certeza qual deles desenvolverá mais os eventos. Agora, no momento de romper a curva de novos casos, estamos em um ponto em que é duplamente problemático prever algo.

Script dos EUA

O mundo hegemônico estava em uma posição extremamente invejável. Atrasado na tomada de decisões importantes que retardariam o crescimento da epidemia no início, ele ainda não conseguiu lidar com o aumento natural de novos casos.

O modelo SEIRD modificado, treinado nos primeiros 33 pontos, a partir de 2 de março, prevê mais ou menos realisticamente o curso da epidemia em abril.

Como você pode ver, o crescimento em abril é quase perfeitamente linear. O modelo exagera um pouco a mortalidade em abril, mas o quadro geral está correto.

Esta imagem mostra o aumento diário de novos casos e mortes nos Estados Unidos. É muito semelhante ao que o modelo previu para a Rússia.

Cenário da Alemanha

Alemães disciplinados conseguiram virar a curva a seu favor, e seu crescimento é mais lento que linear. Além disso, para tornar o modelo relevante, tive que adicionar manualmente um aumento na taxa de recuperação em 6 de abril

$\ gama$ 1,7 vezes, caso contrário, uma queda acentuada no número de casos em termos do modelo SEIRD não pode ser explicada.

O modelo foi treinado nos primeiros 27 pontos, a partir de 10 de março. Eu também mudei a função não linear. Para a Alemanha, o expoente dependente do tempo é melhor:

$\ mathcal {\ varphi} (t) = K_0 e ^ {- \ alpha t}.$

Esse tipo de função indica um aumento cumulativo no número de laços sociais interrompidos e, consequentemente, na disseminação da infecção. Aqui você tem uma ilustração clara dos benefícios do auto-isolamento.

O exposto acima mostra o aumento diário de novos casos e mortes. Como no caso dos Estados Unidos, os dados reais contêm flutuações pronunciadas com um período de 7 dias. Isso significa que, nos fins de semana, o número de contatos aumenta e, consequentemente, o número de contatos infectados aumenta. [UPD: pelo contrário, está diminuindo, como o índice de auto-isolamento Yandex pode testemunhar indiretamente . ]

Conclusão

Fazer previsões - tanto a curto quanto a longo prazo - não é apenas um tributo à curiosidade. No caso de uma epidemia, você precisa saber quantas camas devem ser preparadas, quantos ventiladores produzir, por quantos meses estocar equipamentos de proteção individual para os médicos. Os funcionários devem entender se as medidas adotadas são suficientes ou se novas proibições e restrições devem ser introduzidas. Idealmente, o modelo deve refletir a realidade tão bem que seria possível ver a força de ação de cada medida recém-adotada e, em seguida, seria possível fortalecer medidas úteis e cancelar decisões que se mostraram inúteis.

Embora com algumas reservas, todos os que ainda não adquiriram imunidade do COVID-19 podem realmente ser descritos com a letra

$inline$ , toda a variedade de pessoas doentes - a carta

$inline$ etc. Além disso, o modelo SEIRD pode até explicar algo. Mas ela pode prever algo no futuro distante de maneira muito aproximada.

Eu citei deliberadamente no artigo apenas um cenário negativo quando repetimos o destino dos Estados Unidos em termos da dinâmica da epidemia. Se esse cenário for verdadeiro, até o final de junho teremos mais de 300 mil casos registrados da doença e mais de 10 mil mortes. Embora existam pré-requisitos para o fato de que esse cenário não se transforme em realidade, aconselho que você se relacione com ele de acordo com o princípio: "espere o melhor, prepare-se para o pior". Como diz o ditado, se as melhores mentes da NASA enfrentaram a luta contra a epidemia nos Estados Unidos , isso é realmente uma coisa ruim.

Até agora, tudo o que resta para nós é visitar menos lugares públicos, usarcom respiradores adequados , lave as mãos, não se esqueça de limpar o smartphone com álcool depois de tirá-lo da rua e siga outras recomendações simples.

Ainda assim, é possível fazer previsões mais precisas? Sim, claro. Mas sobre isso em outro momento.

Se você deseja brincar com o código-fonte e sugerir um cenário de desenvolvimento, então aqui está o link para o github .

Pandemia de COVID-19 através dos olhos de um matemático, ou por que o modelo SEIRD clássico não funciona