Biology 2020, 9(5), 97
Lionel Roques, Etienne K. Klein, Julien Papaïx, Antoine Sar, Samuel Soubeyrand
L'épidémie de COVID-19 a commencé en décembre 2019 dans la province du Hubei, en Chine. Depuis lors, la maladie s'est propagée dans le monde entier pour atteindre le stade de la pandémie, selon l'OMS [1], le 11 mars. Les premiers cas ont été détectés en France le 24 janvier. Le taux de mortalité par infection (IFR), défini comme le nombre de décès divisé par le nombre de cas infectés, est une quantité importante qui nous renseigne sur le nombre attendu de victimes à la fin d'une épidémie, lorsqu'une proportion donnée de la population a été infecté. Bien que les données sur le nombre de décès dus au COVID-19 soient probablement exactes, le nombre réel de personnes infectées dans la population n'est pas connu. Ainsi, en raison du nombre relativement faible de tests de dépistage qui ont été réalisés en France (environ cinq personnes sur 10000 en France contre 50 sur 10000 en Corée du Sud jusqu'au 15 mars 2020; sources: Santé Publique France et Korean Centre pour Disease Control), le calcul direct de l'IFR n'est pas possible. Sur la base des cas confirmés par PCR chez des résidents internationaux rapatriés de Chine en janvier 2020, Verity et al. [2] ont obtenu une estimation du taux de mortalité par infection (IFR) de 0,66% en Chine et, en ajustant les taux d'attaque non uniformes par âge, un IFR de 0,9% a été obtenu au Royaume-Uni [3]. En utilisant les données du bateau de croisière Diamond Princess mis en quarantaine au Japon et en corrigeant les retards entre la confirmation et le décès, Russel et al. [4] ont obtenu un IFR de 1,3%.
En utilisant les premières données (jusqu'au 17 mars) disponibles en France, nos objectifs sont: (1) de calculer l'IFR en France, (2) d'estimer le nombre de personnes infectées par le COVID-19 en France, et (3) de calculer un nombre de reproduction de base R0.
Nous avons obtenu le nombre de cas positifs et de décès en France, jour après jour, auprès du Johns Hopkins University Center for Systems Science and Engineering [5]. Les données sur le nombre de tests effectués ont été obtenues auprès de Santé Publique France [6]. Certaines données (cas positifs, décès, nombre de tests) n'étant pas totalement fiables (exemple: 0 nouveau cas détecté en France le 12 mars 2020), nous avons lissé les données avec une moyenne mobile sur 5 jours. Les données officielles sur le nombre de décès par COVID-19 en France ne prennent en compte que les personnes hospitalisées. Environ 728 000 personnes en France vivent dans des maisons de retraite médicalisées (EHPAD, source: DREES [7]). Des données récentes de la région Grand Est (source: Agence Régionale de Santé Grand Est [8]), rapportent un total de 570 décès dans ces maisons de retraite, qu'il faut ajouter au décompte officiel (1015 décès au 31 mars).
Le formalisme mécaniste-statistique, qui devient un standard en écologie [9 , 10 , 11] permet à l'analyste de coupler un modèle mécaniste décrivant une variable latente, ici un modèle d'équation différentielle ordinaire (EDO) de type SIR, et des données incertaines non exhaustives. Pour combler le fossé entre le modèle mécaniste et les données, l'approche utilise un modèle probabiliste décrivant le processus de collecte de données. Une méthode statistique est ensuite utilisée pour l'estimation des paramètres du modèle mécaniste.
Modèle mécaniste. La dynamique de l'épidémie est décrite par le modèle compartimental SIR suivant: \begin{equation}\tag{1} S'(t)=−\frac{\alpha}{N} S(t)I(t),\quad I'(t)=\frac{\alpha}{N}S(t)I(t)−\beta I(t),\quad R'(t)=\beta I(t), \end{equation} avec S la population sensible, I la population infectée, R la population rétablie (individus immunisés) et N=S+I+R la population totale, censée être constante. Le paramètre α est le taux d'infection (à estimer) et 1/β est le temps moyen jusqu'à ce qu'un infecté se rétablisse. Sur la base des résultats de [12], la période médiane d'excrétion virale est de 20 jours, mais l'infectiosité a tendance à se disparaître avant la fin de cette période: les résultats de [13] montrent que l'infectiosité commence à partir de 2,5 jours avant l'apparition des symptômes et diminue dans les 7 jours suivant le début de la maladie. Sur la base de ces observations, nous supposons ici que1/β=10 journées.
Les conditions initiales sont S(t0)=N-1, I(t0)=1 et R(t0)=0, où N=67 106 correspond à la taille de la population. Le modèle SIR est démarré à un moment donné t=t0, qui sera estimée et devrait se rapprocher de la date d'introduction du virus en France (ce point est brièvement discuté à la fin de cet article). Le système EDO (1) est résolu grâce à un algorithme numérique standard, utilisant le solveur Matlab ® ode45 .
Ensuite, nous désignons par D(t) le nombre de décès dus à l'épidémie. Notez que l'impact du compartiment D(t) sur la dynamique du système SIR et sur la population totale est ici négligée. La dynamique de D(t) dépend de I(t) par l'équation différentielle: \begin{equation}\tag{2} D'(t)=\gamma I(t), \end{equation} avec γ(t) le taux de mortalité des personnes infectées.
Modèle d'observation. On suppose que le nombre de cas testés positifs au jour t, noté \(\hat{\delta}_t\), suit des lois binomiales indépendantes, conditionnellement au nombre de tests nt effectuée le jour t et à pt la probabilité d'être testé positif dans cet échantillon: \begin{equation}\tag{3} \hat{\delta}_t \sim Bi (n_t,p_t) . \end{equation} La population testée comprend une fraction des cas infectés et une fraction des susceptibles: nt = τ1(t) I(t) + τ2(t) S(t). Donc, \begin{equation}\tag{4} p_t=\frac{σ \tau_1(t) I(t)}{\tau_1(t) I(t) +\tau_2(t) S(t)}=\frac{σ I(t)}{I(t) +\kappa_t S(t)} \end{equation} avec κt:=τ2(t)/τ1(t), la probabilité relative de subir un test de dépistage pour un individu de type S par rapport à un individu de type I (probabilité d'être testé conditionnellement à être S / probabilité d'être testé conditionnellement à être I). Nous supposons que le ratio κ ne dépend pas de t au début de l'épidémie (c'est-à-dire sur la période que nous utilisons pour estimer les paramètres du modèle). Le coefficient σ correspond à la sensibilité du test. Dans la plupart des cas, des tests RT-PCR ont été utilisés et les données existantes indiquent que la sensibilité de ce test utilisant des écouvillons pharyngés et nasaux est d'environ 63 à 72% [14]. Nous prenons ici σ= 0,7 (70% de sensibilité).
Les paramètres inconnus sont α, t0 et κ. Le paramètre γ(t) est calculé indirectement, en utilisant la valeur estimée de I(t), les données sur D(t) (supposé exact) et la relation (2). La probabilité L est définie comme la probabilité des observations (ici, les incréments {\(\hat{\delta}_t\)}) conditionnellement aux paramètres. En utilisant le modèle d'observation (3), et en supposant que les incréments \(\hat{\delta}_t\) sont indépendants conditionnellement au processus SIR sous-jacent et que le nombre de tests nt est connu, on obtient: \begin{equation}\tag{5} L (\alpha , t_0, \kappa) := P ({\hat{\delta}_t} | \alpha , t_0, \kappa) = \prod_{t =t_i}^{t_f} \frac{n_t!}{(\hat{\delta}_t)! (n_t-\hat{\delta}_t)!} p_t^{\hat{\delta}_t} (1 -p_t)^{n_t-\hat{\delta}_t}, \end{equation} avec ti la date de la première observation et tf la date de la dernière observation. Dans cette expression, L(α ,t0, κ) dépend de α ,t0, κ à travers pt.
L'estimateur du maximum de vraisemblance (MLE, c'est-à-dire les paramètres qui maximisent L), est calculé à l'aide de l'algorithme de minimisation contraint BFGS, appliqué à -ln(L), via la fonction Matlab ® fmincon . Afin de trouver un maximum global de L, nous appliquons cette méthode à partir de valeurs initiales aléatoires pour α ,t0, κ choisies uniformément dans les intervalles suivants: α ∈ (0 , 1), t0∈ (1 , 31), (1er-31 janvier) et κ ∈ (0 , 1). L'algorithme de minimisation est appliqué à 10 000 valeurs initiales aléatoires des paramètres.
La distribution postérieure des paramètres (α ,t0, κ) est calculée avec une méthode bayésienne, en utilisant des distributions a priori uniformes dans les intervalles donnés ci-dessus (une distribution a priori plus informative a également été testée, voir l' annexe A). Cette distribution postérieure correspond à la distribution des paramètres conditionnellement aux observations: \begin{equation}\tag{6} P (\alpha , t_0, \kappa | {\hat{\delta}_t} ) = \frac{L (\alpha , t_0, \kappa) \pi(\alpha , t_0, \kappa)}{ C}, \end{equation} où π(α ,t0, κ) correspond à la distribution a priori des paramètres (donc uniforme) et C est une constante de normalisation indépendante des paramètres. Le calcul numérique de la distribution postérieure est effectué avec un algorithme Metropolis – Hastings (MCMC), en utilisant quatre chaînes indépendantes, chacune avec 106 itérations, à partir de valeurs aléatoires proches de la MLE (seule la seconde moitié des itérations est utilisée pour générer le postérieur). Les codes Matlab ® sont disponibles en tant que matériaux supplémentaires .
Les données \(\hat{\delta}_t\) utilisés pour calculer la MLE et la distribution postérieure sont ceux correspondant à la période du 29 février au 17 mars.
De plus, la distribution du nombre cumulé de cas infectés (I(t) + R(t)) à travers le temps est présentée dans la figure 2. On observe qu'il est beaucoup plus élevé que le nombre total de cas observés (comparer avec la figure 1). Le rapport moyen estimé entre le nombre réel d'individus infectés et les cas observés(I(t) + R(t)) /Σt est de huit (IC à 95%: 5–12) sur la période considérée.
Nombre de reproduction de base. Avec les systèmes SIR de la forme (1), le nombre de reproduction de base R0 peut être calculé directement, sur la base de la formule R0= α / β [15]. Quand R0< 1, l'épidémie ne peut pas se propager dans la population. Quand R0> 1, le compartiment infecté I augmente tant que R0S > N = S+I+R. Nous avons calculé la distribution postérieure marginale du nombre de reproduction de base R0. Cela conduit à une valeur moyenne de R0 de 3,2 (IC à 95%: 3,1–3,3). La distribution complète est disponible dans l' annexe A (figure A4).
Sensibilité des résultats par rapport aux paramètres fixes du modèle. Nous avons calculé la MLE avec une période infectieuse plus longue (1 / β) de 20 jours estimé par [12]. Cela conduit à un nombre reproductif de base beaucoup plus grand R0=4,8 et un facteur ×15 entre les cas signalés et le nombre réel de cas. Cependant, la valeur de l'IFR reste inchangée (0,5%). Nous avons également vérifié si la largeur de la fenêtre du lissage (moyenne mobile sur 5 jours) avait un impact sur nos résultats. Des calculs du MLE avec une largeur de fenêtre de 3 jours (et β= 1 / 10) conduit aux mêmes résultats que ceux présentés ci-dessus, à savoir une R0= 3,2 et un IFR de 0,5%.
Sur l'IFR et le nombre de cas infectés. Le nombre réel d'individus infectés en France est probablement beaucoup plus élevé que les observations (on retrouve ici un facteur × 8), ce qui conduit à un taux de mortalité inférieur à celui calculé sur la base des cas observés: on retrouve ici un IFR de 0,5 % basé sur les données de comptage des décès à l'hôpital, à comparer avec un taux de létalité (TFC, nombre de décès sur nombre de cas diagnostiqués) de 2% le 17 mars. En ajustant le nombre de décès dans les maisons de soins infirmiers, nous avons obtenu un IFR de 0,8%. Ces valeurs pour l'IFR sont cohérentes avec les résultats de [2] (0,66% en Chine) et [3] (0,9% au Royaume-Uni). La valeur de 1,3% estimée sur le bateau de croisière Diamond Princess [4] se situe au-dessus de l'extrémité supérieure de notre IC à 95%. Cela reflète la répartition par âge du navire, qui était biaisée vers les individus plus âgés (âge moyen: 58 ans), parmi lesquels l'IFR est plus élevé [3 , 4].
L'objectif de notre étude était d'estimer l'IFR sur la base des premières données, avant que des enquêtes à grande échelle ne deviennent disponibles. Fin avril, de nouvelles données et études préliminaires sont disponibles et peuvent être comparées à nos résultats. Une étude sur les anticorps à New York publiée le 24 avril 2020 montre qu'environ 14% ont été testés positifs, correspondant à 2,7 millions de cas, à comparer avec les 271000 cas confirmés et un total de 15500 décès dans tout l'État. Cela correspond à un IFR de 0,6%. En France, une autre étude préliminaire menée par l'Institut Pasteur [16], et sur la base d'une estimation conjointe des données françaises jusqu'au 14 avril (données de comptage des décès dans les hôpitaux) et des données des navires de croisière Diamond Princess trouve un IFR de 0,5% confirmant ainsi notre résultat. Au 28 avril, le nombre de décès dus au COVID-19 en Lombardie (Italie) est de 13575 (source: Ministero della Salute), pour une population de 10 millions de personnes, ce qui montre que l'IFR est d'au moins 0,14%. En revanche, en Corée du Sud où le nombre de cas détectés a rapidement atteint un plateau, suggérant une faible proportion de cas non détectés, le rapport entre le nombre de décès et le nombre de cas positifs est de 244/10 752 ≈ 2,3% (source: Johns Hopkins University Center for Systems Science and Engineering [5]), qui peut être considérée comme une limite supérieure de l'IFR, bien que surestimée.
Si le virus conduisait à contaminer 80% de la population française [3], le nombre total de décès à déplorer en l'absence de variation du taux de mortalité (augmentation induite par exemple par la saturation des structures hospitalières, ou diminution liée à de meilleur soins) serait de 336 000 (IC 95%: 192 000–537 000), en excluant le nombre de décès dans les maisons de retraite. Cette estimation pourrait être corroborée ou invalidée lorsque 80% de la population sera infectée, éventuellement sur plusieurs années, en supposant qu'un individu infecté est définitivement immunisé. Il faut noter que les mesures de confinement ou de distanciation sociale peuvent diminuer à la fois le pourcentage d'individus infectés dans la population et le degré de saturation des structures hospitalières.
Sur la valeur de R0. La distribution estimée en France est élevée par rapport aux estimations récentes (2,0–2,6, voir [3]) mais cohérente avec les résultats de [17] (2,24–3,58). Une estimation directe, par une méthode non mécaniste, des paramètres (ρ ,t0) d'un modèle de la forme \(\hat{\delta}_t =e^{ρ(t -t_0)}\) donne t0= 36 (5 février) et ρ =0.22. Avec le modèle SIR, I'(t) ≈ I(α - β) pour les petits temps (S≈ N), ce qui conduit à un taux de croissance égal à ρ ≈ α - β, et une valeur de α ≈ 0,32, c'est-à-dire R0= 3,2 , ce qui est cohérent avec notre distribution de R0. Notez que nous avons supposé ici une période d'infectiosité de 10 jours. Une période plus courte conduirait à une valeur inférieure de R0.
Sur l'incertitude liée aux données. L'incertitude sur le nombre réel d'infectés et donc l'IFR est très élevée. Il faut donc interpréter avec prudence les inférences qui peuvent être faites à partir des données dont nous disposons actuellement en France. De plus, nous ne faisons pas ici de prévisions: la dynamique future sera fortement influencée par les mesures de confinement qui seront prises et devra être modélisée en conséquence.
Sur la sensibilité des résultats par rapport aux paramètres fixes du modèle. Nous avons délibérément choisi un modèle parcimonieux avec quelques paramètres pour éviter les problèmes d'identifiabilité. Cependant, nous devions corriger certaines valeurs de paramètres. En particulier, nous avons supposé une durée moyenne de la période infectieuse (1 / β) de 10 jours. Une période infectieuse beaucoup plus longue de 20 jours (correspondant à la période médiane d'excrétion virale trouvée dans [12]) conduirait à un nombre de reproduction de base beaucoup plus grand R0= 4,8 (mais toujours compris entre 1,4 et 6,49 décrit dans [18]) et un facteur ×15 entre les cas signalés et le nombre réel de cas. Cependant, notre principal résultat sur l'IFR resterait inchangé (0,5%). Nous avons également évalué la sensibilité de l'inférence par rapport aux connaissances antérieures, en proposant un ensemble de distributions antérieures uniformes plus informatives que l'ensemble spécifié dans le texte principal. Dans l'ensemble, cette modification préalable n'influence pas significativement les distributions postérieures; voir l' annexe A .
Sur les hypothèses sous-jacentes au modèle. Les données utilisées ici contiennent peu d'informations, d'autant plus que la période d'observation considérée est courte et correspond à la phase initiale de la dynamique épidémique, qui peut être fortement influencée par des événements discrets. Cette limite nous a conduit à utiliser un modèle particulièrement parcimonieux afin d'éviter des problèmes d'identifiabilité des paramètres. Les hypothèses qui sous-tendent le modèle sont donc relativement simples et les résultats doivent être interprétés au regard de ces hypothèses. Par exemple, la date de l'introduction t0 doit être considérée comme une date d'introduction efficace pour une dynamique où une seule introduction serait déterminante pour l'éclosion et les autres introductions (antérieure et postérieure) auraient un effet insignifiant sur la dynamique.
Un modèle épidémiologique plus complexe de l'épidémie de COVID-19 en Chine a été proposé dans [19], avec une classe infectieuse divisée en plusieurs compartiments (individus asymptomatiques, infectieux symptomatiques non observés et infectieux symptomatiques observés). Les auteurs utilisent ce modèle dans [20] pour faire des prévisions sur le nombre cumulé de cas en Chine, tout en tenant compte des stratégies de gestion. Dans ces deux études, les auteurs soulignent l'importance de pouvoir estimer la fraction de cas infectieux non observés afin de prévoir la dynamique de l'épidémie. Notre étude, bien que basée sur un modèle SIR plus simple, montre que cette fraction peut être estimée sur la base de données précoces.
