David Madore's WebLog: Sur le modèle SIR avec susceptibilité hétérogène

Je continue dans ce billet de blog une série sur l'épidémiologie mathématique que j'avais commencée avec cette entrée sur le modèle SIR classique, celle-ci sur une variante de SIR où le rétablissement se fait en temps constant, accessoirement celle-ci sur la différence entre seuil d'immunité collective et taux d'attaque, et plus indirectement celle-ci sur des modèles d'hétérogénéité basés sur les graphes aléatoires ; je ne présuppose pas la lecture des billets en question, même si celle du premier a un intérêt, mais je vais en résumer rapidement le contenu.

Mon but aujourd'hui est d'expliquer un peu en détails, mathématiquement, comment on peut modifier le modèle SIR classique (dont le vais rappeler les grandes lignes dans un instant), lequel décrit l'évolution d'une épidémie dans laquelle tout le monde est également susceptible à l'infection, pour le cas d'une susceptibilité hétérogène, c'est-à-dire que certains individus sont plus ou moins susceptibles d'être infectés (= ont plus ou moins de chances d'être infectés dans des circonstances identiques), et on va voir que ces hétérogénéités de susceptibilité ont un impact important. (Je ne me prononce pas sur la cause de ces différences de susceptibilité : elles pourraient être dues à des différences biologiques — certaines personnes s'infectent plus facilement que d'autres — ou sociales — certaines personnes sont plus fréquemment exposées à des conditions infectieuses. Néanmoins, comme le modèle que je vais développer ici suppose que la variation de susceptibilité n'est pas corrélée à une variation d'infectiosité, c'est-à-dire que les personnes plus susceptibles ne sont pas spécialement plus infectieuses — si c'était le cas l'effet que je décris ici serait encore plus accentué — il vaut peut-être mieux imaginer le cas d'une origine biologique, parce qu'une hétérogénéité sociale a plus de chances d'être symétrique.)

Ce qui est assez surprenant, c'est que cette idée, qui peut paraître compliquée à traiter, complique en fait extrêmement peu le modèle SIR, et qu'on peut trouver des réponses exactes à essentiellement les mêmes questions que pour SIR classique (du genre quel est le nombre maximal d'infectés ?) dans ce cadre plus complexe, donnée la distribution (initiale) de susceptibilité dans la population. En général les réponses feront intervenir la transformée de Laplace de la distribution de susceptibilité (je vais expliquer ce que c'est plus bas), mais dans un cas particulier assez naturel (celui d'une distribution Γ, par exemple la distribution exponentielle), on peut tout traiter complètement.

[Un résumé de ce post de blog est contenu dans ce fil Twitter (17 tweets ; ici sur ThreadReaderApp), pour ceux qui préfèrent ce format ou qui veulent surtout les points importants (noter que tweet 11/17 il y a une typo, il faut lire φ′(0)=−1 et pas φ′(0)=1). ※ Une version anglaise (un petit peu plus longue) est contenu dans ce fil Twitter (25 tweets ; ici sur ThreadReaderApp).]

Ajout (2021-11-06) : Comme je commence à avoir écrit pas mal de choses au sujet d'épidémiologie sur Twitter, voici un méta-fil (ici sur ThreadReaderApp) rassemblant les plus importants fils que j'ai pondus.

❦

Je commence par rappeler les grandes lignes du modèle SIR classique.

Le modèle SIR classique, donc, étudie l'évolution d'une épidémie dans une population en distinguant trois classes d'individus : les Susceptibles, les Infectieux (qui dans ce modèle sont les mêmes que les infectés) et les Rétablis (qui sont immuns — ou, en fait, morts). Parmi les nombreuses hypothèses simplificatrices faites par ce modèle, il y a les suivantes (j'en oublie certainement) : l'immunité acquise par l'infection est parfaite et permanente, les individus sont infectieux dès qu'ils sont infectés, et ils vont donc soit rester dans l'état S, soit passer succesivement par les étapes S,I,R ; la population est homogène, c'est-à-dire que tous les individus sont également susceptibles et également infectieux une fois infectés, ils ont les mêmes probabilités de se faire infecter, la taille de la population est constante, et elle est assez grande pour être traitée de façon continue déterministe, et les contacts obéissent à une hypothèse de mélange parfait (au sens où tous les contacts sont également plausibles) ; le comportement de la population est constant dans le temps et notamment indépendant de l'évolution de l'épidémie ; les contaminations et le rétablissement obéissent à une cinétique du premier ordre I+S → I+I et I → R respectivement, avec des constantes β (d'infectiosité) et γ (de rétablissement) respectivement, c'est-à-dire le nombre de nouveau infectés par unité de temps est simplement proportionnel au produit du nombre d'infectieux par le nombre de susceptibles, et que le nombre de nouvellement rétablis est simplement proportionnel au nombre d'infectieux.

