© The authors. Published by EDP Sciences, 2020

1 Le modèle

La figure 1(a) montre le nombre cumulé de cas confirmés de coronavirus en France entre le 25 février et le 29 mars 2020 ; ces données incluent à la fois les laboratoires de biologie médicale de ville et les patients hospitalisés [10]. Il faut distinguer la date du 15 mars à partir de laquelle des mesures drastiques ont été prises subitement pour arrêter l’épidémie (fermeture des écoles, des restaurants, etc.). Pour ces trois dates, le nombre cumulé de cas confirmés est passé de 13 à 5423 puis à 40 174. La figure 1(b) montre les mêmes données en coordonnées logarithmiques ainsi que des droites de régression linéaire. On observe trois périodes : dans la première, qui va jusqu’au 6 mars, la croissance est rapide mais assez irrégulière ; dans la deuxième, qui va jusqu’au 15 mars, la croissance est un peu moins rapide mais régulière ; dans la troisième, à partir du 16 mars, la croissance est ralentie mais toujours régulière. Si l’on ajuste une droite sur l’ensemble des deux premières périodes, qui va du 25 février au 15 mars, on trouve que le nombre de cas croît comme e^λt avec un taux λ ≃0,31 par jour (droite rouge). Le temps de doublement est (ln2)∕λ ≃2,2 jours. Si en revanche on se limite à la deuxième période, avec des données qui sont particulièrement bien alignées en échelle logarithmique, on obtient λ ≃0,225 par jour et un temps de doublement de 3,1 jours (droite bleue). Comme les données du début de l’épidémie sont perturbées par une part importante de nouveaux cas importés et par des effets stochastiques, c’est probablement la deuxième estimation qui est la plus fiable. Pour la troisième période, après la mise en place de mesures drastiques, le temps de doublement passe à 4,9 jours.

On va étudier un modèle mathématique inspiré de cette épidémie. Divisons la population en cinq compartiments, selon une variante du modèle classique S-E-I-R (voir par exemple [3], p. 61) : susceptible d’être infecté (S), infecté en phase latente autrement dit non encore infectieux (E), infectieux sans protection (I), et retiré de la chaîne de transmission en étant comptabilisé parmi les cas confirmés (R₁) ou en ne l’étant pas (R₂). Chacun de ces deux derniers compartiments regroupe donc à la fois ceux qui sont encore infectieux mais isolés et ceux qui ne sont plusinfectieux car guéris ou décédés. Certains malades ont des symptômes de faible gravité et restent chez eux sans se faire tester, d’autres habitent dans des maisons de retraites et n’ont pas non plus été testés malgré les complications voire le décès ; ce sont ces catégories que l’on retrouve dans le compartiment R₂. On peut évidemment raffiner à l’infini ce modèle pour le rendre plus réaliste mais on a essayé ici de limiter au maximum le nombre de paramètres inconnus ; on a aussi pour but principal d’obtenir un résultat théorique relatif à la taille finale de l’épidémie.

Notons N la population totale (supposée grande) de sorte que N = S(t) + E(t) + I(t) + R₁(t) + R₂(t). Notons a le taux de contact effectif, b le taux auquel les personnes infectées en phase latente deviennent infectieuses, et c le taux moyen auquel les personnes infectieuses sont isolées et donc retirées de la chaîne de transmission. Notons f la fraction d’individus infectieux qui sont comptabilisés parmi les cas confirmés au moment de l’isolement (0 ≤ f ≤ 1) ; cette fraction peut varier au fil du temps mais on supposera pour simplifier qu’elle est constante. Alors

Pour faire le lien avec les données de la Figure 1, R₁ (t) correspond au nombre cumulé de cas confirmés à l’instant t. Si l’on note R(t) = R₁(t) + R₂(t), on remarque que (1.6)

avec R₁(0) = R₂(0) = 0, on en déduit que R₁(t) = f R(t) et R₂ (t) = (1 − f)R(t) pour tout t ≥ 0.

Au début de l’épidémie, le nombre de cas reste très petit par rapport à la population totale de sorte que S(t) ≃ N, ce qui conduit à la linéarisation

Les compartiments E et I mais aussi les compartiments R₁ et R₂ tendent donc à croître exponentiellement comme e^λt, où λ est la plus grande valeur propre de la matrice (1.7)

Le polynôme caractéristique est (1.8)

