Comment modéliser la propagation d’une épidémie

Introduction

A la différence des maladies chroniques non transmissibles (obésité, diabète, cancer…), les maladies infectieuses sont causées par un agent pathogène (bactérie, virus, parasite, champignon…).

La propagation de cet agent est un phénomène dynamique. Le nombre d’individus sains, malades, immunisés évolue au cours du temps selon qu’un individu change de statut de non malade à infecté et infectant d’autres individus.

Tout cela peut être modélisé par des équations différentielles qui caractérisent l’évolution.

Le modèle compartimental SIR

Le modèle le plus simple est le modèle SIR développé par Kermack et McKendrick en 1927. C’est un modèle déterministe. Il est caractérisé par 3 populations d’individus :

S pour Sains/susceptibles d’être infectés
I pour Infectés
R les personnes Rétablies (Recovered = guéries). Les personnes rétablies ne peuvent plus être ré-infectées (elles sont immunisées).

d’où S(t) + I(t) + R(t) = P

Ce système peut être représenté par des boites (compartiments) traverseés par des flux entrants et sortants.

modèle SIR Par exemple dans le compartiment des individus infectés, on a des nouveaux entrants venant des individus sains qui deviennent infectés à un taux d’incidence β (infection) et on a des individus sortant qui sont guéris au bout d’un temps λ (temps de guérison) dans une proportion I/λ.

SIR <- function(time, state, parameters) {
  par <- as.list(c(state, parameters))
  with(par, {
    dS <- -beta * I * S / N
    dI <- beta * I * S / N - gamma * I
    dR <- gamma * I
    list(c(dS, dI, dR))
  })
}

Adapter le modèle à nos données.

Données

La source des données:coronavirus, R package developé par Rami Krispin.

devtools::install_github("RamiKrispin/coronavirus")
library(coronavirus)
data(coronavirus)
`%>%` <- magrittr::`%>%`

# extract the cumulative incidence
df <- coronavirus %>%
  dplyr::filter(Country.Region == "France")%>%
  dplyr::group_by(date, type) %>%
  dplyr::summarise(total = sum(cases, na.rm = TRUE)) %>%
  tidyr::pivot_wider(
    names_from = type,
    values_from = total
  ) %>%
  dplyr::arrange(date) %>%
  dplyr::ungroup() %>%
  dplyr::mutate(active = confirmed - death - recovered) %>%
  dplyr::mutate(
    confirmed_cum = cumsum(confirmed),
    death_cum = cumsum(death),
    recovered_cum = cumsum(recovered),
    active_cum = cumsum(active)
  )
df

## # A tibble: 87 x 9
##    date       confirmed death recovered active confirmed_cum death_cum
##    <date>         <int> <int>     <int>  <int>         <int>     <int>
##  1 2020-01-22         0     0         0      0             0         0
##  2 2020-01-23         0     0         0      0             0         0
##  3 2020-01-24         2     0         0      2             2         0
##  4 2020-01-25         1     0         0      1             3         0
##  5 2020-01-26         0     0         0      0             3         0
##  6 2020-01-27         0     0         0      0             3         0
##  7 2020-01-28         1     0         0      1             4         0
##  8 2020-01-29         1     0         0      1             5         0
##  9 2020-01-30         0     0         0      0             5         0
## 10 2020-01-31         0     0         0      0             5         0
## # ... with 77 more rows, and 2 more variables: recovered_cum <int>,
## #   active_cum <int>

Visualisation des donneés

library(tidyverse)
library(tidyr)
library(highcharter)
df %>% 
  select(c(date, confirmed_cum, death_cum, recovered_cum))%>% 
  pivot_longer(-date, 
               names_to = "cases",
               values_to = "cumul")%>%
  hchart(type="line",hcaes(x=date,y=cumul,group=cases)) %>% 
  hc_add_theme(hc_theme_ffx()) %>% hc_title(text="Cumul des cas COVID 19 par jour en France") %>% 
  hc_subtitle(text="source:Rami Krispin")

Pour adapter le modèle aux données, deux fonctions sont utilisées:

La fonction ode () du package deSolve qui facilite la résolution des équations différentielles ordinaires et la fonction optim () pour trouver les valeurs optimales pour nos deux paramètres inconnus β, γ(paramètres que nous souhaitons estimer).

