Petit cours d'épidémiologie mathématique
Le modèle de Kermack et McKendrick

Julien Arino

Department of Mathematics & Data Science Nexus
University of Manitoba*

Centre canadien de modélisation des maladies (CCMM)
NSERC-PHAC EID Modelling Consortium (CANMOD, MfPH, OMNI/RÉUNIS)

Plan de ce cours

• Modèle SIR épidémique de Kermack et McKendrick
• Taille finale d'une épidémie
• Dernières remarques

Modèle SIR épidémique de Kermack et McKendrick

Ceci est un cas particulier !

La série de papiers vaut le coup

• Kermack & McKendrick. A contribution to the mathematical theory of epidemics (1927)

Suivi de "Contributions to the mathematical theory of epidemics."

• II. The problem of endemicity (1932)
• III. Further studies of the problem of endemicity (1933)
• IV. Analysis of experimental epidemics of the virus disease mouse ectromelia (1937)
• V. Analysis of experimental epidemics of mouse-typhoid; a bacterial disease conferring incomplete immunity (1939)

Quelle est la taille d'une épidémie?

• Si l'on s'intéresse à la possibilité d'un pic épidémique
  • A-t-il toujours lieu?
  • S'il a lieu, quelle est son intensité?
• Si une épidémie traverse une population, est-ce que tout le monde est infecté?

Le modèle SIR sans démographie

• La période de temps que l'on considère est suffisamment courte que l'on peut négliger la démographie (on parle aussi de modèle sans dynamique vitale)
• Différent du modèle précédent, qui incorpore la démographie mais la considère constante
• Les individus dans la populations sont soit susceptibles, soit infectieux avec la maladie. Après guérison ou mort, on les retire du compartiment infectieux ()
• Incidence du type action de masse

Le modèle que nous considérons maintenant est typiquement appelé modèle SIR de Kermack-McKendrick (KMK)

Reduction du problème

3 compartiments, mais quand on inspecte en détail, on remarque que les retirés n'ont pas d'influence directe sur la dynamique de ou , dans le sens où n'apparaît pas dans ou

De plus, la population totale (incluant les morts qui sont aussi classés dans ) satisfait

Par conséquent, est constant et la dynamique de peut être déduite de
Donc on considère à présent

Équilibres

Considérons les équilibres de

De

• soit
• ou

Substituant dans

• dans le premier cas,
• dans le second cas, n'importe quel est un équilibre

Le second cas est un problème: la linéarisation usuelle ne fonctionne pas lorsque l'on a un continuum d'équilibres, puisque les équilibres ne sont pas isolés!

Une autre approche - On étudie

Quelle est la dynamique de ?

pourvu que

Attention! Il convient de se souvenir que et sont et .. l'équation décrit ainsi la relation entre et le long de solutions de l'EDO d'origine -

On intègre sans aucun mal et on obtient des trajectoires dans l'espace d'état

avec

CI , et la solution de - est, en tant que fonction de

(puisque )

Trajectoires du système normalisé dans l'espace des phases avec CI et

Étudions

On a

