Introduction

Ce travail rentre dans le cadre d’un projet d’UE : modélisation en épidémiologie (MEPI).

Dans le cadre de ce projet, nous allons construire un modèle de dynamique épidémiologique afin de répondre à la question suivante : Quel est l’impact de la distanciation physique (post 11 mai) relative au confinement sur la dynamique épidémique en France ?

Les réponses à cette question influent directement sur les mesures mises en place par le gouvernement. La distanciation physique comme le confinement de la population, sont des mesures dont le but est d’agir sur un paramètre précis de l’épidémie de COVID-19 : le taux de transmission de la maladie. Le taux de transmission dépend de la transmissibilité de la maladie et du taux de contact entre individus. C’est en particulier sur ce dernier paramètre que les mesures agissent.

Le but de notre travail est donc de trouver un modèle qui décrit l’épidémie sur deux périodes différentes : confinement et post-confinement. Une fois ce modèle calibré, nous pourrons faire varier le taux de contact pour voir l’impact de la distanciation physique sur la dynamique de l’épidémie, mais aussi d’autres mesures sanitaires.

Importations

Le modèle que nous allons mettre en place s’appuie sur des équations différentielles qui ne peuvent être résolues analytiquement. Nous allons donc utiliser le module “deSolve” pour résoudre numériquement ces équations. Notre modèle sera calibré sur les données d’hospitalisation publiées par santé publique France.

library(deSolve)
options(warn=-1)  # desactive les messages d'erreur
rm(list=ls())     # nettoie l'environnement de travail
graphics.off()    # ferme les graphes précédents

# donnees
dataA = read.table("Donnees-nouveaux-hopital-1Dec.csv", 
                                     header=TRUE, sep=";") # Admissions
dataC = read.table("Donnees-hopital-1Dec-classes_age.csv", 
                                     header=TRUE, sep=";") # Hospitalisés par classe

Selection des données

Les fichiers .csv importés contiennent de nombreuses données, dont une partie ne nous est pas utile. Nous allons donc manipuler ces tableaux afin de garder uniquement : le nombre de personnes hospitalisées par classe d’âge et le nombre total d’hospitalisations par jour en France métropolitaine. Nous utiliserons le même regroupement des classes d’âge que dans cette publication de Di Domenico et al.

# donnees
dataC = dataC[complete.cases(dataC), ] # enleve les NA

dataA = dataA[,c("dep", "jour", "incid_hosp")]
dataC = dataC[,c("reg","cl_age90","jour","hosp")]

n = 101                     # nombre total de departements fr
LA = length(dataA[,1])/n    # nombre de jours dans les donnees A
m = 18 * 11                 # nombre de region (18) * nombre de classes d'age (11)
LC = length(dataC[,1])/m    # nombre de jours dans les donnees C

# verification de la detection du format date
dataA[,"jour"] = as.Date(dataA[,"jour"], format="%Y-%m-%d")
dataC[,"jour"] = as.Date(dataC[,"jour"], format="%d/%m/%Y")

# creation des tableaux pour les donnees finales
rdyA = data.frame("jour" = unique(dataA[,"jour"]), 
                                    "admi" = rep(0,LA), stringsAsFactors=FALSE)
rdyC = data.frame("jour" = unique(dataC[,"jour"]), 
                                    "0" = rep(0,LC), 
                                    "9" = rep(0,LC), 
                                    "19" = rep(0,LC), 
                                    "29" = rep(0,LC), 
                                    "39" = rep(0,LC), 
                                    "49" = rep(0,LC), 
                                    "59" = rep(0,LC), 
                                    "69" = rep(0,LC),
                                    "79" = rep(0,LC), 
                                    "89" = rep(0,LC),
                                    "90 plus" = rep(0,LC), 
                                    stringsAsFactors=FALSE)

# selection des 96 departements de france metropolitaine (corse incluse)
for (i in 1:LA)
{
    Ajour = dataA[which(dataA[,"jour"]==rdyA[i,1]),]
    rdyA[i,2] = sum(Ajour[1:96,3])
}

# on retire les dept d'outre mer (de 01 a 06)
dataC = dataC[which(is.element(dataC[,"reg"], (10:100))),]

