David Madore's WebLog: Un tout petit peu d'épidémiologie mathématique

[Index of all entries / Index de toutes les entréesLatest entries / Dernières entréesXML (RSS 1.0) • Recent comments / Commentaires récents]

↓Entry #2639 [older| permalink|newer] / ↓Entrée #2639 [précédente| permalien|suivante] ↓

(samedi)

Un tout petit peu d'épidémiologie mathématique

Dans l'entrée précédente, je soulevais entre autres la question de comment calculer (et de comment appeler !) le nombre, que j'y appelais r, de personnes qui sont finalement infectés par une épidémie (quelle que soit l'issue de cette infection) puisque c'est un des facteurs du produit f·r qui donnera le taux de mortalité due à l'infection (l'autre étant la proportion f des cas qui conduisent à un décès) ou de tout autre calcul analogue (comme g·r pour le nombre de cas graves où g est la proportion correspondante). Dans plusieurs mises à jour ultérieures de cette entrée, j'ai signalé que j'ai fini par apprendre que r s'appelle le taux d'attaque et un raisonnement simpliste pour l'estimer, que je reproduis ici parce que je vais vouloir le comparer à une estimation donnée par un modèle différent :

[Essentiellement recopié de ce fil Twitter :] Une amie m'a expliqué le rapport que je cherchais à comprendre entre le taux de reproduction de base R₀ (= nombre de personnes que chaque personne infectée infecte à son tour) et le taux d'attaque final r (= proportion de la population qui sera infectée à terme pendant l'épidémie) : dans le modèle le plus simpliste, c'est r = 1 − 1/R₀ ; en effet, tant que le taux de reproduction est >1, l'épidémie croît exponentiellement ; mais si une proportion r a déjà été infectée, le taux effectif de reproduction est ramené à R₀·(1−r) parce que, en supposant que les personnes déjà infectées sont immunisées et sont également réparties dans la population (j'ai bien dit, modèle simpliste !), seule une proportion 1−r est encore susceptible d'être contaminée ; donc l'épidémie cesse de progresser lorsque R₀·(1−r) redescend à 1, c'est-à-dire r = 1 − 1/R₀. C'est probablement la raison pour laquelle certains ont prédit r ~ 70% en l'absence de contre-mesures efficaces pour réduire R₀ qui a été initialement mesuré à R₀ ~ 3. Encore une fois, ceci est un modèle extrêmement simpliste.

Dans la suite, je vais noter plutôt κ que R₀ ce nombre de reproduction, parce que même si R₀ est la notation standard elle serait source de confusion dans le modèle SIR où la lettre R désigne les cas rétablis (guéris, recovered en anglais ; enfin, avec une drôle de définition de rétablis puisque dans le modèle qui va suivre on ne cherche pas à compter les décès et on les compte avec les guérisons). Par ailleurs, plutôt que le taux d'attaque final noté r ci-dessus (ce qui, par chance, colle bien, à la limite, avec l'usage de la lettre R que je viens d'évoquer), je vais m'intéresser plutôt à la proportion complémentaire s = 1−r, i.e., la proportion de la population qui échappe à l'épidémie, et dont le raisonnement simpliste que je viens de recopier prédit donc qu'il s'agit de 1/κ.

Maintenant, en suivant de près ce fil Twitter (ou ici sur Thread Reader), que je développe un peu un peu, je vais essayer d'expliquer la prédiction que fait un modèle basique en épidémiologie, le modèle SIR :

Le modèle SIR modélise une infection en traduisant l'évolution dans le temps de trois variables : s (susceptible) la proportion de la population qui n'a pas encore contracté l'infection (et qui est donc susceptible de l'attraper), i (infectée) la poportion de la population qui est actuellement infectée, et r (rétablie) la proportion de la population qui n'est plus infectée, que ce soit suite à une guérison ou un décès (cf. ci-dessus : on ne s'intéresse pas à la différence ici). On a s + i + r = 1 puisqu'il s'agit de trois parties exclusives et exhaustives : il y a donc seulement deux variables indépendantes. Le modèle fait toutes sortes d'hypothèses simplificatrices : notamment, que la population est constante (puisque les décès comptent parmi les guéris, ce n'est pas idiot), et surtout, que les personnes ayant contracté l'infection ne peuvent pas la contracter une seconde fois (soit parce qu'elles sont immunisées soit parce qu'elles sont décédées).

