Les modèles épidémiologiques

De gauche à droite : Daniel Bernoulli, Ronald Ross, William Ogilvy Kermack, Anderson Gray McKendrik et Richard Doll

\bullet Historique :

\begin{tabular}{|l|c|} \hline Bloc & 12/02/2024-V1 \\ \hline Historique & Cr\'eation \\ \hline Sommaire & Cr\'eation \\ \hline Pr\'esentation & Cr\'eation \\ \hline Les mod\`eles \'epid\'emiologiques & Cr\'eation \\ \hline Annexe th\'eorique & Cr\'eation \\ \hline Exemple & Cr\'eation \\ \hline Application sous R & Cr\'eation \\ \hline Application sous SAS & Cr\'eation \\ \hline Bibliographie & Cr\'eation \\ \hline \end{tabular}

\bullet Sommaire :

  • Présentation
  • Les modèles épidémiologiques
    • Du modèle SI au modèle SIRD
    • Le compartiment E et le modèle SEIRD
    • Le modèle métapopulationnel
    • Les effets naissances, migrations et décès naturels
    • Quelques autres compartiments classiques
    • Le nombre de reproduction effectif
  • Annexe théorique
    • Les Equations Différentielles Ordinaires
    • L’Algorithme d’Euler
    • L’Algorithme de Runge-Kutta
  • Exemple
    • Le modèle SEIRD
    • Le modèle SEIRD métapopulationnel
  • Application sous R
  • Application sous SAS
  • Bibliographie

\bullet Présentation :

Les modèles compartimentaux épidémiologiques sont une famille d’outils d’anticipation et de compréhension des dynamiques de propagation des maladies, jouant un rôle central dans le domaine de la Santé publique. Ils ont fait énormément parler d’eux au cours de la crise Covid-19 mondiale traversée de 2020 à 2022. Leur principe consiste à s’appuyer sur une liste de paramétrages contextuels pour ensuite proposer une série de projections sur un nombre fixé de temps.

Les modèles épidémiologiques ont évolué au fil du temps grâce aux contributions de nombreux chercheurs, mathématiciens et épidémiologistes. Néanmoins, certaines figures clés ont marqué des étapes cruciales dans leur développement. On pourra citer Daniel Bernoulli qui, au XVIIIème siècle, a jeté les bases de la modélisation des phénomènes aléatoires, élément clé des modèles épidémiologiques. Bases qui ont été ensuite développées par Ronald Ross en contribuant à la modélisation mathématique des maladies infectieuses, notamment le paludisme, éclairant la compréhension de sa transmission par les vecteurs. Parmi les autres célébrités à évoquer, et cet article pêchera dans l’exhaustivité de cette liste, on peut également parler de William Ogilvy Kermack et d’Anderson Gray McKendrick qui, en 1927, publient un article fondateur introduisant le modèle SIR précurseur des modèles épidémiologiques complexes modernes. Enfin, deux dernières personnalités emblématiques : Richard Doll et Austin Bradford Hill, qui ont joué un rôle important dans la formalisation des études épidémiologiques et notamment établi les principes cruciaux sur la causalité entre des facteurs spécifiques et des maladies.

Ces outils se basent sur la notion de système d’Equations Différentiels Ordinaires (EDO) et de leur résolution afin de déterminer à chaque itération les différents volumes associés aux compartiments considérés, permettant alors à tracer la courbe d’évolution de l’épidémie selon autant d’angles qu’il y a de compartiments.

\bullet Les modèles épidémiologiques :

Hypothèse préliminaire : Un jeu de paramètres en adéquation avec les différents compartiments utilisés.

Du modèle SI au modèle SIRD

Il faut définir au préalable une unité de temps : \Delta, et son maximum : T. De cette manière, on fixe déjà la temporalité  \forall t \in \lbrace 1, 1 + 1 \cdot \Delta, 1 + 2 \cdot \Delta, \cdots, 1 + T \cdot \Delta \rbrace sur laquelle la modélisation sera réalisée. Enfin, on note N la somme des effectifs de tous les compartiments utilisés, soit la population total étudiée.

Il y a deux compartiments classiques dans les modèles épidémiologiques : le S pour « Susceptibles », soit les susceptibles d’être contaminés car sains, et il ne faut pas confondre un individu sain avec un individu immunisé. Le I pour « Infectés », soit les individus sains qui ont fini par être contaminés. On fixe alors le paramètre \beta \in [0,1], probabilité qu’un individu fasse la transition : S \rightarrow I à chaque instant t. Ces deux premiers compartiments forment à eux-seuls le modèle SI, regroupant de façon simpliste l’ensemble des infections virales chroniques n’ayant aucune incidence sur la mortalité (on pourra citer par exemple l’herpès simplex dont on ne meurt et ne guérit pas). On peut exprimer les volumes/effectifs à l’instant t au sein de chacun de ces deux compartiments sous la forme généralisée :

S_t = S(t) = S_0 - \int_0 ^t \beta \cdot \frac{S(x) \cdot I(x)}{N} dx

I_t = I(t) = I_0 + \int_0 ^t (\beta \cdot \frac{S(x) \cdot I(x)}{N} - (\gamma + \delta) \cdot I(x)) dx

, à ce stade les paramètres \gamma, \delta sont supposés nuls, on y reviendra plus tard. On définit alors les états initiaux : S_0, I_0 \in \mathbb{R} ^+, et on érige le système SI sous la forme du système d’Equations Différentielles à résoudre suivant :

\forall t, \begin{cases} S_t = S_{t-1} + \partial S = S_{t-1} - \beta \cdot \frac{S_{t-1} \cdot I_{t-1}}{N} \cdot \Delta \\ I_t = I_{t-1} + \partial I = I_{t-1} + \beta \cdot \frac{S_{t-1} \cdot I_{t-1}}{N} \cdot \Delta \end{cases}

A noter que l’un des indicateurs phares utilisés pour le suivi des épidémies est celui du nombre ou incidence de nouveaux cas. On peut l’obtenir alors très facilement en l’extrayant de la formule, et de forme :

\beta \cdot \frac{S_{t-1} \cdot I_{t-1}}{N}

On enrichit le modèle en introduisant maintenant la notion de guérison grâce au compartiment R. Dès lors, on passe du phénomène assez rare d’infections virales chroniques sans réelle conséquence à n’importe quelles maladies dont on guérit mais ne meurt pas, spectre particulièrement large et dont on pourrait citer le rhume en guise d’exemple. L’ajout de ce compartiment se fait avec celui du paramètre \gamma \in [0,1] traduisant mathématiquement le taux de guérison. Le volume/effectif associé au compartiment R à l’instant t est alors :

R_t = R(t) = R_0 + \int_0 ^t \gamma \cdot I(x) dx

A noter que l’insertion de ce nouveau compartiment a également une incidence sur le compartiment I comme on peut le voir sur l’intégrale généralisée décrite précédemment. On définit alors les états initiaux : S_0, I_0, R_0 \in \mathbb{R} ^+, et on érige le système SIR sous la forme du système d’Equations Différentielles à résoudre suivant :

\forall t, \begin{cases} S_t = S_{t-1} + \partial S = S_{t-1} - \beta \cdot \frac{S_{t-1} \cdot I_{t-1}}{N} \cdot \Delta \\ I_t = I_{t-1} + \partial I = I_t + (\beta \cdot \frac{S_{t-1} \cdot I_{t-1}}{N} - \gamma \cdot I_{t - 1}) \cdot \Delta \\ R_t = R_{t-1} + \partial R = R_{t-1} + \gamma \cdot I_{t-1} \cdot \Delta \end{cases}

Enfin, on rajoute maintenant un dernier compartiment relatif aux décès : D. Il s’agit des personnes infectées par la maladie et qui en sont décédées. Ces modèles SIRD permettent ainsi d’intégrer toute maladie à risque endémique avec un potentiel risque de mortalité. Sa mise en place s’accompagne de l’ajout du paramètre \delta retranscrivant le taux de mortalité. Le volume/effectif associé à ce compartiment et à l’instant t est alors :