# Attribution des valeurs totales de chaque classe d'age (somme des valeurs des départements)
ages = c(0, 9, 19, 29, 39, 49, 59, 69, 79, 89, 90)
for (i in 1:LC)
{
    Cjour = dataC[which(dataC[,"jour"]==rdyC[i,1]),]
    for (j in 1:11)
    {
        rdyC[i,j+1] = sum(Cjour[which(Cjour[,"cl_age90"] == ages[j]),4])
    }
    
}

# Fusion des classes d'ages en groupes 
jeune = rdyC[,c("X9")]
ado = rdyC[,c("X19")]
adulte = rdyC[,"X29"] + rdyC[,"X39"] + rdyC[,"X49"] + rdyC[,"X59"]
senior = rdyC[,"X69"] + rdyC[,"X79"] + rdyC[,"X89"] + rdyC[,"X90.plus"]

# Représentations graphiques des données brutes
par(mfrow=c(1,1))
plot(1:LA, rdyA[,2],
         main="Admissions à l'hôpital depuis le 18 mars 2020",
         xlab="Jours écoulés depuis le 18 mars 2020",
         ylab="Nombre d'admis par jour")

plot(1:LC, rdyC[,2], 
         main="Individus hospitalisés depuis le 18 mars 2020",
         xlab="Jours écoulés depuis le 18 mars 2020",
         ylab="Nombre d'hospitalisés par jour")

plot(1:LC, rdyC[,3], col = 3*9, type = 'l', lty = 1, 
         ylim = c(0,max(rdyC[,3:12])),
         main ="Données pour chaque classe d'âge",
         xlab="Jours écoulés depuis le 18 mars 2020",
         ylab="Nombre d'hospitalisés par jour")
for (k in 4:12)
{
    points(1:LC, rdyC[,k], col = k*9, type = 'l', lty = 1)
}
legend("topright", lty = 1,
       legend = names(rdyC)[3:12],
             col = c(seq(3,12)*9))

Modélisation

Hypothèses de travail

Nous avons considéré les hypothèses suivantes :

  • Les temps de guérison et d’hopitalisation sont connus et fixés
  • Les mortalités et natalités naturelles sont négligées
  • Les effets de migrations sont négligés
  • La dynamique est homogène sur tout le pays
  • Les hommes et les femmes sont égaux vis-à-vis du virus
  • Les individus hospitalisés ne transmettent plus l’infection
  • La guérison confère une immunité permanente à la maladie

Modèle

Pour notre modèle, nous avons considéré une dynamique de type SIHR divisée en 4 classes d’âge (le modèle est décrit en détail ci-après). En effet, nous disposons des données d’hospitalisations et d’admissions à l’hôpital. Nous pourrons donc estimer I et H aisément. Les classes d’âge permettent de modéliser des comportements différents selon l’âge, et donc des taux de contact différents (inter et intra-classes). On définit les états suivants :

  • S : l’état sensible, population pouvant être infectée
  • I : l’état infectieux, population qui transmet la maladie
  • H : l’état hospitalisé, population infectieuse nécessitant des soins
  • R : l’état retiré, population guérie ou immunisée à la maladie
  • N : population totale

Pour les différentes classes, on définit les limites suivantes :

  • Jeunes (J) : 0-9 ans
  • Adolescents (T) : 10-19 ans
  • Adultes (A) : 20-59 ans
  • Seniors (S) : 60 ans et plus

Chaque état est indicé pour désigner la classe d’âge auquelle il correspond. Un état non indicé désigne la population totale de cet état, sans distinction de classe d’âge.

La dynamique de l’épidémie est modélisée suivant les équations suivantes (exemple avec la classe ‘Jeunes’) :