Donc (1.9)

[9] indique que la durée d’incubation, c’est-à-dire la période avant apparition des symptômes, est de 5 à 6 jours. La période de latence peutêtre un peu plus courte puisqu’on peut devenir infectieux avant de montrer des symptômes. On fixe la durée moyenne 1∕b dans la phase latente E à 4 jours ; donc b = 0,25 par jour.

La durée moyenne dans le compartiment I avant isolement, qui vaut 1∕c, est plus difficile à estimer car elle dépend de nombreux facteurs. Elle dépend des caractéristiques biologiques du virus, des caractéristiques des individus telles que leur âge, mais aussi de la promptitude avec laquelle les cas sont isolés, ce qui varie d’un pays à l’autre. L’épidémie en France a lieu alors que les habitants sont déjà bien au courant de l’existence de la pandémie ; les malades ne tardent pas trop à être isolés. Certains ne seront pas du tout infectieux, d’autre le seront plusieurs jours avant d’être isolés. Supposons que la moyenne soit de l’ordre de 1 jour, la forme du modèle sous-entendant que la distribution est exponentielle. On aurait une moyenne de cet ordre dans un modèle plus raffiné si par exemple 80 % des infectés restaient 0 jour infectieux et si 20 % restaient 5 jours infectieux avant d’être isolés. En résumé, on a choisi c = 1 par jour.

On déduit de la formule (1.9) que (1.10)

ce qui permettrait de calculer numériquement le taux de contact effectif a à partir du taux de croissance observé λ.

Imaginons que des mesures de santé publique puissent diviser le taux de contact effectif par un nombre k qui soit supérieur à 1. Combien doit valoir au minimum k pour arrêter l’épidémie ? Cette valeur de k, traditionnellement notée à la suite de Lotka et appelée par lui ≪ reproductivité ≫ ([6], p. 102), s’obtient simplement en remarquant que lorsque a est remplacé par , le nouveau taux de croissance de l’épidémie λ′ doit être nul, ce qui d’après l’équation (1.8) conduit à et à

si l’on utilise la valeur numérique λ ≃0,225 par jour suggérée par la courbe épidémique de la Figure 1. Vues les incertitudes sur les paramètres b et c, ceci ne peut être qu’une valeur approchée¹.

Revenons au modèle S-E-I-R non linéaire (1.1)–(1.6). Rappelons comment déterminer la taille finale de l’épidémie en l’absence complète d’intervention ; c’est une adaptation facile et d’ailleurs bien connue de la méthode utilisée pour le modèle S-I-R (voir par exemple [4], p. 76). L’équation (1.1) montre que

Donc en intégrant de t = 0 à t = +∞, (1.11)

où S(∞) désigne la limite quand t → +∞ de la fonction S(t) qui est décroissante et positive. Comme R(0) = 0, l’équation (1.6) montre que (1.12)

Par ailleurs, on a à tout instant S(t) + E(t) + I(t) + R(t) = N. Quand t → +∞, l’épidémie finit par s’arrêter de sorte que E(t) et I(t) tendent vers 0. À la limite, il ne reste donc que les personnes qui ont échappé à l’épidémie et celles qui ont été infectées mais qui sont passées dans les compartiments R : (1.13)

En combinant (1.11)–(1.13), on voit que

Au début de l’épidémie, il n’y a que quelques personnes infectées dans la population, donc S(0) ≃ N. L’équation implicite pour la taille finale de l’épidémie peut s’écrire comme (1.14)

qui se trouve avoir la même forme que pour le modèle S-I-R [4]. Avec , on trouve numériquement R(∞)∕N ≃ 87 %. Seule une fraction f de ces cas est recensée.

Figure 1 (a) Nombre cumulé de cas confirmés en France entre le 25 février et le 29 mars 2020, d’après Santé publique France et [11]; (b) Logarithme népérien de ce nombre et droites de régression linéaire.

2 Deuxième phase avec intervention drastique

Imaginons qu’à une certaine date T, des mesures drastiques soient prises de sorte que le nouveau taux de contact effectif soit réduit à 0 alors qu’il y a R₁ (T) cas confirmés cumulés. Par exemple, il y avait en France 5423 cas confirmés cumulés au 15 mars, date à laquelle sont entrées en vigueur les mesures concernant les écoles et les lieux publics. Peut-on alors prévoir quelle aurait été sous ces hypothèses idéales la nouvelle taille finale de l’épidémie R(∞), ou du moins celle confirmée R₁(∞) ?