Il s'agit d'écrire une équation différentielle (non-linéaire, du premier ordre) portant sur ces variables. L'idée est d'écrire le type d'équations utilisées en cinétique chimique : imaginez qu'on aurait deux réactions chimiques, la réaction d'infection S + I → I + I (une personne infectée en infecte une autre) et la réaction de rétablissement, I → R (les personnes infectées se rétablissent toutes seules avec le temps, je rappelle une fois de plus que rétablir ici compte les décès, tout ce qui m'intéresse est que ces personnes ne puissent plus en contaminer d'autres). Ce qu'on fait en cinétique chimie (de façon ultra-simplifiée…) pour modéliser des réactions de type X + Y → Z est qu'on va écrire que l'occurrence d'une telle réaction, i.e., la variation de concentration due à cette réaction (qui va compter positivement dans la concentration de Z et négativement pour X et Y) est proportionnelle à une certaine constante cinétique (positive) fois le produit des concentrations de X et de Y à des puissances appelées l'ordre de la cinétique dans chacun de ces réactifs, typiquement 1. Dans le modèle épidémiologique SIR, les deux réactions d'infection et de rétablissement seront supposées d'ordre 1. On va appeler β et γ leurs constantes cinétiques respectives : les termes de vitesse de l'infection et du rétablissement seront donc β·i·s et γ·i respectivement. Autrement dit :

Si je note x′ la dérivée dx/dt par rapport au temps (t) de la variable x, les équations du modèle SIR seront :

  • s′ = −β·i·s
  • i′ = β·i·sγ·i
  • r′ = γ·i

(La somme de ces trois quantités fait évidemment zéro, comme il se doit puisqu'on doit conserver s+i+r=1 : comme en chimie, rien ne se crée, rien ne se perd, mais tout se transforme.) La première équation, donc, modélise le fait que la population non encore infectée décroît par la vitesse infection dans le temps β·i·s qui est proportionnelle à une constante β fois les proportions de personnes infectées i et susceptibles de l'être s : si l'on préfère, cela signifie qu'une personne susceptible a une probabilité β·i de devenir infectée par unité de temps (très petite) ; la troisième modélise le fait que les personnes infectées deviennent rétablies avec la vitesse γ·i : si l'on préfère, cela signifie qu'une personne infectée a une probabilité γ de devenir rétablie par unité de temps (très petite) ; et l'équation du milieu, donc, assure l'équilibre s+i+r=1.

La nouvelle (et énorme !) hypothèse simplificatice qu'on a faite en écrivant ces équations, c'est de supposer que le comportement « local » de l'épidémie et de la population ne change ni avec le temps ni avec le progrès de l'épidémie : la probabilité d'infection par rencontre S+I, ou de guérison, ne changent pas : ceci exclut, par exemple, le fait que la population changerait ses habitudes avec la progression de l'épidémie (prendrait des mesures prophylactique), que le système de santé soit débordé (ce qui jouerait possiblement sur le temps de guérison), que le pathogène mute pour devenir plus ou moins virulent, et toutes sortes d'autres scénarios sortant de notre modèle extrêmement basique.

Les constantes cinétiques β et γ ont pour grandeur l'inverse d'un temps : il s'agit essentiellement de l'inverse du temps espéré d'infection si toute la population est infectée et du temps espéré de guérison. Remarquons donc qu'en changeant l'échelle de temps on multiplie β et γ par la même constante : le seul paramètre sans dimension dans le modèle est le rapport κ := β/γ, qu'on interprète comme le nombre de personnes qu'une personne infectée infectera en moyenne dans une population entièrement susceptible avant d'être elle-même rétablie. Comme il s'agit du seul paramètre sans dimension, toute discussion doit se faire sur κ. C'est ce κ = β/γ qu'on appelle nombre de reproduction et qui est souvent noté R₀, mais que je préfère noter κ ici pour éviter la confusion avec la variable r.

Pour ce qui est des paramètres initiaux, on peut prendre ce qu'on veut, mais on prendra normalement r=0 (en fait, si r ne valait pas initialement 0, on pourrait diviser s et i par 1−r et multiplier β par cette quantité, pour se ramener à ce cas), et i=i₀ extrêmement petit (l'idée étant d'amorcer l'épidémie par un tout petit nombre de cas), ce qui permet de considérer que la valeur initiale s₀ de s est quasiment 1.

Les équations étant posées, que peut-on dire du comportement décrit par ce système ?

