⏯️ 🚨 🌜 Pandémie COVID-19 vue par un mathématicien, ou pourquoi le modèle SEIRD classique ne fonctionne pas 🛀 🤵🏼 👨🏾‍🏫

Abstrait, ou sur les loisirs des jeunes scientifiques

Au cours des dernières semaines, mes collègues et moi avons terminé la journée de travail en rivalisant dans l'exactitude des prévisions pour le développement de l'épidémie de COVID-19 en Russie en utilisant diverses méthodes de régression non linéaire. Et si les prévisions pour demain s'avèrent inévitablement bonnes, alors les prévisions pour une période de plus d'une semaine ne reflètent la réalité qu'en termes généraux. Il semblerait que tout soit clair: il y a des modèles épidémiologiques , il y a des méthodes d'optimisation , il y a des données suffisamment détaillées - il suffit de les combiner et d'obtenir à l'avance des prévisions précises pour un mois, voire six mois. Dans cet article, je partagerai mes réflexions sur ce qui ne va pas avec le modèle SEIRD classique et comment y remédier. Et, bien sûr, j'ouvrirai le voile du secret qui enveloppe notre avenir avec vous.

Asseyez-vous, un matan furieux nous attend pour ceux qui savent ce que sont les équations différentielles (pour le reste, de belles photos sont jointes).

La figure ci-dessus montre le nombre total de cas confirmés de COVID-19 sur une échelle logarithmique pour la Russie et les trois pays européens dans le top 5 en termes de nombre de personnes infectées. L'explication se trouve plus loin dans le texte.

UFO Care Minute

La pandémie COVID-19, une infection respiratoire aiguë potentiellement grave causée par le coronavirus SARS-CoV-2 (2019-nCoV), a été officiellement annoncée dans le monde. Il y a beaucoup d'informations sur Habré sur ce sujet - rappelez-vous toujours qu'il peut être à la fois fiable / utile, et vice versa.

Nous vous invitons à critiquer toute information publiée.

Sources officielles
C
C
()

, .

Lavez-vous les mains, prenez soin de vos proches, restez à la maison dans la mesure du possible et travaillez à distance.

Lire les publications sur: coronavirus | travail à distance

Modèle SEIRD

Le modèle épidémique SEIRD appartient à la classe des soi-disant modèles compartimentaux , dont l'essence est de diviser la population en plusieurs groupes ( compartiments anglais ), dans notre cas:

$inline$ ( Sensible à l' anglais ) - sensible,

$inline$ ( Eng. Exposés) - ceux qui ont la maladie pendant la période d'incubation,

$inline$ ( Anglais infectieux) - malade,

$inline$ ( Anglais récupéré) - récupéré,

$inline$ ( Anglais mort) - les morts. Ensuite, la taille de chaque groupe est comparée à une variable du système d'équations différentielles, en résolvant laquelle, vous pouvez prédire la dynamique de l'épidémie. Il y a beaucoup de modifications du modèle SEIRD, par exemple, SEIR est un modèle simplifié qui ne considère pas séparément les cas de récupération et de décès. Pour vous familiariser avec d'autres modèles, je peux vous recommander un bon article sur ce sujet.

Un peu de théorie

Pour la première fois, un modèle épidémique sous la forme d'un système de trois équations différentielles pour les variables

$inline$ est apparu dans les travaux de W. Kermak et A. McKendrick en 1927.
Ces équations différentielles ont la forme:

$\ 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}$

où, en plus des variables qui nous sont familières, les constantes suivantes apparaissent:

$inline$ - taille totale de la population,

$\ beta$ - taux de transmission des infections,

$\ gamma$ - vitesse de récupération.

La signification de l'équation de Kermak et McKendrick est la suivante: le nombre de sujets sensibles diminue proportionnellement à leur nombre multiplié par la proportion moyenne de personnes infectées dans la population

$inline$ , le nombre d'infectés croît au même rythme, ajusté du fait que certains d'entre eux

$\ gamma I$ en convalescence, et le nombre de convalescents augmente en raison d'une diminution du nombre de personnes infectées. Il est à noter que le modèle SIR contient une non-linéarité

