TL;DR : Ne lisez pas ce billet, lisez plutôt celui ci : SEAIR et cette page générale du écrite par le groupe de modélisation de l’équipe ETE (Laboratoire MIVEGEC, CNRS, IRD, Université de Montpellier).

Dès qu’il y a des données, l’envie de comprendre et représenter fait jour. D’où ce travail de modélisation.

L’idée m’est venue en suivant le travail de modélisation d’Yves Peysson, physicien au CEA. Je me suis inspiré du travail de Kai Sasaki, ainsi que de ce billet. Attention, ô lecteur: ce modèle naïf présente à coup sûr de nombreuses limites et d’éventuelles horreurs qu’un véritable épidémiologiste n’aurait jamais commises. Moi, si. Il y a une fortes variation des prédictions, le modèle est très sensible aux paramètres. C’est une limite. L’exercice est cependant intéressant à mes yeux, même s’il conduit à un nombre de décès assez affolant.

Le code est disponible ici. Contact : thomas # alt-gr.tech.

[update 02/04/2020]
Un solveur basé sur l’algorithme de Basin-hopping a été essayé, afin de mieux explorer l’espace paramétrique, de dimension cinq : deux beta, un gamma, un temps d’installation  du confinement et le taux de décès. Il permet d’effectuer des sauts non contraints dans l’espace des paramètres, et d’être plus exhaustif dans sa couverture. La partie d’optimisation locale est confiée à un des algorithmes de la routine minimize avec contraintes (L-BFGS-B, TNC, SLSQP).

[update 06/04/2020]
Le modèle ne permet pas de retrouver les valeurs observées de manière satisfaisante. Un outil de recherche exhaustif permet d’explorer l’espace des solutions de manière systématique (optimize.brute), ce qui autorise à être beaucoup plus large, notamment dans la date de début de l’épidémie. Ces études sont chronophages en temps machine et en temps humain, je les mets donc en suspens.


Modèle

Modèle compartimenté SIR+D avec paramètres libres (β, γ, taux de décès).
Modèle amélioré SEIRD avec taux d’incubation pour le passage exposed vers infected.
Prise en compte du confinement de la population par un paramétrage temporel de β par une sigmoïde avec deux valeurs pour β.
Autres paramètres : date de début de l’épidémie, taux de contamination général (~30% ou ~60% de la population), pourcentage des cas graves, dynamique temporelle d’installation du régime de confinement.
Nombre total de paramètres : neuf.

Résolution numérique du système d’équations linéaire (solve_ivp) et optimisation des paramètres.
Full python : numpy, scipy, pandas.
Optimisation numérique (bibliothèque optimize de scipy) :
•  routine minimize avec un des algorithme d’optimisation sous contrainte ;
•  solveur basé sur l’algorithme de Basin-hopping, sorte de recuit simulé, afin de mieux explorer l’espace paramétrique.
•  solveur idiot utilisant la force brute (optimize.brute), histoire de.

Observations utilisées : uniquement les nombres de décès, les cas graves (=hospitalisations constatées) sont jugés trop peu fiables.

Source : ECDC
Archive du code : source

[résultats partiels présentés ci-dessous à partir des données du 27 mars]


France

Hypothèse de population touchée : 66%

Nombre de décès total : 692k.
Pics prévus :
• sur cas graves : 01/05/2020, journée 76, avec ~52k cas et ~15k décès journaliers ;
• sur décès : 08/05/2020, journée 83, avec 43k cas et ~ 17k décès.

Les paramètres exhibés par le processus d’optimisation ne semblent pas trop déconnants :
• beta : 0.42/jour/personne avant confinement, 0.24 après ;
• gamma : 0.10/jour, soit une durée de maladie de 9.5 jours ;
• R0 : 4.1 avant confinement, 2.3 après.
R0 représente le taux de contamination par personne infectée.

Résultats France : nouveaux cas à gauche, prédictions générales cumulées à droite.

 

Résultats France : zoom observations vs prédictions.


Italie

Hypothèse de population touchée : 66%


Nombre de décès total : 358k.
Pics prévus :
• sur cas graves : 14/04/2020, journée 73, avec 53k cas et ~7.9k décès journaliers ;
• sur décès : 21/04/2020, journée 96, avec 41k cas et ~9.2k décès.