\(dS_J/dt = - S_J(t) \cdot \left(\frac{\beta_{JJ} \cdot I_J(t)}{N_J} + \frac{\beta_{JT} \cdot I_T(t)}{N_T} + \frac{\beta_{JA} \cdot I_A(t)}{N_A} + \frac{\beta_{JS} \cdot I_S(t)}{N_S}\right)\)
\(dI_J/dt = S_J(t) \cdot \left(\frac{\beta_{JJ} \cdot I_J(t)}{N_J} + \frac{\beta_{JT} \cdot I_T(t)}{N_T} + \frac{\beta_{JA} \cdot I_A(t)}{N_A} + \frac{\beta_{JS} \cdot I_S(t)}{N_S}\right) - (\alpha_J + \rho) \cdot I_J(t)\)
\(dH_J/dt = \alpha_J \cdot I_J(t) - \gamma \cdot H_J(t)\)
\(dR(t) = N - S(t) - I(t) - H(t)\)

Tous les paramètres seront déduits des données extraites au début de ce script. Les mesures pour H0 seront considérées comme exactes pour réduire le nombre de conditions initiales à estimer.

Paramètres

Les coefficients des équations ci-dessus sont décrits ici :

  • \(\alpha\) est le taux d’hospitalisations des infectés (indice \(J\) signifie : taux d’hospitalisation dans la classe \(J\))
  • \(\beta\) est le taux de tansmission (indice \(JA\) signifie : transmission de \(J\) vers \(A\))
  • \(\gamma\) est l’inverse du nombre de jours moyen passés à l’hopital
  • \(\rho\) est l’inverse du nombre de jours moyens avant guérison

\(\beta\) est donc le taux de transmission, directement influencé par le taux de contact. Nous cherchons à le connaitre avec précision afin de le manipuler pour répondre à notre problématique. \(\alpha\) est estimé pour chaque classe afin de nous approcher des données importées. La durée passées à l’hopital et la durée de contagion sont fixées, égales pour chaque classe, et issues de https://www.santepubliquefrance.fr/ et https://www.pasteur.fr/. Nous ne cherchons pas à les étudier pour répondre à notre problématique.

# Parametres fixes
rho = 1/11                    # temps moyen avant guerison
gamma=1/21                    # Duree moyenne de l'hospitalisation

# Variables (valeurs d'initialisation)
alpha = c(0.1, 0.1, 0.1, 0.1) # taux d'hospitalisation des infectieux
beta = rep(0.01,16)           # transmission  

P0=c(alpha,beta)              # Vecteur des parametres à estimer

Conditions initiales

Pour déterminer les conditions initiales, c’est-à-dire le nombre d’individus initialement dans chacun des états du modèle, nous allons nous servir des données que nous avons importées précédemment.

Périodes d’étude

Nous allons considérer deux périodes :

  • le premier confinement : du 17 Mars au 11 Mai
  • la période de post-confinement : à partir du 11 Mai au 11 juillet (nous n’étudierons pas les vacances d’été)

Le “Lag”, ou période d’incubation, initialisé ci-dessous, permet de nous assurer que tous les individus infectés avant le confinement présentent des symptômes de la maladie au début de notre simulation.

Lag = 12                 # Periode d'incubation maximale supposee

# Confinement
t0c = 1                  # Debut de la periode
tfc = t0c + 54           # Fin de la peride
tc = (t0c+Lag):(tfc+Lag) # Vecteur temps de la période

# post confinement
t0p = tfc
tfp = (t0p + 61)
tp = (t0p+Lag):(tfp+Lag)

A présent, nous pouvons estimer les conditions initiales à partir des données. \(H_0\) est fixé par les données d’hospitalisations. \(R_0\) est fixé selon l’hypothèse qu’en début de confinement, 1% de la population française était immunisée au virus. Pour \(I_0\) et \(S_0\), on initialise seulement les valeurs mais elles seront éstimées lors de l’optimisation du modèle. Ainsi, les données ne nous permettent uniquement de fixer le point de départ que pour \(H_0\). Nous utilisons la répartition prédite par l’INSEE de la popultation française par tranche d’âge pour déduire certaines valeurs.

N=64.8e6                       # Population française en metropole (approx.)
dist <- c(11.9, 11.9, 49.8, 26.4)*(1e-2) # Distribution pop par tranche d'age

H0 = c(jeune[tc[1]],           # Nombre d'hospitalises dans chaque classe
             ado[tc[1]],
             adulte[tc[1]],
             senior[tc[1]])

I0 = dist*rdyA[tc[1],2]/alpha  # Nombre d'admissions dans chaque classe (approx.)

