David Madore's WebLog: Sur une variante à temps de rétablissement constant du modèle épidémiologique SIR

Il y a quatre éternités semaines, quand nous n'étions pas encore maintenus prisonniers chez nous, j'ai parlé ici du modèle épidémiologique SIR, le plus basique qui soit. Je rappelle brièvement les principes qui le définissent :

• l'immunité acquise est permanente, les individus sont successivement S (susceptibles, c'est-à-dire jamais infectés donc susceptibles de l'être), I (infectés et infectieux) et R (rétablis, c'est-à-dire guéris ou morts) (il existe toutes sortes de variantes, par exemple le modèle SEIR ajoutant un état E (exposé) pour les individus infectés mais non encore infectieux) ;
• la population est homogène (fongible) avec mélange parfait dans les contacts (j'ai parlé ici de l'effet de modifier cette hypothèse) ;
• la contamination et le rétablissement se font selon une cinétique d'ordre 1, c'est-à-dire que la contamination se fait proportionnellement aux proportions d'infectés et de susceptibles (avec une constante cinétique β), et que le rétablissement se fait proportionnellement à la proportion d'infectés (avec une constante cinétique γ).

Rappelons brièvement ce que j'ai exposé la dernière fois. Les équations de ce modèle SIR basique, que j'appellerai (*) pour m'y référer plus tard, sont les suivantes (il s'agit d'un système d'équations différentielles ordinaires non-linéaire, du premier ordre et autonomes) :

• s′ = −β·i·s
• i′ = β·i·s − γ·i
• r′ = γ·i
• (s+i+r=1)

où s,i,r≥0 sont les proportions de susceptibles, d'infectieux et de rétablis dans la population ; les solutions de ces équations ne semblent pas pouvoir s'exprimer en forme close, mais on peut exprimer s en fonction de r (à savoir s = exp(−κ·r) dans les conditions exposées ci-dessous).

Je rappelle les principales conclusions que j'avais exposées dans mon entrée sur ce modèle (*), en supposant qu'on parte d'une population presque entièrement susceptible avec une proportion infinitésimale d'infectés (plus exactement, on s'intéresse aux solutions pour lesquelles s→1 quand t→−∞) ; on notera κ := β/γ le nombre de reproduction, que je suppose >1 :

• tant que s reste très proche de 1 (si on veut, t→−∞), les proportions i et r croissent comme des exponentielles de pente logarithmique β−γ = β·((κ−1)/κ), avec un rapport 1/(κ−1) entre les deux, autrement dit comme i = c·exp((β−γ)·t) = c·exp(β·((κ−1)/κ)·t) et r = c·(γ/(β−γ))·exp((β−γ)·t) = c·(1/(κ−1))·exp(β·((κ−1)/κ)·t) (ergotage : dans l'entrée sur le sujet, j'avais mis un −1 aux exponentielles pour r, parce que je voulais partir de r=0, mais je me rends compte maintenant qu'il est plus logique de partir d'une solution où i/r tend vers une constante en −∞, cette constante étant κ−1) ;
• au moment du pic épidémique (maximum de la proportion i d'infectés), on a s = 1/κ et i = (κ−log(κ)−1)/κ et r = log(κ)/κ ;
• quand t→+∞, la proportion i tend vers 0 (bien sûr) et s tend vers Γ := −W(−κ·exp(−κ))/κ (en notant W la fonction de Lambert) l'unique solution strictement comprise entre 0 et 1 de l'équation Γ = exp(−κ·(1−Γ)) (qui vaut 1 − 2·(κ−1) + O((κ−1)²) pour κ proche de 1, et exp(−κ) + O(κ·exp(−2κ)) pour κ grand), tandis qu'évidemment r, lui, tend vers 1−Γ.

Je veux ici explorer la modification d'une hypothèse de ce modèle (*), celle qui concerne le rétablissement. Quand j'écris ci-dessus que le rétablissement se fait proportionnellement à la proportion d'infectés (avec une constante cinétique γ), au niveau individuel, cela signifie la chose suivante :

Pendant chaque intervalle de temps de longueur (durée) dt très courte, la probabilité qu'un individu infecté (I) se rétablisse (I→R) vaut γ·dt et ce, indépendamment d'un individu à l'autre et d'un instant à l'autre.

Autrement dit, le temps de rétablissement d'un individu infecté donné suit une distribution de probabilité exponentielle d'espérance 1/γ.

Autant l'hypothèse analogue sur la cinétique de la contamination est relativement plausible (si on admet le principe éminemment discutable d'une population homogène et du mélange parfait !), autant l'hypothèse sur le temps de rétablissement est médicalement insensé : on est en train de dire que si vous êtes malade, votre probabilité de guérir (ou d'ailleurs, de mourir) ne dépend pas de l'avancement de votre maladie mais est la même pendant la première heure que pendant la 1729e (si tant est que vous soyez encore malade à ce stade-là). Une maladie ne se comporte pas comme ça !