Le cas κ<1 étant peu intéressant (l'épidémie disparaît tout de suite, au moins si on est parti avec i₀ très petit), je vais prendre κ>1, c'est-à-dire β>γ.

Le départ de l'épidémie est très simple : tant que s est suffisamment proche de 1 pour être identifié à lui, on a i′ = (βγi qui vaut donc i₀·exp((βγt). Autrement dit, i croît exponentiellement avec une dérivée logarithmique βγ (ou si on veut, une constante de temps 1/(βγ), ou encore un temps de doublement ln(2)/(βγ)). Quant à r, il est donné en intégrant r′ = γ·i, soit r = i₀·(γ/(βγ))·(exp((βγt)−1) = i₀·(1/(κ−1))·(exp((βγt)−1) soit, quitte à négliger le −1, la même croissance exponentielle mais avec une valeur κ−1 fois plus petite.

On peut utiliser ce comportement pour tenter d'estimer les paramètres : βγ comme la pente logarithmique de la croissance exponentielle observée, et 1/γ comme la durée moyenne des symptômes (jusqu'à guérison ou décès) ; on devrait aussi pouvoir estimer κ−1 comme le rapport entre nombre de cas en cours et nombre de guéris+décès, mais dans la pratique ça marche très mal. Très grossièrement, dans l'épidémie de Covid-19 hors de Chine, on est en période de croissance exponentielle, avec une pente logarithmique de l'ordre de 0.2/j, et 1/γ de l'ordre de deux semaines, donc κ serait de l'ordre de 4 avec ces calculs. En fait, il est quand même plutôt de l'ordre de 3, parce que les guérisons ne se produisent vraiment pas selon un processus poisson mais plutôt en temps constant, ce qui justifie que la croissance exponentielle se fasse avec une pente logarithmique β plutôt que βγ initialement. Bref.

[Graphes des courbes s(t),i(t),r(t) pour β=3 et γ=1]Le modèle SIR prédit que la dérivée logarithmique de i vaut β·sγ, donc, si elle part de βγ quand toute la population est susceptible, elle devrait diminuer progressivement, s'annuler pour s = 1/κ (maximum de l'infection, donc) et ensuite devenir négative sans jamais pouvoir passer en-dessous de −γ. Le graphique ci-contre (tracé avec Sage en utilisant une méthode de Runge-Kutta d'ordre 4) montre l'allure des courbes, pour β=3 et γ=1 (donc κ=3), avec s tracé en vert, i en rouge et r en bleu.

Pour pouvoir calculer des quantités au pic de l'épidémie ou à sa limite quand t→+∞, il faut éliminer le temps dans les équations. Ceci se fait en divisant par r′ = γ·i les deux premières équations, soit

  • s′/r′ = −κ·s
  • i′/r′ = κ·s − 1

Or s′/r′ et i′/r′ s'interpètent comme ds/dr et di/dr respectivement (la fonction r est C¹ de dérivée strictement positive, on peut la prendre comme variable). Mais l'équation ds/dr = −κ·s se résout en s = s₀·exp(−κ·r), soit, puisqu'on part avec s₀ très proche de 1, s = exp(−κ·r), ce qui nous donne un lien reliant r et s.

Notamment, si on cherche la valeur du pic de l'épidémie (le moment où i est maximal, i.e., i′ s'annule), il est donné par s = 1/κ (d'après la valeur de i′) et donc r = ln(κ)/κ (en utilisant la relation qu'on vient de trouver) et ainsi i = 1 − (ln(κ)+1)/κ.

Maintenant, que dire des valeurs de s, i, r quand t→+∞ ? Manifestement, i→0 (démonstration : s et r ont des limites donc i en a une ; mais en regardant la troisième équation, r′ a une limite, qui ne peut pas être autre que 0 puisque r lui-même a une limite, donc i tend vers 0). Donc il y a deux quantités à déterminer, que je vais noter s et r (et, rapidement, s et r tout court). On a vu s = s₀·exp(−κ·r), donc notamment s = s₀·exp(−κ·r), mais par ailleurs r = 1 − s. En remplaçant l'une dans l'autre, on a l'équation s = s₀·exp(−κ·(1−s)). Avec s₀ très proche de 1, ceci devient s = exp(−κ·(1−s)).