R0=N*0.01                      # Nombre d'immunises naturels (hypothese)

S0=(0.99*N*dist)-I0-H0
X0=c(S0,I0)                    # Conditions initiales a estimer

Fonctions du modèle

Une première fonction prend en arguments nos paramètres et états au temps t, elle calcule l’évolution de chaque état avec leurs dérivées. Les solutions seront estimées par la fonction ode de la librairie deSolve. Cette fonction permet d’obtenir les vecteurs correspondants à la trajectoire de chaque état SIH. L’état R peut être déduit en soustrayant ces états à la population totale mais nous ne nous y intéresserons pas. Cette fonction aurait pu être optimisée par l’utilisation de calculs matriciels mais par soucis de clarté, nous avons opté pour une écriture linéaire.

Ensuite, il nous faut une fonction qui calcule l’écart de nos prédictions aux données. Pour cela nous créons la fonction Erreur. Elle prend en arguments les prédictions du modèle et les données importées puis calcule la somme des carrés des écarts pour plusieurs courbes. Nous avons préféré cette méthode à celle de la vraissemblance afin de privilégier une représentation graphique proche des données. En effet, ce qui nous intéresse n’est pas d’estimer de façon absolue les paramètres et états du modèle, mais de comparer de façon relative les deux périodes (confinement et post-confinement). Les erreurs sont normalisées pour limiter que les totaux aient plus d’influence que les classes d’âge.
Nous avons décidé d’ajuster les paramètres \(\alpha\) pour chaque classe à l’ensemble des deux périodes considérées. Cependant, les paramètres \(\beta\) sont différents d’une période à l’autre. En effet, nous avons émis l’hypothèse qu’entre le confinement et le post-confinement, seul le taux de tansmission varie; le taux d’hospitalisation si infecté restant le même.

Nous utiliserons la fonction optim implantée dans R afin de minimiser l’erreur entre les distributions des effectifs modélisées et les données.

# modele
SIHR = function(t, X, P)
{
  # parametres
  alphaJ = P[1]  # alpha (valeur seule)
  alphaT = P[2]
  alphaA = P[3]
  alphaS = P[4]
  
  betaJJ = P[5]  # beta (matrice)
  betaJT = P[6]
  betaJA = P[7]
  betaJS = P[8]
  
  betaTJ = P[9]  # beta (matrice)
  betaTT = P[10]
  betaTA = P[11]
  betaTS = P[12]
  
  betaAJ = P[13]  # beta (matrice)
  betaAT = P[14]
  betaAA = P[15]
  betaAS = P[16]
  
  betaSJ = P[17]  # beta (matrice)
  betaST = P[18]
  betaSA = P[19]
  betaSS = P[20]
  
  gamma = gamma   # gamma (fixé)
  rho = rho       # tho (fixé)
  
  # variables externes
  # distribution de la population
  NJ = N*dist[1]
  NT = N*dist[2]
  NAd = N*dist[3]
  NS = N*dist[4]
  
  # conditions initiales
  SJ = X[1]     # vecteur des S (un par classe)
  ST = X[2]
  SA = X[3]
  SS = X[4]
  
  IJ = X[5]     # vecteur des I (un par classe)
  IT = X[6]
  IA = X[7]
  IS = X[8]
  
  HJ = X[9]     # vecteur des H (un par classe)
  HT = X[10]
  HA = X[11]
  HS = X[12]
  
  # Calcul de simplification
  yJ = (betaJJ*SJ*IJ/NJ)+(betaJT*SJ*IT/NT)+(betaJA*SJ*IA/NAd)+(betaJS*SJ*IS/NS)
  yT = (betaTJ*ST*IJ/NJ)+(betaTT*ST*IT/NT)+(betaTA*ST*IA/NAd)+(betaTS*ST*IS/NS)
  yA = (betaAJ*SA*IJ/NJ)+(betaAT*SA*IT/NT)+(betaAA*SA*IA/NAd)+(betaAS*SA*IS/NS)
  yS = (betaSJ*SS*IJ/NJ)+(betaST*SS*IT/NT)+(betaSA*SS*IA/NAd)+(betaSS*SS*IS/NS)
  
  # Dérivées des états sensibles
  dSJ = - yJ
  dST = - yT
  dSA = - yA
  dSS = - yS
  
  # Dérivées des états infectieux
  dIJ = yJ - (alphaJ + rho)*IJ
  dIT = yT - (alphaT + rho)*IT
  dIA = yA - (alphaA + rho)*IA
  dIS = yS - (alphaS + rho)*IS
  
  # Dérivées des états hospitalisés
  dHJ = alphaJ*IJ - gamma*HJ
  dHT = alphaT*IT - gamma*HT
  dHA = alphaA*IA - gamma*HA
  dHS = alphaS*IS - gamma*HS
  
  # Renvoie une liste utilisable par ode()
  return(list(c(dSJ,dST,dSA,dSS,dIJ,dIT,dIA,dIS,dHJ,dHT,dHA,dHS)))
}