Cherchons donc à remplacer cette hypothèse par une autre, tout aussi simpliste, mais néanmoins un peu plus proche de la réalité médicale, celle du rétablissement en temps constant.

Un individu infecté (I) se rétablit toujours au bout du même temps T après son moment d'infection.

Autrement dit, le temps de rétablissement d'un individu infecté donné suit une distribution de Dirac concentrée en T (qui est, du coup, son espérance).

Je vais appeler (†) (que je dois encore expliciter) le modèle SIR construit sur l'hypothèse que je viens de dire au lieu de celle qui précède. Je répète que les deux sont extrêmement simplificatrices et invraisemblables dans la réalité, mais l'hypothèse de temps constant est nettement moins fausse ou irréaliste que celle de temps à distribution exponentielle (surtout que ce qui compte vraiment est le temps pendant lequel la personne est infectieuse en pratique) : on pourrait chercher des hypothèses plus fines et moins simplistes (par exemple la somme d'une constante et d'une variable à distribution exponentielle ; ou une distribution de Weibull ou que sais-je), mais mon propos est juste de comparer (*) et (†) pour voir ce qui change et avoir au moins une idée du sens dans lequel la modélisation est affectée par le changement d'une distribution exponentielle vers une constante.

Alors pour commencer, quelles sont les nouvelles équations ? Le terme d'infection β·i·s ne va pas être modifié, mais le terme de guérison, lui, va l'être : au lieu d'avoir une proportion γ·i des infectés qui guérissent par unité de temps, on va faire guérir tous ceux qui étaient nouvellement infectés il y a T unités de temps. Comment exprimer ceci ? Eh bien, si je note f_T la fonction f translatée de T, c'est-à-dire f_T(t) := f(t−T), le terme représentant les guérisons à l'instant t est donné par les celui représentant les infection à l'instant t−T, c'est-à-dire β·(i·s)_T = β·i_T·s_T. Autrement dit, mon système devient :

• s′ = −β·i·s
• i′ = β·i·s − β·i_T·s_T
• r′ = β·i_T·s_T
• (s+i+r=1)

Ce n'est plus un système d'équations différentielles mais d'équations différentielles à retard, c'est-à-dire que la dérivée imposée des fonctions à un instant donné ne dépend pas seulement de leurs valeurs à cet instant mais aussi de celles à un ou plusieurs instants dans le passé. Je ne sais essentiellement rien sur ces équations ; en particulier, j'ignore ce qu'il y a comme résultat d'existence et d'unicité de solutions : on peut toujours essayer d'imposer les valeurs de s et i sur un intervalle de largeur T, disons [−T;0] et utiliser les équations pour en déduire ces valeurs en [0;T], puis en [T;2T] et ainsi de suite, mais ce procédé n'a aucun intérêt si on ne trouve pas moyen de faire en sorte que les valeurs se recollent correctement aux extrémités communes de ces intervalles, c'est-à-dire les multiples de T ; or mon but est de trouver des solutions qui soient régulières partout (au moins C^∞, si possible analytiques), pas quelque chose de recollé sur les multiples de T. Le même problème fait qu'il n'est pas évident de simuler numériquement une telle équation, faute de savoir comment la démarrer : la notion même de valeur initiale n'a pas de sens clair comme elle en a pour les équations différentielles.

