Lors d'une épidémie, un paramètre fondamental à mesurer est le
taux de reproduction de base, noté $R_0$, défini par
Dans les modèles épidémiologiques compartimentaux, les individus sont divisés en classes épidémiologiques, comme les individus susceptibles, ceux infectés ou ceux rétablis.
Le modèle SIR comporte trois classes: Suceptibles, Infectés, Rétablis. Les susceptibles ($S$) peuvent être infectés par contact avec un individu infecté ($I$) et peut guérir et être immunisé ($R$). Les équations suivantes peuvent décrire ce modèle $$S(t+1) = S(t) - \beta S(t) \frac{I(t)}{N},$$ $$I(t+1) = I(t) + \beta S(t) \frac{I(t)}{N} - \gamma I(t),$$ $$R(t+1) = R(t) + \gamma I(t).$$ La probabilité pour un susceptible d'être infecté durant un pas de temps est égale à la probabilité $c$ d'entrer en contact avec un individu fois la probabilité $f$ que l'individu soit infecté fois la probabilité $p$ que l'infection se propage lors de ce contact. (On néglige la probabilité d'avoir plusieurs contact durant un pas de temps. Pour un pas de temps suffisament petit, cette probabilité devient effectivement négligeable.) Le nombre moyen de nouveau cas d'infection est donc $$S c f p = \beta S \frac{I}{N},$$ où $f=I/N$ est la probabilité d'interagir avec un infecté et $\beta = cp$ est la probabilité de contact et d'infection. Le paramètre $\gamma$ est la probabilité de guérison d'un infecté durant un pas de temps. Dans ce modèle, le nombre total d'individus $S+I+R = N$ est constant. Le modèle possède un équilibre sans infection $\bar S = N, \bar I = 0, \bar R = 0$. (L'équilibre est dégénéré, il existe pour tout $N$ positif.)
Quel est le taux de reproduction de base $R_0$ pour ce modèle ? Pour le calculer, on se place dans le voisinage de l'équilibre sans infection et on regarde le système linéaire autour de cet équilibre. La partie linéaire du modèle SIR est donné par la matrice jacobienne à l'équilibre $$ J = \left( \begin{array}{ccc} 1 & -\beta & 0 \\ 0 & 1 + \beta - \gamma & 0 \\ 0 & \gamma & 1 \end{array} \right) $$
(26/10) Le terme $J_{1,2}$ avait une erreur de signe. Elle est maintenant corrigée. Ce terme n'a pas d'effet sur le reste de l'analyse.
La stabilité de l'équilibre $(N,0,0)^*$ est donnée par la plus grande valeur propre en valeur absolue de $J$, le rayon spectral $\rho(J)$. Si le rayon spectral est strictement inférieur à 1, l'équilibre sera localement stable, c.-à-d. qu'il existe un voisinage de l'équilibre tel que pour toute condition initiale dans ce voisinage, les solutions du système convergent vers l'équilibre. Le polynôme caractéristique de $J$ est $$p_J(z) = (z - 1) (z - 1 - \beta + \gamma) (z - 1).$$ Deux des trois valeurs propres valent 1. La troisième détermine la stabilité: si $\lambda = 1 + \beta - \gamma > 1$, l'équilibre sans infection est instable, sinon il est stable. La condition sur les paramètres est donc $ \beta > \gamma$ ou encore $\beta/\gamma > 1$. Il se trouve que $\beta$ est la probabilité d'un contact menant à une infection par pas de temps et $1/\gamma$ est le nombre moyen de pas de temps que dure l'infection. Le premier individu infecté aura donc en moyenne $\beta/\gamma$ contacts menant à une infection, ou encore infectera $\beta/\gamma$ individus (tous les autres individus sont susceptibles). C'est exactement la définition de $R_0$, et donc $R_0 = \beta/\gamma$.
Dans les modèles plus complexes, le nombre de classes d'infection peut être élevé, par exemple on peut introduire une classe d'exposés (infectés mais pas encore contagieux). Dans ce cas le calcul de $R_0$ est moins évident. La méthode de la matrice de la génération suivante (next-generation matrix) peut être utilisée. Dans ce DM nous allons montrer comment calculer $R_0$ avec cette méthode et l'appliquer à un modèle épidémiologique compartimental
Soit un vecteur d'état $X \in \mathbb{R}^n$, où les $m$ premiers coefficients sont les nombres (ou densités) d'individus dans les classes infectées et les $n-m$ derniers coefficients sont les nombres d'individus dans les classes non-infectées. Le vecteur d'état s'exprime comme $X = (X_0^*, X_1^*)^*$, avec $X_0 \in \mathbb{R}^m$ les individus infectés et $X_1 \in \mathbb{R}^{n-m}$ les individus non-infectés. Le vecteur d'état $X$ évolue en fonction du temps, $X(t)$ où le temps $t$ est discret $t \in \mathbb{N}$. On suppose que les populations au temps $t+1$ ne dépendent que de l'état au temps $t$ $$X_0(t+1) = G_0\bigl(X(t)\bigr),$$ $$X_1(t+1) = G_1\bigl(X(t)\bigr).$$
Les fonctions $G_0: \mathbb{R}^n \to \mathbb{R}^m$ et $G_1: \mathbb{R}^n \to \mathbb{R}^{n-m}$ sont continument différentiable. (On a $X(t+1) = G(X(t))$ avec $G(X)=(G_0(X)^*,G_1(X)^*)^*$.) On suppose qu'il existe un unique état d'équilibre sans infection. Un état $\bar X$ est un état d'équilibre si $\bar X = G(\bar X)$. L'état sans infection est $\bar X = (0^*, \bar X_1^*)^*$, où $\bar X_1$ est non-nul.
Dans le voisinage de l'équilibre sans infection, la partie linéaire du système détermine le comportement des solutions. La partie linéaire est donnée par la matrice jacobienne $J(X)$ du système évaluée à l'équilibre. Les coefficients de $J$ sont $$ J_{ij} = \frac{\partial G^{(i)}}{\partial X^{(j)}}, $$ pour $i,j \in \{1,...,n\}$. A l'équilibre, $J(\bar X)$ s'exprime comme une matrice triangulaire par bloc. $$ J = \left( \begin{array}{cc} F + T & O \\ A & C \end{array} \right). $$ La matrice $O$ est la matrice de zéros, $C$ est une matrice $n-m \times n-m$ et les matrices $F$ et $T$ sont des matrices $m \times m$. La matrice $F$ est la matrice de fertilité et $T$ est la matrice de transition. La matrice de fertilité inclut les nouveaux cas infectés qui restent infectés durant un pas de temps. La matrice de transition décrit tous les autres changements (guérison, passage d'exposé à contagieux, etc). La matrice $C$ décrit l'évolution des non-infectés, en l'absence d'infectés.
On suppose qu'en l'absence d'infection, l'état d'équilibre sans infection est stable, ce qui implique que $\rho(C)<1$. Par hypothèse de modélisation, on peut aussi supposer que le rayon spectral de $T$ est inférieur à 1. On suppose enfin que la matrice $F+T$ possède une valeur propre dominante positive, c.-à-d. que le rayon spectral de $F+T$ est la plus grande valeur propre de $F+T$.
Question 1 Montrez que le rayon spectral de $J$ à l'équilibre est inférieur à 1 si et seulement si le rayon spectral de $F+T$ est inférieur à 1.
Il s'agit ici de noter que les que les valeurs propres de $J$ sont données par les valeurs propres de $F+T$ et les valeurs propres de $C$ (pourquoi?). Utilisez $\rho(C)$ pour conclure.
(26/10) Les valeurs propres de la matrice jacobienne $J$ ne se trouvent pas nécessairement sur la diagonale de $J$. Pour le calcul des valeurs propres, il faut utiliser le fait que la matrice jacobienne pour les modèle épidémiologiques, évaluée à l’équilibre sans infection, possède un bloc de zéro dans la partie supérieure droite. La forme générale de la matrice est donc triangulaire par bloc. Ceci n’implique pas que les valeurs propres soient sur la diagonale.
Question 2 Montrez que la matrice de la génération suivante est bien définie, c.-à-d. que la matrice $I-T$ est inversible.
Question 3 Montrez que la matrice de la génération suivante peut aussi s'écrire $$Q = F ( I + T + T^2 + T^3 + ...).$$
Il faut montrer que $(I-T)^{-1} = \sum_{k=0}^{\infty} T^k$ si le rayon spectral $\rho(T) < 1$. La démonstration suit la démonstration de la somme d'une série géométrique. Utilisez le fait que $T^k \to 0$ lorsque $k \to \infty$ si $\rho(T) < 1$.
Le taux de reproduction de base $R_0$ est défini comme le rayon spectral de la matrice de la génération suivante: $R_0 = \rho(Q)$.
Question 4 Soit $r=\rho(F+T)$ et $R_0 = \rho(Q)$. Montrez les inégalités suivantes:
Cette question est difficile.
(28/10) Quelques précisions sur les hypothèses à utiliser pour cette question. Vous n'avez pas à démontrer ces hypothèses.
Il existe donc deux vecteurs non-négatifs $x$ et $y$, tels que $$(F+T)x = rx,$$ $$Qy = R_0 y.$$
Notez que si $y$ est vecteur propre de $Q$, associé à la valeur propre $R_0$, on a $Qy = F(I-T)^{-1} y = R_0 y$. Posez $u = (I-T)^{-1}y$ pour obtenir $F u = R_0 (I-T) u$, ou encore $(F + R_0 T) u = R_0 u$. La valeur $R_0$ est donc valeur propre de la matrice $F + R_0 T$. Que se passe-t-il quand $R_0 = 1$? Quand $R_0 < 1$, $R_0 > 1$?
L’important ici est de considérer la matrice $F + \lambda T$, pour $\lambda$ réel, comme une famille de matrices indicée sur $\lambda$. Intuitivement, vous devrez démontrer que la norme de $(F+\lambda T)x$ est croissante, mais avec pente inférieure à 1, par rapport à $\lambda$. On cherche un point fixe $||(F+\lambda T)u|| = \lambda||u||$ en fonction de $\lambda$, avec norme du vecteur propre $||u||=1$.
Quand $\lambda = 1$, la matrice devient $F+T$. C’est exactement la matrice dont on connait la plus grande valeur propre. Vous devez montrer en 1. que $r = 1$ si et seulement si la matrice $F + \lambda T$ possède aussi une valeur propre égale à 1 lorsque $\lambda = 1$. Maintenant, que se passe-t-il lorsque $\lambda > 1$? Autrement dit, comment la valeur propre maximale de $F + \lambda T$ varie-t-elle lorsque $\lambda$ varie? Si $x$ est le vecteur propre de norme unité de $F + T$, associé à la valeur propre $r$, on a par définition $(F+T)x = rx$. Comparez les vecteurs $(F+T)x = rx$ et $(F+\lambda T)x = (F + T)x + (\lambda-1) Tx = rx + (\lambda -1)Tx$. Prenez leurs normes et utilisez l’inégalité du triangle pour montrer que quand $\lambda$ augmente, la norme est croissante et qu’il existe une valeur critique $\lambda = R_0 > r$ telle que $(F+R_0 T) u = R_0 u$.
La démonstration pour $r < 1$ est la même mais dans l’autre sens.
Notez qu'il n'est pas nécessaire d'utiliser directement les notions de rayon spectral, ou de norme de matrice. Les hypothèses de travail sont suffisantes pour n'avoir à considérer que la valeur propre maximale et des normes de vecteurs.
On considère le modèle SEIR suivant: $$S(t+1) = S(t) - \beta S(t) \frac{I(t)}{N} + \mu N - \mu S(t),$$ $$E(t+1) = E(t) + \beta S(t) \frac{I(t)}{N} - (\mu + a) E(t),$$ $$I(t+1) = I(t) + a E(t) - (\gamma + \mu) I(t),$$ $$R(t+1) = R(t) + \gamma I(t) - \mu R(t).$$ En plus des paramètres et variables du modèle SIR, on a $a$ la probabilité de passer de la classe exposée (E) à infectieux (I), et $\mu$ une probabilité de mort naturel dans la population. Un nombre de naissance $\mu N$ par pas de temps est rajouté chez les susceptibles.
Question 5 Vérifiez que l'état d'équilibre sans infection $(\bar S, 0, 0, 0)^*$ existe. Que vaut $\bar S$?
Question 6 Déterminez la matrice jacobienne $J$ à l'état d'équilibre sans infection et déterminez les matrices $F, T, A$ et $C$. (La matrice $F$ ne tient compte que des nouveaux cas exposés (E) ou contagieux (I).) Calculez $r = \rho(F+T)$.
La matrice $F+T$ est une matrice 2 par 2, le rayon spectral est la plus grande des deux valeurs propres
Question 7 Déterminez la matrice de la génération suivante et calculez $R_0$. Vérifiez que $R_0 = 1$ si et seulement si $r=1$.
Question 8 Déterminez l'existence d'un état d'équilibre avec infection, c.-à-d. un état endémique: $(\tilde S, \tilde E, \tilde I, \tilde R)^*$, avec coefficients positifs. Exprimez cet état en terme de $R_0$.
Question 9 - Simulations numériques En utilisant Python, ou un autre outil de votre choix, écrivez une routine qui prend en entrée un vecteur d'état initial $X(0)$ et un temps $T$ final et qui retourne $X(t)$ pour $t = 1,...,T$. Vous pouvez par exemple utiliser le prototype suivant.
def seir(X0,T,beta,mu,a,gamma,N):
"""Solution d'un modèle SEIR discret
Usage:
X = seir(X0,T,beta,mu,a,gamma,N)
Arguments:
X0: vecteur 4x1 (S,E,I,R)^*
T: temps final
beta: probabilité d'un contact menant à une infection
mu: probabilité de mort naturelle
a: probabilité pour un exposé de devenir infecté
gamma: probabilité pour un infecté de se rétablir
N: population totale
Retour:
X: matrice 4x(T+1) solution
R0: taux de reproduction de base
"""
code...
return X, R0
Question 10 - Simulations numériques Pour des valeurs de paramètres adéquatement choisies, effectuez des simulations et représentez graphiquement les solutions.
Références: Compartmental models in epidemiology.