# Calcul de l'ecart aux donnees
Erreur = function(theta)
{
    ### Confinement
    t = tc                        # vecteur temps
    X = c(theta[1:8], H0)         # conditions initiales 
    P = c(theta[9:28], gamma, rho)# parametres
    Y = ode(X,t,SIHR,P)           # resolution du systeme d'EDO
    
    # Solutions estimees
    # Totaux
    h = rowSums(tail(Y[,10:13],-1)) # Hospitalisations modelisees : H(t)
    a = rowSums(P[1:4]*Y[,6:9])     # Admissions modelisees : alpha*I(t)
    
    # Hospitalisations de chaque classe
    hJ = tail(Y[,10],-1)
    hT = tail(Y[,11],-1)
    hA = tail(Y[,12],-1)
    hS = tail(Y[,13],-1)
    
    # Difference realite - modele (normalisée)
    SSC=((rdyC[t,"X0"]-h)**2)/(max(h)-min(h))
    SSA=((rdyA[t,"admi"]-a)**2)/(max(a)-min(a))
    
    SSJ=((jeune[t]-hJ)**2)/(max(hJ)-min(hJ))
    SST=((ado[t]-hT)**2)/(max(hT)-min(hT))
    SSA=((adulte[t]-hA)**2)/(max(hA)-min(hA))
    SSS=((senior[t]-hS)**2)/(max(hS)-min(hS))
    SS_conf = sum(c(SSC,SSA,SSJ,SST,SSA,SSS))
    
    ### Post-confinement
    t = tp
    X = c(Y[tfc,2],Y[tfc,3],Y[tfc,4],Y[tfc,5],
          Y[tfc,6],Y[tfc,7],Y[tfc,8],Y[tfc,9],
          Y[tfc,10],Y[tfc,11],Y[tfc,12],Y[tfc,13])
    P = c(theta[9:12],theta[29:44], gamma, rho)
    Z = ode(X,t,SIHR,P)
    
    # Solutions estimees
    # Totaux
    h = rowSums(tail(Z[,10:13],-1)) # Hospitalisations modelisees : H(t)
    a = rowSums(P[1:4]*Z[,6:9])     # Admissions modelisees : alpha*I(t)
    
    # Hospitalisations de chaque classe
    hJ = tail(Z[,10],-1)
    hT = tail(Z[,11],-1)
    hA = tail(Z[,12],-1)
    hS = tail(Z[,13],-1)
    
    # Difference realite - modele (normalisées)
    SSC=((rdyC[t,"X0"]-h)**2)/(max(h)-min(h))
    SSA=((rdyA[t,"admi"]-a)**2)/(max(a)-min(a))
    
    SSJ=((jeune[t]-hJ)**2)/(max(hJ)-min(hJ))
    SST=((ado[t]-hT)**2)/(max(hT)-min(hT))
    SSA=((adulte[t]-hA)**2)/(max(hA)-min(hA))
    SSS=((senior[t]-hS)**2)/(max(hS)-min(hS))
    SS_post = sum(c(SSC,SSA,SSJ,SST,SSA,SSS))
    
    # Total (avec vérification des résultats)
    SS_tot = sum(SS_conf,SS_post)
    if (SS_tot == Inf) {SS_tot = 1e12}
    
    return(SS_tot)  # Renvoie l'erreur totale
}