J'étais très pessimiste quant à la possibilité de résoudre explicitement ce système (†) qui a l'air plutôt plus compliqué que (*) (lequel n'admet apparemment déjà pas de solution explicite). Un peu miraculeusement, et à ma grande surprise, j'y suis quand même arrivé (après un temps assez invraisemblable passé à explorer toutes sortes de techniques de calcul, et en étant à chaque fois à deux doigts d'abandonner).

Pour ceux qui se demandent comment j'ai trouvé les solutions qui suivent, j'ai commencé par me pencher sur le cas où s vaut constamment 1 (c'est-à-dire que i et r sont infinitésimaux), ce qui linéarise le système, et il est alors raisonnable d'en chercher des solutions à croissance exponentielle, ce qui est assez facile (on obtient ce que je décris ci-dessous comme comportement à i et r petits) ; j'ai ensuite eu l'idée de chercher en général des séries formelles en cette exponentielle, j'ai calculé les premiers coefficients de cette série pour savoir si la série avait au moins des chances de converger, et j'ai constaté que miraculeusement ces coefficients avaient une forme extrêmement simple que je pouvais sommer explicitement, d'où les expressions rationnelles-exponentielles qui suivent.

(Ce qui suit est essentiellement un développement de cette question+réponse sur MathOverflow et ce petit fil Twitter. La réponse MathOverflow est peut-être plus commode pour lire les formules, d'ailleurs, que mon HTML basique dans ce qui suit.)

Voici la solution explicite que j'ai trouvée : pour l'exprimer, il faut que j'introduise les notations

• κ := β·T le nombre de reproduction, que je suppose toujours >1 ;
• Γ := −W(−κ·exp(−κ))/κ, comme évoqué plus haut, l'unique solution strictement comprise entre 0 et 1 de l'équation Γ = exp(−κ·(1−Γ)), qui sera la limite de s quand t→+∞ (aussi bien dans le modèle (†) que (*)) ;
• X := exp(β·(1−Γ)·t), fonction exponentielle du temps dans laquelle les solutions seront exprimées (changement de variable ; remarquer que soustraire T à t revient à multiplier X par Γ, ce qui est bien l'intérêt du changement de variable en question) ;

par ailleurs, c sera un paramètre réel strictement positif arbitraire (qui sert simplement à translater les solutions dans le temps).

Les solutions que j'ai trouvées sont alors données par :

• s = ((1−Γ)² + Γ·c·X) / ((1−Γ)² + c·X)
• i = ((1−Γ)⁴·c·X) / (((1−Γ)² + c·X) · ((1−Γ)² + Γ·c·X))
• r = (Γ·(1−Γ)·c·X) / ((1−Γ)² + Γ·c·X)

(Comme pour (*), je m'intéresse uniquement aux solutions pour lesquelles s→1, c'est-à-dire que i et r tendent vers 0, quand t→−∞. Si on veut des solutions plus générales, on peut multiplier i et r par une même constante, quitte à multiplier β par l'inverse de cette constante.)

On a exprimé s,i,r comme des fonctions rationnelles d'une fonction exponentielle X du temps t. Les dénominateurs sont strictement positifs puisque (1−Γ)² + c·X et (1−Γ)² + Γ·c·X sont des sommes de quantités strictement positives ; les numérateurs sont eux aussi strictement positifs ; et comme s+i+r=1 comme on peut vérifier avec un petit calcul, on a bien affaire à des proportions comme souhaitées. Pour vérifier qu'elles satisfont aux équations demandées, il s'agit simplement de dériver chacune par rapport à t, ce qui revient à dériver par rapport à X et multiplier par β·(1−Γ)·X, et vérifier qu'on trouve bien le terme de droite souhaité (en utilisant le fait que soustraire T à t revient à multiplier X par Γ).

Sous cette forme, cette solution n'est évidemment pas très parlante ! Mais pour y voir plus clair, on peut mener la même étude que pour (*) en regardant le comportement des solutions trouvées à (†) pour t→−∞, au niveau du pic épidémique, et enfin pour t→+∞ :

• tant que s reste très proche de 1 (si on veut, t→−∞), les proportions i et r croissent comme des exponentielles de pente logarithmique β·(1−Γ), avec un rapport Γ/(1−Γ) entre les deux, autrement dit comme i = c·exp(β·(1−Γ)·t) et r = c·(Γ/(1−Γ))·exp(β·(1−Γ)·t) ; qualitativement, la phase de croissance exponentielle est plus rapide (à infectiosité β et nombre de reproduction κ donnés !) dans le modèle (†) que dans le modèle (*), et le rapport rétablis/infectés pendant cette phase est plus faible ;
• au moment du pic épidémique (maximum de la proportion i d'infectés), qui se produit pour c·X = (1−Γ)²/√Γ (mais peu importe puisque c'est un temps complètement arbitraire), on a s = √Γ et i = (1−√Γ)² et r = (√Γ)·(1−√Γ) ; qualitativement, le pic épidémique est plus haut et plus resserré (à infectiosité β et nombre de reproduction κ donnés !) dans le modèle (†) que dans le modèle (*) ;
• quand t→+∞, la proportion i tend vers 0 (bien sûr) et s tend vers Γ, tandis qu'évidemment r, lui, tend vers 1−Γ : ce sont exactement les mêmes valeurs que dans le modèle (*) (à nombre de reproduction κ donné).

Voici des courbes illustrant deux épidémies modélisées par les modèles (*) et (†) avec les mêmes paramètres β et κ (en l'occurrence β=3, un temps de rétablissement espéré de 1, donc nombre de reproduction κ=3, et ceci conduit à Γ≈0.034, c'est-à-dire un taux d'attaque final de 96.6%, je ne réexplique pas pourquoi l'hypothèse de mélange parfait, commune à mes deux modèles, conduit à des taux d'attaques irréalistement élevés, ce n'est pas le propos ici). Les deux courbes du haut utilisent le modèle SIR classique (*) avec rétablissement selon un processus exponentiel (ici γ=1), tandis que les deux du bas utilisent le modèle SIR modifié (†) que je viens d'exposer avec rétablissement en temps constant (ici T=1). Les deux courbes de gauche sont en échelle linéaire, les deux de droite en échelle logarithmique :

(Le code Sage que j'ai utilisé pour les générer, et qui ne présente aucun intérêt particulier, est ici.)

Puisque l'origine des temps est arbitraire, je l'ai placée au niveau du pic épidémique, ce qui aide à la comparaison. Rappelons par ailleurs que le seul vrai paramètre est le nombre de reproduction κ (le paramètre β peut s'absorber en faisant un changement d'échelle du temps).

La mauvaise nouvelle, donc, si l'objectif est de flatten the curve, est que l'hypothèse (†) (plus réaliste !) de guérison en temps constant fait tout le contraire, elle resserre et amplifie le pic (tout ceci étant dit à infectiosité, nombre de reproduction et taux d'attaque constants). On peut néanmoins y voir une sorte de contrepartie : l'infection prédite par le modèle (†) paraît plus terrifiante que ce qu'elle est réellement, au sens où sa croissance initiale laisse présager (si on l'interprète selon le modèle (*)) un nombre de reproduction et donc un taux d'attaque plus importants que ce qu'ils sont réellement.

Alors que dans le modèle (*) le pic épidémique a d'abord une croissance exponentielle de pente logarithmique β−γ = β·(1−1/κ) puis une décroissance exponentielle de pente logarithmique β·Γ−γ = β·(Γ−1/κ) (sur mon graphique, 2.000 puis −0.821), dans le modèle (†) le pic est parfaitement symétrique, avec une croissance exponentielle de pente logarithmique β·(1−Γ) suivie d'une décroissance exponentielle de pente opposée (sur mon graphique, 2.821 puis −2.821). Une morale de l'histoire est qu'il est difficile de lire le nombre de reproduction à partir de la pente logarithmique de la partie exponentielle de l'épidémie ! Il faut une information assez fine non seulement sur le taux de rétablissement moyen mais aussi sur sa distribution.

Le rapport entre le nombre d'infectieux et le nombre de rétablis, pendant la phase de croissance exponentielle, est beaucoup plus important dans le modèle (†) (à savoir (1−Γ)/Γ) que dans le modèle (*) (à savoir κ−1) : sur mon graphique (κ=3), ce rapport est de 15.8 dans le modèle (†) et 2 dans le modèle (*). Mais on peut aussi exprimer ce fait différemment en considérant le retard des rétablis sur les infectés+rétablis, c'est-à-dire combien de temps il faut attendre pour avoir un nombre de rétablis égal au nombre d'infectés+rétablis à un moment donné (toujours pendant la phase de croissance exponentielle) : ce retard est de κ·log(κ)/(β·(κ−1)) dans le cas (*) et de T = κ/β dans le cas (†) ; il serait intéressant de comparer ces deux fonctions de κ, mais j'avoue manquer un peu de patience, là ; sur mes graphiques, ce retard est de 0.549 pour le modèle (*) et évidemment 1.000 pour le modèle (†).

Reste à dire un mot sur le nombre de reproduction κ. D'abord pour justifier la manière dont je l'ai défini (β/γ dans le cas (*) et β·T dans le cas (†)) : la définition classique est le nombre de personnes moyen qu'une personne donnée infecte avant d'être rétablie, dans une population entièrement susceptible. Dans mes deux modèles, une personne infectée au milieu d'une population entièrement susceptible infecte β·dt autres personnes par unité de temps dt (c'est la définition de β) ; ceci étant simplement proportionnel au temps qu'elle restera infectieuse, l'espérance de ce nombre vaut simplement β fois l'espérance du temps qu'elle restera infectieuse, c'est-à-dire 1/γ dans le cas (*) et T dans le cas (†). Maintenant, pouvait-on prévoir a priori que le taux d'attaque final en fonction de κ serait le même, 1−Γ = 1 + W(−κ·exp(−κ))/κ, dans les modèles (*) et (†) ? Et si oui, quelle hypothèse sur la distribution de probabilité du processus de rétablissement assure ce fait ? Là, j'avoue ne pas du tout savoir, je ne connais sans doute pas assez de probas pour m'exprimer clairement. (Je suis tenté de chercher à exprimer le processus de contamination de manière indépendante du temps, en fonction d'une sorte de génération comme je l'avais fait dans mes modèles stochastiques, pour espérer que le comportement soit identique dans les deux cas, mais je n'arrive pas à exprimer ni même définir les choses correctement.)