Ergebnis für URL: http://www.ummisco.ird.fr/perso/bacaer/2014JMB2.html Sur les processus linéaires de naissance et de mort dans un environnement aléatoire
J. Math. Biol. 69 (2014) 73-90
Nicolas Bacaër
Institut de recherche pour le développement
Bondy, France
nicolas.bacaer@ird.fr
Abdelkarim Ed-Darraz
Université Cadi Ayyad, Département de mathématiques, Marrakech, Maroc
* [de] Um den Artikel auf Deutsch zu lesen, benutzen Sie den in Ihrem Browser
integrierten automatischen Übersetzer.
* [en] To read the article in English, use your browser's machine translator.
[es] Para leer el artículo en español, utilice el traductor automático de su
navegador.
[it] Per leggere l'articolo in italiano, utilizzare il traduttore automatico
del browser.
[ja]
日本語de記事wo読muniha,_BuRaU6Zani内
;蔵saretei5ru自動翻訳機能wogo࠷
3;用kudasai5._
[pt] Para ler o artigo em português, utilize o tradutor automático do seu
navegador.
* [ru] CHtoby prochitat' stat'yu na russkom yazyke, vospol'zujtes' vstroennym
avtomaticheskim perevodchikom Vashego brauzera.
* [zh]
要阅读中文文本,请使
992;你的浏览器内置的自&
#21160;翻译器._
Résumé
On étudie la probabilité d'extinction pour des processus linéaires de naissance
et de mort à un ou plusieurs types dans un environnement évoluant suivant une
chaîne de Markov. La probabilité d'extinction est presque sûrement égale à 1 si
et seulement si la reproductivité est inférieure ou égale à 1. Le point clé est
de trouver la définition convenable de la reproductivité pour que ce résultat
soit vrai.
1. Introduction
Un article récent de Bacaër et Ait Dads (2012) étudie la probabilité
d'extinction w pour un processus linéaire de naissance et de mort avec plusieurs
types dans un environnement périodique. Dans ce cas w=1 si la reproductivité
\(\,R_0\,\) (Dublin et Lotka, 1925) est inférieure ou égale à 1, tandis que w < 1
si \(R_0 > 1\). L'emploi de la reproductivité est motivé par des applications en
épidémiologie. La preuve utilise la technique standard basée sur l'équation aux
dérivées partielles linéaire du premier ordre satisfaite par une fonction
génératrice (Kendall, 1948). Pour les modèles périodiques de période T à un seul
type, avec un taux de naissance \(a(t)\) et avec une mortalité \(b(t)\,\),
\[R_0=\frac{\int_0^T a(t)\, dt}{\int_0^T b(t)\, dt}.\] La même reproductivité
sert aussi de seuil pour les modèles de population périodiques sans stochasticité
démographique (Bacaër et Guernaoui, 2006, Sect. 5).
L'objectif ici est d'étudier les processus linéaires de naissance et de mort
à plusieurs types dans un environnement aléatoire. Essayons de résumer la
littérature sur ce sujet. Pour des modèles de population en temps discret dans un
environnement aléatoire mais sans stochasticité démographique, l'espérance de la
population peut croître à l'infini même si l'extinction se produit presque
sûrement (Lewontin et Cohen, 1969). Voir aussi (Haccou et coll., 2005). Athreya
et Karlin (1971) ont étudié les processus de branchement en temps discret dans un
environnement aléatoire, à la fois dans le cas d'un seul type et celui de
plusieurs types. Cogburn et Torrez (1981) ont étudié les processus linéaires de
naissance et de mort à un seul type dans un environnement aléatoire. S'il y a un
nombre fini d'environnements,
* \(a_k\,\) est le taux de naissance
* \(b_k\,\) est la mortalité
* \(u_k\,\) est la proportion moyenne du temps passée dans l'environnement k.
Leur corollaire 3.2 montre que w=1 si et seulement si \begin{equation}\tag{1}
\sum_k (a_k-b_k)u_k\leq 0\, . \end{equation} Leur preuve est basée sur des
résultats dus à (Kaplan, 1973).
Britton et Lindholm (2009) ont étudié le même processus, avec un seul type
d'individu, dans le cas particulier de deux environnements. Ils ont suggéré que
w=1 si et seulement si \(\,R_\star\leq 1\,\), avec \[R_\star=m_1 m_2, \quad
m_k=\int_0^\infty q_k\, e^{-q_k\tau} e^{(a_k-b_k)\tau}\, d\tau.\] \(q_k\,\) est
la taux auquel l'environnement quitte l'état k. Ils ont aussi montré que
\(\,R_\star\) avait la même position par rapport à 1 que le rayon spectral de la
« matrice de prochaine génération » \begin{equation}\tag{2} \left
(\begin{array}{cc} a_1 & 0\\ 0 & a_2 \end{array}\right ) \left (\begin{array}{cc}
b_1+q_1 & -q_2\\ -q_1 & b_2+q_2 \end{array}\right )^{-1}. \end{equation} Ils ont
appelé le rayon spectral de cette matrice \(\,R_0\).
(Gray et coll., 2012) ont étudié un modèle de population en temps continu et
environnement aléatoire mais sans stochasticité démographique. Ils ont remarqué
que la position de \begin{equation}\tag{3} \frac{\sum_k a_k \, u_k}{\sum_k b_k\,
u_k} \end{equation} par rapport à 1 sert de seuil pour l'extinction. Ils ont
nommé ce nombre \(\,T_0\). Ils ont gardé la notation \(\,R_0\,\) pour le rayon
spectral de (2). Remarquer que le seuil donné par (3) est le même que (1).
Bacaër et Khaladi (2012) ont montré que (3) était le rayon spectral d'un
opérateur « de prochaine géneration » en dimension infinie et ont suggéré de
garder la notation \(\,R_0\,\). (Hernandez-Suarez et coll., 2012) et (Artalejo et
coll., 2012) ont également abordé des problèmes reliés concernant la
reproductivité. Comme on le verra ci-dessous, à part le problème des notations,
(3) et \(\,R_\star\,\) (ou le rayon spectral de (2)) peuvent ne pas avoir la même
position par rapport à 1.
Dans la section 2, on présente deux preuves alternatives du résultat de
Cogburn et Torrez (1981) pour les modèles à un seul type. La première preuve
utilise une formule de Kendall (1948) pour la probabilité d'extinction dans un
environnement variable. La seconde preuve ramène le problème en temps continu au
cadre en temps discret d'Athreya et Karlin (1971). On discute aussi en détail un
exemple simple avec juste deux environnements, qui espérons-le clarifiera les
problèmes concernant la définition de la reproductivit&é mentionnés ci-dessus. La
section 3 étudie la probabilité d'extinction pour des populations de plusieurs
types dans un cadre avec temps continu. On réduit une nouvelle fois le problème
au cas en temps discret d'Athreya et Karlin (1971). Ceci ne semble possible que
pour un nombre fini d'environnements. Il se peut que l'approche suivant Kendall
(1948) puisse aussi être généralisée au cas de plusieurs types en utilisant les
résultats de (Chueshov, 2002) ou de (Benaïm et Schreiber, 2009), éventuellement
sans la restriction d'un nombre fini d'environnements. On présente aussi des
simulations numériques. La conclusion suggère d'autres possibilités de
généralisation.
2. Le modèle à un seul type
2.1 La probabilité d'extinction
Soit un processus linéaire de naissance et de mort à un seul type dans un
environnement variable. \(a(t)\,\) est le taux de naissance et \(b(t)\) la
mortalité au temps t. On suppose \(a(t)=a_{\theta(t)}\) et
\(b(t)=b_{\theta(t)}\). \(\,\theta(t)\) est un processus stochastique à valeurs
dans \(\{1,2,\ldots,K\}\), qui représentent différents environnements. On suppose
\(a_k > 0\) et \(b_k > 0\,\) pour tout k. On suppose que les basculements entre
les environnements suivent une chaîne de Markov homogène en temps continu.
\(Q_{k,\ell}\geq 0\,\) est le taux auquel l'environnement bascule de \(\ell\,\)
vers k (\(k\neq \ell\)). Q est la matrice correspondante avec
\[Q_{\ell,\ell}=-q_\ell,\quad q_\ell=\sum_{k\neq \ell} Q_{k,\ell}.\] On suppose
que la matrice Q est irréductible. On a donc \(\,q_k > 0\,\) pour tout k. Par
conséquent, il existe une unique distribution de probabilités stationnaire et
strictement positive u avec \[Qu=0,\quad \sum_k u_k=1, \quad u_k > 0\quad \forall
\, k\] (Pardoux, 2008, p. 147). \(R_0\,\) est donné par la formule (3). Soit w la
probabilité d'extinction partant d'un individu au temps 0 dans l'environnement
\(\,k_0\). La proposition suivante, bien qu'un cas spécial des résultats de
Cogburn et Torrez (1981), sera démontrée d'une manière un peu plus simple et de
deux façons différentes.
Proposition 1.
* w=1 presque sûrement si \(\,R_0\leq 1\),
* w < 1 presque sûrement si \(\,R_0 > 1\).
Preuve. On sait d'après (Kendall, 1948, équation (18)) que
\begin{equation}\tag{4} \omega=1-\frac{1}{1+\int_{0}^\infty b(s)\, \exp \bigl
[\int_{0}^s (b(v)-a(v))\, dv\bigr ]\, ds}\, , \end{equation} que l'intégrale au
dénominateur soit finie ou infinie. D'après le théorème ergodique (Pardoux, 2008,
p. 150), \[\lim_{s\to +\infty} \frac{1}{s} \int_{0}^s (b(v)-a(v))\, dv = \sum_k
(b_k-a_k)\, u_k \] presque sûrement. Dans le cas \(\,R_0 < 1\,\), on a \(\,\sum_k
(b_k-a_k)\, u_k > 0\,\) et l'intégrale au dénominateur de (4) diverge. On a donc
w=1. Dans le cas \(\,R_0 > 1\,\), on a \(\,\sum_k (b_k-a_k)\, u_k < 0\,\) et
l'intégrale au dénominateur de (4) converge. On a donc w < 1.
On n'a pas considéré le cas critique \(R_0=1\). Cependant, voici une seconde
preuve de la proposition 1, qui couvre aussi le cas critique.
Preuve. Le processus en temps continu à K états est lié à une chaîne de
Markov en temps discret, dont l'espace d'états est
\[\mathcal{X}=\{1,2,\ldots,K\}\times \mathbb{R}_+.\] Chaque pas de temps
correspond au temps entre deux basculements de l'environnement. Au lieu de dire
par exemple que l'environnement est dans l'état k pendant t unités de temps puis
dans l'état \(k'\) pendant \(t'\) unités de temps (avec \(k'\neq k\)), on dit
qu'il y a la transition suivante \[(k,t)\in \mathcal{X}\longrightarrow (k',t')\in
\mathcal{X}.\] La probabilité pour que l'environnement \(\,k'\,\) suive
l'environnement k est donnée par \[\Pi_{k',k}=\left \{\begin{array}{lll}
\frac{Q_{k',k}}{q_k},& &k'\neq k,\\ 0, & &k'= k. \end{array}\right.\]
L'environnement \(k'\) dure entre \(t'\) et \(t'+dt'\) unités de temps (\(dt'\)
infiniment petit) avec une probabilité \(q_{k'}\, e^{-q_{k'}t'}\, dt'\). La
probabilité pour que l'environnement \(\,k'\,\) dure entre \(t'\) et \(t'+dt'\,\)
unités de temps et suive l'environnement k qui a duré t unités de temps, est
\[P_{(k,t)\to(k',t')}\, dt' = \Pi_{k',k} \ q_{k'}\, e^{-q_{k'}t'}\, dt'\, .\] On
a \[\int_0^\infty \sum_{k'} P_{(k,t)\to(k',t')} dt'=1\, .\] Plus généralement,
pour \(z\in \{1,2,\ldots\}\,\), la probabilité de transition en z étapes est
donnée par \[P^{(z)}_{(k,t)\to(k',t')}\, dt' = (\Pi^z)_{k',k} \ q_{k'}\,
e^{-q_{k'}t'}\, dt'\, .\] \(\Pi^z\) est la produit itéré de la matrice
\(\Pi=(\Pi_{k',k})\). \(Qu=0\) est équivalent à \[\sum_{k\neq k'} Q_{k',k}
u_k=q_{k'} u_{k'},\quad \forall k'.\] Avec \begin{equation}\tag{5}
w_{k,t}=\frac{q_k u_k}{\sum_\ell q_\ell u_\ell}\, q_k\, e^{-q_k t} \end{equation}
on a donc \[ \sum_k \int_0^\infty w_{k,t}\, dt =1,\quad \quad \sum_k
\int_0^\infty w_{k,t}\, P_{(k,t)\to(k',t')} \, dt = w_{k',t'},\quad \forall
(k',t')\in \mathcal{X}.\] \((w_{k,t})\) est donc une distribution de probabilités
stationnaire et la chaîne de Markov dans \(\,\mathcal{X}\,\) est récurrente
positive (Meyn et Tweedie, 1993). Athreya et Karlin (1971) mentionnent dans une
remarque suivant leur théorème 4 que leurs résultats restent valables non
seulement pour les chaînes de Markov irréductibles récurrentes positives sur un
espace d'états dénombrable mais aussi pour un processus stationnaire ergodique
sur un espace d'états général (non dénombrable) tel que \(\,\mathcal{X}\).
Entre les basculements de l'environnement, on a un processus linéaire de
naissance et de mort dans un environnement constant. Dans un environnement k
durant t unités de temps, soit \(\,\phi_{k,t}(x)\) la fonction génératrice de la
population au bout de l'intervalle de temps, partant d'un individu au début de
cet intervalle de temps: \[\phi_{k,t}(x)=\frac{b_k(1-x) e^{t(a_k-b_k)} + a_k
x-b_k}{a_k(1-x) e^{t(a_k-b_k)}+a_k x-b_k}\] si \(a_k\neq b_k\) tandis que
\[\phi_{k,t}(x)=\frac{x+(1-x)a_k t}{1+(1-x) a_k t}\] si \(a_k=b_k\,\) (Hillion,
1986, p. 118). L'espérance de la population au bout de l'intervalle de temps est
égale à \[\phi_{k,t}'(1)=e^{(a_k-b_k)t}.\] La probabilité que la population soit
éteinte au bout de l'intervalle de temps est
\[\phi_{k,t}(0)=\frac{1-e^{(a_k-b_k)t}}{1-e^{(a_k-b_k)t}\, a_k/b_k}\, \] si
\(a_k\neq b_k\,\), tandis que \(\phi_{k,t}(0)=a_k\, t/(1+a_k\, t)\) si
\(a_k=b_k\). Avec \(\,a_k < b_k\,\), on a \[1-\phi_{k,t}(0)\sim_{t\to +\infty}
(1-a_k/b_k)\, e^{(a_k-b_k)t}.\] Avec \(a_k=b_k\,\), on a
\[1-\phi_{k,t}(0)\sim_{t\to +\infty} 1/(a_k\, t).\] Donc on peut facilement
vérifier que \begin{align*} &\mathbb{E}(|\log (1-\phi(0))|)=\sum_k \int_0^\infty
w_{k,t} | \log (1- \phi_{k,t}(0)) | \, dt < +\infty\, ,\\ &\mathbb{E}([\log
\phi'(1)]^+)=\sum_k \int_0^\infty w_{k,t} [ \log \phi_{k,t}'(1)]^+ \, dt <
+\infty\, , \end{align*} à cause de la décroissance exponentielle de
\(w_{k,t}\,\) par rapport à t. Avec (Athreya et Karlin, 1971), on conclut que w=1
si et seulement si \[ \mathbb{E}(\log \phi'(1)) = \sum_k \int_0^\infty w_{k,t}
\log \phi'_{k,t}(1)\, dt \leq 0.\] Avec \(\int_0^\infty t\, e^{-q_kt}\, dt =
1/(q_k)^2\), on obtient \[\mathbb{E}(\log \phi'(1)) = \sum_k \int_0^\infty
\frac{q_k u_k}{\sum_\ell q_\ell u_\ell}\, q_k\, e^{-q_k t} [(a_k-b_k)t]\,
dt=\frac{\sum_k (a_k-b_k)u_k}{\sum_\ell q_\ell u_\ell}\, .\] Donc w=1 si et
seulement si \(R_0\leq 1\).
2.2 Un exemple et quelques remarques
Comme dans l'exemple de Britton et Lindholm (2009, section 3), on suppose
qu'il y a deux environnements:
* \(a_1= \mbox{2,7}\) et \(b_1=2\)
* \(a_2=\mbox{0,8}\) et \(b_2=2\).
On définit \[Q=\left (\begin{array}{cc} -q_1 & q_2\\ q_1 & -q_2 \end{array}\right
)\] avec \(q_1=q_2=1\). On a alors \(\,u_1=q_2/(q_1+q_2)=\mbox{0,5}\),
\(u_2=q_1/(q_1+q_2)=\mbox{0,5}\) et \(R_0= \mbox{0,875} < 1\). On a donc
\(\,\omega=1\). Des simulations numériques tendent à confirmer cette conclusion
(figure 1a). [2014JMB2Figure1a.png] [2014JMB2Figure1b.png] Figure 1. À gauche :
100 simulations de la population en fonction du temps t dans le cas où \(a_1=
\mbox{2,7}\,\), en partant d'un individu dans l'environnement 1. Toutes les
simulations conduisent à l'extinction. À droite : en partant d'un individu dans
l'environnement 1 mais avec \(a_1=\mbox{5,4}\,\). On a simulé 1000 histoires pour
l'environnement et calculé la probabilité d'extinction w par la formule (4). La
figure montre un histogramme des valeurs de w (\(0\leq \omega\leq 1\)).
Comme autre exemple, considérons les mêmes valeurs des paramètres sauf que
\(a_1\) est doublé: \(a_1=\mbox{5,4}\). Dans ce cas on a \(R_0= \mbox{1,55} >
1\). La figure 1b montre un histogramme de la probabilité d'extinction w en
partant d'une personne dans l'environnement 1. La moyenne de la probabilité
d'extinction est environ 0,61. Cette moyenne serait 0,85 en partant de
l'environnement 2. L'histogramme a été obtenu en approchant la chaîne de Markov
en temps continu, qui gouverne l'environnement \(\,(0 < t < 100)\,\), par une
chaîne de Markov en temps discret avec un pas de temps \(\,\varepsilon=0,00005\).
L'ordinateur choisit aléatoirement 1000 réalisations de cette chaîne de Markov,
formant 1000 histoires environnementales. La formule (4) pour w est ensuite
estimée numériquement. Avec \(\,b_1=b_2\), \(a_2\leq a_1\) et \(a_1\geq b_1\,\),
on s'aperçoit facilement avec (4) que \(\omega \geq b_1/a_1\). Cette borne
inférieure correspond à la probabilité d'extinction si l'environnement est
toujours 1. Avec les valeurs numériques ci-dessus, on obtient \(\omega\geq 2/
\mbox{5,4} \simeq \mbox{0,37}\,\), en accord avec la figure 1b.
Sous-criticalité du modèle discrétisé. Pour dessiner la figure 1, on a
discrétisé le processus en temps continu en utilisant un pas de temps e>0. Pour
simplifier, supposons comme dans l'exemple qu'il n'y ait que K=2 environnements.
La matrice de transition de la chaîne de Markov en temps discret est
\[\mathcal{P}=\left (\begin{array}{cc} 1-q_1\varepsilon & q_2\varepsilon \\
q_1\varepsilon & 1-q_2\varepsilon \end{array}\right ).\] Sa distribution
stationnaire \(\varpi\) est telle que \(\mathcal{P}\varpi=\varpi\) avec
\(\sum\varpi_k=1\). On obtient \(\,\varpi_1=q_2/(q_1+q_2)\) et
\(\varpi_2=q_1/(q_1+q_2)\), qui sont indépendants de e et coïncident avec les
distributions stationnaires \(u_1\) et \(u_2\,\) du processus en temps continu.
Si l'environnement est de type k (\(k=1\) ou 2), on a supposé que pendant un pas
de temps chaque individu a une probabilité de donner naissance \(a_k\varepsilon\)
et une probabilité de mourir \(b_k\varepsilon\). Au pas de temps suivant, chaque
individu conduit
* à aucun individu avec une probabilité \((1-a_k\varepsilon)b_k\varepsilon\)
[pas de naissance, une mort],
* à 1 individu avec une probabilité
\((1-a_k\varepsilon)(1-b_k\varepsilon)+a_k\varepsilon b_k\varepsilon\) [soit
pas de naissance et pas de mort, soit une naissance et une mort],
* à 2 individus avec une probabilité \(a_k\varepsilon (1-b_k\varepsilon)\) [une
naissance, pas de mort].
La moyenne de cette distribution est \(1+a_k\varepsilon-b_k\varepsilon\). Donc
suivant la théorie des processus de branchement en temps discret en environnement
aléatoire (Athreya et Karlin, 1971), le processus est sous-critique et conduit à
l'extinction presque sûrement si et seulement si \[T(\varepsilon)=\sum_k \varpi_k
\log (1+(a_k-b_k)\varepsilon)\leq 0.\] Rappelons que \(\,\varpi_k=u_k\). On a
\[T(\varepsilon)\mathop{\sim}_{\varepsilon\to 0} \varepsilon \sum_k u_k
(a_k-b_k).\] Cette expression a le même signe que \(\,R_0-1\).
Une autre manière d'utiliser les résultats d'Athreya et Karlin. Une autre
manière de voir cet exemple avec seulement deux environnements est de le
considérer comme un processus de branchement dans une suite « d'environnements »
aléatoires indépendants et identiquement distribués. En effet la séquence
d'environnements \(\,1\to 2\,\) se répète à l'identique, le temps passé dans
chaque environnement étant aléatoire. La probabilité pour que l'environnement k
(\(k=1\) ou 2) dure entre \(t_k\) et \(t_k+dt_k\) unités de temps est \(q_k\,
e^{-q_k t_k}\, dt_k\). Ces probabilités sont indépendantes. Le nouvel espace des
« environnements » est donc \(\{(t_1,t_2)\in (\mathbb{R}_+)^2\}\). Si chaque
environnement dure \(\,t_k\) unités de temps, alors la croissance moyenne durant
une séquence est \[\mathcal{M}=e^{(a_1-b_1)t_1} e^{(a_2-b_2)t_2}.\] La théorie
des processus de branchement dans les environnements indépendants et
identiquement distribués (Athreya et Karlin, 1971) montre que \(\,\omega=1\) si
et seulement si \[\mathbb{E}(\log \mathcal{M}) = \int_0^\infty \!\!\int_0^\infty
q_1\, e^{-q_1 t_1} q_2\, e^{-q_2 t_2} \log \left ( e^{(a_1-b_1)t_1}
e^{(a_2-b_2)t_2} \right ) dt_1\, dt_2 \leq 0.\] Mais on voit que \begin{align*}
\mathbb{E}(\log \mathcal{M}) &= \int_0^\infty \!\!\int_0^\infty q_1\, e^{-q_1
t_1} q_2\, e^{-q_2 t_2} \bigl [(a_1-b_1)t_1+(a_2-b_2)t_2 \bigr ] dt_1\, dt_2\\ &=
\int_0^\infty q_1\, e^{-q_1 t_1} (a_1-b_1)t_1\, dt_1 + \int_0^\infty q_2\,
e^{-q_2 t_2} (a_2-b_2)t_2\, dt_2. \end{align*} On a donc \[ \mathbb{E}(\log
\mathcal{M}) =\frac{(a_1-b_1)q_2+(a_2-b_2)q_1}{q_1 q_2}\, .\] Le signe de cette
quantité est le même que celui de \(R_0-1\). Il y a extinction presque sûrement
si et seulement si \(\,R_0\leq 1\).
Une autre manière de calculer la probabilité d'extinction. \(\,p_{k,n}(t)\,\)
est la probabilité que la population soit dans l'environnement k au temps t avec
n individus. Le processus est considéré comme une chaîne de Markov homogène en
temps continu. L'espace d'états est \(\,\{1,2,\ldots,K\}\times \mathbb{N}\).
Comme dans l'équation (2) de (Yechiali, 1973), on obtient \begin{align}
\frac{dp_{k,n}}{dt} =& -(a_k+b_k)n \,p_{k,n} + b_k (n+1) p_{k,n+1} + a_k (n-1)
p_{k,n-1} \nonumber\\ &+ \sum_{\ell\neq k} \left (Q_{k,\ell} p_{\ell,n} -
Q_{\ell,k} p_{k,n} \right )\, . \tag{6} \end{align} Considérons la chaîne de
Markov induite en temps discret, obtenue en ne considérant que les sauts du
processus en temps continu :
* \(a_k n /(a_k n + b_k n +q_k)\) est la probabilité de la transition \((k,n)
\to (k,n+1)\),
* \(b_k n /(a_k n + b_k n +q_k)\) est la probabilité de la transition \((k,n)
\to (k,n-1)\)
* \(Q_{\ell,k} /(a_k n + b_k n +q_k)\) est la probabilité de la transition
\((k,n) \to (\ell,n)\) si \(\ell\neq k\).
On range les couples \((k,n)\) dans l'ordre suivant : \[(1,0), \ldots, (K,0),
(1,1), \ldots, (K,1),\] etc. On définit
* \(\pi_{k,n}(j)\,\) la probabilité que le système soit après j transitions
dans la position \((k,n)\)
* \(\pi(j)=(\pi_{k,n}(j))\) le vecteur formé de ces probabilités, avec les
indices \((k,n)\) rangés comme ci-dessus
* \(\delta_{k,\ell}=1\) si \(k=\ell\) et \(\delta_{k,\ell}=0\) si \(k\neq
\ell\). \(\delta_{k,\ell}\) est le symbole de Kronecker.
Pour tout n (avec n>=0 ou n>=1) et tout \(1\leq k,\ell\leq K\), on définit
\[L_{k,\ell}^{(n)} = \frac{(n-1)a_\ell\,
\delta_{k,\ell}}{(n-1)(a_\ell+b_\ell)+q_\ell},\quad M^{(n)}_{k,\ell} =
\frac{Q_{k,\ell} (1-\delta_{k,\ell})}{n(a_\ell+b_\ell)+q_\ell},\]
\[N_{k,\ell}^{(n)} = \frac{(n+1)b_\ell\,
\delta_{k,\ell}}{(n+1)(a_{\ell}+b_{\ell})+q_\ell}\, .\] On a alors
\[\pi_{k,n}(j+1)=L_{k,k}^{(n)} \pi_{k,n-1}(j) + \sum_{\ell\neq k}
M_{k,\ell}^{(n)} \pi_{\ell,n}(j) + N_{k,k}^{(n)} \pi_{k,n+1}(j)\, .\] On a donc
\[\pi(j+1)=H \pi(j),\] où la matrice H a la structure triangulaire par blocs
\[H=\left (\begin{array}{cccc} M^{(0)} & N^{(0)} & 0 & \cdots\\ L^{(1)} & M^{(1)}
& N^{(1)} & \ddots \\ 0 & L^{(2)} & M^{(2)} & \ddots \\ \vdots & \ddots & \ddots
& \ddots \end{array}\right ),\] comme dans l'équation (1.1) de (Gaver et coll.,
1984). Remarquer au passage que \(L^{(1)}=0\). On a \[\pi(j)=H^j \pi(0) \quad
\forall \, j\geq 0.\] L'ensemble des positions \((k,n)\,\) avec n=0 est un
sous-ensemble absorbant. Donc la probabilité d'extinction \(\,\Omega_{k,n}\,\),
en partant de n personnes dans l'environnement k, est la plus petite solution du
système \[\Omega=\Omega H,\quad \Omega_{k,0}=1\ \forall k,\] où W* est le vecteur
ligne \((\Omega_{k,n})\,\) avec des indices ordonnés comme avant (Bouleau, 1988,
p. 76). Plus explicitement, on a \[\Omega_{k,n}= \frac{\Omega_{k,n-1} n
b_k}{n(a_k+b_k)+q_k} + \sum_{\ell\neq k} \frac{\Omega_{\ell,n}
Q_{\ell,k}}{n(a_k+b_k)+q_k}+\frac{\Omega_{k,n+1} n a_k}{n(a_k+b_k)+q_k}\, .\]
C'est équivalent à l'équation (4.1) de Cogburn et Torrez (1981). \(\Omega\) peut
être calculé numériquement en tronquant les matrices à un ordre suffisamment
grand et en prenant la limite quand \(i\to +\infty\) de \(\Omega^{(i)}\) avec
\(\Omega^{(i+1)}=\Omega^{(i)} H\) et \(\Omega_{k,n}^{(0)}=\delta_{n,0}\) pour
tout k et tout \(n\). Pour les exemples numériques, on a tronqué à n=500 et itéré
jusqu'à i=20000. Avec \(a_1= \mbox{2,7}\) on obtient \(\Omega_{1,1}\simeq
\Omega_{2,1}\simeq \mbox{1,0}\). Avec \(\,a_1= \mbox{5,4}\) on obtient
\(\Omega_{1,1}\simeq \mbox{0,61}\) et \(\Omega_{2,1}\simeq\mbox{0,84}\,\), en
accord très proche avec la probabilité moyenne d'extinction calculée auparavant.
2.3 Lien avec l'espérance de la population
\(p(t,n)\,\) est la probabilité que la population soit de taille n au temps
t. Le seuil pour w est le même que pour l'espérance de la population
\[\mathcal{E}(t)=\sum_{n\geq 1} n\, p(t,n)\] pour une histoire environnementale
choisie au hasard. En effet, \[\frac{d\mathcal{E}}{dt} =
(a(t)-b(t))\mathcal{E}(t).\] Presque sûrement, on a \[\frac{1}{t} \log
\mathcal{E}(t)\mathop{\longrightarrow}_{t\to +\infty}
(a_1-b_1)u_1+(a_2-b_2)u_2.\] Cette limite a le même signe que \(\,R_0-1\).
Une autre manière de voir cela est de considérer la suite d'environnements
\[1\to 2\to 1\to 2\cdots\] \(\,\tau_n^{(k)}\,\) est la durée aléatoire passée
dans l'environnement k (k=1 ou 2) lors du cycle n (n>=1). En d'autres termes,
l'environnement est
* dans l'état 1 pendant \(\tau_1^{(1)}\) unités de temps,
* puis dans l'état 2 pendant \(\tau_1^{(2)}\) unités de temps,
* puis dans l'état 1 pendant \(\tau_2^{(1)}\) unités de temps, etc.
Après N périodes \(\,1\to 2\,\), l'espérance de la population engendrée par un
individu est \[M_N = \exp \left (\sum_{n=1}^N (a_1-b_1) \tau_n^{(1)} + (a_2-b_2)
\tau_n^{(2)} \right ).\] On a donc \begin{align*} \frac{\log M_N}{N} &= (a_1-b_1)
\frac{\sum_{n=1}^N \tau_n^{(1)}}{N}+(a_2-b_2) \frac{\sum_{n=1}^N
\tau_n^{(2)}}{N}\\ &\mathop{\longrightarrow}_{N\to +\infty}
\frac{a_1-b_1}{q_1}+\frac{a_2-b_2}{q_2}=\frac{(a_1-b_1)q_2+(a_2-b_2)q_1}{q_1q_2}\
, . \end{align*} La limite a le même signe que \(R_0-1\) parce que
\(u_1=q_2/(q_1+q_2)\) et \(u_2=q_1/(q_1+q_2)\). On a donc \(\,M_N\to 0\) si \(R_0
< 1\) et \(\,M_N\to +\infty\) si \(R_0 > 1\).
Noter cependant que ces remarques ne sont pas directement liées à la
proposition 1 car elle donne des informations sur l'espérance de la population et
non sur la probabilité que cette population s'éteigne. En réalité, pour des
processus de branchement en temps discret dans un environnement aléatoire, la
population peut très bien être sous-critique et tendre vers l'extinction presque
sûrement, même si l'espérance de la population converge vers l'infini (Haccou et
coll., 2005, p. 51).
2.4 Autres remarques
Un autre paramètre lié à la croissance d'une espérance mathématique.
(Gray et coll., 2012) ont montré que la position de (3) par rapport à 1 sert
de seuil entre l'extinction et la persistance pour un modèle épidémique de type
SIS constitué d'équations différentielles dans un environnement aléatoire. Il n'y
a pas de stochasticité démographique. Gray et coll. (2012) utilisent pour (3) la
notation \(\,T_0\). Britton et Lindholm (2009) utilisent la notation
\(\tilde{R}_0\,\).
Ces deux références utilisent la notation \(\,R_0\,\) pour un nombre
différent, à savoir le rayon spectral de (2) dans le cas de deux environnements.
Nommons-le \(R^*\,\) pour éviter la confusion. Dans le cas des modèles en temps
discret, (Bacaër et Khaladi, 2012) nomme ce nombre \(\,R_*\). Comme expliqué
ci-dessous, la position de ce nombre par rapport à 1 décide si une certaine
espérance croît ou décroît. \(R^*\,\) s'obtient en général de la manière
suivante. Pour un calcul informel dans le cas particulier K=2, voir la section 3
de (Gray et coll., 2012). Considérons une fois encore la chaîne de Markov (6) en
temps continu sur l'espace d'états \(\{1,2,\ldots,K\}\times \mathbb{N}\). On
définit \[E_k(t)=\sum_{n\geq 1} n\, p_{k,n}(t).\] Alors on voit facilement que
\[\frac{dE_k}{dt} = (a_k-b_k) E_k + \sum_{\ell\neq k} \left (Q_{k,\ell} E_{\ell}
- Q_{\ell,k} E_{k} \right )\, .\] On définit
\[A=\mathrm{diag}(a_1,\ldots,a_K),\quad B=\mathrm{diag}(b_1,\ldots,b_K)\] et
\(E(t)=(E_1(t),\ldots,E_K(t))\), le vecteur des espérances. On a alors
\[\frac{dE}{dt} = (A-B+Q)E.\] D'après (Diekmann et coll., 2013), le vecteur des
espérances converge vers l'infini si et seulement si le rayon spectral \(R^*\) de
\(A(B-Q)^{-1}\) est strictement supérieur à 1. Pour l'exemple numérique ci-dessus
avec \(a_1= \mbox{2,7}\), on obtient \[R^*=\rho(A(B-Q)^{-1})=\rho \left (
\begin{array}{cc} \frac{a_1(b_2+q_2)}{b_1b_2+q_1b_2+q_2b_1} & \frac{a_1
q_2}{b_1b_2+q_1b_2+q_2b_1} \\ \frac{a_2 q_1}{b_1b_2+q_1b_2+q_2b_1} &
\frac{a_2(b_1+q_1)}{b_1b_2+q_1b_2+q_2b_1} \end{array}\right ) \simeq \mbox{1,057}
> 1.\] Rappelons dans ce cas que \(R_0= \mbox{0,875} < 1\,\) et w=1. Donc la
position de \(\,R^*\) par rapport à 1 décide de la croissance de l'espérance mais
ne donne pas le bon seuil pour l'extinction.
Encore un autre paramètre. Comme déjà noté ci-dessus, la suite des
environnements est périodique quand il n'y a que deux environnements possibles:
\(\,1\to 2\to 1\to 2\cdots\). Si l'environnement est k, alors l'espérance de la
population engendrée par un seul individu est \(\,e^{(a_k-b_k)t}\) après t unités
de temps. L'environnement k dure entre t et \(\,t+dt\) unités de temps avec une
probabilité \(q_k\, e^{-q_k t} dt\). Donc l'espérance de la population engendrée
par un individu dans l'environnement k est \begin{equation} m_k=\int_0^\infty
q_k\, e^{-q_k t} e^{(a_k-b_k)t} dt = \frac{q_k}{b_k+q_k-a_k} =
\frac{1}{1-\frac{a_k-b_k}{q_k}}, \quad k\in \{1,2\} \end{equation} pourvu que
\(b_k+q_k > a_k\). Cette condition est vérifiée si \(\,a_1=\mbox{2,7}\). On a
\(\,b_1+q_1-a_1= \mbox{0,3}\) et \(b_2+q_2-a_2= \mbox{2,2}\). Noter que
\(m_k=+\infty\) lorsque \(b_k+q_k\leq a_k\). On définit \(\,R_\star=m_1 m_2\) (à
ne pas confondre avec le \(R^*\,\) de la précédente remarque). On a alors
\(\,R_\star=1/ \mbox{0,66} \simeq \mbox{1,52} > 1\). Britton et Lindholm (2009,
Section 3) ont suggéré que \(\,R_\star > 1\) était équivalent à \(\omega < 1\).
Mais ici on a \(R_\star > 1\) tandis que \(\omega=1\). Donc on voit que la
position de \(\,R_\star\,\) par rapport à 1 ne donne pas le bon seuil pour
l'extinction. Si K=2, le §5.1 de Britton et Lindholm (2009) montre (avec nos
notations) que \(R_\star > 1\) si et seulement si \(R^* > 1\).
3. Les modèles à plusieurs types
On étudie les processus linéaires de naissance et de mort &avec plusieurs
types d'individus dans un environnement variable dans le temps.
\(p(t,n_1,\ldots,n_m)\) est la probabilité d'avoir \(n_i\,\) personnes de type i
au temps t (\(1\leq i\leq m\)). La fonction génératrice
\[g(t,x_1,\ldots,x_m)=\sum_{n_1,\ldots,n_m\geq 0} p(t,n_1,\ldots,n_m)\
x_1^{n_1}\ldots x_m^{n_m}\] est solution de l'équation \begin{equation}\tag{7}
\frac{\partial g}{\partial t}=\sum_{i,j} \bigl [A_{i,j}(t) x_j -
B_{i,j}(t)](x_i-1) \, \frac{\partial g}{\partial x_j} \end{equation} (Bacaër et
Ait Dads, 2012). La matrice des naissances \(A(t)=(A_{i,j}(t))_{1\leq i,j\leq
m}\,\) est à coefficients positifs ou nuls. La matrice des mortalités
\(B(t)=(B_{i,j}(t))\) est de la forme \begin{equation}
B_{i,j}(t)=-b_{i,j}(t)\quad \forall i\neq j,\quad
B_{j,j}(t)=b_{j,j}(t)+\sum_{i\neq j} b_{i,j}(t)\quad \forall j,\quad
b_{i,j}(t)\geq 0 \quad \forall i,j. \end{equation} On suppose que les matrices
\((A(t),B(t))\) appartiennent à une liste finie d'environnements
\(((A^{(k)},B^{(k)}))_{1\leq k\leq K}\), c'est-à-dire \(A(t)=A^{(\theta(t))}\) et
\(B(t)=B^{(\theta(t))}\) avec \(\theta(t)\) un processus stochastique à valeurs
dans \(\{1,2,\ldots,K\}\). On suppose encore une fois que les transitions entre
les environnements suivent une chaîne de Markov homogène en temps continu.
\(Q_{k,\ell}\) est le taux auquel l'environnement peut basculer de (\(\ell\to
k\), \(k\neq \ell\)). Soit Q la matrice de transition correspondante avec
\[Q_{\ell,\ell}=-q_\ell,\quad q_\ell=\sum_{k\neq \ell} Q_{k,\ell}.\] On suppose
que la matrice Q soit irréductible. Par conséquent, il y a une unique
distribution stationnaire u telle que \(\,Qu=0\) et \(\sum_k u_k=1\). Enfin on
suppose que \(b_{j,j}^{(k)} > 0\,\) pour tout k et j. Cette hypothèse implique
que le plus grand exposant de Lyapunov du système différentiel aléatoire
\(\,dZ/dt = -B(t) Z(t)\) est strictement négatif.
Au temps t=0, on suppose que
* il y a \(\nu_i\,\) personnes de type i, avec \(\nu_i\in \mathbb{N}\)
* \(\exists\, i,\ \nu_i > 0\).
On a alors \begin{equation} g(0,x_1,\ldots,x_m)=x_1^{\nu_1}\cdots x_m^{\nu_m}.
\end{equation} L'objectif est de calculer la probabilité d'extinction \[\omega=
\lim_{t\to +\infty} p(t,0,\ldots,0) = \lim_{t\to +\infty} g(t,0,\ldots,0).\]
C'est une variable aléatoire qui dépend de l'histoire environnementale.
Comme l'expliquent par exemple Bacaër et Ait Dads (2012), w peut être calculé
en utilisant les courbes caractéristiques de (7). Avec \(\,\tau\geq 0\,\),
\(Y^{(\tau)}\) est la solution unique du système \begin{equation}\tag{8}
\frac{dY_j^{(\tau)}}{ds}(s) = \sum_i \left [A_{i,j}(-s)\, \left
(1-Y^{(\tau)}_j(s)\right ) - B_{i,j}(-s) \right ] Y^{(\tau)}_i(s) \end{equation}
avec la condition initiale \[Y_j^{(\tau)}(-\tau)=1\quad \forall \, j.\] On a
alors \[\omega= (\omega_1)^{\nu_1} \cdots (\omega_m)^{\nu_m},\quad \quad
\omega_j=1-\lim_{\tau \to +\infty} Y_j^{(\tau)}(0)\, .\] La question est de
savoir si \(\omega=1\) ou \(\omega < 1\). Le résultat dépend de la stabilité du
système suivant d'équations différentielles aléatoires (Arnold, 1998, Sect. 2.2)
\begin{equation}\tag{9} \frac{dX}{dt}=(A(t)-B(t))X(t). \end{equation} C'est
l'équation satisfaite par le vecteur des espérances des populations au temps t.
Cette stabilité dépend du signe du plus grand exposant de Lyapunov de (9)
\[\lambda_1(A,B);\] En suivant Bacaër et Khaladi (2012), la stabilité peut
alternativement être formulée en termes de reproductivité \(\,R_0\), qui est
l'unique solution de \begin{equation} \lambda_1(A/R_0,B)=0. \end{equation}
Une manière d'étudier la probabilité d'extinction w serait d'adapter la
méthode utilisée par Bacaër et Ait Dads (2012) pour des environnements
périodiques au cas des environnements aléatoires. On prendrait avantage du fait
que le système (8) est coopératif et sous-homogène comme dans les travaux de
(Chueshov, 2002) ou (Benaïm et Schreiber, 2009). Ceci conduirait à des
difficultés techniques telles que le lien entre le plus grand exposant de
Lyapunov de (9) et le plus grand exposant de Lyapunov de la linéarisation près de
zéro de (8), qui est le système transposé de (9). Voir (Arnold et Wihstutz, 1986)
ou (Barreira et Valls, 2008). La preuve de la persistance de (8) quand \(\,R_0 >
1\,\) peut aussi être difficile. Pour éviter ces difficultés, on va utiliser la
même idée que dans la seconde preuve de la proposition 1. Pour un nombre fini
d'environnements markoviens, le problème en temps continu peut être ramené à un
processus de branchement avec plusieurs types, en temps discret, dans un
environnement aléatoire. Les résultats d'Athreya et Karlin (1971) peuvent alors
être appliqués.
Proposition 2. On suppose que la matrice suivante est irréductible pour tout
k \[C^{(k)}:=A^{(k)}-B^{(k)}.\] On suppose en plus \[ \forall \, k, \exists \,
(i,j), \quad A^{(k)}_{i,j} > 0.\]
* On a w=1 presque sûrement si \({R}_0\leq 1\)
* On a w 1\) .
Preuve. On prend \(\,t_0=0\). \((t_n)_{n\geq 1}\) avec \(0 < t_1 < t_2 <
\cdots\,\) sont les temps auxquels l'environnement bascule. \(k_n\,\) est
l'environnement dans l'intervalle de temps \(]t_n,t_{n+1}[\,\). Dans
l'environnement k, un individu de type h au départ engendre une population t
unités de temps plus tard dont la fonction génératrice
\(\phi^{(k,h)}(t,x_1,\ldots,x_m)\) est solution de \begin{equation}\tag{10}
\frac{\partial \phi^{(k,h)}}{\partial t}=\sum_{i,j} \bigl [A^{(k)}_{i,j} x_j -
B^{(k)}_{i,j}](x_i-1) \, \frac{\partial \phi^{(k,h)}}{\partial x_j}, \quad
\forall t > 0,\ \forall (x_1,\ldots,x_m)\in ]0,1[^m. \end{equation} avec
\[\phi^{(k,h)}(0,x_1,\ldots,x_m)=x_h.\] On définit \begin{equation}\tag{11}
M_i^{(k,h)}(t)=\frac{\partial \phi^{(k,h)}}{\partial x_i}(t,1,\ldots,1),\quad
M^{(k)}(t)=\left (M_i^{(k,h)}(t)\right )_{i,h}. \end{equation} \(M_i^{(k,h)}(t)\)
est l'espérance de la population de type i. En partant de (10) ou en se référant
à (Athreya et Ney, 1972), on voit (cf. appendice) que \begin{equation}\tag{12}
\frac{dM^{(k)}}{dt}(t) = C^{(k)} M^{(k)}(t), \quad \forall t > 0, \end{equation}
et \(M^{(k)}(0)=I\,\) (la matrice identité). Par conséquent, \begin{equation}
\mathcal{M}_n:=M^{(k_n)}(t_{n+1}-t_n) = \exp\left [C^{(k_n)}(t_{n+1}-t_{n})\right
]\, . \end{equation} Noter avec (9) et (12) que \[\lambda_1(A,B) =\lim_{n\to
+\infty} \frac{1}{t_n} \log \| \mathcal{M}_{n-1} \mathcal{M}_{n-2}\cdots
\mathcal{M}_0\| \] presque sûrement. D'après Athreya et Karlin (1971, Théorème
12), le signe de cette limite décide s'il y a extinction presque sûrement ou pas.
Mais d'abord il faut vérifier les trois hypothèses de ce théorème.
L'irréductibilité de \(\,C^{(k_n)}\) implique que tous les coefficients de
\(\mathcal{M}_n\) sont strictement positifs (Berman et Plemmons, 1994, Théorème
6.3.12) : la première hypothèse est satisfaite. On définit maintenant
\begin{equation}\tag{13} S^{(k,h)}_{i,j}(t)=\frac{\partial^2
\phi^{(k,h)}}{\partial x_i \partial x_j}(t,1,\ldots,1)\ ,\quad S^{(k,h)}(t)=\left
(S^{(k,h)}_{i,j}(t)\right )_{i,j} . \end{equation} On peut montrer que tous les
coefficients de la matrice \(S^{(k_n,h)}(t_{n+1}-t_n)\,\) sont aussi strictement
positifs (cf. appendice): c'est la seconde hypothèse. Enfin on a aussi \[ -
\sum_k \int_0^\infty \!\! w_{k,t} \log \left [\sum_{h=1}^m \left (
1-\phi^{(k,h)}(t,0,\ldots,0) \right ) \right ] dt < +\infty\, ,\] avec
\(w_{k,t}\) donné par (5). En effet, d'une part, \(w_{k,t}\) décroît
exponentiellement si \(t\to +\infty\). Dautre part,
\(1-\phi^{(k,h)}(t,0,\ldots,0)\) ne peut s'approcher de 0 plus vite que
\(e^{-ct}\,\) avec c>0. c est donné par le taux auquel la solution de (8) peut
s'approcher de 0 dans un environnement k qui est sous-critique. Donc la troisième
condition est aussi satisfaite.
Dans le cas \(R_0\leq 1\), on a \(\lambda_1(A,B)\leq 0\). On conclut que w=1
presque sûrement grâce à
* (Athreya et Karlin, 1971, Théorème 12(i)) si \(\,\lambda_1(A,B) < 0\)
* (Kaplan, 1974, Théorème 2) si \(\,\lambda_1(A,B)=0\).
Dans le cas \(R_0 > 1\), on a \(\lambda_1(A,B) > 0\). On conclut que w 0\)
et \(q_2 > 0\). La distribution stationnaire est telle que
\(\,u_1=q_2/(q_1+q_2)\) et \(u_2=q_1/(q_1+q_2)\). On suppose \[A(t)=\left
(\begin{array}{ccc} 0 & &\beta(t)\\ 0 & & 0 \end{array} \right ),\quad B(t)=\left
(\begin{array}{cc} \alpha+\mu & 0\\ -\alpha & \gamma+\mu \end{array}\right ), \]
avec \(\beta(t)=\beta_1 > 0\) ou \(\beta_2 > 0\) selon l'environnement. \(\alpha
> 0\) est le taux auquel les gens infectés mais pas encore infectieux deviennent
infectieux, \(\mu > 0\) est la mortalité et \(\gamma > 0\,\) est le taux de
guérison. La reproductivité est l'unique nombre positif pour lequel le plus grand
exposant de Lyapunov du système suivant est égal à 0 \[\frac{dX}{dt} =
(A(t)/R_0-B(t))X(t).\] Noter que si \(\beta(t)\) était constant et égal à sa
moyenne temporelle \(u_1\beta_1+u_2\beta_2\), on aurait
\(R_0=(u_1\beta_1+u_2\beta_2)\alpha/((\alpha+\mu)(\gamma+\mu))\). Des formules
analytiques approchées pour \(\,R_0\) dans un environnement aléatoire peuvent
être déduites des formules de (Arnold et Kloeden, 1989) pour le plus grand
exposant de Lyapunov des systèmes bidimensionnels.
La probabilité pour que le processus soit éteint au temps t>0, en partant de
\((E_0,I_0)\) personnes au temps 0, est \[(1-Y_1^{(\tau)}(0))^{E_0}
(1-Y_2^{(\tau)}(0))^{I_0},\] avec pour \(-\tau < s < 0\) \begin{align}
\frac{dY_1^{(\tau)}}{ds}(s) &= -(\alpha+\mu) Y_1^{(\tau)}(s) + \alpha\,
Y_2^{(\tau)}(s),\tag{14}\\ \frac{dY_2^{(\tau)}}{ds}(s) &= \beta(-s)\,
Y_1^{(\tau)}(s)\, (1-Y_2^{(\tau)}(s)) -(\gamma+\mu) Y_2^{(\tau)}(s),\tag{15}
\end{align} tandis que \(Y_1^{(\tau)}(-\tau)=1\) et \(Y_2^{(\tau)}(-\tau)=1\). Il
y a des erreurs de signes dans les équations correspondantes données par Bacaër
et Ait Dads (2012) ; les figures 3 et 4 dans cette référence sont néanmoins
correctes. \(\omega_1\) et \(\omega_2\) sont les probabilités d'extinction
ultimes, en partant soit d'une personne infectée mais non infectieuse, soit d'une
personne infectieuse: \[\omega_j=\lim_{\tau\to +\infty} 1-Y_j^{(\tau)}(0),\quad
j=1,2.\] La proposition 2 montre \(\omega_1=\omega_2=1\) presque sûrement si
\(R_0\leq 1\) et que \(\omega_1 < 1\) et \(\omega_2 < 1\) presque sûrement si
\(R_0 > 1\).
Dans le cas \(\beta_2 \leq \beta_1\,\), un théorème de comparaison pour le
système (14)-(15) montre que \(\omega_1\) et \(\omega_2\) sont supérieures aux
probabilités correspondantes pour le processus de naissance et de mort où
l'environnement est toujours 1. Si ce dernier processus est surcritique alors ces
probabilités \(\xi_1\) et \(\xi_2\) se calculent facilement :
* soit on détermine l'état stationnaire de (14)-(15) avec \(\beta(-s)\)
remplacé par \(\beta_1\),
* soit on considère le processus de Bienaymé-Galton-Watson induit, avec
plusieurs types :
\[\xi_1=\frac{\mu}{\alpha+\mu}+\frac{\gamma+\mu}{\beta_1},\quad
\xi_2=\frac{\alpha+\mu}{\alpha}\, \frac{\gamma+\mu}{\beta_1}.\]
Prenons \(q_1=q_2=1\), \(\beta_1=2\), \(\beta_2=1\), \(\alpha=1\), \(\mu=
\mbox{0,01}\) et \(\gamma=1\). On obtient \(\,R_0\simeq \mbox{1,45} > 1\). Noter
que pour le système moyenné en temps, on a \(\,R_0\simeq \mbox{1,47}\). De plus,
on obtient \(\xi_1\simeq \xi_2\simeq \mbox{0,51}\). La figure 2 montre
l'histogramme pour la probabilité d'extinction partant d'une personne infectée
mais non infectieuse dans l'environnement 1. Il a été obtenu avec 1000
simulations de l'histoire environnementale. La moyenne est proche de 0,69 (elle
serait 0,66 partant de l'environnement 2). On a pris t=100 et un pas de temps de
0,001. La figure confirme que \(\omega_1 < 1\) (et \(\omega_2 < 1\)) presque
sûrement si \(R_0 > 1\).
[2014JMB2Figure2.png] Figure 2. Histogramme de la probabilité d'extinction
\(\omega_1\) partant d'une personne infectée mais non infectieuse dans
l'environnement 1.
4. Conclusion
Certaines questions restent ouvertes. On peut se demander, comme Britton et
Lindholm (2009), ce qui arrive si la survie n'est pas distribuée
exponentiellement, c'est-à-dire pour des processus de Crump-Mode-Jagers tels que
dans le travail de (Ball et Donnelly, 1995). Il est aussi clair que la plupart de
résultats restent vrais pas seulement pour un nombre fini d'environnements
markoviens mais aussi pour des environnements ergodiques plus généraux.
Pour les applications biologiques, il serait plus réaliste de supposer que la
matrice de transition est de la forme \[Q(t)=\left (\begin{array}{cc} -q_1(t) &
q_2(t)\\ q_1(t) & -q_2(t) \end{array}\right ),\] où par exemple
\(q_1(t)=k_1(1+\varepsilon_1 \cos \omega t)\), \(q_2(t)=k_2(1+\varepsilon_2 \sin
\omega t)\), \(k_1 > 0\), \(k_2 > 0\), \(\varepsilon_1\in (0,1)\) et
\(\varepsilon_2 \in (0,1)\). De cette manière, l'année est plus ou moins divisée
en deux saisons (disons l'été et l'hiver), l'une qui serait favorable à la
croissance et l'autre moins favorable. Cela pourrait être le cas pour les
maladies infectieuses où le taux effectif de contact dépend de la température.
Une matrice périodique \(\,Q(t)\) comme ci-dessus est un modèle de saisonnalité
plus réaliste qu'une matrice Q constante. Dans ce dernier cas avec seulement deux
environnements, les deux saisons alternent mais les longueurs des saisons sont
indépendantes. Pour être réaliste il faut qu'un été particulièrement court soit
suivi d'un hiver particulièrement long pour garder plus ou moins la périodicité
annuelle.
Complément
En prenant la dérivée de (10) par rapport à \(x_I\), on obtient \begin{align}
\frac{\partial^2 \phi^{(k,h)}}{\partial t \, \partial x_I} =& \sum_{i,j} \left [
A^{(k)}_{i,j} \delta_{j,I} (x_i-1) + \left (A^{(k)}_{i,j} x_j -
B^{(k)}_{i,j}\right ) \delta_{i,I} \right ] \frac{\partial \phi^{(k,h)}}{\partial
x_j} \nonumber\\ &+ \sum_{i,j} \left (A^{(k_n)}_{i,j} x_j - B^{(k_n)}_{i,j}\right
) (x_i-1) \frac{\partial^2 \phi^{(k,h)}}{\partial x_I \partial x_j}\, .\tag{16}
\end{align} \(\delta\,\) est le symbole de Kronecker. On prend
\(\,x_1=\cdots=x_m=1\). On obtient \[\frac{\partial}{\partial t} \left
[\frac{\partial \phi^{(k,h)}}{\partial x_I}(t,1,\ldots,1)\right ]=\sum_{j} \left
(A^{(k)}_{I,j} - B^{(k)}_{I,j}\right ) \frac{\partial \phi^{(k,h)}}{\partial
x_j}(t,1,\ldots,1).\] Ceci est identique à (12). On prend la dérivée par rapport
à \(\,x_J\,\) de (16). On prend une nouvelle fois \(\,x_1=\cdots=x_m=1\). On
obtient \begin{align} \frac{\partial}{\partial t}\left [\frac{\partial^2
\phi^{(k,h)}}{\partial x_I \partial x_J} (t,1,\ldots,1)\right ] = & A^{(k)}_{J,I}
\frac{\partial \phi^{(k,h)}}{\partial x_I}(t,1,\ldots,1) \nonumber\\ &+
A^{(k)}_{I,J} \frac{\partial \phi^{(k,h)}}{\partial x_J}(t,1,\ldots,1)
\nonumber\\ &+ \sum_{j} \left (A^{(k)}_{I,j} - B^{(k)}_{I,j}\right )
\frac{\partial^2 \phi^{(k,h)}}{\partial x_j \partial
x_J}(t,1,\ldots,1)\nonumber\\ &+ \sum_{j} \left (A^{(k)}_{J,j} -
B^{(k)}_{J,j}\right ) \frac{\partial^2 \phi^{(k,h)}}{\partial x_I \partial
x_j}(t,1,\ldots,1)\, . \tag{17} \end{align} Un calcul semblable se trouve par
exemple dans le § V.7.3 du livre d'Athreya et Ney (1972). Rappelons les notations
(11) et (13) pour les dérivées premières et secondes. On définit
\[G_{i,j}^{(k,h)}(t)=A_{i,j}^{(k)} M^{(k,h)}_j(t)+A_{j,i}^{(k)} M^{(k,h)}_i(t)\ ,
\quad G^{(k,h)}(t)=\left (G_{i,j}^{(k,h)}(t)\right )_{i,j} .\] Alors (17) est
identique à \[\frac{dS^{(k,h)}_{I,J}}{dt}(t) = \sum_j \left [C_{I,j}^{(k)}
S^{(k,h)}_{j,J}(t) + C_{J,j}^{(k)} S^{(k,h)}_{j,I}(t)\right ] +
G_{I,J}^{(k,h)}(t),\quad t > 0.\] On utilise la symétrie de la matrice
\(\,S^{(k,h)}(t)\,\). Cette équation devient \[\frac{d S^{(k,h)}}{dt}(t) =
C^{(k)} S^{(k,h)}(t) + S^{(k,h)}(t) \left (C^{(k)}\right )^* + G^{(k,h)}(t)\, .\]
\(^*\,\) désigne la matrice transposée. Mais
\(\,\phi^{(k,h)}(0,x_1,\ldots,x_m)=x_h.\) On a donc \(S^{(k,h)}(0)=0\). Comme
dans le livre d'Athreya et Ney (1972, formule (15), p. 203), on obtient
\begin{equation}\tag{18} S^{(k,h)}(t)= \int_{0}^{t} e^{(t-u)C^{(k)}} G^{(k,h)}(u)
\left (e^{(t-u)C^{(k)}} \right )^* du, \quad \forall \, t\geq 0. \end{equation}
On remarque
* \(M^{(k,h)}_i(t) > 0\,\) pour tout i et tout t>0,
* par hypothèse, il existe \((i,j)\) avec \(A^{(k)}_{i,j} > 0\).
La matrice à coefficients positifs ou nuls \(G^{(k,h)}(t)\,\) a donc au moins un
coefficient strictement positif pour tout t>0. Tous les coefficients de la
matrice \(\,e^{(t-u)C^{(k)}}\) et de sa transposée sont strictement positifs si
\(\,u < t\). Tous les coefficients de la matrice sous l'integrale de (18) sont
donc strictement positifs pour \(u < t\). Tous les coefficients de la matrice
\(\,S^{(k,h)}(t)\) sont donc aussi strictement positifs pour \(t > 0\).
Remerciements
On remercie le professeur Khaladi pour ses encouragements. On remercie aussi
les professeurs Ball et Bansaye et particulièrement le professeur Jagers pour
leurs commentaires sur certaines parties de ce travail.
Références bibliographiques
* Arnold L (1998) Random dynamical systems. Springer, Berlin.
* Arnold L, Kloeden P (1989) Lyapunov exponents and rotation number of
two-dimensional systems with telegraphic noise. SIAM J Appl Math 49:1242-1274
* Arnold L, Wihstutz V (1986) Lyapunov exponents: a survey. In: Arnold L,
Wihstutz V (eds) Lyapunov exponents, Lecture Notes in Mathematics 1186,
Springer, Berlin, p. 1-26
* Artalejo JR, Economou A, Lopez-Herrero MJ (2012) Stochastic epidemic models
with random environment: quasi-stationarity, extinction and final size. J
Math Biol, doi:10.1007/s00285-012-0570-5
* Athreya KB, Karlin S (1971) On branching processes with random environments:
I Extinction probabilities. Ann Math Stat 42:1499-1520
* Athreya KB, Ney PE (1972) Branching processes. Springer, Berlin
* Bacaër N, Ait Dads E (2012) On the probability of extinction in a periodic
environment. J Math Biol, doi:10.1007/s00285-012-0623-9
* Bacaër N, Guernaoui S (2006) The epidemic threshold of vector-borne diseases
with seasonality. J Math Biol 53:421-436
* Bacaër N, Khaladi M (2012) On the basic reproduction number in a random
environment. J Math Biol, doi:10.1007/s00285-012-0611-0
* Ball F, Donnelly P (1995) Strong approximations for epidemic models. Stoch
Proc Applic 55:1-21
* Barreira L, Valls C (2008) Stability of nonautonomous differential equations.
Springer, Berlin
* Benaïm M, Schreiber SJ (2009) Persistence of structured populations in random
environments. Theoret Popul Biol 76:19-34
* Berman A, Plemmons RJ (1994) Nonnegative matrices in the mathematical
sciences. SIAM, Philadelphia
* Bouleau N (1988) Processus stochastiques et applications. Hermann, Paris
* Britton T, Lindholm M (2009) The early stage behaviour of a stochastic SIR
epidemic with term-time forcing. J Appl Probab 46:975-992
* Chueshov I (2002) Monotone random systems. Springer, Berlin
* Cogburn R, Torrez WC (1981) Birth and death processes with random
environments in continuous time. J Appl Probab 18:19-30
* Diekmann O, Heesterbeek H, Britton T (2013) Mathematical tools for
understanding infectious disease dynamics. Princeton University Press,
Princeton
* Dublin LI, Lotka AJ (1925) On the true rate of natural increase. J Am Stat
Assoc 20(151):305-339
* Gaver DP, Jacobs PA, Latouche G (1984) Birth-and-death models in randomly
changing environments. Adv Appl Probab 16:715-731
* Gray A, Greenhalgh D, Mao X, Pan J (2012) The SIS epidemic model with
Markovian switching. J Math Anal Appl 394:496-516
* Haccou P, Jagers P, Vatutin VA (2005) Branching processes: variation, growth,
and extinction of populations. Cambridge University Press, Cambridge
* Hernandez-Suarez C, Rabinovich J, Hernandez K (2012) The long-run
distribution of births across environments under environmental stochasticity
and its use in the calculation of unconditional life-history parameters.
Theor Popul Biol 82:264-274
* Hillion A (1986) Les théories mathématiques des populations. Presses
Universitaires de France, Paris
* Kaplan N (1973) A continuous time Markov branching model with random
environments. Adv Appl Probab 5:37-54
* Kaplan N (1974) Some results about multidimensional branching processes with
random environments. Ann Probab 2:441-455
* Kendall DG (1948) On the generalised 'birth-and-death' process. Ann Math
Statist 19:1-15
* Lewontin RC, Cohen D (1969) On population growth in a randomly varying
environment. Proc Natl Acad Sci USA 62:1056-1060
* Meyn S, Tweedie R (1993) Markov chains and stochastic stability. Springer,
Berlin
* Pardoux E (2008) Markov processes and applications. Wiley, Chichester
* Yechiali U (1973) A queuing-type birth-and-death process defined on a
continuous-time Markov chain. Operations Research 21:604-609
Usage: http://www.kk-software.de/kklynxview/get/URL
e.g. http://www.kk-software.de/kklynxview/get/http://www.kk-software.de
Errormessages are in German, sorry ;-)