Plus précisément, ce que nous souhaitons faire est de minimiser la somme des différences au carré RSS entre I(t), qui est le nombre de personnes dans le compartiment I au temps t, et le nombre correspondant prédit (estimé) par notre modèle ^I(t).

#vecteur: Infected
#Periode choisie; 20 février au 20 mars 
library(lubridate)
Infected <- subset(df, date >= ymd("2020-02-20") & date <= ymd("2020-03-20"))$active_cum
# vecteur; Jour de la même longueur que notre
# vecteur de cas
Day <- 1:(length(Infected))
### Valeurs  des initiales
N <- 6200000
init <- c(
  S = N - Infected[1],
  I = Infected[1],
  R = 0
)

Ajustement du modèle

RSS <- function(parameters) {
  names(parameters) <- c("beta", "gamma")
  out <- ode(y = init, times = Day, func = SIR, parms = parameters)
  fit <- out[, 3]
  sum((Infected - fit)^2)
}

Enfin, le modèle SIR est prêt à être testé à nos données pour trouver les valeurs β, γ qui minimisent la somme résiduelle des carrés entre le nombre de cas observé et le nombre de cas estimé. Nous devons également vérifier que notre modèle a bien convergé.

## [1] "ERROR: ABNORMAL_TERMINATION_IN_LNSRCH"

Calcul des paramétres β, γ :

Opt_par <- setNames(Opt$par, c("beta", "gamma"))
Opt_par

##      beta     gamma 
## 0.6314203 0.3685817

N’oubliez pas que β contrôle la transition entre S et I, et γ contrôle la transition entre I et R. Ces valeurs ne signifient pas grand-chose, mais elles nous sont utiles pour calculer le nombre de cas ajusté dans les compartiments de notre modèle SIR.

sir_start_date <- "2020-02-20"
t <- 1:as.integer(ymd("2020-03-21") - ymd(sir_start_date))
# obtenir les valeurs ajustées de notre modèle 
fitted_cumulative_incidence <- data.frame(ode(
  y = init, times = t,
  func = SIR, parms = Opt_par
))

# ajouter une colonne date et les données observées
library(dplyr)
fitted_cumulative_incidence <- fitted_cumulative_incidence %>%
  mutate(
    Date = ymd(sir_start_date) + days(t - 1),
    Country = "France",
    cumulative_incident_cases = Infected
  )
head(fitted_cumulative_incidence)

##   time       S         I         R       Date Country cumulative_incident_cases
## 1    1 6199993  7.000000  0.000000 2020-02-20  France                         7
## 2    2 6199988  9.104309  2.950908 2020-02-21  France                         7
## 3    3 6199981 11.841201  6.788905 2020-02-22  France                         7
## 4    4 6199973 15.400829 11.780658 2020-02-23  France                         7
## 5    5 6199962 20.030513 18.272999 2020-02-24  France                         7
## 6    6 6199947 26.051917 26.717033 2020-02-25  France                         2

# Visualisation
library(ggplot2)
 fitted_cumulative_incidence%>%
  ggplot(aes(x = Date)) +
  geom_line(aes(y = I), colour = "red") +
  geom_point(aes(y = cumulative_incident_cases), colour = "blue") +
  labs(
    title = "Estimé vs Observé, France",
    y = "Nombre de cas ") +
  theme_minimal()

Les deux tendances se chevauchent; cela indique que la pandémie est clairement dans une phase exponentielle.

fitted_cumulative_incidence %>%
  ggplot(aes(x = Date)) +
  geom_line(aes(y = I), colour = "red") +
  geom_point(aes(y = cumulative_incident_cases), colour = "blue") +
  labs(
    y = "Nombre de cas",
    title = "Estimé vs Observé, France"
   
  ) +
  theme_minimal() +
  scale_y_log10(labels = scales::comma)

Le graphique suivant est similaire au précédent, sauf que l’axe des y est mesuré sur une échelle logarithmique.

jusqu’au début mars, le nombre de cas observé est resté inférieur par rapport au nombre éstimé dans une phase exponentielle.
le nombre de cas observé est resté constant (moins de 10 cas) du début février jusqu’au 24 février.
du début mars au 30 mars, le nombre de cas observé a continué d’augmenter à un rythme proche d’un taux exponentiel.

Le taux de reproduction de base R0

Le taux de reproduction de base R0 constitue le nombre de cas secondaires produits par un individu infectieux moyen au cours de sa période d’infectiosité dans une population entièrement constituée de personnes sensibles (non immunisées).

