Licence CC 0

Modéliser la propagation du SARS-CoV-2

Combinons un modèle mathématique simple et des données de terrain pour prédire l'évolution de l'épidémie en France

Publié :
Auteur :
Catégories :
Temps de lecture estimé : 13 minutes

Dans ce billet nous allons voir comment élaborer un modèle simple pour étudier la propagation d’un virus dans une population, puis nous utiliserons les données sur la propagation du SARS-CoV-2 pour effectuer une regression et prédire le nombre de cas futurs.

Le but de ce billet n’est pas de donner une évaluation précise mais plutot de présenter des méthodes et outils qui pourront être utilisés ailleurs et que j’utilise moi-même, je ne suis pas du tout un expert des épidémies et je suis donc incapable d’émettre un quelconque avis sur le sujet ^^

Pré-requis: maths niveau L1 pour comprendre toutes les équations.

Modèle mathématique

Croissance exponentielle

Un virus se propage d’individus en individus: plus il y aura de gens infectés, plus il y aura de gens pour transmettre le virus.

Si l’on note P(t)P(t) le nombre d’infectés à un instant tt donné, alors à un instant t+dtt+\mathrm{d}t, le nombre de nouveaux infectés sera proportionel à P(t)P(t). Le coefficient de proportionalité dépendra bien sûr de plusieurs facteurs:

  • La durée d’exposition dt\mathrm{d}t
  • La probabilité de transmettre le virus à quelqu’un lors d’un contact
  • Le nombre de personnes que l’on peut croiser pendant dt\mathrm{d}t
  • Le temps pendant lequel un individu est contagieux
    Pour le moment nous allons confondre les trois derniers paramètres en un seul produit que nous noterons aa (les épidémiologistes parlent de R0R_0).
    Le nombre d’infectés à t+dtt+\mathrm{d}t est donné par la relation suivante: P(t+dt)=nouveaux infecteˊs+anciens infecteˊs=a×P(t)×dt+P(t)P(t+\mathrm{d}t) = \text{nouveaux infectés} + \text{anciens infectés} = a\times P(t) \times \mathrm{d}t + P(t)

D’ou l’on tire l’équation différentielle linéaire de premier ordre:
P=aPP' = aP
Qui a pour solutions:
P(t)=Kexp(at),KRP(t) = K\exp(at), \qquad K\in\mathbb{R}

Une exponentielle ça grandit vite !
Une exponentielle ça grandit vite !

On a donc affaire à une croissance exponentielle sans fin, dans quelques mois l’humanité sera donc complètement infectée o_O ! Ou pas …

Plafonnement logistique

Notre modèle ne colle pas à la réalité, en effet le nombre d’infectés ne peut pas grandir indéfiniment, déjà parce que nous ne sommes qu’un nombre fini de terriens. On pourrait donc supposer que la population mondiale PWP_W est une borne supérieure de PP. Dans l’équation différentielle cela se traduit par le fait que l’on ne peut infecter que des gens qui n’ont pas déjà été infectés. Plus il y a de contaminés, moins il y a de gens à contaminer Le coefficient aa dépend donc de PP et a fortiori du temps:

a(P)=α×Proportion de la population non infecteˊe=PWPPWa(P) = \alpha \times \text{Proportion de la population non infectée} = \frac{P_W-P}{P_W}

On a donc:

P=α×(1PPW)×PP' = \alpha \times \left(1-\frac{P}{P_W}\right) \times P

En réalité PP_{\infty} sera inférieur à PWP_W, ici on supposait un cas idéal (pour le virus, pas pour nous :p ) où l’ensemble de la population est mélangé aléatoirement à chaque étape du procesus, en pratique ce n’est pas le cas certains pays ne jouent pas le jeu et ferment leurs frontières, se confinent bref… des barrières sont mises :pirate:

Posons donc que le nombre final d’infectés est une inconnue que l’on nomme NN. Il nous faut donc résoudre l’équation différentielle de premier ordre:

P=α×(1PN)×PP' = \alpha \times \left(1-\frac{P}{N}\right) \times P

Si vous êtes allergique aux maths, sautez les trois lignes suivantes :D

En séparant les variables on obtient:
dP(1P/N)P=αdt\int \frac{\mathrm{d}P}{(1-P/N)P} = \int \alpha \mathrm{d} t
On décompose en éléments simples à gauche, on intègre les fractions et finalement:

P(t)=Kexp(αt)+K/N,KRP(t) = \frac{K}{\exp(-\alpha t) + K/N}\qquad, \qquad K \in \mathbb{R}

Toute courbe exponentielle cache une courbe logistique
Toute courbe exponentielle cache une courbe logistique