$inline$ , à cause de laquelle la solution analytique du système d'équations devient généralement impossible, mais, heureusement, les méthodes de différenciation numérique peuvent facilement faire face à cette tâche.

Ajout d'une autre variable ici

$inline$ (le nombre de personnes atteintes de la maladie pendant la période d'incubation), nous obtenons le modèle 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}$

où une autre constante apparaît

$\ kappa$ - le taux de transition de la maladie du stade d'incubation à l'ouvert. Figure tirée de l' article .

Il est immédiatement clair que le modèle SEIR n'est pas très approprié pour décrire COVID-19, ne serait-ce que parce que dans ce modèle il y a des porteurs cachés d'infection

$inline$ pas contagieux. Cette carence peut être corrigée en introduisant, suivant Pengpeng et al. paramètre facultatif

$\ theta$ , caractérisant le degré d'infectiosité des porteurs d'infection latents par rapport aux malades. Le modèle SEIR modifié, que nous allons essayer d'appliquer à l'épidémie actuelle, ressemblera à:

$\ 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}$

À première vue, le modèle résultant promet d'être complètement crédible.

Expérience numérique avec le modèle SEIR

Pour la modélisation, nous essaierons de prendre les paramètres suivants, en nous concentrant sur les données ouvertes . En supposant que la maladie dure en moyenne 14 jours (au moins combien de temps dure la forme bénigne, qui représente jusqu'à 80% des cas), nous trouvons la valeur

$\ gamma = 1/14 = 0,0714$ . Va accepter

$\ beta = 3/14 = 0,2143$ . Ordre de grandeur

$\ theta = 0,6$ emprunté à Pengpeng et al . Compte tenu de la période d'incubation moyenne de 3 jours, prendre

$\ kappa = 1/3 = 0,33$ . Nous acceptons la population de la Russie égale

$N = 144,5 \ cdot10 ^ 6$ la personne.

Comme conditions initiales, nous utilisons les données pour la Russie le 2 avril, lorsque les mesures introduites fin mars pour limiter la propagation de l'infection devaient avoir leur effet, à savoir:

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

Évaluation

$inline$ nous l'avons pris de manière relativement arbitraire, et cela n'a pas d'importance, car, comme vous le savez, quelque chose va mal se passer.

Suite à la modélisation par la méthode d'Euler avec un pas de 1 jour du 2 avril au 24 avril inclus, nous obtenons les graphiques ci-dessous: à gauche dans l'échelle linéaire, à droite dans le logarithmique.

Des marqueurs ronds marquent des données réelles sur le nombre total de cas en Russie, des marqueurs carrés sur le nombre de patients. À première vue, les résultats semblent bons, sauf un: avec les paramètres du modèle, nous n'avons clairement pas deviné. Et ici, les méthodes d'optimisation viennent à notre aide.

L'optimiser

Les méthodes d'optimisation, si le lecteur ne les connaît pas, sont des algorithmes qui permettent de trouver le minimum d'une fonction objective. Dans notre cas, nous avons devant nous le problème de la régression non linéaire: comment choisir le vecteur de paramètres de l'équation différentielle

$\mathbf{x} = (\beta, \gamma, \kappa, E_0)^\top$ de sorte que l'ensemble des points de solution de l'équation différentielle

$inline$ était aussi proche que possible de l'ensemble des points d'observation

$\mathcal{F}$ .

Nous utilisons l'écart type comme mesure de l'erreur de modèle. La fonction objectif prendra la forme

$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},$

Où

$inline$ - le nombre de points

$inline$ - le nombre total de cas d'infection donnés par le modèle,

$\mathcal{F}$ - le nombre total réel de cas,

$inline$ - le nombre de patients actuellement, ce qui donne le modèle,

$\mathcal{G}$ - le total réel des cas actifs en cours.

À l'aide de la boîte à outils d'optimisation dans MATLAB, nous ajusterons les paramètres du modèle aux données d'observation. En conséquence, nous obtenons la solution indiquée dans la figure ci-dessous.