D_t = D(t) = D_0 + \int_0 ^t \delta \cdot I(x) dx

Tout comme pour le compartiment R, l’insertion du compartiment D a également une incidence sur le compartiment I comme on peut le voir sur l’intégrale généralisée décrite précédemment. On définit alors les états initiaux : S_0, I_0, R_0, D_0 \in \mathbb{R} ^+, et on érige le système SIR sous la forme du système d’Equations Différentielles à résoudre suivant :

\forall t, \begin{cases} S_t = S_{t-1} + \partial S = S_{t-1} - \beta \cdot \frac{S_{t-1} \cdot I_{t-1}}{N} \cdot \Delta \\ I_t = I_{t-1} + \partial I = I_t + (\beta \cdot \frac{S_{t-1} \cdot I_{t-1}}{N} - (\gamma + \delta) \cdot I_{t-1}) \cdot \Delta \\ R_t = R_{t-1} + \partial R = R_{t-1} + \gamma \cdot I_{t-1} \cdot \Delta \\ D_t = D_{t-1} + \partial D = D_t + I_{t-1} \cdot \delta \cdot \Delta \end{cases}

– Le compartiment E et le modèle SEIRD

Le principal défaut que l’on peut trouver au modèle SIRD, outre sa simplicité, est de ne pas prendre en compte le principe d’incubation d’une maladie. Il s’agit plus particulièrement du temps écoulé entre la contamination et le début des symptômes, soit, d’une certaine manière, la longueur de la flèche : S \rightarrow I. Ce qui vaut à ce principe un compartiment à part entière et noté : E, c’est qu’en plus de créer un décalage dans le temps de la propagation du virus étudié, il permet de prendre en compte le stade de progression de la maladie chez un individu infecté qui, méconnaissant son statut réel pourra en contaminer d’autres qui sont sains et sans s’en rendre compte. Dès lors, le terme d’Exposés est utilisé. Cet enrichissement permet aux modèles épidémiologiques de s’approcher toujours un peu plus des mécanismes des maladies à risque endémique et modifie sensiblement la forme généralisée du volume/effectif associée au compartiment I du fait de liaison imbriquée :

E_t = E(t) = E_0 + \int_0 ^t (\beta \cdot \frac{S(x) \cdot I(x)}{N} - \sigma \cdot E(x)) dx

I_t = I(t) = I_0 + \int_0 ^t (\sigma \cdot E(x) - (\gamma + \delta) \cdot I(x)) dx