On a donc trois grandes inconnues qui sont KK, NN et α\alpha. Nous allons tenter de déterminer leurs valeurs pour effectuer une prédiction du nombre de cas total en France.

Mais avant intéressons nous sur le sens physique de ces paramètres. C’est bien beau de faire des modèles mais encore faut-il qu’ils aient un sens !

NN a été posé comme étant le nombre de cas limite, et en effet on obtient bien:

limt+P(t)=N\lim\limits_{t \rightarrow +\infty} P(t) = N

aa est notre coefficient de transmission et on peut voir qu’il ait bien homogène à l’inverse du temps comme nous l’avions défini au début.

Et KK alors késako ? à vous de trouver !

KK est homogène à une population comme NN et représente en fait la proportion de gens infectés par rapport aux gens non infectés au début de l’épidémie. En effet on peut vérifier que K=NP0/(NP0)K= NP_0/(N-P_0). Ce qui est un peu flou c’est la définition que l’on prend de début de l’épidémie ici c’est quand t=0t=0, en pratique cela ne correspond à rien de physique, il n’y a pas un moment précis ou l’épidémie commence…

Modèle informatique

Ici le but est donc d’effectuer une regression sur un jeu de données pour tenter de deviner ces paramètres, notre modèle étant non linéaire et non linéarisable on choisit d’utiliser l’algorithme de Levenberg–Marquardt pour minimiser le χ2\chi^2. (code disponible ici )

Je rappelle que le χ2\chi^2 désigne ici une mesure de l’erreur quadratique:

χ2=i(yifp(xi))2\chi^2 = \sum_i (y_i - f_p(x_i))^2

fpf_p représente la fonction que l’on suppose être solution du fit et pp ses paramètres, dans notre pas fpf_p est la fonction logistique et pp est le vecteur (NKα)(N K \alpha)

Si vous ne comprenez pas bien ce qui vient d’être dit, cela n’a aucune importance pour comprendre la suite.

Nous allons commencer par faire un premier fit en nous basant sur les données du site worldometers et nous utiliserons le logiciel gnuplot. L’implémentation de LM y est très puissante et permet de trouver facilement les paramètres de beaucoup de courbes ! Souvenez vous cependant que l’on cherche à minimiser le χ2\chi^2 et que pour cela l’algorithme va chercher à tâtons une solution, il faut donc avoir une bonne idée de la solution pour que cela fonctionne. En effet si vous démarrez avec un pp initial trop éloigné de la réalité, vous risquez de ne jamais converger et de vous perdre dans les limbes des algorithmes de regression … :-°

On obtient alors la courbe suivante: Regression sur les données du 8 mars aux 2 avril
Et les résultats du fit sont:

Final set of parameters            Asymptotic Standard Error
=======================            ==========================
K           = 474.653              +/- 49.27        (10.38%)
N           = 105988               +/- 6271         (5.916%)
a           = 0.177523             +/- 0.00514      (2.895%)

correlation matrix of the fit parameters:
            K      N      a  
K           1.000 
N           0.852  1.000 
a          -0.986 -0.921  1.000                                                                                            

À première vue le fit semble correct et on obtient moins de 10%
d’erreurs pour chaque paramètre.

Mais pourquoi tu ignores les jours avant le 8 mars ?

Parceque qu’il faut bien commencer quelque part ! En traçant les données du tout début de l’épidémie on peut voir qu’elles sont très fluctuantes et semblent donc peu pertinentes. Cela ne se voit pas forcément en échelle linéaire mais en traçant en échelle logarithmique comme on doit toujours le faire d’ailleurs ;)

Échelle logarithmique: la vraie, l'unique !
Échelle logarithmique: la vraie, l'unique !

Ok mais tu te moques de nous ! Tu prédis 106 000 cas et là on est le 13 avril on a dépassé les 130 000 ! Ça marche pas ton truc !

Remarque judicieuse ! Et d’ailleurs si vous y avez prêté attention α\alpha vaut environ 0.18 s1^{-1} alors que dans la réalité les épidémiologistes s’accordent pour dire que sa valeur est supérieur à deux pour une population non confinée. Où est l’arnaque ?

Alors bien sûr je n’ai aucun moyen d’en être sûr mais pour moi le problème vient de la sous-estimation des cas: en effet du à un manque de test en France il est impossible de savoir au jour le jour combien de personnes sont touchées. Nos données sont biaisés, notre résultats l’est donc aussi !

Et d’ailleurs si on s’intéresse aux données après le 2 avril on voit que tout change:

Changement radical de direction
Changement radical de direction
Final set of parameters            Asymptotic Standard Error
=======================            ==========================
K               = 429.453          +/- 83.48        (19.44%)
N               = 170329           +/- 6874         (4.036%)
a               = 0.174692         +/- 0.007079     (4.052%)

La faute à quoi ? Et bien là encore difficile d’être sûr mais apparemment cette augmentation brusque correspond à la prise en compte de nombreux cas comme ceux détectés dans les EHPAD ou dépistés par les médecins généralistes.

Critique du modèle / de la regression

Tentons de prendre un peu de recul pour analyser ce qui représente la réalité et ce qui est à côté de la plaque:

  • Notre fit converge, il y a donc cohérence entre les données de terrain et un modèle purement mathématiques
  • Le R0R_0 est très sous-évalué, probablement du à une mauvaise estimation des cas contaminés. On pourrait essayer de régler le problème en gonflant artificiellement les données (on peut par exemple inférer le nombre de cas à partir du nombde de morts et du taux de mortalité) mais on risque d’introduire de nouveaux paramètres peu ou mal connus donc pas sûr que ça soit beaucoup mieux … Voir par exemple cet article pour en savoir plus sur la difficulté de déterminer le taux de mortalité d’une épidémie en cours.
  • Il en résulte que le NN est probablement lui aussi sous-évalué
  • Notre modèle reste très simpliste par rapport au modèle SIR par exemple qui lui prend en compte l’évolution du nombre de décés, de cas rétablis etc … mais cela est voulu: dans le cas du modèle SIR faire des regression sur un jeu de données en cours devient beaucoup plus compliqué étant donné que l’on a plus de solution explicite. Voir ce lien pour plus de détail.
  • J’ai aussi tenté de créer un modèle qui prenait en compte le confinement pour obtenir deux R0R_0 différents, un avant et un après le confinement mais la différence n’était pas parlante alors je n’en ai pas parlé.

En conclusion: on a vu comment à parti d’outils simples et facilement accessible on peut modéliser une situation physique et effectuer des prédictions plus ou moins réalistes ^^ Ici le problème est que nous n’avons pas de données très fiables pour le moment sur lesquelles baser notre esimation, comme c’est souvent le cas en science: il faut faire avec ce que l’on a ! On a rarement un chemin tout tracé vers une solution merveilleuse … J’espère que ma démarche pourra vous être utile !

3 commentaires

Oui ! Enfin un truc là dessus sur ZdS :)

Par contre, je suis étonné. Ok, certes données pas complètes et pas très précises, mais N170000N\approx 170\, 000, ce qui est au final très peu au regard de la population française. Est ce qu’on doit y voir un effet du confinement malgré tout ?

PS: oui, doublement oui et même quadruplement oui pour l’échelle logarithmique pour le nombre d’infection. Non mais. Par contre, pour le temps, je trouve ça plus douteux.

Édité par pierre_24

#JeSuisToujoursArius • Doctorant et assistant en chimiedev' à temps partiel (co-réalisateur ZEP-12, recherche et template LaTeX)

+0 -0

Ok, certes données pas complètes et pas très précises, mais N170000N\approx 170\, 000, ce qui est au final très peu au regard de la population française. Est ce qu’on doit y voir un effet du confinement malgré tout ?

La magnitude de ce nombre est complètement masquée par le fait que la majorité des cas ne sont même pas testés, par le fait que les données sont biaisés par des reports parfois irréguliers (comme les 20000 cas des EHPAD concentrés en 2 jours), et par le fait que modèle est simpliste (typiquement, la mise en application des mesures barrières puis le confinement affectent la valeur de α\alpha qui devrait alors dépendre du temps, et les fortes restrictions sur les voyages empêchent complètement de pouvoir voir la population comme un seul gros groupe où une loi logistique avec des paramètres uniques a du sens). On peut vraiment pas dire grand chose du NN prévu par ce modèle. Ça reste intéressant en tant que modèle simple d’une épidémie, mais dans le cas réel présent, on peut pas en faire grand chose.

Édité par adri1

I don’t mind that you think slowly, but I do mind that you are publishing faster. — W. Pauli

+0 -0

Merci pour vos retours, en effet comme signalé par adri1 l’incertitude sur le nombre de cas est beaucoup trop grande pour que la valeur de 170000 ait une réelle signification, même les analyses des meilleurs spécialistes sont donnés avec des incertitudes énormes à ce jour.

+1 -0
Connectez-vous pour pouvoir poster un message.
Connexion

Pas encore membre ?

Créez un compte en une minute pour profiter pleinement de toutes les fonctionnalités de Zeste de Savoir. Ici, tout est gratuit et sans publicité.
Créer un compte