Donc dans les courbes précédentes, le max de a lieu lorsque ( dans l'exemple). En ce point,

rhs_SIR_KMK <- function(t, x, p) {
  with(as.list(c(x, p)), {
    dS = - beta * S * I
    dI = beta * S * I - gamma * I
    dR = gamma * I
    return(list(c(dS, dI, dR)))
  })
}
# Condition initiale pour S (pour calculer R_0)
S0 = 1000
gamma = 1/14
# On règle beta pour que R_0 = 1.5
beta = 1.5 * gamma / S0 
params = list(gamma = gamma, beta = beta)
IC = c(S = S0, I = 1, R = 0)
times = seq(0, 365, 1)
sol <- ode(IC, times, rhs_SIR_KMK, params)

plot(sol[, "time"], sol[, "I"], type = "l",
     main = TeX("SIR de Kermack et McKendrick, $R_0=1.5$"),
     xlab = "Temps (jours)", ylab = "Prévalence")

Taille finale d'une épidémie

Calcul de la taille finale de l'épidémie

Pour la fonction à valeurs positives et intégrable , notons

On note aussi . Dans le sous-système

calculons la somme de et , en prenant garde de bien faire apparaître la dépendance en le temps

Intégrons de 0 à :

Le terme de gauche donne

car

Le terme de droite prend la forme

On a donc

Considérons maintenant :

Divisons les deux cotés par :

Intégrons de 0 a :

Exprimons et en termes de et égalisons

On a donc

L'équation (transcendante) de la taille finale

que l'on écrit

On cherche par conséquent les zéros de la fonction

On cherche dans t.q. pour

Notons pour commencer que lorsque

Différentiant par rapport à , on obtient

Lorsque , , donc décroit jusqu'à

Donc si , la fonction est décroissante sur , tandis qu'elle a un minimum si

Cas

• On a vu que est décroissante sur

• Par ailleurs, (le cas est exclu car trivial)

• Du reste, est continue

il existe un unique dans t.q.

Cas

• On a vu que est décroissante sur

• Pour ,

• Comme précédemment,

• est continue

il existe un unique dans t.q. . Plus précisemment, dans ce cas,

On va résoudre numériquement. Il faut une fonction

final_size_eq = function(S_inf, S0 = 999, I0 = 1, R_0 = 2.5) {
  OUT = S0*(log(S0)-log(S_inf)) - (S0+I0-S_inf)*R_0
  return(OUT)
}

et on résoud facilement en utilisant uniroot, ici avec les valeurs par défaut qu'on a mis dans la fonction:

uniroot(f = final_size_eq, interval = c(0.05, 999))
$root
[1] 106.8819
$f.root
[1] -2.649285e-07
$iter
[1] 10
$init.it
[1] NA
$estim.prec
[1] 6.103516e-05

Pour passer des valeurs différentes, on peut par exemple faire

N0 = 1000
I0 = 1
S0 = N0-I0
R_0 = 2.4
uniroot(
  f = function(x) 
    final_size_eq(S_inf = x, 
                  S0 = S0, I0 = I0, 
                  R_0 = R_0),
  interval = c(0.05, S0))

Une figure avec toutes les informations

S = seq(0.1, S0, by = 0.1)
fs = final_size(S, S0 = S0, I0 = I0, R_0 = R_0)
S_inf = uniroot(f = function(x) final_size_eq(S_inf = x, 
                                              S0 = S0, I0 = I0, 
                                              R_0 = R_0),
                interval = c(0.05, S0))
plot(S, fs, type = "l", ylab = "Valeur de l'équation (8)")
abline(h = 0)
points(x = S_inf$root, y = 0, pch = 19)
text(x = S_inf$root, y = 0, labels = "S_inf", adj = c(-0.25,-1))

Un peu plus joli

valeurs = expand.grid(
  R_0 = seq(0.01, 3, by = 0.01),
  I0 = 1:100
)
valeurs$S0 = N0-valeurs$I0
L = split(valeurs, 1:nrow(valeurs))

valeurs$S_inf = sapply(X = L, FUN = final_size)

valeurs$taille_finale = valeurs$S0-valeurs$S_inf+valeurs$I0
valeurs$taux_attaque = (valeurs$taille_finale / N0)*100

levelplot(taux_attaque ~ R_0*I0, data = valeurs, 
          xlab="R_0", ylab = "I0",
          col.regions = viridis(100))

(requiert les librairies lattice et viridis)

Taux d'attaque (en %)

Dernières remarques

Normaliser ou pas?

• Dans le SIS du Cours 04 et ici, puisque la population totale est constante, on peut normaliser à
• Cela simplifie pas mal certains calculs
• Toutefois, je ne suis pas très fan: il est important de toujours garder en tête les grandeurs biologiques
• Si vous normalisez, en tout cas, pour un papier à visée biomathématique, pensez à exprimer vos résultats en grandeur réelle

Là où nous en sommes

• Un modèle SIS endémique dans lequel le seuil est t.q. quand , la maladie s'éteint, tandis que lorsque , la maladie s'établit dans la population
• Un modèle SIR épidémique (le KMK SIR) dans lequel la présence ou absence de vague épidémique est caractérisée par la valeur de
• Le SIS et le KMK SIR sont intégrables dans un certain sens. C'est une exception!!!