Petit cours d'épidémiologie mathématique
Modèles dans des métapopulations

Julien Arino

Department of Mathematics & Data Science Nexus
University of Manitoba*

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

* The University of Manitoba campuses are located on original lands of Anishinaabeg, Cree, Oji-Cree, Dakota and Dene peoples, and on the homeland of the Métis Nation.

Plan de ce cours

  • Formulation des modèles en métapopulation
  • Analyse mathématique élémentaire
  • n'est pas la panacée - Un centre urbain et des villes satellites
  • Problèmes spécifiques aux métapopulations
  • Investigation computationnelle des grands systèmes

Formulation des modèles en métapopulation

Propagation des maladies dans un monde juridictionnel

Principes généraux (1)

  • lieux géographiques (patchs) dans un ensemble (ville, région, pays..)
  • Les patchs sont les noeuds d'un graphe
  • Each patch contient des compartiments
    • individus susceptibles à la maladie
      • individus infectés par la maladie
      • différentes espèces affectées par la maladie
      • etc.

Principes généraux (2)

  • Les compartiments peuvent se déplacer entre patchs, vaec taux de mouvement des individus du compartiment depuis le patch vers le patch
  • Mouvement instantané et pas de mortalité pendant le mouvement
  • , définit un digraphe avec arcs
  • Arc de à si , absent sinon
  • compartiments, donc tout peut avoir au plus flèches multi-digraphe

Le modèle de mobilité sous-jacent

population du compartiment dans le patch

Supposons pas de naissance ni de mort. Balançons les flux entrants et sortants

si l'on écrit

Le modèle SLIRS dans les patchs

est le taux de naissance (typiquement ou )

= infectés de façon latente ( exposés, bien que ce terme soit ambigu)

Modèle -SLIRS

avec incidence

-SLIRS (espèces multiples)

un ensemble d'espèces

avec incidence

-SLIRS (résidents-voyageurs)

avec incidence

Modèles généraux en metapopulation

compartiments non infectés et compartiments infectés, et

Pour , et ,

où et

Analyse mathématique élémentaire

Analyse - Système illustratif

Considérons le -SLIRS

et l'incidence

Système de équations

La taille n'est pas si gênante ..

Système de équations !!!

Toutefois, beaucoup de structure:

  • copies d'unités individuelles comprenant chacune 4 équations
  • Dynamique des unités constitutives bien comprise
  • Couplage linéaire

Un bon cas de système de grande taille (l'analyse matricielle est très utile ici)

Notations dans ce qui suit

  • une matrice carrée

  • si pour tout (pourrait être la matrice nulle); si et avec ; si . Même notations pour les vecteurs

  • spectre of

  • rayon spectral

  • abscisse spectrale (ou module de stabilité)

  • est une M-matrice si c'est une Z-matrice ( pour ) et que , avec et

Comportement de la population totale

Considérons le comportement de . On a

Donc

Forme vectorielle / matricielle de l'équation

On a

En forme vectorielle

où -matrices avec

La matrice de mouvement

Soit un compartiment . Alors les propriétés suivantes sont vraies:

  1. et correspond au vecteur propre à gauche
  2. M-matrice singulière
  4. irréductible de multiplicité 1

Le cas sympa

On rappelle que

Supposons les taux de mouvement égaux pour tous les compartiments, i.e.,

Alors

Équilibres

si, bien sur, (ou, de façon équivalente, ) est inversible.. L'est-elle ?

une matrice de mouvement et une matrice diagonale. Les résultats suivants sont vrais:

  1. pour tout
  2. et est associé à un vecteur propre . Si, de plus, est irréductible, alors a multiplicité 1 et est associé à
  3. M-matrice non singulière et
  4. irréductible et M-matrice non singulière irréductible et

Non-singularité de

Utilisons une translation du spectre

Ceci donne une contrainte: pour que la population totale se comporte bien (en général, on veut ça), il faut que tous les taux de mortalité soient strictement positifs

Supposons qu'ils le sont (en d'autres termes, supposons que non-singulière). Alors non-singulière et unique

Comportement de la population totale
Cas des mouvements égaux

attire les solutions de

En effet, on a

Puisqu'on suppose à présent que non-singulière, on a (translation du spectre & propriétés de )

irréductible (pourvu que , bien sur)

Comportement de la population totale
Mouvement réductible