À première vue, tout va bien. L'écart s'est avéré égal

$f(\mathbf{x}) = 131,98$ , et «à l'œil de la mer convexe» l'ajustement de la solution s'est bien passé. Regardons les paramètres reçus:

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

La valeur de près de 8 millions de patients latents, avec environ 60 000 cas enregistrés
le 24 avril, est quelque chose de douteux. Nous avons également constaté que le temps de transition moyen vers la phase active de la maladie est

$1/\kappa = 2079$ journées.

Pourquoi est-ce arrivé? Tout deviendra clair si nous analysons la forme de la courbe sur une longue échelle de temps. Pour ce faire, prenez notre modèle SEIR avec des paramètres "plausibles" et simulez-le sur une longue période de temps (dans cette expérience j'ai pris un nouveau sens

$\beta=0,186$ ):

La courbe correspondant au nombre total de cas a une forme en S caractéristique sur une échelle linéaire. Le programme d'optimisation a essayé de donner à cette forme une courbe. De plus, la prévision elle-même avec des paramètres «plausibles» est horrible - selon elle, près de 90% de la population du pays sera malade d'ici septembre - elle est évidemment irréaliste si vous regardez les résultats pour d'autres pays (la même image qu'au début de l'article, uniquement sur une échelle linéaire ):

Ici, je compare les trois pays européens dans le top 5 du nombre de cas, et la Russie. On peut voir que nous accusons un retard d'environ un mois sur le rythme de l'épidémie et que pendant un mois, la croissance du nombre total de cas dans les trois pays est presque linéaire (et même plus lente que linéaire), contrairement aux résultats obtenus dans le modèle SEIR. Cela soulève trois questions:

?
SEIR, ?
, , ?

Je vais commencer par répondre à la troisième question. Lorsque nous prédisons quelque chose, nous sommes confrontés à une tâche plutôt désagréable: les données sur lesquelles nous construisons le modèle ne sont pas idéales - elles contiennent des erreurs, du bruit, et le modèle construit sur leur base contient également des erreurs. Lorsque nous continuons la série chronologique avec nos points de modèle, l'erreur s'accumule - et assez rapidement, si nous prédisons une fonction qui augmente dans le temps. Et c'est le cas dans notre cas. De plus, le modèle est un modèle qui reflète la situation réelle est très limitée. Le développement soudain de l'épidémie dans une nouvelle grande ville, l'utilisation d'une méthode de traitement plus efficace, un changement dans la méthode de collecte des informations - tout cela peut faire tellement d'erreurs dans les données réelles qu'une prévision à long terme sera complètement loin de la réalité.

Béquilles et vélos: nous modifions le modèle SEIR

Essayons de répondre à la question de savoir pourquoi la croissance de l'épidémie ralentit pour devenir linéaire. Avec le nombre de personnes infectées que nous avons actuellement, les économies d'échelle associées à la vitesse de communication limitée entre les personnes jouent un rôle important.

Plus précisément, rappelons-nous: le nombre de cas dans le modèle SEIR est directement proportionnel au nombre moyen de cas dans la population

$inline$ . Cette règle fonctionne bien dans les petites populations, où tout le monde peut communiquer avec tout le monde, et les patients sont répartis uniformément. En réalité, en particulier à l'échelle de dizaines et de centaines de milliers de personnes, si vous prenez deux personnes malades au hasard, il s'avère qu'elles non seulement ne communiquaient jamais entre elles et ne se voyaient pas, elles ne montaient même pas dans la même voiture de métro. Et en effet, ils vivent dans des villes différentes. Tout ce qui les unit, c'est une chaîne de liens sociaux, ce qui a conduit au fait que le virus leur a été transmis.

À titre d'exemple, j'ai construit un modèle épidémique sous la forme d'un automate cellulaire, où chaque cellule interagit avec seulement 4 cellules voisines. Cela équivaut au fait que chaque individu dans la population a 4 contacts sociaux - c'est un très petit nombre pour la population humaine, mais plus vite l'effet de restriction des liens sociaux se manifestera. A chaque itération avec une probabilité de 0,1, chacun des 4 voisins de la cellule infectée peut être infecté. La maladie dure en moyenne 14 jours. Les résultats de simulation pour un pool de cellules 200x200 sont présentés dans la figure ci-dessous, où