J'écris maintenant simplement s et r pour s et r puisque je m'intéresse à ces limites en l'infini. Il s'agit donc de résoudre l'équation s = exp(−κ·(1−s)) pour trouver, en fonction du nombre de reproduction κ, la proportion s de la population qui reste non-infectée jusqu'au bout. (J'ai longtemps cru que le modèle SIR convergeait vers s=0, mais comme on va le voir, et comme c'est déjà clair sur cette équation, ce n'est pas le cas.) Encore une fois, si κ<1, ce n'est pas intéressant (l'épidémie s'arrête tout de suite, s vaut 1), donc prenons κ>1.

L'équation s = exp(−κ·(1−s)) n'a pas de solution avec des fonctions usuelles, mais elle en a une, à savoir s = −W(−κ·exp(−κ))/κ, en utilisant une fonction spéciale appelée fonction transcendante W de Lambert. Essentiellement, cette fonction est définie par le fait que w := W(z) vérifie w·exp(w) = z (il faut prendre la bonne branche complexe de cette fonction : disons qu'on parle ici, pour z≥−1/e réel, de l'unique solution réelle ≥−1 de l'équation w·exp(w) = z). La transformation de s = exp(−κ·(1−s)) en la bonne forme se fait simplement en notant que −κ·s·exp(−κ·s) = −κ·exp(−κ) d'où −κ·s est l'autre solution (à part −κ) de w·exp(w) = −κ·exp(−κ).

[Graphes des courbes −W(−κ·exp(−κ))/κ et 1/κ]Bref, le modèle SIR prédit que la proportion s de la population qui reste non-infectée jusqu'au bout vaut −W(−κ·exp(−κ))/κκ>1 est le nombre de reproduction. Comment cette fonction se comporte-t-elle ?

Et aussi, comment se compare-t-elle à l'estimation s = 1/κ citée au début de ce post ? En gros : c'est bien pire. Essentiellement, parce que le modèle SIR suit l'épidémie au-delà de son pic où κ·s=1 : elle a beau être en phase décroissante, elle continue à faire de nouvelles infections. Le graphique ci-contre montre les courbes donnant les valeurs s = −W(−κ·exp(−κ))/κ (en bleu) et s = 1/κ (en rouge) en fonction du nombre de reproduction κ.

Si le nombre de reproduction κ est juste un peu au-delà de 1, disons 1 + h avec h petit, alors s = −W(−κ·exp(−κ))/κ vaut 1 − 2·h + (8/3)·h² + O(h³). C'est en gros deux fois plus de cas que n'en donne l'estimation simpliste 1/κ soit 1 − h + h² + O(h³).

Mais surtout, si κ est grand, le développement W(z) = zz² + O(z³) donne s = −W(−κ·exp(−κ))/κ valant exp(−κ) + κ·exp(−2κ) + O(κ²·exp(−3κ)). Il y a bien une fraction s>0 qui reste non contaminée, mais elle décroît exponentiellement avec le nombre de reproduction κ.

Application numérique : pour κ≈3, le modèle SIR donne s = −W(−κ·exp(−κ))/κ ≈ 6% restant non-infectés, soit un taux d'attaque final de r ≈ 94% (et au moment du pic de l'épidémie, on a i ≈ 30% infectés contre s ≈ 33% encore susceptibles et r ≈ 37% rétablis). Je ne suis pas en train de dire que Covid-19 infectera 94% de la population mondiale (ou française), ni qu'au pic de l'épidémie il y aura 30% de malades ; outre que je suis mathématicien et pas épidémiologue, c'est ce que ferait une infection contre laquelle on ne prendrait aucune mesure et aucun changement de comportement, dans un modèle extrêmement simpliste, avec le paramètre κ≈3 qu'on a approximativement mesuré pour une phase non-contrôlée de l'épidémie. Néanmoins, cela donne une idée de ce que peut faire un nombre de reproduction de l'ordre de 3 (et je ne parle pas de la rougeole pour laquelle c'est plutôt de l'ordre de 15).

En un certain sens, en fait, la valeur simpliste 1/κ pour s est peut-être plus informative : c'est la valeur à laquelle on peut stabiliser l'épidémie : tant qu'on a une proportion r ≪ 1 − 1/κ de personnes immunisées (soit parce qu'elles ont déjà été contaminées ; mais cela pourrait être parce qu'on les a vaccinées si on a un vaccin), la situation est instable, l'épidémie peut (re)démarrer à tout moment, alors que si r ≫ 1 − 1/κ, la situation est stable, l'épidémie est contenue.

Il serait peut-être intéressant de voir ce que devient le modèle si on remplace le terme β·i·s par Β·i·(1−is partout, modélisant de façon simpliste le phénomène que les gens prennent peur et limitent leurs contacts au fur et à mesure que l'épidémie prend de l'ampleur (il y a peut-être une meilleure façon de modéliser ce phénomène, c'est juste ça qui me vient à l'esprit) ; on pourra encore moins calculer les limites faute d'avoir une forme sympathique pour la solution de ds/dr = −κ·(s+rsκ := Β/γ (Sage et Mathematica pensent que c'est exp(−½κ·r²) / (1 + (√(½πκ))·erf((√(κ/2))·r)), je vais dire que je les crois mais ce n'est pas très utilisable comme expression), mais on peut peut-être quand même en dire quelque chose.

Suite : voir cette entrée ultérieure pour une variation sur le modèle SIR où le rétablissement se fait en temps constant plutôt que selon un processus exponentiel.

↑Entry #2639 [older| permalink|newer] / ↑Entrée #2639 [précédente| permalien|suivante] ↑

[Index of all entries / Index de toutes les entréesLatest entries / Dernières entréesXML (RSS 1.0) • Recent comments / Commentaires récents]