Supposons réductible. Soit le nombre d'ensembles minimaux absorbants dans le graphe de connection associé . Alors

  1. L'abscisse spectrale a multiplicité
  2. Associaté à est un vecteur propre positif t.q.
    • si est un nœud dans un ensemble minimal absorbant
    • si est un nœud transient (au sens Markov)

De Foster & Jacquez, Multiple zeros for eigenvalues and the multiplicity of traps of a linear compartmental system, Mathematical Biosciences (1975)

Le cas pas si sympa

Rappelons que

Supposons que les taux de mouvement sont similaires pour tous les compartiments, i.e., la structure de zéros/non-zéros dans toutes les matrices est la même, mais les entrées non-nulles peuvent différer
Soient

et

Cool, non ? Non !

On a alors

Moi, tous les 6 mois: Oooh, coooool, une inclusion différentielle linéaire !

Moi, 10 minutes plus tard: Quel con !

En effet, et ne sont pas des matrices de mouvement (en particulier, leurs sommes de colonnes ne sont pas toutes nulles)

Pas d'ouverture de ce côté là, donc..

Toutefois, non lasciate ogne speranza, on peut quand même faire des choses !

L'équilibre sans maladie (ÉSM)

Supposons le système à l'équilibre et pour . Alors et

On veut résoudre pour . Ici, il est crucial de se souvenir d'un peu d'algèbre linéaire. Écrivons le système en forme vectorielle:

où , -matrices, diagonales

en l'ÉSM

Rappelons la seconde équation:

Donc on a une solution unique si inversible. L'est-elle ?

Cela me rappelle quelque chose !

Par translation du spectre,

Donc, étant donné , est l'unique équilibre et

en l'ÉSM,

en l'ÉSM

ÉSM a et , i.e.,

Rappel: M-matrice singulière. En raisonnant comme avant, a son module d'instabilité translaté vers la droite de . Donc:

  • inversible
  • M-matrice non-singulière

Second point (on aurait si est irréductible)

Donc l'ÉSM a un sens, avec

Calcul de

On utilise la matrice de prochaine génération avec ,

Différentiation p.r. à :

Notons que

dès lors que , donc

Évaluons en l'ÉSM

Si , alors

Si , alors

  • en l'ÉSM
  • en l'ÉSM

Dans les 2 cas, le bloc est zéro, donc

Calculons et évaluons en l'ÉSM

où . Inverse de facile (bloc triangulaire inférieure ):

où

comme

Matrice de prochaine génération

où est le bloc dans . Donc

i.e.,

Stabilité asymptotique locale de l'ÉSM

Définissons pour le -SLIRS par

Alors l'ÉSM

est localement asymptotiquement stable si et instable si

De PvdD & Watmough, Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission, Bulletin of Mathematical Biology 180(1-2): 29-48 (2002)

n'est pas la panacée

Un centre urbain et des villes satellites

JA & S Portet. Epidemiological implications of mobility between a large urban centre and smaller satellite cities. Journal of Mathematical Biology 71(5):1243-1265 (2015)

Contexte de l'étude

Winnipeg centre urbain, 3 villes satellites plus petites: Portage la Prairie, Selkirk et Steinbach

  • densité de population basse à très basse hors de Winnipeg
  • Réseau routier du MB bien étudié par MB Infrastructure Traffic Engineering Branch

Quantités connues et estimées

Ville Pop. (2014) Pop. (2022) Dist. # moyen voyages/jour
Winnipeg (W) 663,617 749,607 - -
Portage la Prairie (1) 12,996 13,270 88 4,115
Selkirk (2) 9,834 10,504 34 7,983
Steinbach (3) 13,524 17,806 66 7,505

Estimation des taux de mouvement

Supposons que est taux de mouvement de la ville à la ville . Ceteris paribus, , donc . Donc, après 1 jour, , i.e.,

Maintenant, , où le nombre d'individus allant de à / jour. Donc

On calcule ça pour toutes les paires et de villes

Sensitivité de aux variations de


with disease: ; without disease: . Chaque boite et moustaches correspondantes sont 10,000 simulations

Une connectivité basse peut contrôler

PLP et Steinbach ont des populations comparables mais avec les paramètres utilisés, seul PLP peut induire des valeurs du général supérieures à 1 quand

Ceci est dû aux taux de mouvement: si , alors

puisque est alors bloc diagonale

