Analyse phase space COVID-19





1 Introduction

L’évolution du nombre de décès dus au Covid-19 est analysée et comparée pour 6 pays à la population a peu près comparable : France, Italie, Espagne, Allemagne, Grande Bretagne (GB) ainsi que la province chinoise de Hubei. Il faut préciser d’emblée que la définition du nombre de décès dus au Covid-19 varie d’un pays à l’autre en fonction de ses réglementations sanitaires et administratives. Il le ne faut donc pas à chercher à comparer les dynamiques pour porter un jugement sur l’efficacité des mesures prises dans un pays. En revanche, suivre l’évolution du nombre de décès au cours du temps peut apporter des renseignements sur la dynamique de l’épidémie, et permet – éventuellement – de faire des prévisions pour son futur.

Notre hypothèse est qu’au-delà des différences propres à chaque pays, il existerait un modèle commun ne dépendant que d’un faible nombre de paramètres, capable de décrire la dynamique de la première vague de l’épidémie. Nous comparons deux approches. La première consiste à modéliser la dynamique des décès par une courbe logistique. Même si il est sans doute simpliste, le modèle logistique semble être suffisamment robuste et général pour modéliser la dynamique de l’épidémie dans sa phase aiguë, et cela dans les 6 pays étudiés. Dans une deuxième partie, nous modélisons la dynamique à l’aide d’un système d’équations différentielles (EDO) classique en épidémiologie, le modèle SIR. Celui-ci propose une modélisation de la dynamique de l’épidémie en considérant trois classes dans la population; les susceptibles, les infectés et les immunisés. Les deux approches sont liées, et ces liens sont décrits.

Nous montrons comment réaliser l’estimation des paramètres pour ces deux approches, et nous montrons qu’elles produisent toutes deux des ajustements crédibles. Ils permettent également de faire une estimation du nombre de décès à l’issue de la phase la plus aigüe. En revanche, on s’attend à ce que ce modèle ne soit pas nécessairement opérant dans la période qui va s’ouvrir à l’issue du confinement initial.

2 Les données

Les données proviennent du dépôt Github maintenu par le CSSE de Johns Hopkins University, à l’URL CSSE Johns Hopkins. La définition exacte d’un décès dû au Covid-19 dépend d’un pays à l’autre, ce qui fait que la comparaison directe entre les chiffres bruts des différents pays n’a pas beaucoup de sens.

En France, la définition du nombre de décès a changé le 1er avril 2020. Avant cette date, il fallait que le décès soit constaté en hôpital et que le test Covid-19 soit positif; à partir du premier avril les décès en EHPAD (testés positifs) ont été également comptabilisés. Johns Hopkins a fait le choix de refléter ce changement de définition. Les données françaises complètes, qui différencient les deux comptages, sont téléchargées à partir du dépôt maintenu par Santé Publique France qui distingue les différents types de décès, à l’URL OpenCovid19. Nous ferons deux analyses: la première concerne les décès à l’hôpital; la second utilisera le nombre de décès total. Dans ce cas, les données antérieures au 1er avril seront corrigées par un facteur multiplicatif ajusté sur les données acquises après cette date.

3 Portrait de phase par un modèle logistique

3.1 La fonction logistique

La fonction logistique s’écrit \(x(t) = K/(1 + \exp( -r(t-t_0)) )\), où \(t\) est une date (un jour dans l’année), et \(x(t)\) est le nombre de décès cumulés à la date \(t\). Ce modèle a trois paramètres: \(K\) est un paramètre de volume (le nombre total de décès lorsque \(t \to \infty\)), \(r\) est un paramètre de vitesse et \(t_0\) le paramètre de calage temporel qui correspond à la date où \(K/2\) décès ont été enregistrés.

