David Madore's WebLog: Quelques considérations de graphes aléatoires pour l'épidémiologie

J'en viens à ce dont je voulais vraiment parler : un modèle basé sur la percolation dans des graphes aléatoires et permettant de modéliser (de façon simpliste !) la manière dont la variance du nombre de contacts infectieux modifie le taux d'attaque d'une épidémie à nombre de reproduction R₀ donné. C'est ce que représentent les courbes ci-contre, en l'occurrence pour R₀=2.5 (contacts infectieux par individu en moyenne), avec l'écart-type σ du nombre de contacts infectieux en abscisse, et en ordonnée le taux d'attaque prédit (en bleu par un modèle basé sur un graphe orienté, en rouge par un modèle symétrique) : je veux expliquer un peu comment lire ces courbes et comment elles ont été calculées.

L'idée est qu'on va modéliser l'épidémie non plus comme un processus dynamique traduit par une équation différentielle (comme dans le modèle SIR classique ou une équation différentielle à délai comme dans la variante que j'avais étudiée plus tard) mais comme un parcours (un phénomène de percolation) dans un « graphe infectieux ». L'avantage est de pouvoir modéliser des phénomènes un peu plus fins, ici l'hétérogénéité des contacts ; l'inconvénient est qu'on perd complètement la notion du temps, toute la dynamique de l'épidémie, il faut faire l'hypothèse (évidemment irréaliste en vérité — mais forcément pessimiste) que cette dynamique ne change pas au cours du temps (notamment, le nombre de reproduction est une constante), et on ne peut que calculer une seule chose, à savoir le taux d'attaque final, c'est-à-dire la proportion de la population qui se fait infecter (sans qu'on puisse dire quand). Comme j'aime bien le rappeler, quand on a affaire à un phénomène complexe, il faut multiplier les approches, étudier une complexité à la fois, essayer des modèles simples, voire simplistes, qui permettent de bien cerner l'effet de chacune de ces complications, avant d'envisager un modèle plus complexe et réaliste. Donc ici on supprime la complication « dynamique » et on rajoute la complication « hétérogénéité » (mais on va heureusement trouver la même prédiction lorsque les deux complications sont absentes ! à savoir le taux d'attaque final prédit par le modèle SIR au point que j'appelle le point de Poisson dans ce qui suit). Dans tous les cas, on travaille sur des hypothèses simplistes sur le fait que l'immunité acquise est permanente : une fois un individu rétabli, il cesse d'être infectieux et ne pourra plus le redevenir.

Je rappelle, donc, que le modèle SIR classique prédit, pour un nombre de reproduction R₀=κ (j'utiliserai de façon quasi interchangeables les notations R₀ et κ ici), un taux d'attaque (= proportion de la population qui est infectée à la fin de l'épidémie) donné par la formule exacte 1 + W(−κ·exp(−κ))/κ où W est une fonction spéciale (la fonction transcendante W de Lambert, ou logarithme-produit), et approximativement 1 − exp(−κ) pour κ grand (et 2(κ−1) si κ est juste un peu supérieur à 1). Il prédit aussi un pic épidémique avec un taux d'attaque à ce moment (qu'on appelle le seuil d'immunité grégaire) à 1 − 1/κ, mais ce n'est pas celui qui va m'intéresser ici (il va juste recevoir une mention au passage), c'est vraiment le taux d'attaque final que je veux : entre les deux il y a un overshoot, simplement parce que l'épidémie continue à progresser même si elle est en phase de ralentissement.

Ici on va faire complètement autrement : on définit un graphe de contamination (commençons par la variante orientée), dont les sommets sont les individus de la population considérée, et avec la propriété que les personnes infectées par l'épidémie sont toutes celles qui sont atteignables depuis le « patient zéro » x₀, i.e., tous les z pour lesquels il existe un chemin orienté, représentant une chaîne de contamination possible, x₀ → x₁ → ⋯ → z, dans le graphe. La définition exacte de ce que sont les arêtes x→y du graphe est un petit peu épineuse : ce ne sont pas juste les contacts potentiels, ni même les contacts réels, ni les contacts réels potentiellement infectieux, mais les contacts potentiellement infectieux pendant une période d'infectiosité de x, je vais y revenir. Le fait est que le nombre moyen d'arêtes x→y sortant d'un sommet x, i.e., le degré sortant moyen (moyenné sur l'ensemble des sommets), et qui est aussi le nombre moyen d'arêtes x→y aboutissant à un sommet y, i.e., le degré entrant moyen, est ce qu'on appelle le nombre de reproduction κ=R₀, de l'épidémie. L'idée est donc de chercher des moyens de calculer le taux d'attaque, c'est-à-dire la proportion des sommes accessibles depuis x₀, en fonction du degré moyen et d'informations sur la distribution des degrés (entrants et/ou sortants) : pour ça, bien sûr, il faut faire des hypothèses sur la manière dont le graphe est construit, selon un modèle de graphe aléatoire. Sauf qu'en fait, on ne construit pas vraiment le graphe (qui de toute façon est aléatoire, donc pas très intéressant), on utilise des résultats venus des probabilités, de la théorie des graphes aléatoires, qui donnent des informations sur la taille des composantes de graphes aléatoires, et notamment des formules pour calculer (l'espérance de) la taille de la « composante géante », en fonction de la distribution des degrés utilisée en entrée du modèle de graphe aléatoire. Au final, il y a une formule assez simple (que je vais expliquer) permettant de calculer cette taille d'après la fonction génératrice de la distribution de probabilité des degrés ; puis il faut faire une hypothèse sur la distribution utilisée permettant de jouer sur la variance.

Une distribution de probabilités qui va jouer un rôle proéminent ici est la distribution de Poisson qui correspond à la variance σ²=κ (sur mes courbes ci-dessus, c'est le point (1.581, 0.893) où les courbes rouge et bleue se croisent : σ=1.581 est la racine carrée de κ=2.5), et dont il va donc falloir que je reparle. Je souligne d'ores et déjà son importance, d'une part parce que c'est le point où on retrouve le taux d'attaque (89.3% pour un nombre de reproduction de 2.5) prédit par le modèle classique SIR, et d'autre part parce que ça a été une source de confusion fréquente quand j'ai parlé de cette histoire, les gens s'imaginent que si tout le monde se comporte de la même manière les degrés ont une variance nulle (et me demandent pourquoi mes courbes donnent 1 en σ=0) alors que la variance « de base », si j'ose dire, c'est celle de la distribution de Poisson (on peut aller en-dessous mais c'est bizarre). Un autre point qui va jouer un rôle, plus secondaire, c'est celui qui correspond à la distribution géométrique, donnant la variance σ²=κ(κ+1) (soit σ = √(κ(κ+1)) ≈ κ, ce qui peut entraîner une confusion avec la distribution de Poisson !), (2.958, 0.600) sur ma courbe bleue, pour l'instant je me contente de le mentionner pour qu'on ne le confonde pas avec l'autre.

Je précise que je vais commettre des approximations mathématiques dans ce qui suit, qui me vaudraient certainement les foudres de Nicolas Bourbaki. À titre d'exemple, comme je considère un graphe avec un nombre très important de sommets (auquel je ne vais même pas me fatiguer à donner un nom), je vais confondre allègrement le degré moyen d'un sommet et l'espérance de la distribution selon laquelle les degrés ont été tirés, ou bien la variance du degré avec la variance de la distribution de probabilités selon laquelle ils sont tirés : cette confusion dans l'exposition est motivée par la loi des grands nombres, mais à la limite ce n'est même pas le problème, ici (de toute façon, ce qui compte est bien l'espérance et la variance de la distribution à partir de laquelle on a tiré les degrés, le graphe lui-même n'est pas intéressant, il n'importe que pour l'intuition) ; et je commettrai des approximations bien plus graves que ça et des « démonstrations » en agitant les mains : je fournirai cependant des pointeurs vers des théorèmes mathématiques plus précis.

Reprenons. On veut ignorer complètement le temps qui passe (mais bien sûr obtenir quand même une description correcte de l'état final !, en faisant l'hypothèse que la dynamique reste constante dans le temps) et modéliser l'épidémie par un graphe de contamination. (Considérer ce graphe suppose qu'on laisse, conceptuellement, l'épidémie se dérouler jusqu'à son terme, et on peut ensuite définir le graphe et calculer des choses dessus.) Quel est ce graphe, au juste ? (Je commence par le cas orienté, parce qu'il me semble plus naturel et plus simple à décrire.)

Les sommets du graphe (je dirai parfois aussi les nœuds) sont les individus. Mais la définition d'une arête x→y est un peu subtile. Si x est infecté au cours de l'épidémie, alors les arêtes x→y sortant de x pointent vers les individus y qui ont été contaminés par x ou qui seraient contaminés s'ils sont effectivement susceptibles (si y est déjà infecté, ou déjà immun, on fait quand même pointer une arête x→y s'ils ont un contact tel que y aurait été infecté s'il n'était pas déjà infecté ou immun). Autrement dit, il s'agit des contacts potentiellement infectieux qu'a x au cours de sa période d'infectiosité. Mais si x n'est jamais infecté, on va quand même définir des arêtes x→y : pour ça, on choisit une période d'infectiosité typique (c'est-à-dire selon la loi de probabilité que suit effectivement l'infectiosité des personnes infectées) pour x, et on trace les contacts qui seraient infectieux pendant cette période si x était infectieux. Ça peut sembler vraiment bizarre de faire ça, mais je vais expliquer la raison.

Bref, une arête du graphe, et je dirai simplement un contact infectieux x→y, est un contact entre x et y qui ① a lieu pendant une période d'infectiosité typique pour x, et, si x est effectivement infecté à un moment donné, sa période d'infectiosité réelle, ② peu importe que y soit susceptible (=infectable) ou non, et ③ qui rendrait y infecté si y était susceptible.

Pourquoi cette définition tarabiscotée (que j'ai mis beaucoup de temps à éclaircir dans ma tête) ? Si on s'intéresse à la propagation de l'épidémie, pourquoi ne pas simplement mettre une arête x→y si x infecte y et basta ? La raison est que dans ce cas le graphe ne ressemblerait pas à un graphe aléatoire facilement modélisable : les personnes infectées tard dans l'histoire de l'épidémie en infectent moins autour d'eux parce que plus de gens autour d'eux sont immuns, donc le degré des sommets loin de x₀ serait plus faible que celui des sommets proches de x₀ et cela rendrait le calcul plus difficile à modéliser (comme chaque y aurait au plus une arête entrante, le degré entrant moyen serait exactement le taux d'attaque qu'on cherche à calculer ! au lieu que ce soit le nombre de reproduction à partir duquel on cherche à mener le calcul) : on préfère avoir une arête x→y même si y n'est plus infectable, de façon à ce que chaque x infecté ait le même degré sortant typique (aux hétérogénéités intrinsèques près, qu'on cherche justement à étudier !), et il reste vrai que les sommets accessibles à partir de x₀ sont les personnes qui se font infecter ; ceci explique qu'on ne se limite pas aux y qui sont susceptibles. Quant qu fait qu'on ne se limite pas aux x qui sont effectivement infectés, la raison est analogue : on veut pouvoir raisonner sur le degré moyen sans avoir à distinguer les individus infectés et ceux qui ne le sont pas (puisque le but est justement de calculer la proportion qui se fait infecter, on ne veut pas avoir à la connaître a priori). Comme on a fait l'hypothèse d'indépendance du temps qui assure que le nombre de personnes potentiellement infectées par x ne dépend pas de la période infectieuse choisie, seulement de sa durée, on peut prendre une période essentiellement quelconque (mais si x est effectivement infecté, on prendra la vraie pour pouvoir relier ce graphe à la réalité).

Bref, la définition du graphe étant posée, les personnes effectivement infectées sont celles qui sont accessibles depuis le patient zéro x₀ (ou, en fait, depuis un sommet « typique ») par un chemin menant jusqu'à elle.

Je rappelle que le degré sortant d'un sommet x dans un graphe orienté est le nombre de ses voisins sortants, c'est-à-dire le nombre d'arêtes x→y partant de x (en disant ça je suppose implicitement que mon graphe n'a pas d'arêtes multiples : je n'ai pas été terriblement clair à ce sujet dans la définition des arêtes, mais ce n'est pas très important parce que mon modèle de graphe aléatoire ne va pas en produire pour un nombre important de sommets : je travaille toujours dans des conditions de « mélange parfait » des individus). Symétriquement, le degré entrant d'un sommet u est le nombre de ses voisins entrants, c'est-à-dire le nombre d'arêtes x→y aboutissant à y.

Le degré sortant de x est donc le nombre de personnes contaminées par x, ou qui le seraient si elles étaient effectivement susceptibles, dans la mesure où x est effectivement infectieux à un moment donné (et si ce n'est pas le cas, on fait comme s'il l'était, et on compte sur une période arbitraire). Le degré entrant de y est le nombre de personnes qui contamineraient y si elles sont effectivement infectieuses et que y est encore susceptible à ce moment, mais en comptant aussi les personnes non infectieuses (pendant une période d'infectiosité arbitrairement choisies pour elles).

Voici un fait évident mais essentiel : le degré sortant moyen égale le degré entrant moyen. (C'est un argument que les matheux appellent un double comptage : on compte les arêtes soit en les regroupant par leur source, soit en les regroupant par leur cible.) Ce nombre moyen s'appelle le nombre de reproduction de l'épidémie (normalement il y a un nombre de reproduction en fonction du temps, mais ici je fais justement l'hypothèse d'une dynamique stationnaire), noté R₀ et que je préfère noter κ.

En revanche, si on ne peut pas étudier les variations en fonction du temps, on peut étudier les variations entre individus, les écarts à la moyenne, c'est-à-dire, la variance. Je rappelle que la variance σ² (d'une quantité, resp. variable aléatoire) est la moyenne (resp. espérance) des carrés des déviations à la moyenne (resp. espérance), et que l'écart-type σ est la racine carrée de la variance, c'est-à-dire la moyenne quadratique des écarts à la moyenne. Plus la variance (ou ce qui revient au même, l'écart-type) est élevée, plus les valeurs tendent à être dispersées autour de la moyenne.

Bref, même si R₀=2.5 signifie que chaque personne est la source (ou la cible !) de 2.5 contacts infectieux en moyenne, certains en feront moins (ou pas du tout) et d'autres en feront plus. Une variance de 0 signifierait que tout le monde en fait exactement autant (ce qui, s'agissant de 2.5, est manifestement impossible).

Les individus ayant un degré sortant nettement plus élevé que la moyenne (et qui sont effectivement infectieux) ont été appelés supercontaminateurs : leur existence augmente la variance des degrés sortants. Bien sûr, avoir des gens qui en contaminent beaucoup d'autres augmente aussi la moyenne ! Mais à moyenne donnée fixée, une variance élevée signifie qu'il y a des supercontaminateurs et, forcément, des subcontaminateurs (voire non-contaminateurs).

Ce qui peut être surprenant est que ⓐ bien que les degrés sortants et entrants aient la même moyenne, leur variance n'a aucune raison de coïncider et leur distribution n'a pas de rapport particulier à part la moyenne, et ⓑ bien que l'épidémie se propage en suivant les arêtes orientées, donc d'un sommet vers ses voisins sortants, c'est pourtant la variance et la distribution des degrés entrants qui importe pour calculer le taux d'attaque (la propagation ultime). Je vais essayer plus bas d'expliquer pourquoi.

(La variance et distribution des degrés sortants est importante pour quelque chose d'autre, à savoir la probabilité de non-extinction de l'épidémie. Il existe, en fait, une symétrie parfaite entre degrés entrants et taux d'attaque d'une part et degrés sortants et probabilité de non-extinction de l'autre : en effet, le taux d'attaque est la proportion des sommets accessibles depuis un ensemble raisonnable de sommets initiaux, alors que la probabilité de non-extinction est la proportion de sommets (pouvant servir de patient zéro) depuis lesquels on peut atteindre un ensemble raisonnable de sommets finaux. Mais on parle ici d'une épidémie qui, de toute évidence, ne s'est pas éteinte rapidement, donc les degrés sortants ne nous intéressent plus.)

On peut donc être amené à s'intéresser, plutôt qu'aux supercontaminateurs (dont le degré sortant est très élevé), aux supersusceptibles(?), c'est-à-dire les personnes dont le degré entrant est très élevé. J'ai attiré l'attention sur ce concept dans un fil Twitter [lien direct Twitter] un peu ancien, mais le concept a été mal compris (il est vrai qu'il est subtil !), on a pensé que je disais que les supersusceptibles étaient un problème, ou quelque chose dont on doit s'inquiéter, ou sur quoi on devrait essayer d'agir : donc rappelons que dans la mesure où on parle de l'effet de la variance à moyenne donnée (on varie un paramètre à la fois, merci !), si on augmente le nombre de supercontaminateurs ou de supersusceptibles, on augmente aussi le nombre de subcontaminateurs ou subsusceptibles.

Le concept de degré sortant d'une personne x est assez clair. Le concept de degré entrant de y est plus subtil, parce que l'existence d'une arête x→y dépend du fait que x (et pas y) est dans sa phase infectieuse, ce qui n'est pas vraiment du ressort de y. Néanmoins, bien sûr, il y a des caractéristiques de y qui vont jouer sur ce degré. D'abord simplement le nombre de contacts de y : si quelqu'un a beaucoup de contacts, son degré entrant comme sortant va avoir tendance à être plus élevé (mais cette symétrie va être une raison de se pencher éventuellement sur des graphes non-orientés où la distinction entre degré entrant et sortant disparaît). Ensuite, il peut y avoir des différences de comportement jouant de façon partiellement asymétrique : certains sont plus attentifs à l'hygiène ou à la distanciation sociale : cela peut jouer sur les degrés entrant et sortant (si ces mesures évitent de contaminer les autres comme d'être soi-même contaminé) ou seulement l'un (et on peut imaginer que beaucoup de gens cherchent avant tout à se protéger : quelqu'un qui porte un masque FFP2 avec valve dès qu'il sort va sans doute avoir un degré entrant beaucoup plus faible que sortant). Enfin, il peut y avoir des raisons médicales : pour certaines infections (et je ne sais pas si c'est le cas du Covid-19), il y a des individus beaucoup plus susceptibles que d'autres, c'est-à-dire facilement infectés, comme il peut y avoir des individus beaucoup plus infectieux : ceci va jouer respectivement sur les degrés entrants et sortants.

Mais il y aussi une forme de variance « basique » (intrinsèque, naturelle), et qui est celle d'un processus de Poisson : il faut que je souligne ça, parce que ça a causé une certaine confusion, plusieurs personnes m'ayant posé des questions au sujet du point σ=0 de mes courbes, pensant spontanément (et c'est naturel) que si tout le monde se comporte de la même manière, on obtient σ=0, alors qu'en fait on obtient la variance de Poisson σ²=κ. Expliquons informellement ce qu'est un processus de Poisson : imaginez que vous ayez un million de personnes et un million de boîtes, et vous demandez aux gens de mettre des boules, autant qu'ils veulent, tirées d'un énorme sac, dans des boîtes, comme ils veulent. On observe que chaque personne place 2.5 boules en moyenne dans l'ensemble des boîtes : c'est le degré sortant moyen ; donc chaque boîte contiendra en moyenne 2.5 boules (degré entrant moyen). Mais même si les boîtes sont absolument identiques et que les participants les choisissent aléatoirement (uniformément), il y aura (évidemment !) des variations aléatoires dans le nombre de boules dans les boîtes. En fait, le nombre de boules dans une boîte suivra (asymptotiquement) une distribution de Poisson de moyenne 2.5 et de variance σ²=2.5. On peut prendre ça comme définition d'un processus de Poisson de moyenne κ : si on place N·κ boules aléatoirement (uniformément et indépendamment) dans N boîtes, le nombre de boules dans une boîte suivra, dans la limite N→∞, une distribution de Poisson de moyenne κ (et sa variance égale sa moyenne). Maintenant, si certaines boîtes sont plus attirantes que d'autres (et donc certaines moins attirantes), cela augmentera la variance, mais pour passer en-dessous de la variance de Poisson, il faut quelque chose d'inhabituel (par exemple, que les gens regardent à l'intérieur des boîtes pour faire leur choix).

Pour une épidémie, si les gens se comportent tous à l'identique et se font contaminer au hasard, on obtient donc une distribution de Poisson sur les contacts infectieux reçus, et σ²=R₀. Faire moins de variance que ça serait bizarre et peu naturel pour une épidémie. (Mais je donne quand même un exemple : l'« infection par la mortalité », où (presque ?) tout le monde se fait infecter par exactement deux personnes, ses parents, donc R₀=2 (dans une population stationnaire), σ=0, et le taux d'attaque est de 100% puisque tout le monde est mortel ; le graphe est simplement l'arbre généalogique de la population. Noter que si la variance des degrés entrants est nulle, celle des degrés sortants est probablement élevée, je la soupçonne d'être nettement au-dessus de celle du processus de Poisson.)

Bref, le point où σ²=R₀ (soit (1.581, 0.893) sur les courbes que j'ai tracées), avec la famille de distributions que je dois encore expliquer, est ce que j'appellerai le point de Poisson : tout le monde se comporte de façon identique, les contacts sont aléatoires, les degrés entrants sont distribués selon un processus de Poisson, et le taux d'attaque qu'on trouve (par le mode de calcul que je dois encore expliquer !) coïncide avec le taux d'attaque 1 + W(−κ·exp(−κ))/κ prédit par le modèle SIR classique (ou sa variante à rétablissement en temps constant), et ce n'est pas un hasard, on peut construire un modèle (SIR sochastique en temps continu sur un graphe complet de taille tendant vers l'infini) qui relie ces deux description (l'espérance de la proportion de sommets contaminés suivra les équations différentielles du SIR classique, tandis que le graphe des contaminations reproduira le modèle que je cherche à décrire ici, ce qui fait le pont entre les deux).

Maintenant, il n'y a pas que la variance des degrés sortants qui importe, il y a la distribution précise des degrés, et il y a en plus, comme toile de fond, le modèle de graphe aléatoire. [Je suis un peu embêté pour décider dans quel ordre raconter tout ça, donc je vais suivre l'ordre de mes fils Twitter, quitte à me répéter, en disant d'abord les choses de façon plus informelle, puis en rentrant dans plus de détails mathématiques.]

J'ai choisi une famille de distributions (sur les entiers naturels) qui permette de faire bouger la variance σ² (ou ce qui revient au même, l'écart-type σ, qui sert d'abscisse pour mes courbes) à espérance κ donnée, et qui inclut la distribution de Poisson qui est manifestement de la plus grande importance ici, comme cas limite, et inclut aussi la distribution géométrique dont j'ai des raisons de penser qu'elle est aussi assez naturelle, comme cas particulier : il s'agit de la famille des distributions binomiales, ou plus exactement, la distribution binomiale pour σ²<κ, et la distribution binomiale négative (ou de Pólya) pour σ²>κ (avec Poisson comme cas limite σ²=κ, et géométrique comme cas particulier σ²=κ(κ+1)). La partie σ²<κ est un petit peu tirée par les cheveux, mais de toute façon j'ai déjà dit qu'elle n'est pas très naturelle et pas très intéressante : elle est juste incluse pour l'élégance mathématique des courbes.

Bref, la courbe bleue représentée plus haut est un taux d'attaque modélisé comme (l'espérance de) la proportion des sommets accessibles, depuis un sommet x₀ pour lequel cette proportion est importante, dans un graphe orienté aléatoire construit (selon un modèle que je n'ai pas encore décrit précisément) avec une distribution des degrés entrants suivant une des distributions décrites au paragraphe précédent, pour la moyenne κ=2.5 et l'écart-type σ en abscisse. Et la courbe rouge ? C'est la même chose mais avec un graphe non-orienté : c'est la taille de la composante connexe géante d'un graphe aléatoire construit (selon le modèle aléatoire dit de la configuration de Bollobás, ou de Molloy-Reed), avec une distribution des degrés de moyenne κ=2.5 et d'écart-type σ en abscisse toujours dans cette même famille. Je décrirai plus bas comment je fais le calcul en pratique (avec des points fixes de fonctions génératrices).

L'explication intuitive de l'allure des deux courbes l'une par rapport à l'autre, c'est que dans la courbe rouge (graphe non-orienté), si x infecte y, comme tout est symétrique, il y a une probabilité plus élevée que y ait un degré élevé, donc il infectera à son tour plus de sommets, tandis que dans la courbe bleue (graphe orienté), si x infecte y, alors y a une probabilité plus élevée d'avoir un degré entrant élevé, mais les degrés entrants et sortants ne sont pas corrélés, donc il n'y a pas de raison particulière pour que y infecte plus de sommets à son tour.

Les deux modèles sont très simplistes (encore une fois, une complexité à la fois !), mais j'ai tendance à penser que le modèle orienté est plus réaliste pour décrire les effets de la variance due à des raisons médicales intrinsèques (différences de susceptibilité), tandis que le modèle non-orienté serait plus réaliste pour des effets de la variance due à des raisons sociales (où si on a plus de contacts entrants on a aussi plus de contacts sortants).

La question épineuse, bien sûr, c'est de savoir combien de variance il faut attendre. Comme je l'ai dit, on attend au moins à σ²≥κ (la variance de Poisson), i.e., σ≥√κ, mais combien plus ? Je n'en sais rien ! Au-delà du point de Poisson σ²=κ, il y a un autre point qui me semble assez naturel, c'est le point σ²=κ(κ+1) correspondant à une distribution géométrique. Dans le cas des degrés sortants, la distribution géométrique est assez naturelle : le modèle SIR classique suppose que le temps de rétablissement suit une distribution exponentielle, et on obtient alors une distribution géométrique sur les degrés sortants (alors qu'elle reste poissonnienne sur les degrés entrants ! c'est quelque chose qui m'avait causé une certaine confusion). Pour les degrés entrants je ne vois pas vraiment de justification à une distribution géométrique, mais j'observe que la valeur obtenue (pour la courbe bleue, i.e., dans le cas du graphe orienté) à ce point σ=√(κ(κ+1)), coïncide avec la valeur 1 − 1/R₀ prédite comme seuil d'immunité grégaire (= taux d'attaque au maximum des infectés) dans le modèle SIR classique : je n'ai pas d'explication directe à cette égalité, mais ce n'est sans doute pas une coïncidence, et je me dis qu'il y a probablement une explication faisant intervenir la distribution géométrique. En tout cas, ceci suggère que σ≈κ n'est pas un ordre de grandeur déraisonnable.

Notons qu'il n'y a rien d'absurde ou logiquement contradictoire à ce qu'on ait, pour une variable aléatoire positive comme un nombre de contacts, un écart-type σ supérieur à la moyenne κ : cela peut sembler bizarre de dire que chaque personne a 2.5 ± 5.0 contacts infectieux alors que le nombre de contacts est positifs, mais cela signifie simplement qu'il y a des valeurs très élevées et assez de valeurs petites (y compris 0 : les personnes recevant 0 contacts infectieux ne seront, de toute évidence, pas infectées) pour compenser la moyenne. La distribution sera très asymétrique, mais cela peut arriver.

Mais bien sûr, le fait que ce soit logiquement possible ne signifie pas que ce soit réaliste. Je ne sais pas si c'est le cas. Je ne suis ni sociologue ni médecin, et même ceux qui le sont n'ont pas de données précises sur les contacts, encore moins les contacts infectieux : mesurer R₀ est déjà assez dur, la variance sortante sera sans doute plus dur, et la variance entrante encore plus ! Je dirais qu'on ne doit pas s'attendre à avoir σ nettement supérieur à κ(=R₀), mais je ne saurais pas en dire plus.

Bref, mes courbes doivent simplement être prises comme une indication du fait que la variance dans le nombre de contacts a de l'importance, et comme un ordre de grandeur de l'importance possible de cet effet, mais guère plus. Ce n'est certainement pas une prédiction.

❧

J'en viens maintenant à des explications mathématiquement plus précises (mais toujours loin d'être rigoureuses !) sur le calcul des courbes.

Je rappelle que le principe est de modéliser la propagation de l'épidémie par un graphe aléatoire sur lequel on va regarder la proportion des sommets accessibles depuis un sommet x₀ (permettant effectivement d'accéder à une proportion importante de sommets, i.e., dans la composante géante pour le cas non-orienté). On cherche donc des formules permettant de calculer (l'espérance de) cette proportion pour un modèle de graphe aléatoire un peu crédible permettant de spécifier la distribution des degrés entrants et/ou sortants.

Dans le cas non-orienté, ce modèle de graphe aléatoire est classique : c'est le modèle de la configuration (de Bollobás), ou de Molloy-Reed, décrit ici sur Wikipédia. La construction est assez naturelle : partant de l'ensemble des sommets (dont je suppose évidemment qu'il est très grand, je vais vouloir le faire tendre vers l'infini), on tire au hasard (et indépendamment) le degré de chaque sommet selon la distribution imposée, on crée autant d'« amorces » d'arêtes, autour du sommet en question, que ce degré qu'on vient de tirer, puis on crée les arêtes en tirant aléatoirement de façon répétée (uniformément et indépendamment) deux amorces, qu'on remplace par une arête reliant les deux sommets en question. (Il faut éventuellement décider quelque chose pour écarter le cas des arêtes multiples ou des arêtes reliant un sommet à lui-même, mais ce n'est pas important pour un nombre de sommets tendant vers l'infini alors que la distribution est de moyenne (et de variance ?) finie.) Je pense que c'est une façon vraiment intuitive de construire un graphe dont on veut fixer la distribution des degrés (tellement naturelle que je l'ai non seulement redécouverte, mais redécouverte au moins deux fois, et certainement plein de gens l'ont fait, d'ailleurs le fait qu'elle ait plusieurs noms le suggère). La formule permettant de calculer la taille de la composante géante dans ce modèle est « bien connue » et je vais la décrire plus bas.

Dans le cas orienté, je ne sais pas si le modèle a un nom, mais il est à peu près aussi naturel et évident que ce que je viens de dire. Si on veut juste spécifier la distribution des degrés sortants, on peut faire la chose suivante : on tire au hasard (et indépendamment) le degré sortant de chaque sommet selon la distribution imposée, on crée autant d'« amorces sortantes » d'arêtes, autour du sommet en question, que ce degré qu'on vient de tirer, puis on crée les arêtes en tirant aléatoirement de façon répétée (uniformément et indépendamment) une amorces sortante comme source et un sommet aléatoire comme cible, qu'on remplace par une arête reliant la source vers la cible. La distribution des degrés sortants sera celle spécifiée, tandis que la distribution des degrés entrantes sera poissonnienne. Il y a une formule permettant de calculer la proportion de sommets accessibles dans ce modèle depuis un sommet x₀ menant effectivement à une proportion importante de sommets, et la probabilité que ce soit le cas. Symétriquement, on peut spécifier la distribution des degrés entrants. Pour spécifier à la fois les distributions entrante et sortante (de même moyenne), la généralisation évidente consiste à faire la même chose avec des amorces entrantes et sortantes (il faudra bidouiller quelque chose, par exemple conditionner, pour que le nombre total d'amorces entrantes et sortantes coïncide, mais ça ne doit pas être important quand le nombre de sommets tend vers l'infini) : je pense que, sous des hypothèses raisonnables, la même formule marche encore dans ce cadre, mais je n'ai pas de preuve, juste une intuition par analogie avec les cas que je viens de décrire.

Essayons d'être plus précis. Je rappelle que si p_i est une distribution de probabilité sur ℕ (c'est-à-dire que p₀,p₁,p₂,… sont des nombres réels positifs de somme 1, où p_i est la probabilité ℙ(X=i) que la variable aléatoire X vaille i selon cette loi), la fonction génératrice associée est la série formelle G(z) := ∑_i p_i·zⁱ dont la valeur en z est l'espérance 𝔼(z^X) de z^X. Manifestement, G est positive, strictement croissante et strictement convexe, vérifie G(1)=1, et G′(1) = ∑_i p_i·i est l'espérance 𝔼(X) de la distribution considérée. (On remarquera aussi, à toutes fins utiles, que la donnée de G permet de retrouver les p_i, puisque p_i est la valeur de la dérivée i-ième de G en 0, divisée par i!.)

Je rappelle aussi ce qu'est un processus de Galton-Watson basé sur une distribution de probabilité p_i comme je viens de dire : on construit un arbre aléatoire, en partant de la racine et, pour chaque nœud qui n'a pas encore été traité, en lui donnant (indépendamment de tous les autres) un nombre de fils tiré au hasard selon la distribution spécifiée. La probabilité d'extinction du processus de Galton-Watson est la probabilité que cet arbre soit fini (le processus s'éteint en un nombre fini de générations). Il y a une formule simple pour la calculer : c'est le plus petit point fixe de G, c'est-à-dire le plus petit u tel que G(u)=u (en se rappelant que G(1)=1 et que, par stricte convexité, il ne peut exister qu'au plus deux points fixes, et il est facile de voir qu'il y en a deux si et seulement si G′(1)<1, c'est-à-dire 𝔼(X)<1). Une explication avec les mains est que, si on appelle u la probabilité d'extinction (en partant d'un seul nœud, donc), la probabilité d'extinction en partant de k nœuds est u^k, or un nœud a une probabilité p_k de donner k nœuds, d'où il résulte (en distinguant selon le nombre de fils de la racine) que u = ∑_k p_k·u^k, c'est-à-dire exactement G(u)=u. Ou, si on préfère, et pour voir que c'est le plus petit, autant faire les choses directement à l'endroit : la probabilité d'extinction en r générations est G^r(0) (où G^r = G∘G∘⋯∘G est l'itérée r-ième) comme on le voit en distinguant tous les cas selon le nombre de descendants possibles de chaque nœud impliqué, et cette itération répétée converge bien vers le plus petit point fixe.

Ceci étant rappelé, pour un nombre de sommets tendant vers l'infini, le taux d'attaque (c'est-à-dire la taille relative de la composante-sortante-géante) dans le graphe orienté aléatoire dont la distribution des degrés entrants est donnée par les p_i coïncide avec la probabilité de non-extinction 1−u du processus de Galton-Watson construit sur la distribution p_i en question. Voici une explication intuitive de ce fait : pour déterminer si un sommet aléatoire z est atteint par l'épidémie, on trace en arrière les sommets qui auraient pu le contaminer, c'est-à-dire les voisins entrants de z, puis les voisins entrants de celui-ci, et ainsi de suite à rebours : comme le nombre de sommets est très grand, il s'agit initialement d'un processus de Galton-Watson (on ne retombe pas deux fois sur le même sommet), et la probabilité que ce processus ne s'éteigne pas est la probabilité qu'il atteigne un ensemble géant de sommets, donc le patient zéro. (Symétriquement, la probabilité de non-extinction de l'épidémie coïncide avec la probabilité de non-extinction du processus de Galton-Watson construit sur la distribution des degrés sortants.)

Tout ça est dit de façon extrêmement informel, mais le théorème précis est le théorème 4 de l'article The strong giant in a random digraph de Mathew Penrose, au moins si une des deux distributions impliquées est de Poisson (je conjecture que la généralisation évidente est valable, mais avec mon niveau en probas je n'ai aucun espoir d'arriver à le prouver).

Bref, pour calculer ce taux d'attaque 1−u on a juste à chercher le plus petit point fixe u de la fonction génératrice G de la distribution des degrés entrants.

Qu'en est-il pour le cas non-orienté ? La formule pour la composante connexe géante d'un graphe aléatoire construit selon le modèle de la configuration est très semblable, mais juste un petit peu plus compliquée : on introduit la fonction G₁(z) := G′(z)/G′(1) où G′(z) = ∑_i i·p_i·zⁱ⁻¹ est la dérivée de G ; c'est-à-dire que G₁ est la fonction génératrice associée à la distribution de probabilités i·p_i / (∑_j j·p_j). Si u est le plus petit point fixe de G₁, la taille recherchée est 1−G(u) — du moins si ceci est >0 c'est-à-dire si u<1. Voici une explication intuitive de ce fait : comme dans le cas orienté, pour déterminer si un sommet aléatoire z est atteint par l'épidémie (appartient à la composante connexe géante), on trace les sommets qui auraient pu l'y mettre, c'est-à-dire les voisins de z, puis les voisins de celui-ci, et ainsi de suite : (pour un nombre de sommets très grand) ceci constitue un processus très semblable à un processus de Galton-Watson, mais avec une différence à la génération racine. En effet, dès lors qu'on a suivi une arête pour arriver à un sommet, ce sommet n'est plus un sommet aléatoire mais un sommet aléatoire pondéré par le fait d'avoir une arête qui y mène, c'est-à-dire choisi avec une distribution de probabilités proportionnelle au degré de chaque sommet, et le degré d'un tel sommet est justement i·p_i / (∑_j j·p_j) dont la fonction génératrice est G₁ : on retrouve donc un processus de Galton-Watsons sur G₁, coiffé par G.

Encore une fois, ceci est très informel, et je n'ai pas parlé de l'hypothèse de criticité (si G₁′(1)<1, il n'y a pas de composante connexe géante). Un théorème précis est énoncé dans l'article The Size of the Giant Component of a Random Graph with a Given Degree Sequence de Molloy & Reed (théorème 1), mais le décodage en ce que je viens de dire n'est pas complètement évident. Pour une discussion plus lisible et plus informelle, je renvoie à Random graphs with arbitrary degree distributions and their applications de Newman, Strogatz & Watts (ils donnent aussi des résultats sur les graphes orientés, mais sans preuve ni référence à une preuve, et je n'ai pas l'impression qu'il y ait le cadre que j'ai évoqué ci-dessus).

Bref, pour calculer ce taux d'attaque 1−G(u) on a juste à chercher le plus petit point fixe u de la fonction G₁ dérivée-renormalisée de la fonction génératrice G de la distribution des degrés, et appliquer G.

À ce stade, on peut complètement oublier les graphes et travailler avec les fonctions génératrices :

• La distribution de Poisson d'espérance κ a pour variance σ²=κ, et est donnée par p_i = exp(−κ)·κⁱ/i!, ce qui correspond à la fonction génératrice G(z) = exp(−κ(1−z)), et on a G₁(z) = G(z) = exp(−κ(1−z)). Pour κ≥1, le plus petit point fixe u de G (ou de G₁, du coup), est −W(−κ·exp(−κ))/κ (avec pour W la fonction de Lambert), ce qui donne le taux d'attaque 1−u = 1 + W(−κ·exp(−κ))/κ, qui est le même que prédit par le modèle SIR classique. (Pour κ≤1, le plus petit point fixe de G est 1, le taux d'attaque prédit est 0 : l'épidémie ne peut pas démarrer.)
• La distribution géométrique d'espérance κ a pour variance σ²=κ(κ+1), et est donnée par p_i = κⁱ/(κ+1)⁽ⁱ⁺¹⁾, ce qui correspond à la fonction génératrice G(z) = 1/(κ+1−κz), et on a G₁(z) = 1/(κ+1−κz)². Pour κ≥1, le plus petit point fixe u de G, est 1/κ, ce qui donne le taux d'attaque 1−u = 1 − 1/κ (dans le cas orienté), qui est le même que le seuil d'immunité grégaire prédit par le modèle SIR classique. Comme je l'ai dit plus haut, je ne sais pas s'il y a une explication à ce fait, mais je suppose que ce n'est pas une coïncidence. Pour le cas orienté, le plus petit point fixe de G₁ vaut (κ+2−√(κ(κ+4)))/(2κ), ce qui n'est pas très intéressant, et le taux d'attaque 1−G(u) prédit est (3−√((κ+4)/κ))/2. Tout ça est un cas particulier (r=1) du point suivant.
• La distribution de Pólya (binomiale négative) d'espérance κ et de variance σ²>κ a pour fonction génératrice G(z) = (r/(r+κ−κ·z))^r = (κ/(σ²−(σ²−κ)·z))^{(κ²/(σ²−κ))} où r = κ²/(σ²−κ). On a G₁(z) = (r/(r+κ−κ·z))^r+1 = (κ/(σ²−(σ²−κ)·z))^{(κ²/(σ²−κ))+1}. Il n'y a pas de formule simple particulière pour le plus petit point fixe de G ou G₁, mais il se calcule aisément.
• La distribution binomiale d'espérance κ et de variance σ²<κ a pour fonction génératrice G(z) = ((n−κ+κ·z)/n)ⁿ = ((σ²−(κ−σ²)·z)/κ)^{(κ²/(κ−σ²))} où n = κ²/(κ−σ²). On a G₁(z) = ((n−κ+κ·z)/n)ⁿ⁺¹ = ((σ²−(κ−σ²)·z)/κ)^{(κ²/(κ−σ²))+1}. Il n'y a pas de formule simple particulière pour le plus petit point fixe de G ou G₁, mais il se calcule aisément. (Il y a une asymétrie agaçante, cependant, entre ce cas et le précédent, à savoir que pour la loi binomiale, n a besoin d'être entier, histoire d'obtenir des probabilités positives, alors que pour la loi binomiale négative, r n'a pas besoin de l'être. Donc mes courbes devraient en fait être en pointillés dans le domaine σ²<κ, les pointillés devenant de plus en plus denses en s'approchant du point de Poisson σ²=κ. Mais ce n'est pas très important, de toute façon ce domaine n'est pas très réaliste et est juste inclus pour la complétude des courbes, alors pourquoi ne pas ajouter des probabilités négatives tant qu'à faire. )

Mes courbes sont donc obtenues simplement en faisant varier σ à κ=2.5 fixé, et à chaque fois, en calculant le plus petit point fixe u de G ou G₁ et en traçant 1−G(u). Le source Sage est ici, il est beaucoup plus court que toutes les explications que je viens de donner.