, en définissant le nouveau paramètre \sigma = \frac{1}{\mbox{Temps d'incubation}}

On fixe les états initiaux : S_0, E_0, I_0, R_0, D_0 \in \mathbb{R} ^+, et on érige le système SEIRD sous la forme du système d’Equations Différentielles à résoudre suivant :

\forall t, \begin{cases} S_t = S_{t-1} + \partial S = S_{t-1} - \beta \cdot \frac{S_{t-1} \cdot I_{t-1}}{N} \cdot \Delta \\ E_t = E_{t-1} + \partial E = E_{t-1} + (\beta \cdot \frac{S_{t-1} \cdot I_{t-1}}{N} - \sigma \cdot E_{t-1}) \cdot \Delta \\ I_t = I_{t-1} + \partial I = I_{t-1} - (\sigma \cdot E_{t-1} - (\gamma + \delta) \cdot I_{t-1}) \cdot \Delta \\ R_t = R_{t-1} + \partial R = R_{t-1} + \gamma \cdot I_{t-1} \cdot \Delta \\ D_t = D_{t-1} + \partial D = D_{t-1} + I_{t-1} \cdot \delta \cdot \Delta \end{cases}

– Le modèle métapopulationnel

Une manière d’améliorer le modèle construit est d’intégrer dans son schéma d’Equations Différentielles l’information sur le taux de contact entre les différentes classes d’âge d’une population. Par exemple, si l’on suppose que l’épidémie que l’on cherche à modéliser se déroule dans un pays avec un moyenne d’âge particulièrement jeune, ajouter un paramètre qui permettra de pondérer le taux de transmission en fonction du taux de contact inter-classe d’âge pourra apporter un boost de précision très intéressant.

On doit donc définir ce l’on appelle une matrice de contacts de forme symétrique suivante :

\mathbf{C} = \begin{pmatrix}  & C_1 & \cdots & C_a  & \cdots & C_A \\ C_1 & c_{1,1} & \cdots & c_{1,a} & \cdots & c_{1,A} \\ \cdots & \cdots & \cdots & \cdots & \cdots & \cdots \\ C_a & c_{a,1} & \cdots & c_{a,a} & \cdots & C_{a,A} \\ \cdots & \cdots & \cdots & \cdots & \cdots & \cdots \\ C_A & c_{A,1} & \cdots & C_{A,a} & \cdots & C_{A,A} \\ \end{pmatrix}

, avec C_1, C_2, \cdots, C_A les différentes classes ou intervalles d’âge paramétrés et \forall a_1, a_2 \in \lbrace 1, \cdots, C_A \rbrace, c_{a_1, a_2} \in [0,1] le taux de contact entre les classes d’âge a_1 et a_2. Le modèle épidémiologique s’en voit fortement enrichi puisqu’il intègre dans son concept les probabilités de transiter d’une classe d’âge à l’autre. Si l’on prend par exemple le modèle SEIRD, on intègre ce concept ainsi :

\forall t, \begin{cases} S_t = S_{t-1} + \partial S = S_{t-1} - (\frac{\beta}{N} \cdot S_{t-1} \cdot (\mathbf{C} \cdot I_{t-1})) \cdot \Delta \\ E_t = E_{t-1} + \partial E = E_{t-1} + (\frac{\beta}{N} \cdot S_{t-1} \cdot (\mathbf{C} \cdot I_{t-1}) + \sigma \cdot E_{t-1}) \cdot \Delta \\ I_t = I_{t-1} + \partial I = I_{t-1} + (\sigma \cdot E_{t-1} - (\gamma + \delta) \cdot I_{t-1}) \cdot \Delta \\ R_t = R_{t-1} + \partial R = R_{t-1} + \gamma \cdot I_{t-1} \cdot \Delta \\ D_t = D_{t-1} + \partial D = D_{t-1} + I_{t-1} \cdot \delta \cdot \Delta \\ \end{cases}

, avec S_t, E_t, I_t, R_t et D_t des vecteurs de taille 1 \times A des volumes ventilés par classe d’âge.

– Les effets naissances, migrations et décès naturels

Les modèles vus précédemment sont plus souvent utilisés pour les épidémies ponctuelles et sont d’ailleurs assez flexibles pour procéder à des réajustements lors d’une propagation en cours ou même établir des scénarios servant à modéliser plusieurs vagues successives. Pour cela il suffit simplement d’adapter les paramètres initiaux à la situation rencontrée. A partir du moment où l’objectif devient la modélisation d’un phénomène s’allongeant dans le temps, il est important d’intégrer dans son modèle les notions de :

  • Nouvelles naissances, qui implique que la population croît régulièrement sur un profil bien spécifique venant ainsi alimenter uniquement le compartiment S. Cet effet sera noté \lambda ;
  • Effets migratoires, avec une perte ou un gain de population. Les notions de cas importés et autochtones prennent tout leur sens, et de fait, venant alimenter tous les compartiments. Cet effet sera noté \eta ;
  • Décès naturels, à ne pas confondre avec les décès du compartiment D directement liés à la maladie modélisée, et donc de diminution de la population, touchant alors tous les compartiments également. Cet effet sera noté \mu.

Si l’on souhaite enrichir à nouveau le modèle SEIRD vu précédemment, cela donne, pour les états initiaux S_0, E_0, I_0, R_0, D_0 \in \mathbb{R} ^+ :

\forall t, \begin{cases} S_t = S_{t-1} + \partial S = S_{t-1} - (\lambda - \beta \frac{\cdot S_{t-1} \cdot I_{t-1})}{N} - (\mu + \eta) \cdot S_{t-1}) \cdot \Delta \\ E_t = E_{t-1} + \partial E = E_{t-1} + (\beta \cdot \frac{S_{t-1} \cdot I_{t-1}}{N} - (\sigma + \mu + \eta) \cdot E_{t-1}) \cdot \Delta \\ I_t = I_{t-1} + \partial I = I_{t-1} + (\sigma \cdot E_{t-1} - (\gamma + \delta + \mu + \eta) \cdot I_{t-1}) \cdot \Delta \\ R_t = R_{t-1} + \partial R = R_{t-1} + (\gamma \cdot I_{t-1} - (\mu + \eta) \cdot R_{t-1}) \cdot \Delta \\ D_t = D_{t-1} + \partial D = D_{t-1} + (\delta \cdot I_{t-1} - (\mu + \eta) \cdot D_{t-1}) \cdot \Delta \\ \end{cases}

Quelques autres compartiments classiques

Enfin, il est possible de rajouter pléthore de compartiments différents et variés dans l’objectif de coller au mieux au contexte : 

  • Le compartiment des vaccinés, noté V, intégrant à la fois le taux de vaccinés de départ et celui de nouveaux vaccinés ;
  • Le compartiment des mis en quarantaine, noté Q, intégrant le taux de mis en quarantaine, de guéris et de décès en quarantaine ;
  • Le compartiment des asymptomatiques, noté A, intégrant le taux d’infectieux asymptomatiques qui présentent un enjeu important lors du contrôle d’une épidémie du fait qu’ils sont capables d’infecter d’autres porteurs sains et quasi-impossible à repérer ;
  • Le compartiment des immunisés après infection et guérison, noté Im, et son taux associé ;
  • Le second compartiment S qui permet de remettre en jeu les individus infectés, guéris et non immunisés dans le premier compartiment S, et son taux associé.

– Le nombre de reproduction de base

Le R_0 ou « Nombre de Reproduction de Base » est une mesure cruciale en épidémiologie, offrant des indications sur la transmission d’une maladie dans une population (à ne pas confondre avec le volume/effectif d’individus du compartiment R, des guéris, à l’instant initial t = 0). Il s’agit du nombre moyen d’infections secondaires résultant d’un individu infecté introduit dans une population totalement susceptible. Autrement dit, il mesure la capacité d’une maladie à se propager à travers une population. Son calcul est dépendant du type de modèle utilisé. En reprenant les notations précédentes, en voici plusieurs de ses formes :

Il est également possible de déterminer le nombre moyen de nouvelles infections à un instant t et noté R_t (nombre de reproduction effectif instantané). Sa formule, quel que soit la forme du modèle est :

R_t = R_0 \cdot \frac{S_{t}}{N}

La différence essentielle entre R_0 et R_t est en fait liée à son objectif d’utilisation. Par exemple, pour qualifier une maladie spécifique, le R_0 apportera une information sur sa vitesse de propagation globale permettant de caractériser son niveau de dangerosité. Alors que le R_t trouve son intérêt dans le suivi d’une épidémie en tant qu’indicateur synthétisant son stade d’évolution, souvent utilisé en Santé publique afin déterminer la dimension des mesures de contrôle à mettre en œuvre.

En termes de lecture, un R_t > 1 correspond à une épidémie en pleine propagation puisqu’un individu infecté en contamine au moins un, = 1 lorsqu’elle stagne avec un infecté qui contamine un seul sain, et R_t < 1 pour son déclin, un infecté en contamine moins d’un.

A noter que les formules du R_0 et du R_t reposent sur des hypothèses simplificatrices, telles que l’homogénéité de la population et la constance des paramètres. Dans des contextes plus complexes, où ces conditions ne sont pas remplies, des ajustements et des modèles plus élaborés peuvent être nécessaires pour les évaluer.

\bullet Annexe théorique :

Les équations différentielles ordinaires

Les modèles épidémiologiques à compartiments se basent sur la notion de système d’Equations Différentielles Ordinaires (EDO). Pour obtenir une prévision à un instant t, on résout ce système donnant les différents volumes associés aux différents compartiments. Une EDO à l’ordre n est de la forme générale : 

F(t, y, \frac{\partial y}{\partial t}, \frac{\partial ^2 y}{\partial t ^2}, \cdots, \frac{\partial ^n y}{\partial t ^n})

, avec y(t) fonction adaptée à la forme des compartiments que l’on veut considérer dans le modèle et de forme : 

Y(t) = \begin{bmatrix} S(t) \\ E(t) \\ I(t) \\ R(t) \\ D(t) \\ etc. \end{bmatrix}

Les EDO disposent de plusieurs propriétés essentielles et offrant aux modèles épidémiologiques compartimentaux des caractéristiques des plus pratiques :

– Le théorème d’existence et d’unicité de Picard-Lindelöf et de Cauchy-Lipschitz permet aux EDO, sous certaines conditions, d’admettre une solution qui existe et est unique, relativement aux conditions initiales données ;

Théorème : Soit \frac{\partial y}{\partial t} = f(t,y) une équation différentielle ordinaire du premier ordre avec la condition initiale y(t_0) = y_0, où f(t,y) est continue dans une région rectangulaire R : | t - t_0 | \leq a et | y - y_0 | \leq b. Si f(t,y) est Lipschitz continue par rapport à y dans R, c’est-à-dire qu’il existe une constante L telle que | f(t,y_1) - f(t,y_2) | \leq L | y_1 - y_2 | pour tous (t,y_1) et (t,y_2) dans R, alors l’EDO a une unique solution dans un intervalle | t - t_0 | \leq h, où h est le plus petit des nombres a et \frac{b}{L}.

– De ce théorème découle également le concept d’EDO autonome, assurant son indépendance en fonction du temps. Cela signifie que les taux de variation des variables dépendent uniquement des valeurs actuelles des variables, et non du temps ;

– Le critère de stabilité de Nyquist et le critère de Hurwitz assure la stabilité linéaire pour les systèmes dynamiques décrits par des EDO. Pour les solutions non linéaires, c’est le théorème de stabilité de Lyapunov qui joue un rôle fondamentale. La notion de stabilité indique comment les solutions réagissent aux petites perturbations des conditions initiales ;

Définition : Soit un système dynamique décrit par l’équation différentielle \frac{\partial y}{\partial t} = f(t,y) alors le critère de stabilité de Nyquist s’applique à un système linéaire sous la forme de fonction de transfert déterminé par,

  1. G(s) avec s variable de Laplace.
  2. Elle admet des valeurs spécifiques de s sur le plan complexe le long de l’axe imaginaire s = j \cdot w, j unité imaginaire et w fréquence angulaire.
  3. Après tracé de la fonction G(s), \forall s, on peut conclure en la stabilité du système si et seulement si la courbe ne contourne pas le point critique (-1,0).

Définition : Soit un système dynamique décrit par l’équation différentielle \frac{\partial y}{\partial t} = f(t,y) alors le critère de stabilité de Hurwitz s’applique à un système polynomial déterminé par,

  1. L’équation : a_n \cdot s ^n + a_{n-1} \cdot s ^{n-1} + \cdots + a_1 \cdot s + a_0 = 0.
  2. On détermine la matrice de Hurwitz associée aux coefficients de cette équation.
  3. On peut conclure en la stabilité du système si et seulement tous les mineurs principaux, ou déterminants associés aux sous-matrices, sont strictement positifs.

Définition : Soit D un domaine contenant l’équilibre x ^* d’un système dynamique décrit par l’équation différentielle : y = f(x), où x est un vecteur d’état dans \mathbb{R} ^n, f : D \rightarrow \mathbb{R} ^n est une fonction continue, et x ^* est un point d’équilibre tel que f(x ^*) = 0.

Supposons qu’il existe une fonction continue V : D \rightarrow \mathbb{R} définie positive, c’est-à-dire que V(x) > 0, \forall x \neq x ^*, avec V(x ^*) = 0. De plus, considérons une fonction définie sur D telle que sa dérivée le long des trajectoires du système soit semi-définie négative, c’est-à-dire \dot{V}(x) \leq 0, \forall x \in D.

Alors, l’équilibre x ^* est stable au sein de Lyapunov. Si, de plus, \dot{x} < 0, \forall x \neq x ^*, alors l’équilibre est asymptotiquement stable au sens de Lyapunov. En d’autres termes, la fonction V est une fonction Lyapunov candidate, et la dérivée \dot{V} joue le rôle d’une mesure de la stabilité du système. Si \dot{V} est toujours négative (ou nulle), cela indique que le système est stable ou asymptotiquement stable, respectivement, autour de l’équilibre x ^*.

Ces caractéristiques des EDO font qu’elles sont couramment utilisées pour modéliser les systèmes dynamiques. Ces systèmes peuvent représenter divers phénomènes, tels que la croissance de la population, les mouvements planétaires, les réactions chimiques, etc. On peut les résoudre selon deux principales méthodes : l’algorithme de Newton et l’algorithme de Runge-Kutta.

L’Algorithme d’Euler

L’algorithme d’Euler est une méthode itérative numérique, proposée par Leonhard Euler au XVIIIème siècle, utilisée pour résoudre des équations différentielles ordinaires. Dans le contexte des modèles épidémiologiques, cette méthode peut être appliquée pour résoudre des systèmes d’équations différentielles non linéaires, généralement associés à des modèles plus complexes. C’est une méthode itérative qui se base sur le principe que, près d’une solution, une fonction peut être approximée par une droite tangente à la courbe. Les itérations successives conduisent à une convergence rapide vers la solution.

Étapes de l’algorithme,

  1. Initialisation : Définir une équation différentielle \frac{\partial y}{\partial t} = f(t,y), un pas de temps h et une condition initiale y(t_0) = y_0 ;
  2. Itérations :
    • y_{n+1} = y_n + h \cdot f(t_n, y_n) ;
    • Avec t_n = t_0 + n \cdot h, y_n une estimation de y(t_n), et f(t_n, y_n) est la dérivée de y par rapport à t évaluée en (t_n, y_n).
  3. Répéter l’étape deux jusqu’à atteindre la fin de l’intervalle temporel spécifié.

L’Algorithme de Runge-Kutta

L’algorithme de Runge-Kutta, développé par les mathématiciens allemands Carl Runge et Martin Kutta au début du XXe siècle, constitue une avancée significative dans la résolution numérique d’EDO. Leur objectif était d’améliorer la méthode d’Euler, en proposant une approche plus précise et stable pour approximer les solutions d’EDO. Les méthodes numériques Runge-Kutta, notamment celle d’ordre 4 (RK4) sont aujourd’hui les plus largement utilisées pour résoudre une variété de problèmes scientifiques, y compris les modèles mathématiques en épidémiologie, en offrant un bon compromis entre précision et complexité. Cet algorithme utilise des coefficients pour pondérer les pentes de la fonction à différents points le long du pas de temps. L’itération est répétée jusqu’à ce que la solution soit calculée sur tout l’intervalle temporel souhaité.

Étapes de l’algorithme,

  1. Initialisation : Définir une équation différentielle \frac{\partial y}{\partial t} = f(t,y), un pas de temps h et une condition initiale y(t_0) = y_0 ;
  2. Itération :
    • Calculer les coefficients intermédiaires k_1, k_2, k_3, et k_4 comme suit,
      • k_1 = h \cdot f(t_n, y_n) ;
      • k_2 = h \cdot f(t_n + \frac{h}{2}, y_n + \frac{k_1}{2}) ;
      • k_3 = h \cdot f(t_n + \frac{h}{2}, y_n + \frac{k_2}{2}) ;
      • k_4 = h \cdot f(t_n + h, y_n + \frac{k_3}{2}).
    • Mettre à jour la solution y_{n+1} = y_n + \frac{1}{6} \cdot (k_1 + 2 \dot k_2 + 2 \dot k_3 + k_4) ;
    • Mettre à jour le temps : t_{n+1} = t_n + h ;
  3. Répéter l’étape deux jusqu’à atteindre la fin de l’intervalle temporel spécifié.

\bullet Exemple :

Le modèle SEIRD

On cherche à modéliser la propagation d’une maladie X sur un territoire Y selon un modèle SEIRD qui correspond parfaitement aux attentes puisque l’on dispose des paramètres associés. Ainsi, on définit le taux de contamination \beta = 5 \% = 0.05, avec une durée d’incubation de 5 « temps » le taux d’incubation \sigma = \frac{1}{5} = 0.2, le taux de guérison \gamma = 1 \% = 0.01 et le taux de décès associé à X, \delta = 0.1 \% = 0.001. On veut que ce modèle donne des prévisions sur 1 000 « temps » avec un pas \Delta = 1.

Enfin, on initialise la modélisation avec S_0 = 9999 individus sains, E_0 = 1 infecté introduit dans la population concernée et dont les symptômes ne sont pas encore apparus (« exposé »), I_0 = 0 infecté avec symptômes, R_0 = 0 guéri et D_0 = 0 décédé par la pathologie X. Soit un effectif total N = 9999 + 1 + 0 + 0 + 0 = 10000

Ainsi, pour t  = 1,

\begin{cases} S_1 = S_0 - 0.05 \times \frac{S_0 \times I_0}{10000} \times 1 \\ E_1 = E_0 + (0.05 \times \frac{S_0 \times I_0}{10000} - 0.2 \times E_0) \times 1 \\ I_1 = I_0 + (0.2 \times E_0 - (0.01 + 0.001) \times I_0) \times 1 \\ R_1 = R_0 + 0.01 \times I_0 \times 1 \\ D_1 = D_0 + 0.001 \times I_0 \times 1 \end{cases} = \begin{cases} S_1 = 9999 - 0.05 \times \frac{9999 \times 0}{10000} \times 1 \\ E_1 = 1 + (0.05 \times \frac{9999 \times 0}{10000} - 0.2 \times 1) \times 1 \\ I_1 = 0 + (0.2 \times 1 - (0.01 + 0.001) \times 0) \times 1 \\ R_1 = 0 + 0.01 \times I_0 \times 1 \\ D_1 = 0 + 0.001 \times I_0 \times 1 \end{cases}

= \begin{cases} S_1 = 9999 - 0.05 \times 0 \times 1 \\ E_1 = 1 + (0.05 \times 0 - 0.2 \times 1) \times 1 \\ I_1 = 0 + (0.2 \times 1  - 0.011 \times 0) \times 1 \\ R_1 = 0 + 0.01 \times 0 \times 1 \\ D_1 = 0 + 0.001 \times 0 \times 1 \end{cases}

= \begin{cases} S_1 = 9999 - 0 \\ E_1 = 1 - 0.2 \\ I_1 = 0 + 0.2 \\ R_1 = 0 + 0 \\ D_1 = 0 + 0 \end{cases}

= \begin{cases} S_1 = 9999 \\ E_1 = 0.8 \\ I_1 = 0.2 \\ R_1 = 0 \\ D_1 = 0 \end{cases}

, avec à ce temps, un nombre de nouveaux cas de :

\beta \times \frac{S_1 \times I_1}{N} = 0.05 \times \frac{9999 \times 0.8}{10000} = 0.039996

Pour t = 2,

\begin{cases} S_2 = S_1 - 0.05 \times \frac{S_1 \times I_1}{10000} \times 1 \\ E_2 = E_1 + (0.05 \times \frac{S_1 \times I_1}{10000} - 0.2 \times E_1) \times 1 \\ I_2 = I_1 + (0.1 \times E_1 - (0.01 + 0.001) \times I_1) \times 1 \\ R_2 = R_1 + 0.01 \times I_1 \times 1 \\ D_2 = D_1 + 0.0001 \times I_1 \times 1 \end{cases} = \begin{cases} S_2 = 9999 - 0.0099 \\ E_2 = 0.8 - 0.150001 \\ I_2 = 0.2 + 0.1578 \\ R_2 = 0 + 0.002 \\ D_2 = 0 + 0.0002 \end{cases}

= \begin{cases} S_2 = 9998.99 \\ E_2 = 0.649999 \\ I_2 = 0.3578 \\ R_2 = 0.002 \\ D_2 = 0.0002 \end{cases}

, avec à ce temps, un nombre de nouveaux cas de :

0.05 \times \frac{9998.99 \times 0.3578}{10000} = 0.01788819

Pour t = 3,

\begin{cases} S_3 = 9998.99 - 0.01788819 \\ E_3 = 0.649999 - 0.1121116 \\ I_3 = 0.3578 + 0.126064 \\ R_3 = 0.002 + 0.003578 \\ D_3 = 0.0002 + 0.0003578 \end{cases} = \begin{cases} S_3 = 9998.972 \\ E_3 = 0.5378874 \\ I_3 = 0.483864 \\ R_3 = 0.005578 \\ D_3 = 0.0005578 \end{cases}

, avec à ce temps, un nombre de nouveaux cas de :

0.05 \times \frac{9998.972 \times 0.483864}{10000} = 0.02419071

On continue les calculs jusqu’à avoir parcouru tous les temps t \in \lbrace 1, 2, 3, 4, 5, \cdots, 1000 \rbrace. Et si l’on affiche le graphe que l’on obtient avec les différents vecteurs associés à chacun des compartiments S, E, I, R, D alors on obtient la figure suivante :

A la fin de cette vague, on extrait les indicateurs de Santé publique associés à l’épidémie de la maladie X suite à son introduction au sein de la population que l’on observe :

9 883.329 \approx 9 883 individus ont été infectés, soit  98.8 \% de la population totale ;

883.9951 \approx 884 individus en sont décédés, soit 8.8 \% de la population totale ;

– Le pic de l’épidémie a eu lieu à t = 285 avec 77.58169 \approx 78 nouveaux cas à ce temps-là.

Le modèle SEIRD métapopulationnel

On cherche à modéliser la propagation d’une maladie X sur un territoire Y selon un modèle SEIRD qui correspond parfaitement aux attentes puisque l’on dispose des paramètres associés. Ainsi, on définit le taux de contamination \beta = 5 \% = 0.05, avec une durée d’incubation de 5 « temps » le taux d’incubation \sigma = \frac{1}{5} = 0.2, le taux de guérison \gamma = 1 \% = 0.01 et le taux de décès associé à X, \delta = 0.1 \% = 0.001. On veut que ce modèle donne des prévisions sur 1 000 « temps » avec un pas \Delta = 1.

On inclut également la matrice de contacts entre les différentes classes d’âge que l’on a à disposition. Ces classes d’âge sont construites selon un découpage en cinq intervalles de même longueur :

C = \begin{pmatrix}  & 0-20 ans & 21-40 ans & 41-60 ans & 61-80 ans & 81 ans+ \\ 0-20 ans & 0.3 & 0.7 & 0.2 & 0.5 & 0.1 \\ 21-40 ans & 0.7 & 0.4 & 0.3 & 0.4 & 0.9 \\ 41-60 ans & 0.2 & 0.3 & 0.1 & 0.8 & 0.1 \\ 61-80 ans & 0.5 & 0.4 & 0.8 & 0.9 & 0.7 \\ 81 ans+ & 0.1 & 0.9 & 0.1 & 0.7 & 0.2 \\ \end{pmatrix}

Enfin, on initialise la modélisation avec S_0 = (4999, 2500, 1500, 500, 500) individus sains pour les différentes classes d’âge, E_0 = (1, 0, 0, 0, 0) infecté, de 0-20 ans, introduit dans la population concernée et dont les symptômes ne sont pas encore apparus (« exposé »), I_0 = (0, 0, 0, 0, 0) infecté avec symptômes, R_0 = (0, 0, 0, 0, 0) guéri et D_0 = (0, 0, 0, 0, 0) décédé par la pathologie X.

Soit un effectif total,

N = \begin{pmatrix} 4999 \\ 2500 \\ 1500 \\ 500 \\ 500 \\ \end{pmatrix} + \begin{pmatrix} 1 \\ 0 \\ 0 \\ 0 \\ 0 \\ \end{pmatrix} + \begin{pmatrix} 0 \\ 0 \\ 0 \\ 0 \\ 0 \\ \end{pmatrix} + \begin{pmatrix} 0 \\ 0 \\ 0 \\ 0 \\ 0 \\ \end{pmatrix} + \begin{pmatrix} 0 \\ 0 \\ 0 \\ 0 \\ 0 \\ \end{pmatrix} = \begin{pmatrix} 5000 \\ 2500 \\ 1500 \\ 500 \\ 500 \\ \end{pmatrix}

Ainsi, pour t  = 1,

\begin{cases} S_1 = S_0 - 0.05 \times S_0 \times (\mathbf{C} \cdot I_0) \times \frac{1}{N} \times 1 \\ E_1 = E_0 + (0.05 \times S_0 \times (\mathbf{C} \cdot I_0) \times \frac{1}{N} - 0.2 \times E_0) \times 1 \\ I_1 = I_0 + (0.2 \times E_0 - (0.01 + 0.001) \times I_0) \times 1 \\ R_1 = R_0 + 0.01 \times I_0 \times 1 \\ D_1 = D_0 + 0.0001 \times I_0 \times 1 \end{cases}

= \begin{cases} S_1 = \begin{pmatrix} 0.49 \\ 0.25 \\ 0.15 \\ 0.05 \\ 0.05 \\ \end{pmatrix} - 0.05 \times \begin{pmatrix} 0.3 & 0.7 & 0.2 & 0.5 & 0.1 \\ 0.7 & 0.4 & 0.3 & 0.4 & 0.9 \\ 0.2 & 0.3 & 0.1 & 0.8 & 0.1 \\ 0.5 & 0.4 & 0.8 & 0.9 & 0.7 \\ 0.1 & 0.9 & 0.1 & 0.7 & 0.2 \\ \end{pmatrix} \cdot \begin{pmatrix} 0 \\ 0 \\ 0 \\ 0 \\ 0 \\ \end{pmatrix} \times \frac{1}{\begin {pmatrix} 4500 \\ 2500 \\ 1500 \\ 500 \\ 500 \\ \end{pmatrix}} \times 1 \\ E_1 = \begin{pmatrix} 0.01 \\ 0 \\ 0 \\ 0 \\ 0 \end{pmatrix} + (0.05 \times 0 \times \begin{pmatrix} 0.3 & 0.7 & 0.2 & 0.5 & 0.1 \\ 0.7 & 0.4 & 0.3 & 0.4 & 0.9 \\ 0.2 & 0.3 & 0.1 & 0.8 & 0.1 \\ 0.5 & 0.4 & 0.8 & 0.9 & 0.7 \\ 0.1 & 0.9 & 0.1 & 0.7 & 0.2 \\ \end{pmatrix} \cdot \begin{pmatrix} 0 \\ 0 \\ 0 \\ 0 \\ 0 \\ \end{pmatrix} \times \frac{1}{\begin{pmatrix} 4500 \\ 2500 \\ 1500 \\ 500 \\ 500 \\ \end{pmatrix}} - 0.2 \times \begin{pmatrix} 0.01 \\ 0 \\ 0 \\ 0 \\ 0 \\ \end{pmatrix}) \times 1 \\ I_1 = \begin{pmatrix} 0 \\ 0 \\ 0 \\ 0 \\ 0 \\ \end{pmatrix} + (0.2 \times \begin{pmatrix} 0.01 \\ 0 \\ 0 \\ 0 \\ 0 \\ \end{pmatrix} - (0.01 + 0.0001) \times \begin{pmatrix} 0 \\ 0 \\ 0 \\ 0 \\ 0 \\ \end{pmatrix} \times 1 \\ R_1 = \begin{pmatrix} 0 \\ 0 \\ 0 \\ 0 \\ 0 \\ \end{pmatrix} + 0.01 \times \begin{pmatrix} 0 \\ 0 \\ 0 \\ 0 \\ 0 \end{pmatrix} \times 1 \\ D_1 = \begin{pmatrix} 0 \\ 0 \\ 0 \\ 0 \\ 0 \\ \end{pmatrix} + 0 \times \begin{pmatrix} 0 \\ 0 \\ 0 \\ 0 \\ 0 \\ \end{pmatrix} \times 1 \end{cases}

= \begin{cases} S_1 = \begin {pmatrix} 4999 \\ 2500 \\ 1500 \\ 500 \\ 500 \\ \end{pmatrix} - \begin{pmatrix} 0 \\ 0 \\ 0 \\ 0 \\ 0 \end{pmatrix} \\ E_1 = \begin{pmatrix} 1 \\ 0 \\ 0 \\ 0 \\ 0 \end{pmatrix} + \begin{pmatrix} -0.2 \\ 0 \\ 0 \\ 0 \\ 0 \end{pmatrix} \\ I_1 = \begin{pmatrix} 0 \\ 0 \\ 0 \\ 0 \\ 0 \end{pmatrix} + \begin{pmatrix} 0.2 \\ 0 \\ 0 \\ 0 \\ 0 \end{pmatrix} \\ R_1 = \begin{pmatrix} 0 \\ 0 \\ 0 \\ 0 \\ 0 \end{pmatrix} + \begin{pmatrix} 0 \\ 0 \\ 0 \\ 0 \\ 0 \end{pmatrix} \\ D_1 = \begin{pmatrix} 0 \\ 0 \\ 0 \\ 0 \\ 0 \end{pmatrix} + \begin{pmatrix} 0 \\ 0 \\ 0 \\ 0 \\ 0 \end{pmatrix} \end{cases} = \begin{cases} S_1 = \begin{pmatrix} 4999 \\ 2500 \\ 1500 \\ 500 \\ 500 \\ \end{pmatrix} \\ E_1 = \begin{pmatrix} 0.8 \\ 0 \\ 0 \\ 0 \\ 0 \end{pmatrix} \\ I_1 = \begin{pmatrix} 0.2 \\ 0 \\ 0 \\ 0 \\ 0 \end{pmatrix} \\ R_1 = \begin{pmatrix} 0 \\ 0 \\ 0 \\ 0 \\ 0 \end{pmatrix} \\ D_1 = \begin{pmatrix} 0 \\ 0 \\ 0 \\ 0 \\ 0 \end{pmatrix} \end{cases}

, avec à ce temps, un nombre de nouveaux cas par classe d’âge de :

\beta \times S_1 \cdot I_1 \times \frac{1}{N} = 0.05 \times \begin{pmatrix} 0.3 & 0.7 & 0.2 & 0.5 & 0.1 \\ 0.7 & 0.4 & 0.3 & 0.4 & 0.9 \\ 0.2 & 0.3 & 0.1 & 0.8 & 0.1 \\ 0.5 & 0.4 & 0.8 & 0.9 & 0.7 \\ 0.1 & 0.9 & 0.1 & 0.7 & 0.2 \\ \end{pmatrix} \cdot \begin{pmatrix} 0.2 \\ 0 \\ 0 \\ 0 \\ 0 \end{pmatrix} \times \frac{1}{\begin{pmatrix} 4500 \\ 2500 \\ 1500 \\ 500 \\ 500 \\ \end{pmatrix}}

= \begin{pmatrix} 0.0029994 \\ 0.007 \\ 0.002 \\ 0.005 \\ 0.001 \\ \end{pmatrix}

Pour t = 2,

\begin{cases} S_2 = S_1 - 0.05 \times S_1 \times \mathbf{C} \cdot I_1 \times \frac{1}{N} \times 1 \\ E_2 = E_1 + (0.05 \times S_1 \times \mathbf{C} \cdot I_1 \times \frac{1}{N} - 0.2 \times E_1) \times 1 \\ I_2 = I_1 + (0.2 \times E_1 - (0.01 + 0.001) \times I_1) \times 1 \\ R_2 = R_1 + 0.01 \times I_1 \times 1 \\ D_2 = D_1 + 0.0001 \times I_1 \times 1 \end{cases}

= \begin{cases} S_2 = \begin{pmatrix} 4999 \\ 2500 \\ 1500 \\ 500 \\ 500 \\ \end{pmatrix} - \begin{pmatrix} 0.0029994 \\ 0.007 \\ 0.002 \\ 0.005 \\ 0.001) \\ \end{pmatrix} \\ E_2 = \begin{pmatrix} 0.8 \\ 0 \\ 0 \\ 0 \\ 0 \\ \end{pmatrix} + \begin{pmatrix} -0.1570006 \\ 0.007 \\ 0.002 \\ 0.005 \\ 0.001 \\ \end{pmatrix} \\ I_2 = \begin{pmatrix} 0.2 \\ 0 \\ 0 \\ 0 \\ 0 \\ \end{pmatrix} + \begin{pmatrix} 0.1578 \\ 0 \\ 0 \\ 0 \\ 0 \\ \end{pmatrix} \\ R_2 = \begin{pmatrix} 0 \\ 0 \\ 0 \\ 0 \\ 0 \\ \end{pmatrix} + \begin{pmatrix} 0.002 \\ 0 \\ 0 \\ 0 \\ 0 \\ \end{pmatrix} \\ D_2 = \begin{pmatrix} 0 \\ 0 \\ 0 \\ 0 \\ 0 \\ \end{pmatrix} + \begin{pmatrix} 0.0002 \\ 0 \\ 0 \\ 0 \\ 0 \\ \end{pmatrix} \end{cases} = \begin{cases} S_2 = \begin{pmatrix} 4998.997 \\ 2499.993 \\ 1499.998 \\ 499.995 \\ 499.999 \\ \end{pmatrix} \\ E_2 = \begin{pmatrix} 0.6429994 \\ 0.007 \\ 0.002 \\ 0.005 \\ 0.001 \\ \end{pmatrix} \\ I_2 = \begin{pmatrix} 0.3578 \\ 0 \\ 0 \\ 0 \\ 0 \\ \end{pmatrix} \\ R_2 = \begin{pmatrix} 0.002 \\ 0 \\ 0 \\ 0 \\ 0 \\ \end{pmatrix} \\ D_2 = \begin{pmatrix} 0.0002 \\ 0 \\ 0 \\ 0 \\ 0 \\ \end{pmatrix} \end{cases}

, avec à ce temps, un nombre de nouveaux cas par classe d’âge de :

0.05 \times \begin{pmatrix} 0.3 & 0.7 & 0.2 & 0.5 & 0.1 \\ 0.7 & 0.4 & 0.3 & 0.4 & 0.9 \\ 0.2 & 0.3 & 0.1 & 0.8 & 0.1 \\ 0.5 & 0.4 & 0.8 & 0.9 & 0.7 \\ 0.1 & 0.9 & 0.1 & 0.7 & 0.2 \\ \end{pmatrix} \cdot \begin{pmatrix} 0.3578 \\ 0 \\ 0 \\ 0 \\ 0 \\ \end{pmatrix} \times \frac{1}{\begin{pmatrix} 4500 \\ 2500 \\ 1500 \\ 500 \\ 500 \\ \end{pmatrix}} = \begin{pmatrix} 0.005365923 \\ 0.012522965 \\ 0.003577995 \\ 0.008944911 \\ 0.001788996 \end{pmatrix}

Pour t = 3,

\begin{cases} S_3 = \begin{pmatrix} 4998.997 \\ 2499.993 \\ 1499.998 \\ 499.995 \\ 499.999 \end{pmatrix} - \begin{pmatrix} 0.005365923 \\ 0.012522965 \\ 0.003577995 \\ 0.008944911 \\ 0.001788996 \\ \end{pmatrix} \\ E_3 = \begin{pmatrix} 0.6429994 \\ 0.007 \\ 0.002 \\ 0.005 \\ 0.001 \\ \end{pmatrix} + \begin{pmatrix} -0.123233957 \\ 0.011122965 \\ 0.003177995 \\ 0.007944911 \\ 0.001588996 \\ \end{pmatrix} \\ I_3 = \begin{pmatrix} 0.3578 \\ 0 \\ 0 \\ 0 \\ 0 \end{pmatrix} + \begin{pmatrix} 0.1246641 \\ 0.0014 \\ 0.0004 \\ 0.001 \\ 0.0002 \\ \end{pmatrix} \\ R_3 = \begin{pmatrix} 0.002 \\ 0 \\ 0 \\ 0 \\ 0 \\ \end{pmatrix} + \begin{pmatrix} 0.003578 \\ 0 \\ 0 \\ 0 \\ 0 \\ \end{pmatrix} \\ D_3 = \begin{pmatrix} 0.0002 \\ 0 \\ 0 \\ 0 \\ 0 \\ \end{pmatrix} + \begin{pmatrix} 0.0003578 \\ 0 \\ 0 \\ 0 \\ 0 \\ \end{pmatrix} \end{cases} = \begin{cases} S_3 = \begin{pmatrix} 4998.9916 \\ 2499.9805 \\ 1499.994 \\ 499.9861 \\ 499.9972 \\ \end{pmatrix} \\ E_3 = \begin{pmatrix} 0.519765443 \\ 0.018122965 \\ 0.005177995 \\ 0.012944911 \\ 0.002588996 \\ \end{pmatrix} \\ I_3 = \begin{pmatrix} 0.4824641 \\ 0.0014 \\ 0.0004 \\ 0.001 \\ 0.000 \end{pmatrix} \\ R_3 = \begin{pmatrix} 0.005578 \\ 0 \\ 0 \\ 0 \\ 0 \\ \end{pmatrix} \\ D_3 = \begin{pmatrix} 0.0005578 \\ 0 \\ 0 \\ 0 \\ 0 \\ \end{pmatrix} \end{cases}

, avec à ce temps, un nombre de nouveaux cas par classe d’âge de :

0.05 \times \begin{pmatrix} 0.3 & 0.7 & 0.2 & 0.5 & 0.1 \\ 0.7 & 0.4 & 0.3 & 0.4 & 0.9 \\ 0.2 & 0.3 & 0.1 & 0.8 & 0.1 \\ 0.5 & 0.4 & 0.8 & 0.9 & 0.7 \\ 0.1 & 0.9 & 0.1 & 0.7 & 0.2 \\ \end{pmatrix} \cdot \begin{pmatrix} 0.4824641 \\ 0.0014 \\ 0.0004 \\ 0.001 \\ 0.0002 \\ \end{pmatrix} \times \frac{1}{\begin{pmatrix} 4500 \\ 2500 \\ 1500 \\ 500 \\ 500 \\ \end{pmatrix}} = \begin{pmatrix} 0.007314486 \\ 0.016949110 \\ 0.004888623 \\ 0.012157263 \\ 0.002514306 \\ \end{pmatrix}

On continue les calculs jusqu’à avoir parcouru tous les temps t \in \lbrace 1, 2, 3, 4, 5, \cdots, 1000 \rbrace. Et si l’on affiche le graphe que l’on obtient avec les différents vecteurs associés à chacun des compartiments  S, E, I, R, D et sommés sur les effectifs ventilés par classe d’âge à chaque instant t, alors on obtient la figure suivante :

A la fin de cette vague, on extrait les indicateurs de Santé publique associés à l’épidémie de la maladie X suite à son introduction au sein de la population que l’on observe :

9 844,474 \approx 9 844 individus ont été infectés, soit  98.4 \% de la population totale ;

894,8136 \approx 895 individus en sont décédés, soit 8.9 \% de la population totale ;

– Le pic de l’épidémie a eu lieu à t = 137 avec 106,1841 \approx 106 nouveaux cas à ce temps-là.

[latex]\bullet[/latex] Application sous R :

On proposera la fonction suivante permettant de produire des modélisations incluant les compartiments S, E, I , R, D :

Model_epidemiologique = function(S0,E0,I0,R0,D0,beta,sigma,gamma,delta,pas,T) {

# La fonction s’adapte en fonction du paramétrage rentré. S0 correspond au volume d’individus susceptibles/sains de départ associé à beta le taux de contamination (si <= 0 alors le fonction ne se lance pas), E0 à celui d’exposés/contaminés n’ayant pas encore développé les symptômes associé à sigma le taux d’incubation (si <= 0 alors le compartiment est exclu), I0 à celui d’infectés, R0 à celui des guéris associé à gamma le temps de guérison (si <= 0 alors le compartiment est exclu) et D0 à celui des décédés suite à la pathologie étudiée associé à delta le taux de mortalité lié à la pathologie (si <= 0 alors le compartiment est exclu). Les paramètres pas et T représentant respectivement l’unité de temps et le nombre de temps total parcouru.  

Model_epidemiologique = function(S0,E0,I0,R0,D0,beta,sigma,gamma,delta,pas,T) {

# Si l’on ne veut pas utiliser le compartiment E, et comme l’erreur est facile, on applique un réajustement ici si on s’est trompé
if (sigma == 0 && I0 == 0 && E0 > 0) {
print(« Attention, le paramètre sigma étant nul, le compartiment E ne sera pas utilisé, or sur votre paramétrage Sigma = 0, E > 0 et I = 0, aussi le programme va intervertir les valeurs de E et I »)
I0 = E0
E0 = 0
}

# Chargement de la bibliothèque pour la résolution d’équations différentielles
library(deSolve)

# Fonction générique du modèle SEIRD avec tous les compartiments voulus pour cette fonction
seird_model <- function(t, y, parameters) {
with(as.list(c(y, parameters)), {
dS = (-beta * S * I / N) * pas
if (sigma > 0) {
dE = (beta * S * I / N – sigma * E) * pas
dI = (sigma * E – (gamma + delta) * I) * pas
}
if (sigma == 0) {
dE = 0
dI = (beta * S * I/N – gamma * I) * pas
}
dR = (gamma * I) * pas
dD = (delta * I) * pas
return(list(c(dS, dE, dI, dR, dD)))
})
}

# Conditions initiales
initial_conditions = c(S = S0, E = E0, I = I0, R = R0, D = D0)

# Création du vecteur des différents temps
times = seq(0, T, by = pas)

# Auto-configuration des paramètres du modèle pour l’utilisation de la fonction ode du package deSolve
NAME = « S »
CHOIX = « S »
initiale_state = 1
LEGEND = « Susceptibles »
COL = « blue »
parameters = c(beta = beta)
if (sigma > 0) {
NAME = paste(NAME, « E »)
CHOIX = c(CHOIX, »E »)
LEGEND = c(LEGEND, « Exposés »)
COL = c(COL, « yellow »)
parameters = c(parameters, sigma = sigma)
}
if (beta > 0) {
NAME = paste(NAME, « I »)
CHOIX = c(CHOIX, »I »)
LEGEND = c(LEGEND, « Infectés »)
COL = c(COL, « red »)
}
if (gamma > 0) {
NAME = paste(NAME, « R »)
CHOIX = c(CHOIX, »R »)
LEGEND = c(LEGEND, « Guéris »)
COL = c(COL, « green »)
parameters = c(parameters, gamma = gamma)
}
if (delta > 0) {
NAME = paste(NAME, « D »)
CHOIX = c(CHOIX, »D »)
LEGEND = c(LEGEND, « Décès »)
COL = c(COL, « purple »)
parameters = c(parameters, delta = delta)
}

# Calcul de la population totale pour application des équations différentielles
N = S0 + E0 + I0 + R0 + D0
parameters = c(parameters, pas = pas, N = N)

# Si le paramètre beta est différent de 0 alors on peut lancer la modélisation, si ce dernier est égal à 0, au regard de la formule des équations différentielles, assez logiquement on ne peut lancer les calculs
if (beta > 0) {
# Résolution des équations différentielles avec la fonction ode
result = ode(y = initial_conditions, times = times, func = seird_model, parms = parameters)
# Sélection des colonnes d’intérêt en fonction du type de modèle décrit par les paramètres entrés
result = result[,c(« time »,CHOIX)]

# Tracé des résultats
plot(times, result[, « S »], type = « l », lty = 1, col = « blue », xlab = « Jours », ylab = « Population », main = NAME, ylim = c(0,N))
for (p in 3:dim(result)[2]) {lines(times, result[, CHOIX[p-1]], lty = 1, col = COL[p – 1])}
legend(« topright », legend = LEGEND, col = COL, lty = 1)
}

if (beta == 0) {print(« A minima, le paramètre Beta doit être > 0 »)}
}

On souhaite lancer la modélisation d’un SEIRD. Pour cela on procède de la manière suivante :

Model_epidemiologique(S0=9999,E0=1,I0=0,R0=0,D0=0,beta=0.05,sigma=0.2,gamma=0.01,delta=0.001,pas=1,T=1000)

Parmi les éléments à insérer les plus importants il faut relever :

–  Le volume de susceptibles de départ : S_0 = 9999 ;

–  Le volume d’exposés de départ : E_0 = 1 ;

–  Le volume d’infectés de départ : I_0 = 0 ;

–  Le volume de guéris de départ : S_0 = 0 ;

–  Le volume de décédés de départ : D_0 = 0 ;

–  Le taux de contamination : \beta = 0.05 ;

–  Le taux d’incubation : \sigma = 0.2 ;

–  Le taux de guérison : \gamma = 0.01 ;

–  Le taux de décès : \delta = 0.001.

On obtient alors les résultats suivants :

, qui sont alors les mêmes que ceux obtenus lors des calculs manuels (cf partie « Exemple »).

\bullet Application sous SAS :

On proposera la macro suivante permettant de produire des modélisations incluant les compartiments S, E, I , R, D :

%macro Model_epidemiologique(S0=, E0=, I0=, R0=, D0=, beta=, sigma=, gamma=, delta=, pas=, T=);
/* La fonction s’adapte en fonction du paramétrage rentré. S0 correspond au volume d’individus susceptibles/sains de départ associé à beta le taux de contamination (si <= 0 alors le fonction ne se lance pas), E0 à celui d’exposés/contaminés n’ayant pas encore développé les symptômes associé à sigma le taux d’incubation (si <= 0 alors le compartiment est exclu), I0 à celui d’infectés, R0 à celui des guéris associé à gamma le temps de guérison (si <= 0 alors le compartiment est exclu) et D0 à celui des décédés suite à la pathologie étudiée associé à delta le taux de mortalité lié à la pathologie (si <= 0 alors le compartiment est exclu). Les paramètres pas et T représentant respectivement l’unité de temps et le nombre de temps total parcouru. */
options nonotes;
%if &beta. > 0 %then %do;
/* Si l’on ne veut pas utiliser le compartiment E, et comme l’erreur est facile, on applique un réajustement ici si on s’est trompé */
%if (&sigma. = 0 and &I0. = 0 and &E0. > 0) %then %do;
%put Attention, le paramètre sigma étant nul, le compartiment E ne sera pas utilisé, or sur votre paramétrage Sigma = 0, E > 0 et I = 0, aussi le programme va intervertir les valeurs de E et I;
%let I0 = &E0.;
%let E0 = 0;
%end;
/* Initialisation des compartiements */
%let N = %sysevalf(&S0. + &E0. + &I0. + &R0. + &D0.);
%let S = &S0.;
%let E = &E0.;
%let I = &I0.;
%let R = &R0.;
%let D = &D0.;
data Modelisation;
t = 0;
S = &S0.; E = &E0.; I = &I0.; R = &R0.; D = &D0.;
dS = 0; dE = 0; dI = 0; dR = 0; dD = 0;
N = &N.;
run;
/* Algorithme de calcul des différents compartiments à chaque instant t */
%do t = 1 %to &T.;
%put &t.;
/* Compartiment Susceptibles */
%let dS = %sysevalf(-(&beta. * &S. * &I. / &N.) * &pas.);
/* Compartiment Exposés & Infectés */
%if &sigma. > 0 %then %do;
%let dE = %sysevalf((&beta. * &S. * &I. / &N. – &sigma. * &E.) * &pas.);
%let dI = %sysevalf((&sigma. * &E. – (&gamma. + &delta.) * &I.) * &pas.);
%end;
/* Version du compartiment Infectés si celui Exposés n’est pas enclenché */
%if &sigma. = 0 %then %do;
%let dE = 0;
%let dI = %sysevalf((&beta. * &S. * &I./&N. – &gamma. * &I.) * &pas.);
%end;
/* Compartiment Guéris */
%let dR = %sysevalf((&gamma. * &I.) * &pas.);
/* Compartiment Décès */
%let dD = %sysevalf((&delta. * &I.) * &pas.);
/* Mise à jour des compartiements */
%let S = %sysevalf(&S. + &dS.);
%let E = %sysevalf(&E. + &dE.);
%let I = %sysevalf(&I. + &dI.);
%let R = %sysevalf(&R. + &dR.);
%let D = %sysevalf(&D. + &dD.);
data ajout;
t = &t.;
S = &S.; E = &E.; I = &I.; R = &R.; D = &D.;
dS = &dS.; dE = &dE.; dI = &dI.; dR = &dR.; dD = &dD.;
N = S + E + I + R + D;
run;
data Modelisation;
set Modelisation ajout;
run;
%end;
/* Tracé des résultats */
proc sgplot data=Modelisation;
series x = t y = S / legendlabel = « Susceptibles » lineattrs = (color = blue);
%if &sigma. > 0 %then %do;
series x = t y = E / legendlabel = « Exposés » lineattrs = (color = yellow);
%end;
series x = t y = I / legendlabel = « Infectés » lineattrs = (color = red);
%if &gamma. > 0 %then %do;
series x = t y = R / legendlabel = « Rétablis » lineattrs = (color = green);
%end;
%if &delta. > 0 %then %do;
series x = t y = D / legendlabel = « Décès » lineattrs = (color = purple);
%end;
xaxis label = « Temps »;
yaxis label = « Population »;
keylegend / location = inside position = topright;
run;
%end;
%if &beta. = 0 %then %do;
%put Beta doit être supérieur à 0;
%end;
options notes;
%mend;

On souhaite lancer la modélisation d’un SEIRD. Pour cela on procède de la manière suivante :

%Model_epidemiologique(S0=9999, E0=1, I0=0, R0=0, D0=0, beta=0.05, sigma=0.2, gamma=0.01, delta=0.001, pas=1, T=100);

Parmi les éléments à insérer les plus importants il faut relever :

–  Le volume de susceptibles de départ : S_0 = 9999 ;

–  Le volume d’exposés de départ : E_0 = 1 ;

–  Le volume d’infectés de départ : I_0 = 0 ;

–  Le volume de guéris de départ : S_0 = 0 ;

–  Le volume de décédés de départ : D_0 = 0 ;

–  Le taux de contamination : \beta = 0.05 ;

–  Le taux d’incubation : \sigma = 0.2 ;

–  Le taux de guérison : \gamma = 0.01 ;

–  Le taux de décès : \delta = 0.001.

On obtient alors les résultats suivants :

, qui sont alors les mêmes que ceux obtenus lors des calculs manuels (cf partie « Exemple »).

\bullet Bibliographie :

– A contribution to the mathematical theory of epidemics de W. O. Kermack et A. G. McKendrick ; 

– Infectious diseases of humans : Dynamics and control de R. M. Anderson et R. M. May ;

–  On the definition and the computation of the basic reproduction ratio R0 in models for infectious diseases in heterogeneous populations de O. Diekmann, J. A. Heesterbeek et J. A. Metz ;

– The mathematics of infectious diseases de H. W. Hethcote ; 

–  Modeling infectious diseases in humans and animals de M. J. .Keeling et P. Rohani ;

– Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission de P. Van den Driessche et J.  Watmough ; 

– Appropriate models for the management of infectious diseasesde H. J. Wearing, P. Rohani et M. J. Keeling ; 

– Le site web : https://images.math.cnrs.fr/Modelisation-d-une-epidemie-partie-1.html ;

– Simulation des Equations Différentielles Stochastiques sous R de A. C. Guidoum.

 

Laisser un commentaire