Le R0 est un seuil * Si R0 > 1, l’épidémie peut s’installer * Si R0 < 1, il y a peu de chance d’avoir une épidémie

Les estimations du R0 sont souvent calculées en fonction de 3 paramètres principaux

la durée de contagiosité après qu’une personne soit infectée
la probabilité d’infection par contact entre une personne sensible et une personne ou un agent infectieux
le taux de contact – avec des paramètres supplémentaires qui peuvent être ajoutés pour décrire des cycles de transmission plus complexes

Le taux de reproduction de base du coronavirus en France

Opt_par

##      beta     gamma 
## 0.6314203 0.3685817

R0 <- as.numeric(Opt_par[1] / Opt_par[2])
R0

## [1] 1.713108

le risque d’être contaminé est de 1.7
la proportion de la population susceptible d’être infectée tout au long de la pandémie = 1-(1/R0)

Application du modèle

t <- 1:120
#obtenir les valeurs ajustées de notre modèle SIR
fitted_cumulative_incidence <- data.frame(ode(
  y = init, times = t,
  func = SIR, parms = Opt_par
))
#ajouter une colonne Date et joindre les données d'incidence observées

fitted_cumulative_incidence <- fitted_cumulative_incidence %>%
  mutate(
    Date = ymd(sir_start_date) + days(t - 1),
    Country = "France",
    cumulative_incident_cases = c(Infected, rep(NA, length(t) - length(Infected)))
  )

## plot the data
fitted_cumulative_incidence %>%
  ggplot(aes(x = Date)) +
  geom_line(aes(y = I), colour = "red") +
  geom_line(aes(y = S), colour = "black") +
  geom_line(aes(y = R), colour = "green") +
  geom_point(aes(y = cumulative_incident_cases),
    colour = "blue"
  ) +
  scale_y_continuous(labels = scales::comma) +
  labs(y = "Nombre de cas", title = "COVID-19 Estimé vs Obsrvé, France") +
  scale_colour_manual(name = "", values = c(
    red = "red", black = "black",
    green = "green", blue = "blue"
  ), labels = c(
    "Susceptible",
    "Recovered", "Observed", "Infectious"
  )) +
  theme_minimal()

## Warning: Removed 90 rows containing missing values (geom_point).

Le même graphique en échelle logarithmique pour l’axe des y.

fitted_cumulative_incidence %>%
  ggplot(aes(x = Date)) +
  geom_line(aes(y = I, colour = "red")) +
  geom_line(aes(y = S, colour = "black")) +
  geom_line(aes(y = R, colour = "green")) +
  geom_point(aes(y = cumulative_incident_cases, colour = "blue")) +
  scale_y_log10(labels = scales::comma) +
  labs(
    y = "Nombre de cas",
    title = "COVID-19 estimé vs observé, france"
  ) +
  scale_colour_manual(
    name = "",
    values = c(red = "red", black = "black", green = "green", blue = "blue"),
    labels = c("Susceptible", "Observed", "Recovered", "Infectious")
  ) +
  theme_minimal()

## Warning: Transformation introduced infinite values in continuous y-axis

## Warning: Removed 90 rows containing missing values (geom_point).

Résumé statistique

D’autres statistiques intéressantes peuvent être calculées à partir de l’ajustement de notre modèle. Par exemple:

la date du pic de la pandémie
le nombre de personnes ayant besoin de soins intensifs
le nombre de décès

fit <- fitted_cumulative_incidence

# pic de la pandémie
fit[fit$I == max(fit$I), c("Date", "I")]

##          Date        I
## 51 2020-04-10 631997.5

# le nombre max des infectés
max_infected <- max(fit$I)

#taux de mortalité supposé de 0,7%
max_infected * 0.007

## [1] 4423.982

Ce cours a pour but d’illustrer la structure d’un modèle épidémiologique compartimenté et de l’appliquer au données Covid 19 France. Nous avons utilisé des données de source chinoises et italiennes. Nous ignorons, cependant, dans quelles conditions et selon quels critères elles ont été recueillies. selon l’analyse de ces données, le pic en France a été atteint 10 Avril avec environ 600 000 personnes infectées, ce qui se traduit par 4 000 décès(en supposant un taux de mortalité de 0,7%).

Comment modéliser la propagation d’une épidémie

Covid 19 en France