Vers la fin de la phase exponentielle où t ≤ T et où le nombre total de cas représente encore une part infime de la population totale, on a

où (u, v) est un vecteur propre associé à la plus grande valeur propre λ de la matrice (1.7). Ainsi, − b u + a v = λ u. Avec l’équation (1.10), on trouve que

Comme dR∕dt ≃ λR pour t < T lorsque t n’est pas trop proche de 0, on a

Mais I(T) ≃ v e^λT, donc

Les contacts étant supposés réduits à zéro, on a pour t > T (2.1)

tandis que les autres équations (1.3)–(1.5) restent identiques. Sans avoir à résoudre ce système, il est clair que la taille finale de l’épidémie sera R(∞) = R(T) + E(T) + I(T), puisqu’il y a E(T) + I(T) individus infectés qui ne sont pas encore dans les compartiments R au temps T. Ainsi

Puisque à tout instant R₁(t) = f R(t), on en déduit aussi que

Ainsi, si les contacts sont réduits à zéro à partir d’une certaine date proche du début de l’épidémie (assez proche pour que l’approximation linéaire soit encore valable mais pas trop proche pour que le système linéarisé ait eu le temps de converger vers le vecteur propre associé à la première valeur propre), alors la taille finale (confirmée ou totale) de l’épidémie est proche de celle que l’on obtient en multipliant le nombre cumulé de cas (confirmés ou au total) à cette date par la reproductivité de l’épidémie. Un résultat semblable s’obtient de la même manière pour un modèle S-I-R. En annexe, on remarque cependant que ce n’est plus qui détermine le rapport R(∞)∕R(T) dans les modèles avec une période infectieuse qui n’est pas distribuée exponentiellement, mais une expression plus compliquée.

Avec R₁(T) = 5423 et , cela donne R₁ (∞) ≃ 12 600. Soulignons encore une fois l’incertitude autour des paramètres b et c, qui se retrouve dans la valeur de R₁(∞).

Notons au passage l’analogie avec le concept de ≪ potentiel d’accroissement ≫ des populations en démographie ([8], p. 176). C’est le rapport entre la population finale stationnaire et la population à un certain instant si la fertilité est divisée subitement à cet instant par la reproductivité , de sorte que la population se retrouve avec un taux asymptotique de croissance nul. Comme dans notre calcul, c’est en supposant que la population à cet instant est ≪ stable ≫ au sens de Lotka (c’est-à-dire donnée par le premier vecteur propre) que Keyfitz a obtenu une formule relativement simple pour le potentiel d’accroissement, formule qui fait aussi intervenir quoique de manière plus compliquée que pour notre modèle S-E-I-R ([8], p. 179).

Notons aussi que l’estimation de R(T) + E(T) + I(T) à partir de la donnée R(T) seule est analogue au problème qui s’était posé aux débuts de l’épidémie de VIH pour estimer le nombre de personnes séropositives à partir du nombre de cas déclarés de sida.

La Figure 2(a) illustre ce modèle à deux phases. On a pris N = 65 × 10⁶ (la population totale de la France) et les conditions initiales (2.2)

Le paramètre a est donné par la formule (1.10) avec λ = 0,225 par jour, comme dans la Figure 1. On dispose de peu d’information au sujet du paramètre f, si ce n’est rétrospectivement qu’un grand nombre de décès dus au virus dans les maisons de retraite n’étaient pas comptabilisés parmi les cas confirmés au début de l’épidémie ; fixons par exemple f =0,5 pour l’illustration. On a pris T = 43,2 jours de sorte que R₁(T) ≃ 5438 soit proche de la donnée 5423 du 15 mars. En poursuivant la simulation un peu plus longtemps que dans la figure, on trouve bien numériquement que .