$inline$ Est le numéro d'itération.

La couleur bleue indique sensible, jaune - malade, vert - récupéré. La chose la plus intéressante est à quoi ressemblent les graphiques du nombre de cas. Et ils ressemblent à quelque chose de prévu: après une courte phase de croissance sous-exponentielle, tout comme dans le modèle SEIR, il y a une phase prolongée de croissance linéaire - comme dans la réalité.

Mon objectif n'était pas d'obtenir une image qui ressemble quantitativement à la réalité. Si vous voulez plus de crédibilité, je peux recommander le projet de Sergey Potekhin, qui a été récemment publié sur Habré . Pour les lecteurs méticuleux, la preuve plus stricte de la croissance linéaire est donnée ci-dessous.

Preuve du théorème d'épidémie de croissance linéaire dans une grande population

$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- .

En ce qui concerne les statistiques mondiales sur l'épidémie, nous verrons la même chose: depuis un mois maintenant, la croissance est linéaire, malgré le fait que les modèles classiques nous promettaient une poursuite de l'augmentation exponentielle du nombre de cas.

Tâches pour les curieux

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

Passons maintenant à la modification du modèle SEIR. La chose la plus simple que nous puissions faire est de multiplier la composante non linéaire par une fonction, en fonction du nombre de patients. Avec un petit nombre de cas, cette fonction doit être proche de 1, avec un grand nombre, elle doit tendre asymptotiquement vers zéro. Le candidat approprié le plus simple est

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

Sélection des paramètres

$\alpha$ et

$inline$ peut compenser la croissance exponentielle dans le modèle d'origine.

Ajouter au modèle, pour plus d'informations, et au composant

$inline$ - le nombre de décès. Obtenez une version modifiée du modèle 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 { aligner}$

Les résultats de la simulation sont présentés dans la figure ci-dessous.

L'erreur standard par rapport au modèle d'origine est presque inchangée. Les valeurs des paramètres sont déjà réalistes. Pour plus de commodité, j'ai désigné le nombre initial de personnes infectées dans la phase active comme

$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}$

Le modèle interpole également très bien les dérivés - les valeurs de l'augmentation quotidienne du nombre de cas et du nombre de décès.

Essayons de faire une prévision. Nous prenons un horizon de prévision de 2 mois et continuons à modéliser les solutions avec les paramètres trouvés par le programme d'optimisation.

À première vue, pas mal, mais vous ne pouvez pas souhaiter une telle prévision à votre patrie bien-aimée: le nombre de nouveaux cas continuera de diminuer, mais le nombre total continuera d'augmenter. Dans ce cas, une épidémie ne peut être arrêtée qu'à l'aide d'un vaccin ou en attendant que la quasi-totalité de la population soit malade. Le nombre de nouveaux décès est fixé à environ 200 par jour. C'est une illustration claire de ce qui se passera si nous ne renforçons pas les mesures de lutte contre l'épidémie. C'est ça qui nous attend? Et pour cet avenir pas très brillant, beaucoup d'entre nous sont assidûment à la maison, après avoir acheté du sarrasin et du papier toilette?

Ci-dessous, je considérerai deux scénarios, et en regardant la distance brumeuse des prochains mois du 28 avril 2020, je ne peux pas dire avec certitude lequel développera davantage les événements. Maintenant, au moment de briser la courbe des nouveaux cas, nous sommes à un point où il est doublement problématique de prévoir quelque chose.

Scénario américain

L'hégémonie mondiale était dans une position extrêmement peu enviable. Retardé dans la prise de décisions clés qui ralentiraient au début la croissance de l'épidémie, il est toujours incapable de faire face à l'augmentation naturelle des nouveaux cas.

Le modèle SEIRD modifié, formé aux 33 premiers points, à partir du 2 mars, plus ou moins prédit de façon réaliste l'évolution de l'épidémie en avril.