On note \(y(t) = x'(t)\) la dérivée première de \(x(t)\) par rapport au temps, qui représente le nombre de décès enregistrés au jour \(t\). Il est assez aisé de montrer que \(x\) et \(y\) sont reliés par l’équation \[y = xr(1-x/K).\] Cette parabole s’annule lorsque \(x=0\) et \(x=K\). Elle atteint son maximum en \(x=K/2\) et vaut alors \(Kr/4\). Il faut également noter que cette parabole ne dépend pas de \(t_0\), c’est à dire qu’elle ne dépend pas de la date à laquelle démarre l’épidémie. Cette représentation ne nécessite donc pas de recalage des données à une date initiale commune.

3.2 Estimation des paramètres de la courbe logistique

3.2.1 Estimation dans l’espace des phases

Les paramètres \(r\) et \(K\) sont estimés par ajustement aux moindres carrés dans la représentation parabolique, que les physiciens appellent également l’espace des phases. Pour un pays \(i\) donné, on note \(Y_i(t)\) les décès journaliers au jour \(t\) et \(X_i(t) = \sum_{k\leq t} Y_i(k)\) les décès cumulés à cette date. On représente ci-dessous ces données ainsi que les courbes \(y = xr(1-x/K)\) obtenues avec les valeurs estimées de \(r\) et \(K\).

Il est intéressant de noter qu’aux fluctuations journalières près, le comportement parabolique semble se confirmer pour l’ensemble des pays, y compris pour la province de Hubei qui est le seul ensemble territorial analysé à avoir parcouru la totalité de la dynamique épidémique. En particulier, il est tout à fait remarquable que l’ajustement soit également d’assez bonne qualité sur la phase finale de la dynamique. Les fluctuations importantes au centre de la courbe sont sans doute dues à la difficulté de maintenir un comptage précis au pic de l’épidémie1.

Paramètres estimés

On obtient les paramètres suivants:

##   F. hosp.   France   Italie  Espagne Allemagne       GB   Hubei
## r     0.16     0.16     0.14     0.18      0.13     0.20    0.15
## K 13279.84 26320.55 24579.74 20772.90   7905.95 17570.95 3186.75

3.2.2 Vérification sur la dynamique temporelle

On vérifie l’ajustement précédent en ré-introduisant la dimension temporelle. Dans une dynamique décrite par un modèle logistique, il existe un lien linéaire entre le temps et le logarithme du ratio entre l’avancée de l’épidémie et ce qui lui reste à parcourir. En termes mathématiques \[\ln \frac{x(t)}{K-x(t)} = r(t-t_0).\] Cette relation offre un moyen simple de contrôler la qualité de l’ajustement logistique. La représentation \[\ln \frac{X}{K-X(t)} \quad \rm{vs.} \quad t\] doit être linéaire et \(t_0\) peut se déduire à partir de l’intercept d’un ajustement linéaire réalisé par moindres carrés.

La Figure suivante montre la grandeur \(\ln \frac{X_i}{K-X(t)}\) en fonction du temps \(t\), compté à partir de la date où on a observé 10 décès cumulés. Globalement, on observe un ajustement assez bon, et même excellent dans certain cas (France hospitalière, Grande-Breagne). Pour la France, l’Italie, l’Espagne et dans une mondre mesure pour l’Allemagne, on observe que la courbe empirique s’infléchit, indiquant un ralentissement de la dynamique logistique. Mathématiquement, cela correspond à un coefficient \(r\) qui décroit au cours du temps. Cela peut être du à un effet du confinement, ou encore à une capacité accrue du système de santé, ou les deux. Les données de la province de Hubei montrent un comportement assez étrange avec une dynamique accélérée au début et à la fin de la période.

3.3 Résultats de l’analyse logistique

L’ajustement linéaire présenté à la Section précédente permet d’estimer le paramètre \(t_0\) qui représente la date du pic épidémique. Maintenant que tous les paramètres ont été estimés, on peut faire une représentation graphique du nombre de décès cumulés et journaliers.

Pour la suite, nous définissons la fin de la phase épidémique aigüe le jour où sera atteint un rythme de décès journaliers inférieur à 100 (c’est à dire moins de 1 mort par département en moyenne). D’autres valeurs pourraient bien entendu être choisies.

  Début Max <100
F. hosp. 2020-03-08 2020-04-08 2020-04-28
France 2020-03-06 2020-04-09 2020-05-04
Italie 2020-02-26 2020-04-03 2020-04-30
Espagne 2020-03-08 2020-04-05 2020-04-27
Allemagne 2020-03-16 2020-04-15 2020-05-02
GB 2020-03-15 2020-04-12 2020-05-01
Hubei 2020-01-23 2020-02-20 2020-02-27

Le tableau ci-dessus indique que le pic est atteint entre le 3 et le 15 avril dans les pays Européen étudiés, une sortie de la phase la plus aigüe (conventionnellement définie par des décès journaliers inférieurs à 100) se faisant à partir du 27 avril (en Espagne) et jusqu’au 4 mai (pour France-total).

4 Portraits de phase via des modèles épidémiologiques de type SIR-D

4.1 Courbe logistique et dynamique du nombre d’infectés

Supposons que le nombre de morts \(D(t)\) dû à l’épidémie de COVID-19 suive une loi logistique, à partir d’un temps \(t=0\). Cela signifie qu’il existe \(r, \, K>0\) tels que \[\begin{equation} \label{eq:Dprime} D'(t)= r\, D(t) \, \left( 1 -\frac{D(t)}{K}\right), \end{equation}\] avec \(D(0)=D_0>0\). Soit \(I(t)\) le nombre d’infectés au temps \(t\). Une part des infectés devient immunisée, avec un taux \(\beta\), et une autre part décède, avec un taux \(\gamma\) (supposé constant). On a donc : \[\begin{equation} \label{eq:D_I} D'(t)=\gamma \, I(t). \end{equation}\] Posons \[G(t):=\int_0^t I(s) ds.\]En intégrant \(D'(t)=\gamma \, I(t)\) entre \(0\) et \(t\), on obtient \(D(t)=\gamma \, G(t)+D_0\), puis \[ \gamma \, I(t)=r \, (\gamma \, G+D_0) -\frac{r}{K} \, (\gamma \, G+D_0)^2, \] ou encore: \[\begin{equation} \label{eq:G1} \gamma \, G'(t)=r \, (\gamma \, G+D_0) -\frac{r}{K} \, (\gamma \, G+D_0)^2. \end{equation}\] Cette équation peut se résoudre explicitement, ce qui conduit à: \[\begin{equation} I(t)=I_0\frac{r^2}{\left( (r-\gamma I_0/D_0)e^{r\, t/2}+(\gamma I_0/D_0)e^{-r\, t/2}\right) ^2}. \end{equation}\]

La trajectoire logistique de \(D(t)\) correspond donc à une croissance des infectés, suivie d’une décroissance vers \(0\). La décroissance du nombre d’infectés débute quand \(D''(t)\) s’annule, donc au point d’inflexion de \(D(t)\) qui correspond au temps \(t\) tel que \(D(t)=K/2\).

4.2 Un modèle SIR-D

Les modèles SIR sont les modèles d’EDO les plus classiques en épidémiologie. Ce sont des modèles dits compartimentaux, qui divisent la population en plusieurs classes. Classiquement, on distingue les trois classes suivantes: les susceptibles, les infectés et les résistants (immunisés, recovered ou removed en anglais), d’où le nom de modèle SIR. Un compartiment correspondant au nombre de décès peut être ajouté, conduisant à un modèle SIR-D: \[\begin{equation} \left\{\begin{array}{l} S'(t)=- \frac{\alpha}{N} \, S(t) \, I(t),\\ I'(t)=\frac{\alpha}{N} \, S(t) \, I(t) - (\beta+\gamma) \, I(t), \\ R'(t)= \beta \, I(t),\\ D'(t)= \gamma \, I(t), \end{array}\right. \end{equation}\] avec \(S(t)\) le nombre de susceptibles, \(I(t)\) le nombre d’infectés, \(R(t)\) le nombre d’immunisés et \(D(t)\) le nombre de morts à l’instant \(t\). Pour simplifier, la population totale \(N\) est supposée constante: on néglige l’impact de la mortalité sur la taille totale de la population. La donnée initiale correspond à une certaine population de susceptibles \(S_0\), \(I_0=1\) (1 infecté), \(R_0=0\) et \(D_0=0\).

Comme dans Roques et al. (2020), nous supposons que la durée de contagion moyenne est de 10 jours~: \(\beta=1/10\). Nous utilisons les valeurs du taux de contact \(\alpha\) et du taux \(\gamma\) estimés dans Roques et al. (2020) à partir des données françaises (jusqu’au 17 mars), soit \(\alpha=0.32\) et \(\gamma=8 \, 10^{-4}\) (médianes des distributions a posteriori des paramètres), correspondant à un taux de létalité \(\gamma/(\gamma+\beta)\) de \(0.8\%\).

4.3 Estimation dans l’espace des phases

Comme ci-dessus, nous faisons un ajustement aux moindres carrés dans l’espace des phases \((D,D')=(D,\gamma \, I)\) de l’unique paramètre inconnu: le nombre initial de susceptibles \(S_0\). Ici, cette population \(S_0\) ne correspond pas à la population totale, mais aux `non confinés’ ou futurs infectés. Les résultats sont présentés sur les Figures ci-dessous. On note un portrait de phase comparable avec ceux obtenus avec le modèle logistique, bien cette fois-ci que légèrement asymétrique.