Les taux de mouvement vers et depuis PLP sont plus bas situation plus proche du découplage et a un impact plus important sur le général

ne raconte pas toute l'histoire !


Graphiques fonctions de à PLP et de la réduction du mouvement entre Winnipeg et PLP. Gauche: général. Droite: taux d'attaque à Winnipeg

Problèmes spécifiques aux métapopulations

Propriétés dynamiques héritées (a.k.a. je suis flemmard)

Soit

avec des propriétés dynamiques connues, que sait-on de des propriétés de

  • Existence et unicité
  • Invariance de sous le flot
  • Bornitude
  • Localisation des individuels et du général ?
  • SAG ?

Un problème d'héritage - Bifurcations à revers

  • Supposons un modèle qui, lorsqu'il est isolé dans un patch unique, voit une bifurcation à revers
  • Donc le modèle admet des équilibres endémiques sous-critiques
  • Que se passe-t-il lorsque l'on couple plusieurs telles unités ensemble dans une métapopulation ?

OUI, coupler ensemble des unités sujettes à une bifurcation à revers peut conduire à un système qui exhibe une bifurcation à revers

JA, Ducrot & Zongo. A metapopulation model for malaria with transmission-blocking partial immunity in hosts. Journal of Mathematical Biology 64(3):423-448 (2012)

Comportements induits par les métapopulations

Problème "complémentaire" du problème d'héritage. Soit

avec des propriétés connues. Est-ce que

peut exhiber des comportements pas observés dans le unités du système non couplé ?

E.g.: les unités ont un comportement

ÉÉÉ

tandis que la métapopulation a des solutions périodiques

Comportements induits - Équilibres mixtes

Peut-il y avoir des situations dans lesquelles certains patchs sont en l'ÉSM tandis que d'autres sont en un ÉÉ ?

Ceci est le problème des équilibres mixtes

Types d'équilibres (niveau patch)

Un patch à l'équilibre est vide si , à l'équilibre sans maladie si , où sont des indices avec et sont strictement positifs, et à un équilibre endémique si

Types d'équilibres (niveau métapopulation)

Un équilibre de métapopulation sans population corresponds à un équilibre où tous les patchs sont vides. Un équilibre de métapopulation sans maladie a tous les patchs en un ÉSM impliquant les mêmes compartiments pour tous les patchs. Un équilibre de métapopulation endémique a tous les patchs à un équilibre endémique

Équilibres mixtes

Un équilibre mixte est un équilibre de métapopulation t.q.

  • tous les patchs sont à un ÉSM mais le système n'est pas à un ÉSM de métapopulation
  • ou il y a au moins 2 patchs qui sont à un équilibre de type différent (vide, ÉSM ou ÉÉ)

E.g.,

est mixte, de même que

Supposons que le mouvement soit similaire pour tous les compartiments (MSTC) et que le système soit à l'équilibre

  • Si le patch est vide, tous les patchs dans (ancêtres de ) sont également vides
  • Si le patch est en un ÉSM, alors le sous-système consistant de tous les patchs dans est en un ÉSM de métapopulation
  • Si le patch est en un ÉÉ, alors tous les patchs dans (descendants de ) sont aussi en un ÉÉ
  • Si est fortement connexe pour un compartiment , alors il n'existe pas d'équilibres mixtes

Notons que MSTC et pour tout

Problèmes intéressants (AMHA)

Il faut faire plus sur le problème de l'héritage, et en particulier la SAG (Li, Shuai, Kamgang, Sallet, et des choses plus anciennes: Michel & Miller, Šiljak)

Incorporer des temps de voyage (retard) et des événements (infection, guérison, mort) pendant le voyage

Quelle est la complexité minimale des fonctions de mouvement ci-dessous

requise pour observer des comportements induits par la structure de métapopulation ?

Investigation computationnelle
des grands systèmes

Pas très compliqué

  • Comme pour l'analyse mathématique (Cours 12), si l'on procède avec attention et que l'on réfléchit un peu, ça n'est pas trop compliqué
  • Enfin, pas plus dur que dans le cas de la basse dimension
  • Important de prendre avantage de la structure vectorielle

Définissons le champ de vecteurs