Bref, si on note s,i,r (quantités réelles entre 0 et 1, fonctions du temps) les proportions de la population formées d'individus susceptibles, infectieux et rétablis respectivement, alors les nouvelles infections par unités de temps se représentent par le terme β·i·s, et les rétablissements par γ·i, du coup le modèle SIR est décrit par le système d'équations différentielles ordinaires (autonomes) du premier ordre suivant :

• ds/dt = −β·i·s
• di/dt = β·i·s − γ·i
• dr/dt = γ·i
• (s+i+r=1)

où on impose en outre généralement les conditions initiales telles que s(−∞)=1, i(−∞)=0 et r(−∞)=0 (je parle bien sûr des limites en −∞), avec i croissant exponentiellement pour t assez proche de −∞ (cf. ci-dessous). La constante β d'infectiosité représente le nombre moyen de personnes qu'une personne infectieuse donnée infecte par unité de temps dans une population entièrement susceptible, tandis que la constante γ de rétablissement représente la proportion moyenne d'infectés qui se rétablissent par unité de temps (donc l'inverse du temps moyen de rétablissement, le temps de rétablissement suivant en fait une loi exponentielle). Notons que β peut aussi, symétriquement, se comprendre comme une constante de susceptibilité, c'est-à-dire comme le nombre moyen de personnes par lesquelles une personne susceptible donnée sera infectée par unité de temps dans une population entièrement infectieuse : c'est la raison pour laquelle je parlerai tantôt de β comme représentant une infectiosité et tantôt une susceptibilité (et comme ici on veut modéliser des variations de susceptibilité, c'est plutôt le deuxième qui va être mis en lumière).

Rappelons quelques uns des points saillants de ce modèle concernant le début, le pic et la fin de l'épidémie, résumé que je recopie de ce billet (plus exactement, comme je viens de le dire, 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'infectieux), on a s = 1/κ et i = (κ−log(κ)−1)/κ et r = log(κ)/κ ; notamment, le moment où l'épidémie commence à régresser correspond à i+r = 1 − 1/κ (seuil d'immunité collective) ;
• 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−Γ (taux d'attaque final).

J'ai parlé dans cette entrée de la différence entre seuil d'immunité collective et taux d'attaque final (qui sont les deux quantités essentielles que le modèle calcule) dans le modèle SIR.

*

Le modèle SIR est simpliste, mais il a ceci de bien qu'il est facile à adapter à toutes sortes de modifications en changeant les « réactions » qui le constituent (je veux dire I+S → I+I et I → R). Je donne à présent quelques exemples (même si ce n'est pas vraiment lié au sujet que je veux aborder) pour illustrer quelques variations possibles sur ce thème (on peut sauter la fin de ce paragraphe, qui servira principalement à motiver la manière dont on introduit des classes de susceptibilité). ❧ Pour commencer, si on veut ajouter un délai (appelons-ça l'état Exposé) entre le moment où on est infecté et le moment où on est infectieux, on remplace la réaction I+S → I+I par I+S → I+E et on introduit E → I (avec une nouvelle constante cinétique, disons α) en gardant I → R, ce qui donnerait les équations ds/dt = −β·i·s, de/dt = β·i·s − α·e, di/dt = α·e − γ·i et dr/dt = γ·i avec s+e+i+r=1. ❧ Pour donner un autre exemple, si on veut une immunité qui ne dure qu'un certain temps (décroît exponentiellement), on introduit une réaction R → S (avec une nouvelle constante cinétique, disons δ), ce qui donnerait les équations ds/dt = −β·i·s + δ·r, di/dt = β·i·s − γ·i et dr/dt = γ·i − δ·r avec s+i+r=1. ❧ Pour donner encore un autre exemple (peut-être plus intéressant de nos jours, mais sur lequel je n'ai pas énormément à dire — en tout cas pas aujourd'hui), si on a deux variants de la maladie, avec des contagiosités différentes mais induisant une parfaite immunité croisée, en appelant I₁ et I₂ les infectés par ces deux variants, on va remplacer la réaction I+S → I+I par I₁+S → I₁+I₁ avec une constante β₁ et I₂+S → I₂+I₂ avec une constante β₂, et bien sûr I → R par I₁ → R et I₂ → R (qui peuvent là aussi avoir deux constantes différentes, γ₁ et γ₂, si les deux variants induisent des temps de rétablissement différents) ; ceci conduit alors aux équations suivantes : ds/dt = −β₁·i₁·s −β₂·i₂·s, di₁/dt = β₁·i₁·s − γ₁·i₁, di₂/dt = β₂·i₂·s − γ₂·i₂ et dr/dt = γ₁·i₁ + γ₂·i₂ avec s+i₁+i₂+r=1. ❧ Bref, on comprend à travers ces différents exemples qu'il est facile de faire toutes sortes de variations de SIR pour décrire des situations du même genre avec différentes complexités additionnelles ou modifications de cet acabit.

❦

Mais je veux dans ce billet évoquer la manière de modéliser la situation suivante : au lieu de faire l'hypothèse (faite dans SIR) que toute la population est également susceptible à la maladie (i.e., que, dans des circonstances données, tout le monde a la même probabilité d'être infecté), que se passe-t-il s'il y a une hétérogénéité de susceptibilité ? Il peut être surprenant, mais pas tant que ça si on a lu cette entrée, d'apprendre que des hétérogénéités de contagiosité seules ne changent absolument rien à la dynamique de l'épidémie (sauf si elles sont, par exemple, corrélées à autre chose) : elles se moyennent simplement ; en revanche, des hétérogénéités de susceptibilité ont un impact énorme, et c'est ce qu'on veut voir ici. (Au cas où la différence entre hétérogénéités d'infectiosité et hétérogénéités de susceptibilité ne serait pas claire, ce que je veux dire c'est que dans la réaction I+S → I+I où un individu infectieux I infecte un individu susceptible S pour donner deux infectieux, s'il y a des différences de cinétique qui dépendent du premier (I) cela ne change rien à la dynamique d'ensemble alors que s'il y en a qui dépendent du second (S), cela change beaucoup les choses et mon but est de l'expliquer ici. Mais du coup je ne comprends décidément pas pourquoi on se focalise tellement sur les superspreaders, donc les hétérogénéités d'infectiosité, mais pas les hétérogénéités de susceptibilité qui semblent avoir été extrêmement peu étudiées.)