La Figure 2(b) montre comment le rapport R₁ (∞)∕R₁(T) varie en fonction de l’instant T où le taux de contact est réduit à zéro. On observe effectivement un plateau où ce rapport est proche de . Quand T → 0, on a R₁(T) → 0 et R₁ (∞) → f(E(0) + I(0)) > 0, donc le rapport R₁ (∞)∕R₁(T) tend vers l’infini. Le rapport se rapproche de lorsque T est de l’ordre de l’inverse de la différence entre les deux valeurs propres de la matrice (1.7). Quand au contraire T →∞, alors l’intervention intervient trop tard ; l’épidémie est déjà passée et R₁ (∞)∕R₁(T) → 1. On s’attend à ce que la largeur du plateau où R₁(∞)∕R₁(T) est proche de croisse comme (lnN)∕λ lorsque N → +∞, puisque tel est par exemple le comportement du temps jusqu’au pic épidémique dans un modèle S-I-R à coefficients constants (voir [2] ou [5], p. 12).

Figure 2 (a) Exemple de simulation du modèle à deux phases. (b) Le rapport R1 (∞)∕R1(T) en fonction de T.

3 Généralisation

Dans la réalité, le taux de contact effectif n’est sûrement pas tout à fait nul pour t > T. La valeur obtenue pour R(∞) peut néanmoins être considérée comme une borne inférieure de la valeur réelle puisqu’il est certain que la taille finale de l’épidémie sera supérieure avec des contacts non nuls qu’avec des contacts nuls pour t > T. Rappelons cependant à ce sujet que les modèles épidémiques de type S-I-R ou S-E-I-R avec un taux de contact variable ne sont pas ≪ monotones ≫, dans le sens qu’une réduction du taux de contact peut parfois conduire à une taille finale de l’épidémie plus grande [1].

Considérons maintenant le cas où le taux de contact n’est pas réduit à 0 mais simplement divisé par un nombre q > 1. La réduction à 0 correspond au cas limite où q tend vers l’infini. On a pour t > T, (3.1)

tandis que les équations (1.3)–(1.5) restent identiques. Par le même raisonnement que dans la section 1, on a pour t > T

En intégrant entre t = T et t = +∞, on en déduit que

où . Comme S(∞) = N − R(∞), on a donc (3.2)

Supposons comme dans la section 2 que le temps T ne soit ni trop petit ni trop grand, c’est-à-dire dans le plateau de la Figure 2(b). En première approximation, S(T) ≃ N et R(T) est encore petit devant N. Deux cas se présentent alors.

Si , alors l’argument graphique classique [4] qui consiste à tracer les membres de gauche et de droite de l’équation (3.2) en fonction de R(∞)∕N montre que la solution R(∞)∕N n’est pas infinitésimale mais proche de la solution strictement positive de l’équation (3.3)

Si au contraire , alors la solution R(∞)∕N de l’équation (3.2) est petite. Comme S(T) = N − E(T) − I(T) − R(T) et , un développement à l’ordre 1 de l’exponentielle dans l’équation (3.2) conduit à

En ne gardant que les termes d’ordre inférieur ou égal à 1, on trouve que

Finalement, (3.4)

Lorsque q → +∞, on retrouve bien que . On remarque aussi que , comme il se doit. Une formule identique à (3.4) lie R₁(∞) et R₁(T).

La Figure 3 montre avec une ligne continue, en fonction du paramètre q, la taille finale de l’épidémie en échelle logarithmique, ln(R(∞)∕N), obtenue parsimulation numérique du système (1.1)–(1.6) pour t < T avec les conditions initiales (2.2) et du système (3.1) pour t > T. Comme dans la Figure 2(a), la population totale est N = 65 × 10⁶ et le paramètre a est donné par la formule (1.10) avec λ = 0,225 par jour ; ona encore pris f = 0,5 et T = 43,2 jours de sorte que R₁(T) ≃ 5 438. La figure montre aussi avec des petits ronds ce que donne la formule (3.4) pour . Elle montre enfin avec de petits losanges la solution strictement positive de l’équation (3.3) pour . On voit que les deux approximations cessent d’être valables au voisinage de .

On notera que la taille finale de l’épidémie varie de plusieurs ordres de grandeurs lorsque le paramètre q est proche de . Comme il est difficile de lequantifier, la prédiction de la taille finale de l’épidémie est également difficile dans cette zone. Il n’y a que si le paramètre q est nettement supérieur à que la prévision avec la formule (3.4) devient moins sensible à la valeur de q.

Figure 3 ln(R(∞)∕N) en fonction de q (ligne continue), comparé avec la formule (3.4) (petits ronds) valable pour et avec la solution de l’équation (3.3) (petits losanges) valable pour .

4 Tentative d’estimation du paramètre q