Résolution des équations

optimisation du modèle

Il ne reste plus qu’à initier les vecteurs des paramètres et des conditions initiales puis à se servir des fonctions créées pour déduire les valeurs optimales de ces vecteurs.

theta = c(X0,P0,beta) # paramètres à faire varier 

# Estimation des coefficients
opt = optim(theta, Erreur, lower = 0, method = "L-BFGS-B")

# Attribution des résultats
X0_conf = c(opt$par[1:8], H0)          # conditions initiales
P0_conf = c(opt$par[9:28], gamma, rho) # parametres
alpha = opt$par[9:12] 

# Mise en forme des beta
beta_conf = matrix(opt$par[13:28], nrow = 4, ncol = 4, byrow = TRUE)
dimnames(beta_conf) = list(c("J","T","A","S"),c("J","T","A","S"))

# resolution
mod_conf = ode(X0_conf, tc, SIHR, P0_conf)

Les valeurs obtenues pour les \(\beta\) sont données dans cette matrice (Se lit "coefficient de transmission de [ligne] à [colonne]):

##             J           T           A          S
## J 0.005188716 0.008105219 0.010839765 0.04490042
## T 0.003141736 0.006365779 0.009599446 0.04961805
## A 0.000000000 0.000000000 0.000000000 0.03524624
## S 0.051951242 0.052860422 0.082299676 0.08263697

Nous pouvons déduire le \(\mathcal R_0\) sur la période du confinement en trouvant la valeur propre dominante de la matrice contenant les coefficients \(\frac{\beta_{ij}}{\alpha_{i} + \rho}\).

M_conf = matrix(data=0, nrow=4, ncol=4)
for (i in 1:4)
{
    M_conf[i,] = beta_conf[i,]/(alpha+rho)
}

R0_conf = Re(eigen(M_conf)$values[1])

Nous obtenons \(\mathcal R_0\) = 0.7851089

Il est bien inférieur à \(1\) ce qui signifie que le nombre d’infectieux diminue. Cette estimation pour le confinement n’est pas si loin de celle de l’INRAE (Roques et al.) et encore moins de celle de l’INSERM (Di Domenico et al.) où le \(\mathcal R_0\) avait été estimé respectivement à 0.47 et 0.66. Les écarts avec notre modélisation viennent probablement du choix d’utiliser les distances aux données plutôt que la vraisemblance et d’autres hypothèses. Mais notre résultat reste totalement plausible et la dynamique du confinement est donc bien modélisée par notre modèle.

A présent, nous effectuons les mêmes étapes pour la période de post-confinement. A ceci près que les conditions initiales sont tirées des conditions finales du confinement : seulement les \(\beta\) varient. Cela nous permet de nous assurer que seul ce paramètre influe sur la dynamique en post-confinement, et nous permettra de répondre à notre problématique.

beta_post = c(opt$par[29:44])             # Initialisation du vecteur

# récupération des états finaux du confinement
X0_post = c(mod_conf[tfc,2],mod_conf[tfc,3],mod_conf[tfc,4],mod_conf[tfc,5],
            mod_conf[tfc,6],mod_conf[tfc,7],mod_conf[tfc,8],mod_conf[tfc,9],
            mod_conf[tfc,10],mod_conf[tfc,11],mod_conf[tfc,12],mod_conf[tfc,13])
P0_post = c(alpha, beta_post, gamma, rho) # vecteur des paramètres

# Resolution
mod_post = ode(X0_post, tp, SIHR, P0_post)

# Mise en forme des beta
beta_post = matrix(beta_post, nrow = 4, ncol = 4, byrow = F)
dimnames(beta_post) = list(c("J","T","A","S"),c("J","T","A","S"))

Les valeurs obtenues pour les \(\beta\) sont données dans cette matrice :

##            J          T           A          S
## J 0.01072529 0.01113679 0.006026447 0.04380151
## T 0.01077240 0.01118614 0.005641177 0.04513690
## A 0.01440965 0.01556748 0.023235999 0.10127002
## S 0.02126602 0.02275430 0.043090514 0.12631315

Elles sont, dans l’ensemble, plus élevées que pendant le confinement. Notre modèle capture donc bien un changement de comportement avec le déconfinement. Nous en déduisons le \(\mathcal R_0\) de la même manière que précédemment.

Nous obtenons \(\mathcal R_0\) = 0.8962448

Nous constatons que le \(\mathcal R_0\) est plus grand en cette période de déconfinement. Cela s’explique par une augmentation des taux de transmission \(\beta\). Naturellement, le virus circule plus lorsque la population n’est plus confinée. La dynamique de l’épidémie est cependant toujours en diminution à ce moment. Le \(\mathcal R_0\) est inférieur à 1 mais de peu.
Si le \(\mathcal R_0\) dépasse 1, cela signifie qu’un individu infectieux transmet la maladie à plus d’une personne en moyenne : l’épidémie repart à la hausse. C’est ce qui a été observé à la fin des vacances d’été, après la rentrée. La frontière entre ces situations est donc très proche.

Représentations

Maintenant quelques graphes : les coefficients permettent de modéliser les données d’hospitalisations avec une bonne précision. Il n’y a que peu d’écarts aux données.

Notre but n’est pas d’analyser les variables une par une. Cependant, nous pouvons rappeler que le taux de transmission est plus grand après le confinement d’où le ralentissement dans la diminution des hospitalisations au déconfinement. Nos modèles s’approchent assez de la réalité pour nous permettre de tester l’influence de leur variation sur la dynamique épidémique. Et notamment, l’influence de la variation du taux de contact.

NB : les points représentent les données tandis que les lignes représentent les prédictions modélisées

Quelques commentaires sur les écarts constatés :

Pour les admissions, les courbes modélisées s’éloignent légèrement des admissions réelles. Cela s’explique par le fait que le taux \(\alpha\) n’est pas stable en réalité. Il évolue dans le temps. Néanmoins, la tendance des données est bien reflétée par les courbes, même si leur valeur absolue s’éloignent des données.

Pour les hospitalisations, les écarts sur la période de déconfinement peuvent s’expliquer par la présence de plusieurs phase de déconfinement (d’après les mesures du gouvernement). Néanmoins, les courbes restent très proches des données et les écarts reflètent surtout la variance des relevés.

Influence du taux de contact

A présent, voyons l’influence du taux de contact sur la dynamique de l’épidémie. Au lieu de considérer Beta, on va décomposer cette variable en deux facteurs : \(\beta = \beta_c \cdot \beta_\mu\)
où \(\mu\) désigne le taux de transmission
et \(c\) les taux de contact entre classes d’âge.

Jusqu’ici, nous ne nous étions pas intéressés à la division de \(\beta\) en deux facteurs. Pour l’introduire, nous allons devoir suivre plusieurs étapes :

Evidemment, nous ne pourrons pas modéliser des changements absolus avec les hypothèses que nous faisons. Mais nous pourrons montrer des évolutions relatives, qui seront donc à nuancer. Ce projet n’est pas à but de prédiction, mais plutôt de comparaison.

# Fonction pour calculer une nouvelle matrice des taux des contacts
matrice_contact = function(mesures){
    ecoles = mesures[1]
    teletravail = mesures[2]
    social = mesures[3]    
    commerce = mesures[4]      
    ephad = mesures[5]   
    distance = mesures[6]
    
    # Initialisation de la matrice
    beta.c = matrix(data=1, nrow=4, ncol=4)
    dimnames(beta.c) = list(c("J","T","A","S"),c("J","T","A","S"))
    
    
    # effet ecoles 
    beta.c["J","J"] = max(0.30*ecoles,0.05) # les enfants restent chez eux
    beta.c["T","T"] = max(0.30*ecoles,0.15) # les ados sortent un peu
    beta.c["A","A"] = 0.90 + (0.10*0.30*ecoles) # Peu d'adultes a l'uni
    
    # effet télétravail
    beta.c["A","A"] = beta.c["A","A"]*(1-teletravail)

    # Ephad (environ 15% de seniors)
    beta.c[,"S"] = beta.c[,"S"]*(1-0.15*(1-ephad))
    beta.c["S",] = beta.c["S",]*(1-0.15*(1-ephad))
    
    # effet du social (hypothèse : loisirs = 1/4 du temps)
    beta.c = beta.c * (1-0.25*(1-social))
    
    # effet des magasins
    beta.c = beta.c * max(commerce,0.95) # représente 10% des interactions (hypothèse)
    
    # effet distanciation sociale (avec coeff pour relativiser l'effet)
    beta.c = beta.c * (1+0.4*(1/distance))

    return(beta.c)
}

# Initialisation de beta.c
ecoles = 0         # Les écoles sont fermées jusqu'à la fin de l'année
teletravail = 0.25 # Petit pourcentage toujours en télétravail
social = 0.75      # Evenements/sports/loisirs réduits
commerce = 1       # Réouverture des commerces non-essentiels
ephad = 0.75       # visites en ephad réduits
distance = 2       # distanciation sociale de 2m 

mesures = c(ecoles, teletravail, social, commerce, ephad, distance)
beta.c = matrice_contact(mesures)

# On déduit beta.mu a partir de beta.c et beta_post
beta.mu = beta_post/beta.c

Les hypothèses que nous avons faites et les variables utilisées sont expliquées ci-dessous :

Rappel : le but est de voir des variations relatives, les calculs de chaque coefficient sont aussi précis que possible mais nous ne prétendons pas à l’exactitude.

Exemples de mesures

A présent, nous pouvons modéliser l’effet relatif des comportements sur la dynamique de l’épidémie en post-confinement. Procédons ici à un exemple ou l’on fait varier seulement la distanciation sociale. Imaginons que la distance de 2m n’est pas été respectée, que se passe-t-il pour la dynamique épidémique ?

mesures = c(0, 0.25, 0.75, 1, 0.75, 1) # distance = 1m
beta.c.ex = matrice_contact(mesures)
beta_post_ex = beta.mu*beta.c.ex
beta_post_ex = c(beta_post_ex)

# Nouveaux paramètres
P0_post_ex = c(alpha, beta_post_ex, gamma, rho)


# resolution
mod_post_ex = ode(X0_post, tp, SIHR, P0_post_ex)
Graphes(mod_conf, mod_post_ex)

Nous constatons que les courbes montent toutes en flêche. Nous constatons également que les taux de contact ont une très grande influence sur la dynamique épidémique. Mais surtout, qu’assouplir les mesures peut vite amener à un emballement de l’épidémie alors que l’inverse (durcir les mesures, réduire les contacts) provoque une amélioration lente.

De plus, on remarque qu’en modifiant seulement un paramètre global : la distanciation physique, toutes les classes d’âge sont affectées mais surtout les Adultes et Seniors, pour lesquelles les hospitalisations doublent (au minimum). Cela montre que les plus vulnérables sont les plus âgés et que des mesures globales, mêmes si elles ne changent pas beaucoup les hospitalisations des plus jeunes, peuvent améliorer la situation des plus âgés.

Un autre exemple à présent : qu’aurait donné un déconfinement abrupt, sans mesures spécifiques mises en place ? Regardons cela :

Ces représentations montrent des dynamiques exponentielles pour la progression de l’épidémie qui sont irréalistes. Cela est dû à nos hypothèses, mais nous voyons bien un effet relatif : les conséquences sont beaucoup plus importantes que pour l’exemple précédent avec une augmentation des admissions et des hospitalisations. Elles justifient un déconfinement “en 3 phases”. Il est facile de faire le parallèle entre cette situation et celle de la remontée des cas avant le second confinement, auquel notre modèle pourrait être appliqué.

Ce travail montre bien pourquoi le confinement total est nécessaire pour enrayer l’épidémie : c’est le seul moyen de limiter assez les contacts pour être sûr que le nombre de malades diminue. D’autres pistes à explorer avec ce modèle sont son application au second confinement et l’ajout d’un état D pour décédés afin de voir les conséquences supplémentaires sur les décès des mesures sanitaires.

Merci d’avoir lu ce rendu de projet. N’oubliez pas regarder l’application shiny qui permet de crééer vos propres scénarios de déconfinement. Alors Sherlock, sauriez-vous enrayer la maladie ?

Références