Le phénomène auquel on s'attend intuitivement est, bien sûr, que les individus les plus susceptibles soient infectés, donc rendus immuns, en premier, donc que l'accumulation d'immunité ne diminue pas seulement le nombre total de susceptible mais aussi la susceptibilité moyenne de ceux qui le sont, et qu'elle soit ainsi plus efficace. Le but de ce qui suit est d'appuyer ce raisonnement intuitif par une modélisation mathématique précise.

(Ce que je vais raconter ici est probablement entièrement contenu dans cet article, mais j'ai préféré retrouver les choses moi-même plutôt que lire ce qu'ils ont écrit, parce que c'est plus instructif, et du coup je n'utilise probablement pas les mêmes notations ni exactement la même approche. Par ailleurs, si on n'est pas intéressé par leur dérivation, on peut sauter directement aux équations voire directement à le cas particulier de celles-ci sur la distribution Γ.)

❦

La façon la plus évidente de modéliser des variations de susceptibilité est, sur le modèle des différentes variations autour de SIR que j'ai évoquées ci-dessus, d'introduire plusieurs classes de susceptibilité, disons deux pour montrer l'exemple, soit S₁ et S₂ (à la place de S), et de remplacer S+I → I+I par les réactions S₁+I → I+I et S₂+I → I+I avec des constantes β₁ et β₂ distinctes (il n'y a pas de raison de distinguer I en I₁ et I₂ si la contagiosité et le temps de rétablissement sont le même). Ceci conduit aux équations suivantes :

• ds₁/dt = −β₁·i·s₁ et ds₂/dt = −β₂·i·s₂
• di/dt = β₁·i·s₁ + β₂·i·s₂ − γ·i ou plutôt di/dt = (β₁·s₁ + β₂·s₂)·i − γ·i
• dr/dt = γ·i
• avec s₁+s₂+i+r=1