SLIAR_metapop_rhs <- function(t, x, p) {
  with(as.list(x), {
    S = x[p$idx_S]
    L = x[p$idx_L]
    I = x[p$idx_I]
    A = x[p$idx_A]
    R = x[p$idx_R]
    N = S + L + I + A + R
    Phi = p$beta * S * (I + p$eta * A)
    dS = S - Phi + p$M %*% S
    dL = Phi - p$epsilon * L + p$M %*% L
    dI = (1 - p$pi) * p$epsilon * L - p$gammaI * I + p$M %*% I
    dA = p$pi * p$epsilon * L - p$gammaA * A + p$M %*% A
    dR = p$gammaI * I + p$gammaA * A + p$M %*% R
    list(c(dS, dL, dI, dA, dR))
  })
}

Mise en place des paramètres

pop = c(34.017, 1348.932, 1224.614, 173.593, 93.261) * 1e+06
countries = c("Canada", "China", "India", "Pakistan", "Philippines")
T = matrix(data = 
             c(0, 1268, 900, 489, 200, 
               1274, 0, 678, 859, 150, 
               985, 703, 0, 148, 58, 
               515, 893, 144, 0, 9, 
               209, 174, 90, 2, 0), 
           nrow = 5, ncol = 5, byrow = TRUE)

Calculons la matrice de mouvement

p = list()
# Use the approximation explained in Arino & Portet (JMB 2015)
p$M = mat.or.vec(nr = dim(T)[1], nc = dim(T)[2])
for (from in 1:5) {
  for (to in 1:5) {
    p$M[to, from] = -log(1 - T[from, to]/pop[from])
  }
  p$M[from, from] = 0
}
p$M = p$M - diag(colSums(p$M))

p$P = dim(p$M)[1]
p$idx_S = 1:p$P
p$idx_L = (p$P + 1):(2 * p$P)
p$idx_I = (2 * p$P + 1):(3 * p$P)
p$idx_A = (3 * p$P + 1):(4 * p$P)
p$idx_R = (4 * p$P + 1):(5 * p$P)
p$eta = rep(0.3, p$P)
p$epsilon = rep((1/1.5), p$P)
p$pi = rep(0.7, p$P)
p$gammaI = rep((1/5), p$P)
p$gammaA = rep((1/3), p$P)
# The desired values for R_0. Here we take something simple
R_0 = rep(1.5, p$P)

Définissons les conditions initiales

# On fixe les conditions initiales
# Par exemple, 2 infectieux au Canada
L0 = mat.or.vec(p$P, 1)
I0 = mat.or.vec(p$P, 1)
A0 = mat.or.vec(p$P, 1)
R0 = mat.or.vec(p$P, 1)
I0[1] = 2
S0 = pop - (L0 + I0 + A0 + R0)
# Vecteur de conditions initiales à passer au solveur
IC = c(S = S0, L = L0, I = I0, A = A0, R = R0)
# Instants auxquels renvoyer une solution 
# Tous les 0.1 jours pour 5 ans ici
tspan = seq(from = 0, to = 5 * 365.25, by = 0.1)

On choisit pour éviter une explosion

On prend pour les patchs quand ils sont isolés. Pour cela, on résoud en fonction de 

for (i in 1:p$P) {
  p$beta[i] = 
    R_0[i] / S0[i] / ((1 - p$pi[i])/p$gammaI[i] + p$pi[i] * p$eta[i]/p$gammaA[i])
}

Et maintenant les problèmes commencent..

# On appelle le solveur
sol <- deSolve::ode(y = IC, times = tspan, 
                    func = SLIAR_metapop_rhs, parms = p)
## DLSODA- At current T (=R1), MXSTEP (=I1) steps
## taken on this call before reaching TOUT
## In above message, I1 = 5000
##
## In above message, R1 = 117.498

La sortie au dessus indique un problème numérique lors de l'intégration. Le problème vient certainement de la différence de taille entre les pays, en particulier le Canada et la Chine

Il faut jouer avec les taux de mouvement et les conditions initiales. Je n'explique pas ici

Épilogue

  • L'espace est une composante fondamentale du processus de propagation épidémique et ne peut pas être ignoré, dans les modèles et dans les décisions de santé publique

  • Une des façons de modéliser l'espace est en utilisant des modèles en métapopulation

  • Les modèles en métapopulation sont faciles à analyser localement et donnent des résultats intéressants au niveau global

  • Les simulations déterministes peuvent être couteuses en RAM et en CPU mais sont faciles

  • Les modèles en métapopulation ne sont pas la seule solution, loin de là !!!