Essayons d’estimer le paramètre q en ajustant une simulation du modèle aux données postérieures au 15 mars, y compris celles jusqu’au 15 avril qui ne figuraient pas dans la Figure 1. [10] avertit néanmoins que ≪ le nombre de cas confirmés en France ne reflète plus de manière satisfaisante la dynamique de l’épidémie ≫, étant donné que ≪ les patients présentant des signes de COVID-19 ne sont plus systématiquement confirmés par un test biologique ≫. Pour cela, on part de la donnée R₁(T) = 5423 et des relations R(T) = R₁(T)∕f et R₂(T) = (1 − f)R(T). Comme les données des 8 jours qui précèdent sont particulièrement bien alignées, on démarre la simulation de notre modèle avec R(T − τ) = e^−λτR(T) où λ = 0,225 par jour et τ = 8 jours, et avec les estimations correspondantes , et S(T − τ) = N − E(T − τ) − I(T − τ) − R(T − τ).

Pour t > T, le taux de contact effectif est a∕q et l’on essaie d’ajuster R₁(t) aux données jusqu’au 15 avril. Le meilleur ajustement se trouve aux alentours de q = 1,7 (voir la Fig. 4). Comme cette valeur est inférieure à , il semblerait que les mesures de confinement soient encore insuffisantes. Mais les tous derniers points de la figure montrent que l’écart avec le modèle grandit dans le sens d’un ralentissement de l’épidémie réelle. Il se peut que la valeur de f choisie ne soit pas appropriée ou qu’elle ait varié au cours de l’épidémie. Ou alors le modèle est peut-être un peu trop simpliste ; on s’attend notamment à ce qu’une distribution non exponentielle des temps passés dans les différents compartiments influence le moment où la courbe commence à s’infléchir.

En conclusion, on a exploré un scénario à deux phases où le taux de contact est réduit à partir d’un certaine date. On a trouvé une formule approchée simple pour la taille finale de l’épidémie en fonction du nombre de cas détectés au moment de la réduction. Il reste néanmoins à énoncer et à démontrer plus rigoureusement ce résultat, probablement en le faisant apparaître comme un résultat asymptotique lorsque N → +∞.

Figure 4 Logarithme népérien du nombre de cas recensés entre le 7 mars et le15 avril (petits ronds, données de Santé publique France [11]) et ln (R1 (t)) en fonction du temps t dans 4 simulations avec de haut en bas q ∈{1,5; 1,7; 2; 2,5}.

Annexe

Considérons un modèle S-I-R avec une période infectieuse qui n’est pas nécessairement distribuée exponentiellement. Soit I(t, x) la densité de personnes infectées depuis x unités de temps au temps t. Soit a(x) le taux de contact effectif et b(x) le taux auquel les personnes infectées cessent de transmettre l’infection. On a au début de l’épidémie

On en déduit, comme dans la théorie des populations stables de Lotka [6] que

où k est une constante et le taux de croissance λ est l’unique solution de l’équation

Si l’on pose , le problème est d’estimer I(T) + R(T) à partir de R(T). Or

On en déduit que

Finalement,

On voit que le membre de droite n’a pas de raison particulière de coïncider avec . Dans le cas spécial où les taux sont constants, avec a(x) ≡ a et b(x) ≡ b, on a cependant λ = a − b et donc .

Remerciements

On remercie Hisashi Inaba, Ali Moussaoui et Frédéric Hamelin pour leurs commentaires sur le manuscrit.

Références

Liste des figures

	Figure 1 (a) Nombre cumulé de cas confirmés en France entre le 25 février et le 29 mars 2020, d’après Santé publique France et [11]; (b) Logarithme népérien de ce nombre et droites de régression linéaire.
Dans le texte

	Figure 2 (a) Exemple de simulation du modèle à deux phases. (b) Le rapport R1 (∞)∕R1(T) en fonction de T.
Dans le texte

	Figure 3 ln(R(∞)∕N) en fonction de q (ligne continue), comparé avec la formule (3.4) (petits ronds) valable pour et avec la solution de l’équation (3.3) (petits losanges) valable pour .
Dans le texte

	Figure 4 Logarithme népérien du nombre de cas recensés entre le 7 mars et le15 avril (petits ronds, données de Santé publique France [11]) et ln (R1 (t)) en fonction du temps t dans 4 simulations avec de haut en bas q ∈{1,5; 1,7; 2; 2,5}.
Dans le texte