Évidemment, ceci peut se faire pour n'importe quelle autre valeur que 2 : si on veut douze classes de susceptibilité distinctes, on voit très bien comment en écrire les équations.

Maintenant, plutôt qu'avoir 12 ou 1729 classes de susceptibilité avec autant de constantes β, on peut préférer l'approche consistant à paramétrer les classes de susceptibilité par une (nouvelle) coordonnée, appelons-la x, proportionnelle à la susceptibilité, cette dernière étant alors β·x où β est une susceptibilité « standard » et x le rapport de la susceptibilité de la personne considérée à cette susceptibilité standard, et maintenant il est logique de passer à la limite continue. Autrement dit, on va avoir pour inconnue dans le système une fonction de deux coordonnées s(x,t) représentant le profil de susceptibilité au temps t, c'est-à-dire la proportion, au temps t, de susceptibles ayant susceptibilité β·x. (Plus exactement, s(x,t) est la limite quand dx tend vers 0 de la proportion de la population formée de susceptibles dont la susceptibilité est comprise entre β·x et β·(x+dx), divisée par dx.) La proportion totale de susceptibles (au temps t) est S(t) := ∫s(x,t)·dx. Et on a toujours des fonctions de la seule coordonnée temps i(t) et r(t). Les équations deviennent (de façon complètement analogue à ci-dessus mais où maintenant s(x,t) joue le rôle des s_j et où β·x joue le rôle des β_j) les suivantes :