Les paramètres exhibés par le processus d’optimisation diffèrent du cas français, et semblent moins bons :
beta : 0.41/jour/personne avant confinement, 0.25 après ;
• gamma : 0.096/jour, soit une durée de maladie de 10.4 jours ;
• R0 : 4.3 avant confinement, 2.6 après.

Résultats Italie : nouveaux cas à gauche, prédictions générales cumulées à droite.

 


Résultats France : zoom observations vs prédictions.


Discussion

La forme globale des courbes semble correcte, avec une tendance nette à privilégier une épidémie de forte ampleur . Les valeurs prédites sont assez sensibles aux conditions de simulation, essentiellement sur la pondération de l’alignement du modèle avec les observations de décès et cas graves. L’effet du confinement est visible, mais son impact est encore limité et peu pris en compte. Il est probable que le nombre de décès soit fortement sur-évalué, avec une prédominance de la dynamique d’immunité collective (approche #yolo).

Les données sont problématiques. En l’état, je fais principalement confiance aux observations de décès (80% du forçage). Les cas graves me semblent globalement sous-estimés. Le modèle fixe un taux de décès (cas / décès) autour d’un 1%. L’estimation est barres d’erreurs sur les prédictions n’a pas été faite : je les soupçonne d’être considérables.

Améliorations du modèle

Beaucoup de choses sont à améliorer. Le paramétrage est fait principalement sur les observations de décès (80% du poids) et sur les cas graves (20%). Les cas graves sont estimés environ à 15% de l’ensemble des cas, avec un forçage total. Il faudrait relâcher cette contrainte en laissant le paramètre évoluer dans un intervalle 1%-15% par exemple, et ne se baser que sur les décès. Sous-pondérer les observations en début d’épidémie, quand la maladie était mal diagnostiquée, est aussi une option.

Est-il possible de prendre en compte plus finement la composition de la population avec un modèle SIR+, ainsi que de nombreux facteurs décrivant les co-morbidités observées … ? Pas sûr. Des médecins m’ont proposé beaucoup de perfectionnement sur ce sujet, que je doute de savoir (pouvoir ?) prendre en compte.

Concernant le décompte des décès, on pourrait largement alourdir le bilan en ajoutant les cas graves qui sont au delà de la capacité des lits de réanimation. Cela revient à ajouter l’aire entre la courbe rouge et la droite = ~8000 pour la France, ~5000 pour l’Italie, ~28000 en Allemagne. Il est d’ailleurs important de prendre ce paramètre en compte dès que les hôpitaux entrent en saturation : ca se verra sur les courbes de décès, qui augmenteront brutalement.

Limites

Le modèle SIR avec régime de confinement est très difficile à faire coller aux observations du terrain.

Aux paramètres classique du modèle (taux de contagion, durée typique de contamination, taux de décès) s’ajoutent les incertitudes majeures sur la proportion d’asymptomatiques ainsi que la fraction de la population qui sera affecté. Trop de paramètres libres. En l’état, mes simulations prédisent des pics courant Avril au mieux, et ensuite des plateaux concernant les nombres de décès. Difficile de trouver des paramètres qui reproduisent le pic actuel qu’on observe en Italie en ce jour. C’est une limite forte, et j’attends donc de disposer de plus de données pour continuer ce travail.

Les quelques publications que j’ai pu lire (ici, par exemple) font appel à des modèles en partie bayésien afin de pouvoir estimer les incertitudes sur les paramètres d’entrée du modèle, dont le R0.

Remarques personnelles

Ce travail tout à fait modeste et sans aucune prétention m’a donné l’occasion de repartir faire de la simulation numérique avec Python. Les bibliothèques sont très fournies et simples d’usage. Intellectuellement, ce travail m’a donné l’occasion de découvrir le monde compliqué de la modélisation épidémiologique et ses incertitudes expérimentales. Pas mal de ressemblances avec mon expérience passée de physicien (2005-2011). Une forme de frustration de ne pas connaître la fin de l’histoire de devoir attendre la fin de l’épidémie pour pouvoir estimer correctement les paramètres.