Comme vous pouvez le voir, la croissance en avril est presque parfaitement linéaire. Le modèle surestime légèrement la mortalité en avril, mais le tableau d'ensemble est correct.

Cette image montre l'augmentation quotidienne des nouveaux cas et des décès aux États-Unis. Il est très similaire à ce que le modèle prédit pour la Russie.

Scénario Allemagne

Des Allemands disciplinés ont réussi à inverser la courbe en leur faveur et sa croissance est plus lente que linéaire. De plus, pour rendre le modèle pertinent, j'ai dû ajouter manuellement une augmentation du taux de récupération le 6 avril

$\ gamma$ 1,7 fois, sinon une telle baisse brutale du nombre de cas en termes de modèle SEIRD ne peut pas être expliquée.

Le modèle a été formé aux 27 premiers points, à partir du 10 mars. J'ai également changé la fonction non linéaire. Pour l'Allemagne, l'exposant dépendant du temps est meilleur:

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

Ce type de fonction indique une augmentation cumulative du nombre de liens sociaux interrompus et, par conséquent, la propagation de l'infection. Ici, vous avez une illustration claire des avantages de l'auto-isolement.

Ce qui précède montre l'augmentation quotidienne des nouveaux cas et des décès. Comme dans le cas des États-Unis, les données réelles contiennent des fluctuations prononcées avec une période de 7 jours. Cela signifie que le week-end, le nombre de contacts augmente et, par conséquent, le nombre de contacts infectés augmente. [UPD: au contraire, elle diminue, comme l' indique l' auto-isolement Yandex peut en témoigner indirectement . ]

Conclusion

Faire des prédictions - à court et à long terme - n'est pas seulement un hommage à la curiosité. En cas d'épidémie, vous devez savoir combien de lits doivent être préparés, combien de ventilateurs à produire, pendant combien de mois pour stocker les équipements de protection individuelle des médecins. Les fonctionnaires doivent comprendre si les mesures prises sont suffisantes ou si de nouvelles interdictions et restrictions doivent être introduites. Idéalement, le modèle devrait refléter si bien la réalité qu'il serait possible de voir la force d'action de chaque mesure nouvellement adoptée, puis il serait possible de renforcer les mesures utiles et d'annuler les décisions qui se sont avérées inutiles.

Bien qu'avec quelques réserves, tous ceux qui n'ont pas encore acquis l'immunité de COVID-19 peuvent vraiment être décrits avec la lettre

$inline$ , toute la variété des malades - la lettre

$inline$ , etc. De plus, le modèle SEIRD peut même expliquer quelque chose. Mais elle peut prédire quelque chose dans un avenir lointain très approximativement.

Je n'ai délibérément cité dans l'article qu'un scénario négatif lorsque nous répétons le sort des États-Unis en termes de dynamique de l'épidémie. Si ce scénario s'avère vrai, nous aurons fin juin plus de 300 000 cas enregistrés de la maladie et plus de 10 000 décès. Bien qu'il y ait des conditions préalables au fait que ce scénario ne se matérialise pas, je vous conseille de vous y rapporter selon le principe: "espérer le meilleur, se préparer au pire". Comme le dit le proverbe, si les meilleurs esprits de la NASA ont pris la lutte contre l'épidémie aux États-Unis , alors c'est vraiment une mauvaise chose.

Jusqu'à présent, il ne nous reste plus qu'à visiter moins de lieux publics, à utiliseravec des respirateurs appropriés , lavez-vous les mains, n'oubliez pas d'essuyer le smartphone avec de l'alcool après l'avoir sorti dans la rue et suivez d'autres recommandations simples.

Mais encore, est-il possible de faire des prédictions plus précises? Oui bien sûr. Mais à ce sujet une autre fois.

Si vous souhaitez jouer avec le code source et suggérer un scénario de développement vous-même, voici le lien vers le github .

Pandémie COVID-19 vue par un mathématicien, ou pourquoi le modèle SEIRD classique ne fonctionne pas