• ∂s/∂t = −β·x·i·s (à lire comme : ∂s(x,t)/∂t = −β·x·i(t)·s(x,t)),
• di/dt = β·∫(x·s·dx)·i − γ·i (où ∫(x·s·dx) dénote l'intégrale de x·s(x,t) par rapport à x, qui est une fonction de t),
• dr/dt = γ·i
• (S+i+r=1, c'est-à-dire ∫s·dx+i+r=1)

Maintenant, sous cette forme, le système est assez peu maniable. Que peut-on en faire ? Regardons la première équation, ∂s/∂t = −β·x·i·s : connaissant i (comme fonction de t), on peut lui trouver une solution sous la forme s(x,t) = s₀(x)·exp(−x·f(t)) : ici, s₀ est une fonction uniquement de la coordonnée x et f une fonction uniquement de la coordonnée t vérifiant df/dt = β·i (on vérifie facilement que c'est ce que devient l'équation ∂s(x,t)/∂t = −β·x·i(t)·s(x,t) appliquée à l'ansatz s(x,t) = s₀(x)·exp(−x·f(t))). Si on impose comme condition initiale (limite en −∞) que, disons, f(−∞)=0, et comme on veut S(−∞)=1 (initialement tout le monde est susceptible), alors s₀(x) est d'intégrale 1, i.e. est une distribution de probabilités sur la coordonnées x.

Ce s₀ se comprend comme le profil de susceptibilité initial, c'est-à-dire la distribution de la variable x (susceptibilité normalisée) dans la population avant toute infection. C'est donc notre donnée fondamentale décrivant l'hétérogénéité de susceptibilité dans la population. Notre but est de comprendre l'évolution de l'épidémie en supposant connu s₀, et de voir comment ce dernier impacte cette évolution. Je vais expliquer dans un instant que ce qui va jouer un rôle clé est surtout la transformée de Laplace φ de s₀. Mais en attendant, notons qu'en plus d'imposer ∫s₀(x)·dx = 1 (i.e., qu'on a affaire à une distribution de probabilités) comme je viens de le dire, on peut de plus imposer ∫x·s₀(x)·dx = 1, c'est-à-dire que cette distribution est d'espérance 1, quitte à modifier la constante β pour que ce soit le cas, c'est-à-dire, en imposant que la susceptibilité « standard » β soit la susceptibilité moyenne de la population avant toute infection. Ce ne sera pas nécessaire, mais ça simplifie un certain nombre de calculs.

La fonction f, quant à elle, est donnée par f(t) = β · ∫_−∞^t i(t)·dt d'après la condition sur sa dérivée et la condition initiale qu'on a choisie : c'est, si on veut, une sorte de compte cumulatif des opportunités d'infections depuis le début de l'épidémie.

Le nombre total S(t) := ∫s(x,t)·dx de susceptibles au temps t est alors égal à ∫s₀(x)·exp(−x·f(t))·dx, c'est-à-dire S(t) = φ(f(t)) où φ(u) := ∫s₀(x)·exp(−u·x)·dx est l'espérance de la quantité exp(−u·x) lorsque x est distribué selon la loi s₀ (le profil initial). Cette fonction φ(u) s'appelle la transformée de Laplace de la fonction s₀ (en probabilités on l'appelle aussi, à un signe près, la fonction génératrice des moments). On peut de même utiliser φ pour exprimer la quantité ∫x·s(x,t)·dx qui intervient dans les équations (est qui est en quelque sorte la susceptibilité totale) : en intégrant par parties, on voit qu'elle vaut −φ′(f(t)), où φ′ est la dérivée de cette transformée de Laplace (par rapport à son paramètre u). Quant à la susceptibilité moyenne, qui est le rapport entre cette susceptibilité totale ∫x·s(x,t)·dx et le nombre S(t) = ∫s(x,t)·dx de susceptibles, on peut remarquer qu'elle veut −φ′/φ (moins la dérivée logarithmique de la transformée de Laplace), toujours évaluée en f(t). Bref, nos équations se réécrivent (de façon un peu redondante) :

• df/dt = β·i
• dS/dt = β·φ′(f)·i
• di/dt = − β·φ′(f)·i − γ·i
• dr/dt = γ·i
• (S+i+r=1 ; S=φ(f))

Ceci ramène donc, une fois connu le profil de susceptibilité s₀, donc sa transformée de Laplace φ, de calculer l'évolution ultérieure de l'épidémie par un système d'équations différentielles ordinaires comme précédemment.

Comme mon but est de me ramener à une présentation aussi proche que possible du SIR initial, je vais maintenant oublier le s de deux variables que j'avais avant, et renommer en s ce qui s'appelle S ci-dessus. L'idée étant que la quantité s (ex-S, donc), nombre total de susceptibles, suffit à déterminer l'évolution de l'épidémie.

En effet, la fonction φ est continue et strictement décroissante, donc injective (i.e., bijective sur son image), ce qui permet légitimement d'écrire f=φ⁻¹(s) à la place de s=φ(f), où φ⁻¹ est la fonction réciproque de φ. Bref, on peut oublier la fonction f et écrire :

• ds/dt = β·i·φ′(φ⁻¹(s))
• di/dt = − β·i·φ′(φ⁻¹(s)) − γ·i
• dr/dt = γ·i
• (s+i+r=1)

où s,i,r sont comme avant les proportions totales de susceptibles, infectieux et rétablis respectivement, et [je répète pour ceux qui selon mes indications auraient sauté directement à ce point] φ est une fonction (strictement décroissante) connue, à savoir la transformée de Laplace φ(u) := ∫s₀(x)·exp(−u·x)·dx du profil s₀ de susceptibilité de la population avant infection (la susceptibilité étant β·x), avec la normalisation que ∫s₀(x)·dx = 1 (i.e., φ(0)=1) et éventuellement ∫x·s₀(x)·dx = 1 (i.e., φ′(0)=−1). (Je donnerai plus loin un exemple de famille de distributions s₀, donc de fonctions φ, qui sont à la fois mathématiquement maniables et biologiquement plausibles.)

Je vais appeler le système ci-dessus le modèle SIR à susceptibilité hétérogène. On a donc le même système que pour le SIR classique (=homogène), mais le terme de nouvelles infections β·i·s est remplacé par −β·i·φ′(φ⁻¹(s)) où φ′∘φ⁻¹ est une fonction connue (remarquons que φ′(φ⁻¹(s)) est la pente du graphe de φ à l'ordonnée s). Le modèle classique (=homogène) correspond au cas où φ(u) = exp(−u) (transformée de Laplace d'une distribution delta de Dirac en 1 puisque tout le monde a la même susceptibilité β, cas qui n'était pas couvert par le système précédent mais qui l'est par le système tel que je l'ai écrit). La quantité −φ′(φ⁻¹(s)) est une « susceptibilité totale » de la population (normalisée par β), et −φ′(φ⁻¹(s))/s s'interprète comme la susceptibilité moyenne restante (c'est-à-dire la susceptibilité moyenne des individus qui sont susceptibles, là aussi normalisée par β).

En effaçant une partie des équations, j'ai cependant perdu quelque chose de précieux permettant de le résoudre partiellement : en effet, φ⁻¹(s) (qui était noté f ci-dessus) a une dérivée par rapport à t valant β·i qui est proportionnelle à celle de r soit γ·i. Donc en notant κ := β/γ le nombre de reproduction et en rappelant qu'on fait l'hypothèse sur les conditions initiales que f et r valent 0 en −∞, on a l'invariant suivant :

☞ s = φ(κ·r)

(autrement dit, non seulement la fonction φ permet d'écrire les équations mais même elle les résout en partie).

J'ai rappelé plus haut trois éléments de l'analyse du comportement du modèle SIR classique : l'exponentielle initiale, le pic épidémique, et le taux d'attaque final. Que deviennent-ils dans le cas hétérogène ?

• Le comportement exponentiel initial n'est pas modifié (si on a normalisé par φ′(0)=−1 comme je le proposais, c'est-à-dire que β est bien la susceptibilité moyenne avant toute infection, alors le comportement est exactement le même : pour s proche de 1, on a −φ′(φ⁻¹(s)) également proche de 1 et tout se passe exactement pareil à l'ordre le plus bas).
• Néanmoins, si on va chercher le terme d'ordre suivant en s de −φ′(φ⁻¹(s)) (pour s≈1), cette quantité vaut −φ′(0) + (φ″(0)/φ′(0))·(1−s) + O((1−s)²). En interprétant −φ′(0) comme l'espérance d'une variable aléatoire distribuée selon la loi s₀ de susceptibilité initiale, et φ″(0) comme l'espérance du carré de cette variable (soit la variance plus le carré de l'espérance), alors ceci nous permet de dire la chose suivante : une petite accumulation d'immunité, dans le cas hétérogène, est 1+v fois plus efficace (sur la diminution du nombre de reproduction effectif) que dans le cas homogène, où v est la variance relative de la susceptibilité initiale, c'est-à-dire le rapport (sans dimension) entre le variance et le carré de l'espérance (1+v = φ″(0)/(φ′(0))²). L'explication intuitive est que le petit nombre (1−s) déduit aux susceptibles s'accompagne d'une baisse de v·(1−s) de la susceptibilité moyenne de ceux qui restent susceptibles.
• Le pic épidémique est atteint lorsque di/dt = 0, soit lorsque φ′(φ⁻¹(s)) = −1/κ, soit φ′(κ·r) = −1/κ en se rappelant que s = φ(κ·r) (ceci permet de retrouver s et r, après quoi i s'en déduit comme 1−s−r). Graphiquement, on cherche le point du graphe de φ où la pente vaut −1/κ (i.e., où la tangente est parallèle à la droite reliant (κ,0) et (0,1)) : son ordonnée vaut s et son abscisse vaut κ·r. Le seuil d'immunité collective (par infection, donc) est 1−s pour ce point. (Cette méthode graphique est illustrée par le premier des deux graphiques ci-contre à droite ; le seuil d'immunité collective est ici environ 0.42.)
• Le taux d'attaque final est obtenu pour i=0, soit s+r=1, et s'obtient donc en résolvant φ(κ·r) + r = 1. Graphiquement, on cherche l'intersection du graphe de φ avec la droite reliant (κ,0) et (0,1) : son ordonnée vaut s et son abscisse vaut κ·r. (Cette méthode graphique est illustrée par le second des deux graphiques ci-contre à droite ; le taux d'attaque est ici environ 0.67.)

Énormément de choses se déduisent donc du graphe de cette transformée de Laplace φ de la distribution de susceptibilité initiale (soit dit en passant, la transformée de Laplace des profils de susceptibilité ultérieurs s'obtient par translation en abscisse de l'initiale).

J'ai déjà dit que le cas du modèle SIR classique (=homogène, tout le monde a la même susceptibilité) correspond au cas où φ(u) = exp(−u) car s₀ est une distribution delta de Dirac en x=1. Y a-t-il d'autres cas à la fois naturels et explicitement traitables ? Une autre distribution de probabilités tout à fait naturelle sur les réels positifs (et qui me semble être un a priori raisonnable si on ne sait rien sur une quantité que son espérance) est la loi exponentielle : si s₀(x) = exp(−x) alors φ(u) = 1/(u+1). Plus généralement, une famille de distributions de probabilités incluant l'exponentielle et ayant la distribution delta comme cas limite, mais qui permet de choisir la variance indépendamment de l'espérance, est la distribution Γ.

Spécifiquement, si le profil initial de susceptibilité s₀ suit une distribution Γ de « forme » a>1, que je peux supposer d'espérance 1 quitte à l'absorber dans β, c'est-à-dire s₀(x) = (a^a/Γ(a))·x^a−1·exp(−a·x) (où Γ(a) = ∫ x^a−1·exp(−x)·dx, servant à normaliser l'intégrale de s₀ à 1, est la fonction gamma d'Euler, qui vaut (a−1)! si a est entier), alors sa transformée de Laplace φ(u) vaut 1/((u/a)+1)^a et la dérivée φ′(u) de celle-ci vaut −1/((u/a)+1)^a+1 = −φ(u)^(a+1)/a : les équations du modèle SIR hétérogène deviennent donc :

• ds/dt = −β·i·s^(a+1)/a
• di/dt = β·i·s^(a+1)/a − γ·i
• dr/dt = γ·i
• (s+i+r=1)

ou encore, si on préfère noter v := 1/a (variance relative de s₀), on est ramené au système suivant

• ds/dt = −β·i·s^1+v
• di/dt = β·i·s^1+v − γ·i
• dr/dt = γ·i
• (s+i+r=1)

— c'est-à-dire exactement le système initial sauf que le terme β·i·s de nouvelles infections a été remplacé par β·i·s^1+v : tout se passe comme si la cinétique I+S → 2I était remplacée par I+(1+v)·S → 2I+v·S. En fait, ce qui se produit est que quand la proportion susceptible diminue par accumulation d'immunité, l'espérance de la susceptibilité relative de ceux qui le sont évolue comme la puissance v-ième de s (et donc la susceptibilité totale comme la puissance (1+v)-ième de s).

Le cas d'une distribution exponentielle est le cas particulier de variance relative v=1 (c'est-à-dire de forme a=1) de la distribution Γ ; le cas homogène (SIR classique : distribution delta) est la limite de variance relative v=0 (c'est-à-dire de forme a→+∞).

Si je reprends dans le cas particulier de la distribution Γ ce que j'ai dit sur le système SIR hétérogène en général sur le comportement en petit temps, le pic épidémique et le taux d'attaque final :

• Le comportement exponentiel initial est le même que pour le SIR classique (=homogène). Mais quand il s'accumule un peu d'immunité, celle-ci est initialement 1+v fois plus efficace (sur la diminution du nombre de reproduction effectif) que dans le cas homogène.
• Le pic épidémique est atteint pour s = κ^−1/(1+v) (soit κ^−a/(a+1)), c'est-à-dire que le seuil d'immunité collective[#] vaut 1 − κ^−1/(1+v) (on a précisément r = (κ^−1/(1+v)−κ⁻¹)/v, et i = 1 − ((1+v)/v)·κ^−1/(1+v) + (1/v)·κ⁻¹ au moment du pic, cette dernière expression donnant donc la proportion maximale d'infectés).
• Le taux d'attaque final est la solution r>0 de (κ·v·r+1)^−1/v + r = 1 ; je ne crois pas qu'on puisse simplifier ça plus que ça, mais pour v=1 (le cas exponentiel) on trouve 1 − 1/κ (c'est-à-dire que le taux d'attaque final dans le cas exponentiel est égal au seuil d'immunité collective dans le cas homogène : je ne sais pas s'il y a une explication non-calculatoire de ce fait).

[#] Pour être bien clair, il s'agit là du seuil d'immunité collective par infection, qui bénéficie des effets d'hétérogénéité que je viens de décrire. Le seuil d'immunité collective par vaccination n'a pas de raison d'être différent de 1 − 1/κ (si on vaccine aléatoirement).

Les graphiques ci-contre (cliquer pour agrandir) illustrent un peu l'allure de ces fonctions : la variance relative v de susceptibilité est en abscisse, entre 0 (correspondant au SIR classique) et 2, avec 1 (le cas exponentielle) au milieu ; l'ordonnée représente une proportion de la population : les courbes rouges sont celles du seuil d'immunité collective par infection, les bleues sont celles du taux d'attaque final pour une épidémie non contrôlée ; le nombre de reproduction vaut 2 pour les courbes pleines, 3 pour les courbes en tirets et 4 pour les courbes en pointillés.

Ajout (2021-03-21) : Au niveau de la dynamique (temporelle, je veux dire) de l'épidémie, l'effet de la variance n'est pas extrêmement frappant sur l'allure des courbes. Les graphes ci-contre (cliquer pour agrandir) montrent l'évolution d'une épidémie décrite par un SIR classique (=homogène, soit v=0) sur le premier jeu de quatre courbes, et par un SIR hétérogène à susceptibilité distribuée selon une loi exponentielle (v=1) sur le second jeu de quatre courbes, dans les deux cas avec un nombre de reproduction de κ=3 : dans chaque image, les courbes sont tracées en fonction du temps compté en temps de rétablissement (1/γ) ; la courbe en haut à gauche montre les valeurs de s (en vert), i (en rouge) et r (en bleu) ; celle en haut à droite est la même courbe i mais à une échelle verticale différente pour plus de lisibilité (j'ai réutilisé du code où j'affichais des choses plus détaillées en haut à droite qui ne sont pas, ici, pertinentes) ; la courbe en bas à gauche est la même qu'en haut à gauche mais en échelle logarithmique ; et la courbe en bas à droite montre le nombre de reproduction effectif en fonction du temps. (Le code Sage est ici, il faut éditer quelques réglages triviaux pour obtenir exactement les courbes ci-contre, mais je suppose que ce sera facile à trouver.) Je suppose qu'on sera d'accord avec moi que la différence qualitative ne saute pas aux yeux : on voit certes que l'épidémie monte moins haut et attaque finalement moins dans le second cas, mais l'allure est très semblable ; on pourrait se dire que le second jeu de courbes est le résultat d'un SIR classique avec un nombre de reproduction plus faible, mais en fait non, parce que la croissance exponentielle des cas quand même bien la même dans les deux cas (ce n'est que quand on accumule un peu d'immunité que l'effet de l'hétérogénéité se fait sentir).

Il faut que j'insiste sur le point suivant : si l'hypothèse que dans une infection réelle la susceptibilité suive une loi Γ est un peu arbitraire, elle est néanmoins naturelle et pas du tout fantaisiste, et en tout cas c'est un modèle approximatif raisonnable d'une situation avec une hétérogénéité non nulle, paramétrée par la variance relative v : même si la distribution n'est pas spécifiquement une Γ, le système ci-dessus devrait être une approximation raisonnable de ce qui se passe avec une variance relative v. Le paramètre (sans dimension) v de variance relative de la susceptibilité doit être considéré comme faisant partie des données épidémique et est aussi essentiel que le nombre de reproduction κ (lui aussi sans dimension) pour modéliser l'épidémie. Postuler que ce paramètre vaut 0 (le modèle SIR classique) est une hypothèse déraisonnable si elle n'est pas appuyée par des observations expérimentales (or dans une épidémie où il est clair que les enfant sont beaucoup moins susceptibles que les adultes, c'est déjà impossible d'avoir v=0, en fait). Cela fait partie de mon slogan général prédire une exponentielle est facile, mais prédire quand cette exponentielle s'arrête est toute la difficulté, or on ne dispose pas des données pour ça. Mais à tout le moins, si on ne sait rien du tout partir sur l'hypothèse que v=1 est plus naturel pour des raisons de simple analyse dimensionnelle (ne sachant rien sur l'écart-type de la susceptibilité, on peut imaginer qu'il est de l'ordre de grandeur de la susceptibilité moyenne elle-même), et de fait, la distribution exponentielle est quelque chose qu'on retrouve assez souvent dans la nature (et le fait qu'il existe des personnes très peu susceptibles plaide en faveur de v≳1). Bref, la formule 1 − 1/κ pour le seuil d'immunité collective est raisonnable pour le seuil d'immunité collective par vaccination, mais par infection il faut considérer que la bonne formule est 1 − κ^−1/(1+v) où, à défaut d'avoir des informations sur v, on prendra v=1, donc 1 − 1/√κ.

Il faudrait que j'explique ce qui se passe quand en plus d'avoir des hétérogénéités de susceptibilité on en a en plus d'infectiosité qui sont corrélées avec elles, et aussi ce qui se passe si on a deux variants qui ont non seulement des nombres de reproduction différents mais même des hétérogénéités différentes (des hétérogénéités d'hétérogénéité, si on veut !), mais je commence vraiment à fatiguer, donc je vais en rester là pour le moment.