A Stochastic SVIR Model for Measles

Moussa Seydou, Ousmane Moussa Tessa

Department of Mathematics and Informatics, Abdou Moumouni University, Niamey, Niger.

DOI: 10.4236/am.2021.123013                

Abstract

In this article, we consider the construction of a SVIR (Susceptible, Vaccinated, Infected, Recovered) stochastic compartmental model of measles. We prove that the deterministic solution is asymptotically the average of the stochastic solution in the case of small population size. The choice of this model takes into account the random fluctuations inherent to the epidemiological characteristics of rural populations of Niger, notably a high prevalence of measles in children under 5, coupled with a very low immunization coverage.

Keywords

Measles, Compartmental Model, SVIR, Basic Reproductive Number, Markov Chains, Lyapunov Function, Stochastic Stability, Stochastic Simulation, Niger

1. Introduction

The measles is caused by a virus belonging to morbillivirus group. It may infect other primates, but is largely specialized on its human host. It is transmitted by direct contact with an infected person or by air [1] [2]. Upon infection, the patient passes through a latent period of 6 to 9 days, followed by 6 to 7 day infective period [3]. The infection results in either death or full recovery of the host. In the last case, the host develops lifelong immunity. However, immunity can also be acquired by vaccination before infection. Before the introduction of measles vaccine in 1963 and widespread vaccination, major epidemics occurred approximately every two or three years and measles caused an estimated 2.6 million deaths each year [1]. In developing countries, like Niger, measles remains one of the main causes of infant mortality because children under 5 remain the most affected, 90% who die have less than 5 years [1] [4] [5] [6]. In sub-Saharan Africa, especially in areas where vaccination coverage is not optimal, the case fatality is one of the highest 5% - 10%, compared to that of high-income countries, where we have 1 death in this age group out of 1000 measles cases [7] [8] [9].

A fundamental concept that has come out of the measles transmission process is that of the basic reproduction number R₀. It is defined as average number of secondary infections produced when one infected individual is introduced into a host population where everyone is susceptible [10] [11]. R₀ is a threshold parameter in the course of the spread of measles disease; indeed, if R₀ < 1, the disease will eventually disappear from the population, while if R₀ > 1, the disease can spread as an epidemic in the absence of health interventions. In a small, isolated population, a measles epidemic cannot persist [12] [13] [14], even if the basic reproduction number is initially greater than 1. Indeed, the spread of the disease subsides at term, due to a progressive immunization of a growing proportion of the population. Thus, in such a context, measles can only be endemic after regular importation of the virus, generally from infected people from large urban centers [15].

Most model used for infections diseases are the compartmental models, originally introduced by Kermack and Mckendrick and their variants [5] [10] [16] [17]. They are based on the partition of the population into distinct classes (or compartments) according to its epidemiological status. Host can move from one class to another (transition). In the case of the SIR (Susceptible, Infected, Recovered) model, an infection is the transition which moves an individual from the susceptible class to the infected class and a recovery leads an infected person to the recovered compartment [5] [10] [11]. In general, the transition rate, which expresses the probability that an individual passes from one class to another per unit of time, depends essentially on the state of the system at a given moment, in particular on the number of individuals in the different compartments and the disease infection force [6] [10].

In our SVIR stochastic model, we consider R_p the effective reproduction number, characterizing the vaccination effort to control the spread of the disease, where p is the proportion of newborns vaccinated and immunized. In the total absence of vaccination (p= 0) against measles, we estimate R₀ the basic reproduction number [3] [5] between 10 and 18.

The rest of the paper is organized as follows: Section 2 describes in detail the deterministic model SVIR and the equilibrium points of the system of differential equations of the model. In Section 3, we formulate our stochastic SVIR model by means of the Kolmogorov Forward equations, precisely by means of a system of differential equations of the mathematical expectations of the number of susceptible, infected and immune (recovered and vaccinated). Section 4 is devoted to the study of the asymptotic behavior of our stochastic model, followed by numerical simulations in the fifth section. Finally, in the last section, we discuss our stochastic approach and scientific conclusions.

2. The Deterministic SVIR Model

In what follows, $S\left(t\right)\mathrm{,}I\left(t\right)\mathrm{,}R\left(t\right)$ denote respectively the number of susceptible, infected and immunized (susceptible vaccinated and recovered patients) at time t.

In this model, the new susceptibles (newborns) are introduced at a constant rate n. A fraction, pn, of newborns has acquired immunity by vaccination. The other fraction $\left(1-p\right)n$ remains susceptible. In addition, we assume that:

· The natural death rate is δ for each compartment.

· Infectious patients recover at the rate of γ.

· Infectious patients have an additional μ death rate from measles.

· We consider the standard incidence $f\left(I,S\right)=\beta SI$, β is the disease transmission coefficient. βis the average probability of an adequate contact (contact sufficient for transmission) between an infected and a susceptible per unit of time.

In Figure 1, a compartmental diagram of the transitions illustrates the relationship between the three classes.

The dynamics of a well-mixed population can be described by the differential equations:

$\{\begin{array}{l}\frac{\text{d}S}{\text{d}t}=n\left(1-p\right)-\beta SI-\delta S\\ \frac{\text{d}I}{\text{d}t}=\beta SI-\left(\delta +\mu +\gamma \right)I\\ \frac{\text{d}R}{\text{d}t}=np+\gamma I-\delta R\end{array}$ (1)

Remark. 1) In the case of equilibrium without disease, the system (1) admits an equilibrium point $\left({\stackrel{\mathrm{<mo>\; *\; </mo>}}{S}}_0\mathrm{,}{\stackrel{\mathrm{<mo>\; *\; </mo>}}{I}}_0\mathrm{,}{\stackrel{\mathrm{<mo>\; *\; </mo>}}{R}}_0\right)$ with

${\stackrel{\mathrm{<mo>\; *\; </mo>}}{S}}_0=\frac{\left(1-p\right)n}{\delta },{\stackrel{\mathrm{<mo>\; *\; </mo>}}{I}}_0=0\text{et}{\stackrel{\mathrm{<mo>\; *\; </mo>}}{R}}_0=\frac{np}{\delta }$ (2)

Setting $R_0=\frac{\beta n}{\delta \left(\delta +\mu +\gamma \right)}$ et $R_p=\left(1-p\right)R_0$, this equilibrium point is asymptotically stable [18] if $R_p<1$. In addition, we have $R_p<R_0$ et $R_p<1$ if and only if $p>1-\frac{1}{R_0}$. We say that $p_c=1-\frac{1}{R_0}$ is the critical vaccination coverage of newborns.

2) If $R_p>1$, an endemic point of equilibrium appears $\left({\stackrel{\mathrm{<mo>\; *\; </mo>}}{S}}_e\mathrm{,}{\stackrel{\mathrm{<mo>\; *\; </mo>}}{I}}_e\mathrm{,}{\stackrel{\mathrm{<mo>\; *\; </mo>}}{R}}_e\right)$ asymptotically stable [18], where

${\stackrel{\mathrm{<mo>\; *\; </mo>}}{S}}_e=\frac{\delta +\mu +\gamma }{\beta },{\stackrel{\mathrm{<mo>\; *\; </mo>}}{I}}_e=\frac{\left(R_p-1\right)\delta }{\beta }\text{et}{\stackrel{\mathrm{<mo>\; *\; </mo>}}{R}}_e=\frac{np\beta +\gamma \delta \left(R_p-1\right)}{\delta \beta }$ (3)

Figure 1. Compartment diagram of model SVIR.

3. The Continuous Stochastic SVIR Model

Let $X_t={\left(S\left(t\right)\mathrm{,}I\left(t\right)\right)}_{t\ge 0}$ be a continuous-time homogeneous Markov chain on the denumerable state space ${\mathbb{N}}^2={\left\{\mathrm{0,1,2,}\cdots \right\}}^2$. First, assume that $\Delta t$ can be chosen sufficiently small such that at most one change in state occurs during the time interval $\Delta t$. In particular, there can be either a new infection, a birth, a death, or a recovery. From of state $\left\{X_t=\left(s,i\right)\right\}$, only the following states are accessible:

$\left(s\mathrm{,}i\right)\mathrm{<mo>\; ;\; </mo>}\left(s+\mathrm{1,}i\right)\mathrm{<mo>\; ;\; </mo>}\left(s\mathrm{,}i-1\right)\mathrm{<mo>\; ;\; </mo>}\left(s-\mathrm{1,}i\right)\mathrm{<mo>\; ;\; </mo>}\left(s-\mathrm{1,}i+1\right)\mathrm{.}$

corresponding to the possible transitions starting from the state $\left(s\mathrm{,}i\right)$. (See Figure 2). $X_t$ has an absorbing set corresponding to disease-free equilibrium states $E_0=\left\{\left(s,i\right),s\ge 0;i=0\right\}$.

Let $V_{\left(s,i\right)}$ be the set of neighbors of state $\left(s,i\right)$ :

$V_{\left(s,i\right)}=\left\{\left(s+\mathrm{1,}i\right)\mathrm{<mo>\; ;\; </mo>}\left(s-\mathrm{1,}i+1\right)\mathrm{<mo>\; ;\; </mo>}\left(s-\mathrm{1,}i\right)\mathrm{<mo>\; ;\; </mo>}\left(s\mathrm{,}i-1\right)\right\}$

Setting ${\tau }_{\left(s,i\right)}=n\left(1-p\right)+\beta is+\delta s+\left(\mu +\delta +\gamma \right)i$, the transition rates are defined by:

${\tau }_{\left(s,i\right),\left(k,l\right)}=\{\begin{array}{ll}n\left(1-p\right)\hfill & \left(k,l\right)=\left(s+1,i\right),s\ge 0,i\ge 0\hfill \\ \beta is\hfill & \left(k,l\right)=\left(s-1,i+1\right),s\ge 1,i\ge 0\hfill \\ \delta s\hfill & \left(k,l\right)=\left(s-1,i\right),s\ge 1,i\ge 0\hfill \\ \left(\mu +\delta +\gamma \right)i\hfill & \left(k,l\right)=\left(s,i-1\right),s\ge 0,i\ge 1\hfill \end{array}$ (4)

The transition probabilities of $X_t=\left(S\left(t\right)\mathrm{,}I\left(t\right)\right)$ are defined by

$P_{\left(s,i\right),\left(k,l\right)}\left(\Delta t\right)=\mathbb{P}\left\{X_{t+\Delta t}=\left(k,l\right)/X_t=\left(s,i\right)\right\}$

We have $\forall s\ge 0$,

$P_{\left(s,i\right),\left(k,l\right)}\left(\Delta t\right)=\{\begin{array}{l}\forall i>0,{\tau }_{\left(s,i\right),\left(k,l\right)}\Delta t+o\left(\Delta t\right)\text{if}\left(k,l\right)\in V_{\left(s,i\right)}\\ 1-{\tau }_{\left(s,i\right)}\Delta t+o\left(\Delta t\right)\text{if}\left(k,l\right)=\left(s,i\right)\\ \forall i=0,P_{\left(s,0\right),\left(s,0\right)}\left(\Delta t\right)=1\end{array}$ (5)

Figure 2. States transition.

The distribution of $X_t$ is $P_{s,i}\left(t\right)=0$ if $s<0$ or $i<0$ and $P_{s,i}\left(t\right)=\mathbb{P}\left\{X_t=\left(s,i\right)\right\}$ if $s\ge 0,i\ge 0$. Therefore, the marginal distributions are given by:

$\mathbb{P}\left\{I\left(t\right)=i\right\}={\displaystyle \underset{s\ge 0}{\sum }}{P}_{s,i}\left(t\right)\text{and}\mathbb{P}\left\{S\left(t\right)=s\right\}={\displaystyle \underset{i\ge 0}{\sum }}{P}_{s,i}(\; t\; )$

From the Equation (5), we obtain the Kolmogorov Forward equations, for all $s\ge 0$ and $i\ge 0$

$\begin{array}{c}\frac{\text{d}{P}_{s,i}}{\text{d}t}=n\left(1-p\right)\left[P_{s-1,i}-P_{s,i}\right]+\beta \left[\left(s+1\right)\left(i-1\right)P_{s+1,i-1}-si{P}_{s,i}\right]\\ +\left(\mu +\gamma +\delta \right)\left[\left(i+1\right)P_{s,i+1}-i{P}_{s,i}\right]+\delta \left[\left(s+1\right)P_{s+1,i}-s{P}_{s,i}\right]\end{array}$ (6)

Hence the system of differential equations verified by the mathematical expectations:

$\{\begin{array}{l}\frac{\text{d}\bar{S}}{\text{d}t}=\left(1-p\right)n-\beta \bar{S}\bar{I}-\delta \bar{S}-\beta co{v}_{SI}\\ \frac{\text{d}\bar{I}}{\text{d}t}=\beta \bar{S}\bar{I}-\left(\mu +\delta +\gamma \right)\bar{I}+\beta co{v}_{SI}\\ \frac{\text{d}\bar{R}}{\text{d}t}=np+\gamma \bar{I}-\delta \bar{R}\end{array}$ (7)

$\begin{array}{l}\bar{S}\left(t\right)={\displaystyle \underset{s=0}{\overset{+\infty }{\sum }}{\displaystyle \underset{i=0}{\overset{+\infty }{\sum }}s{P}_{s\mathrm{,}i}}\left(t\right)}\mathrm{,}\bar{I}\left(t\right)={\displaystyle \underset{s=0}{\overset{+\infty }{\sum }}{\displaystyle \underset{i=0}{\overset{+\infty }{\sum }}i{P}_{s\mathrm{,}i}\left(t\right)}}\\ co{v}_{SI}\left(t\right)={\displaystyle \underset{s=0}{\overset{+\infty }{\sum }}{\displaystyle \underset{i=0}{\overset{+\infty }{\sum }}si{P}_{s\mathrm{,}i}\left(t\right)}}-\bar{S}\left(t\right)\bar{I}\left(t\right)\text{and}\bar{R}\left(t\right)={\displaystyle \underset{r=0}{\overset{+\infty }{\sum }}r\mathbb{P}\left\{R\left(t\right)=r\right\}}\end{array}$

4. Asymptotic Behavior

In this part, we establish that the extinction of the epidemic is done almost surely independently of the number R_p, although this is not a priori guaranteed in infinite dimension.

Let us consider the embedded process ${\left(Y_k\right)}_{k\in \mathbb{N}}$ of ${\left(X_t\right)}_{t\ge 0}$ which is a discrete Markov chain representing the sequence of values taken by ${\left(X_t\right)}_{t\ge 0}$ at transition times.

Setting $\Delta Y_k=Y_{k+1}-Y_k$, we have:

$\mathbb{P}\left\{\Delta Y_k=\left(e_1,e_2\right)/Y_k=\left(s,i\right)\right\}=\{\begin{array}{l}\frac{n\left(1-p\right)}{\tau \left(s,i\right)},\text{if}\left(e_1,e_2\right)=\left(1,0\right)\\ \frac{\beta si}{\tau \left(s,i\right)},\text{if}\left(e_1,e_2\right)=\left(-1,1\right)\\ \frac{\delta s}{\tau \left(s,i\right)},\text{if}\left(e_1,e_2\right)=\left(-1,0\right)\\ \frac{\left(\mu +\gamma +\delta \right)i}{\tau \left(s,i\right)},\text{if}\left(e_1,e_2\right)=\left(0,-1\right)\end{array}$ (8)

note that $\tau \left(s,i\right)=n\left(1-p\right)+\beta si+\delta s+\left(\mu +\gamma +\delta \right)i$.

To establish our results, we need the proposition [1] and the lemmas [2] [3] [4] [5] which are obtained according to the proof of the criterion of ergodicity and recurrence of Markov chains, given by Rosenkrantz [19]. These assertions are essentially based on the Lyapunov-Foster ergodicity criterion, which shows that a Markov chain is recurrent positive. This criterion was subsequently extended by Meyn and Tweedie [20] [21]. The proofs of the lemmas are given in the Appendix.

Proposition 1. Setting $s_0=\max \left(\frac{n\left(1-p\right)}{\delta },\frac{\mu +\gamma +\delta }{\beta }\right)$ ; $D_0=\left\{\left(s,i\right),s>s_0,i>0\right\}$,

$D_1=\left\{\left(s,0\right),s>s_0\right\}$, $D_2=\left\{\left(s_0,i\right),i>0\right\}$ and let $d\left(s\mathrm{,}i\right)=\mathbb{E}\left[\Delta Y_k\mathrm{<mo>\; /\; </mo>}{Y}_k=\left(s\mathrm{,}i\right)\right]$ be the drift vector and $d_j\left(s,i\right)=d\left(s,i\right)$ $\left(s,i\right)\in D_j,0\le j\le 2$. then

$d\left(s,i\right)=\left(\frac{n\left(1-p\right)-\beta si-\delta s}{\tau \left(s,i\right)},\frac{\beta si-\left(\mu +\gamma +\delta \right)i}{\tau \left(s,i\right)}\right)$ (9)

Lemma 2. For all $j\in \left\{\mathrm{0,1,2}\right\}$, we pose $d_j=d\left(s\mathrm{,}i\right)$ where $\left(s\mathrm{,}i\right)\in D_j$. We denote by $\psi =\widehat{\left(n_1\mathrm{,}{d}_0\right)}$ the angle between $n_1$ and $d_0$, ${\psi }_1=\widehat{\left(n_1,d_1\right)}$ the angle between $n_1$ and $d_1$, ${\psi }_2=\widehat{\left(n_2\mathrm{,}{d}_2\right)}$ the angle between $n_2$ and $d_2$, where $n_1=\left(\mathrm{0,1}\right)$ and $n_2=\left(\mathrm{1,0}\right)$.

We have the following results:

1) $0<\psi <{\psi }_1=\frac{\pi }{2}<{\psi }_2\le \pi$.

2) If $R_p\le 1$ : $s_0=\frac{\mu +\gamma +\delta }{\beta }$, ${\psi }_2=\pi$ and $\frac{\pi }{4}<\psi <{\psi }_1=\frac{\pi }{2}$

3) If $R_p>1$ : $s_0=\frac{n\left(1-p\right)}{\delta }$, $\frac{\pi }{2}<{\psi }_2<\pi$.

Proof: See Appendix.

Definition 4.1. Let $\varphi \left(r\mathrm{,}\theta \right)=r^{\alpha }\cos \left(\alpha \theta -{\theta }_1\right)$ where $\alpha =\frac{2\left({\theta }_1+{\theta }_2\right)}{\pi }$, for all reals $r\ge 0$ and $\theta \in \left[\mathrm{0,}\frac{\pi }{2}\right]$ with

· ${\theta }_1\in \left]\mathrm{0,}\frac{\pi }{4}\right[$ and ${\theta }_2\in \left]\frac{\pi }{2}-{\theta }_1\mathrm{,}\frac{\pi }{2}\right[$, in the case where $R_p\le 1$.

· ${\theta }_1\in \left]-\frac{\pi }{2}\mathrm{,}\inf \left(\psi -\frac{\pi }{2}\mathrm{,}-\frac{\pi }{4}\right)\right[$ and ${\theta }_2\in \left]-{\theta }_1\mathrm{,}\frac{\pi }{2}\right[$, in the case where $R_p>1$.

$\psi =\widehat{\left(n_1\mathrm{,}{d}_0\right)}$ is the angle between $n_1$ and $d_0$,

We say that $\varphi$ is the Lyapounov function intervening in the study of the recurrence-transience of $X_t$.

Remark. If $R_p\le 1$, we obtain $1<\alpha <2$, whereas if $R_p>1$, $0<\alpha <1$.

Lemma 3. Let $\varphi$ be the Lyapounov function. For all reals $r\ge 0$ and $\theta \in \left[\mathrm{0,}\frac{\pi }{2}\right]$, we have the following results:

1) $\nabla \varphi \left(r,\theta \right)\cdot d_0<0$, $\nabla \varphi \left(r,0\right)\cdot d_1<0$ and $\nabla \varphi \left(r,\pi /2\right)\cdot d_2<0$.

2) There are real constants C₀ and C₁ such that, uniformly in θ we have:

a) $\underset{r\to +\infty }{\lim \sup }{r}^{1-\alpha }\nabla \varphi \left(r,\theta \right)\cdot d_0\le C_0<0$

b) $\underset{r\to +\infty }{\lim \sup }{r}^{2-\alpha }\left|D_{lj}\varphi \left(r,\theta \right)\right|\le C_1$ and (c) $\underset{r\to +\infty }{\lim \sup }\varphi \left(r,\theta \right)=+\infty$

3) $\underset{r\to +\infty }{\lim \sup }{r}^{1-\alpha }\nabla \varphi \left(r,0\right)\cdot d_1\le C_0<0$ and $\underset{r\to +\infty }{\lim \sup }{r}^{1-\alpha }\nabla \varphi \left(r,\pi /2\right)\cdot d_2\le C_0<0$

$D_{lj}\varphi \left(r\mathrm{,}\theta \right)$ denote the partial derivatives of $\varphi \left(r\mathrm{,}\theta \right)$ with respect to $x_l\left(l=1,2\right)$ and $x_j\left(j=1,2\right)$. r and $\theta$ are the polar coordinates of $x=\left(x_1\mathrm{,}{x}_2\right)$.

Proof: See Appendix.

Remark. Let $x=\left(x_1\mathrm{,}{x}_2\right)\in D_j\mathrm{,0}\le j\le 2$ and $A_j\left(x\right)=\left(\Delta Y_k/Y_k=x\right)$ ; we obtain

$d_j\left(x\right)=\mathbb{E}\left[A_j\left(x\right)\right]\text{and}{A}_j\left(x\right)\in \left\{\left(1,0\right),\left(-1,1\right),\left(-1,0\right),\left(0,-1\right)\right\}$

On $\left\{Y_k=\left(s\mathrm{,}i\right)\right\}$ we have: $\mathbb{P}\left[{\Vert \Delta Y_k\Vert }^2=2\right]=\frac{\beta si}{\tau \left(s,i\right)}=1-\mathbb{P}\left[{\Vert \Delta Y_k\Vert }^2=1\right]$ and

$\mathbb{E}\left[{\Vert A_0\left(s,i\right)\Vert }^2\right]=1+\frac{\beta si}{\tau \left(s,i\right)}<2,\mathbb{E}\left[{\Vert A_1\left(s,i\right)\Vert }^2\right]=1<2,\mathbb{E}\left[{\Vert A_2\left(s,i\right)\Vert }^2\right]=1<2$

An immediate consequence of the lemma 3 is:

Lemma 4. Let $x=\left(x_1,x_2\right),y=\left(y_1,y_2\right)$ two vectors of the plane and $\Vert x\Vert =\sqrt{x_1^2+x_2^2}$, $x\cdot y=x_1{y}_1+x_2{y}_2$.

Then, There are $\epsilon >0$ and $K>0$ such that

1) If $R_p>1$, then $\forall \Vert x\Vert \ge K$, $\mathbb{E}\left[\varphi \left(Y_{k+1}\right)-\varphi \left(Y_k\right)\mathrm{<mo>\; /\; </mo>}{Y}_k=x\right]\le 0$.

2) If $R_p\le 1$, then $\forall \Vert x\Vert \ge K$, $\mathbb{E}\left[\varphi \left(Y_{k+1}\right)-\varphi \left(Y_k\right)\mathrm{<mo>\; /\; </mo>}{Y}_k=x\right]\le -\epsilon$.

Proof: See Appendix.

Lemma 5. Let ${\left(Y_k\right)}_{k\in \mathbb{N}}$ be the embedded process of ${\left(X_t\right)}_{t\ge 0}$, which is a discrete Markov chain representing the sequence of values taken by ${\left(X_t\right)}_{t\ge 0}$ at transition times. Then

1) If $R_p\le 1$, then the Markov chain ${\left(Y_k\right)}_{k\in \mathbb{N}}$ is positive recurrent.

2) If $R_p>1$, then the Markov chain ${\left(Y_k\right)}_{k\in \mathbb{N}}$ is null recurrent.

Proof: See Appendix.

We can state now our main results:

Theorem 6. Let $T_0=\inf \left\{t\ge 0,I\left(t\right)=0\right\}$ with $\inf \varnothing =+\infty$. Then, for all $i\in {\mathbb{N}}^{\mathrm{<mo>\; *\; </mo>}}$, ${\mathbb{P}}_i\left[T_0<+\infty \right]=1$ and ${\lim }_{t\to +\infty }{\mathbb{P}}_i\left[I\left(t\right)=0\right]=1$.

Proof: This result is a consequence of the lemma 5 and the properties of recurrent Markov chains with nonempty absorbing set of states. (see [22], Proposition 5-15). It reflects the absorbent nature of the Markov chain. $\square$

Theorem 7. Let $T_0=\inf \left\{t\ge 0,I\left(t\right)=0\right\}$ with $\inf \varnothing =+\infty$ and $\left({\stackrel{\mathrm{<mo>\; *\; </mo>}}{S}}_0=\frac{\left(1-p\right)n}{\delta },{\stackrel{\mathrm{<mo>\; *\; </mo>}}{I}}_0=0,{\stackrel{\mathrm{<mo>\; *\; </mo>}}{R}}_0=\frac{np}{\delta }\right)$.

If $R_p\le 1$, then (1) $\mathbb{E}\left[T_0\right]=+\infty$ and (2) ${\lim }_{t\to +\infty }\left(\bar{S}\left(t\right)\mathrm{,}\bar{I}\left(t\right)\mathrm{,}\bar{R}\left(t\right)\right)=\left({\stackrel{\mathrm{<mo>\; *\; </mo>}}{S}}_0\mathrm{,}{\stackrel{\mathrm{<mo>\; *\; </mo>}}{I}}_0\mathrm{,}{\stackrel{\mathrm{<mo>\; *\; </mo>}}{R}}_0\right)$.

Proof: The first result reflects the positive recurrence obtained from the lemma 5. The second assertion follows from the fact that the Markov chain is absorbent, and once in the absorbing state, the correlation between $S\left(t\right)$ and $I\left(t\right)$ is identically zero. Therefore, asymptotically the deterministic equations and the mathematical expectation equations have the same equilibrium points. $\square$

Theorem 8. Let $T_0=\inf \left\{t\ge 0,I\left(t\right)=0\right\}$, $\inf \varnothing =+\infty$ and $\left({\stackrel{\mathrm{<mo>\; *\; </mo>}}{S}}_e=\frac{\delta +\mu +\gamma }{\beta },{\stackrel{\mathrm{<mo>\; *\; </mo>}}{I}}_e=\frac{\left(R_p-1\right)\delta }{\beta },{\stackrel{\mathrm{<mo>\; *\; </mo>}}{R}}_e=\frac{np\beta +\gamma \delta \left(R_p-1\right)}{\delta \beta }\right)$

If $R_p>1$, then (1) $\mathbb{E}\left[T_0\right]=+\infty$ and (2) ${\lim }_{t\to +\infty }\left(\bar{S}\left(t\right)\mathrm{,}\bar{I}\left(t\right)\mathrm{,}\bar{R}\left(t\right)\right)=\left({\stackrel{\mathrm{<mo>\; *\; </mo>}}{S}}_e\mathrm{,}{\stackrel{\mathrm{<mo>\; *\; </mo>}}{I}}_e\mathrm{,}{\stackrel{\mathrm{<mo>\; *\; </mo>}}{R}}_e\right)$

Proof: The first assertion is proved by observing that there are asymptotically two distinct equilibrium points, and necessarily $E\left(T\right)=+\infty$ in the case $R_p>1$, otherwise the two equilibrium points would be confused by uniqueness of the stationary measure.

The proof of the second assertion is similar to that of the second assertion of Theorem 7. $\square$

5. Simulation

In what follows, we will denote by $\bar{I}$ and dI numerical solutions of Equations (7) and (1) respectively. The average of the simulated realizations of the number of infected $I\left(t\right)$ is denoted by mI. We used MATLAB software for Monte-Carlo simulations and R software for graphics

Let an initial population of $S_0=100$ susceptibles with an initial number of $I_0=2$ infected for the following values of the parameters:

$\beta =0.69;\delta =0.25;\mu =0.02;\gamma =0.5;n=3.5;p=0.51;R_p=6.15$

In Figure 4, we have the estimate of the covariance $co{v}_{SI}\left(t\right)$ from 50 simulations. Figure 5 give a comparison of $\bar{I}$, dI and mI in the time interval $\left[\mathrm{0,26}\right]$. The time interval is then varied for the same values of the parameters. It appears that for the large values of t, we obtain $\bar{I}\left(t\right)\approx dI\left(t\right)\approx I_e^{\mathrm{<mo>\; *\; </mo>}}$, where $I_e^{\mathrm{<mo>\; *\; </mo>}}$ is the endemic equilibrium of the Equation (1), the expected asymptotic value when $R_p>1$.

For the considered values of the parameters, the endemic equilibrium value is $I_e^{*}=1.8659$. The simulations gave the following values :

$\bar{I}\left(26\right)\approx \mathrm{1.8648,}dI\left(26\right)\approx \mathrm{1.8650,}mI\left(26\right)\approx 0$ (voir Figure 5)

$\bar{I}\left(52\right)\approx \mathrm{1.8650,}dI\left(52\right)\approx \mathrm{1.8650,}mI\left(52\right)\approx 0$ (voir Figure 6)

In Figure 3, two sample paths of $I\left(t\right)$, their mean and the deterministic solution for the following values:

$\begin{array}{l}{S}_0=100;I_0=2;\beta =0.69;\delta =0.25;\mu =0.02;\\ \gamma =0.5;n=3.5;p=0.51;t\in \left[0,26\right]\text{et}{R}_p=6.1\end{array}$

Figure 4 estimated covariance function from 50 sample paths of $X_t=\left(S\left(t\right)\mathrm{,}I\left(t\right)\right)$ and in Figure 5 Deterministic solution (dI), solution of mathematical expectations ( $\bar{I}$ ) and mean of 50 sample paths of $I\left(t\right)$ (mI) for the following values:

$\begin{array}{l}{S}_0=100;I_0=2;\beta =0.69;\delta =0.25;\mu =0.02;\\ \gamma =0.5;n=3.5;p=0.51;t\in \left[0,26\right];R_p=6.15\end{array}$

Figure 3. Two sample paths of I(t), their mean—the average of the simulated values calculated at each instant, estimate of the mathematical expectation of I(t)—and the deterministic solution.

Figure 4. Estimated covariance function from 50 sample paths of $X_t=\left(S\left(t\right)\mathrm{,}I\left(t\right)\right)$.

Figure 5. Deterministic solution (dI), solution of mathematical expectations ( $\bar{I}$ ) and mean of 50 sample paths of I(t) (mI).

In Figure 6 mean of 50 sample paths of $I\left(t\right)$ and solution of mathematical expectations ( $\bar{I}$ ) for previous parameter values are compared to the deterministic solution (dI) over the time interval $\left[\mathrm{0,52}\right]$.

Figure 6. Deterministic solution (dI), solution of mathematical expectations ( $\bar{I}$ ) and mean of 50 sample paths of I(t) (mI) for previous values but $t\in \left[0,52\right]$.

6. Discussions

This paper presents a stochastic compartmental model SVIR of measles. A comparison of our stochastic model with the corresponding deterministic model indicates that the deterministic solution is asymptotically the mean of the stochastic solution. It is well known that mI obtained by random sampling (Monte Carlo methods) before extinction is an estimate of $\bar{I}$. Our result shows that the three trajectories of $\bar{I}$, dI and mI asymptotically coincide. The deterministic solution is the mean of the stochastic solution.

In addition, unlike the deterministic approach, we show that the epidemic is extinguished independently of the threshold R_p with a probability equal to 1. More precisely, if $R_p\le 1$ extinction occurs in a time of finite mean, and if $R_p>1$ the disease eventually disappears in a time of infinite mean.

One of the peculiarities of our model is that the size of the population is not constant and can be quite large. The extinction of the process in this case is not guaranteed unlike in the case where the size of the population is constant. This led us to focus on the probability of absorption of the process.

On the other hand, when $R_0>1$, it is well known for the constant population SIR model [23] that the average duration of the epidemic increases exponentially with the size of the population. This fact is confirmed by the assertion 1. of the theorem 8 and the extinction is done in a time of infinite mean when $R_p>1$.

7. Conclusions

To understand the dynamics of the system before absorption, a commonly used measure is the quasi-stationary distribution [24]. The term quasi-stationarity refers to the distribution of the Markov chain by conditioning on the event that absorption has not occurred yet. It gives a good measure of the behavior before absorption when the absorption time is very long.

If the set of transient states is finite and irreducible, it is well known that the quasi-stationary distribution exists [25]. But if this set is infinite the existence of a quasi-stationary distribution is not guaranteed, and even if it does exist it is practically impossible to determine it explicitly. To elucidate this situation, an extension of our work would be the study of the process in quasi-stationary regime.

The emergence of epidemics often reveals complex dynamic relationships between susceptible individuals, pathogens and their environments. Complex dynamic relationships that result in seasonal epidemic cycles vary over time [26]. In Niger, recent studies [7] reveal two main periodicity of measles, a more accentuated annual periodicity, probably due to seasonal agricultural labour migration and a low and unstable periodicity of 2 to 3 years which is partly explained by heterogeneity in vaccination coverage. To account for this aspect of temporal and environmental variability, it would be necessary to extend our study to the analysis of the time series of cases of infection. The stochastic aspect takes better account of these temporal and environmental fluctuations and may provide a framework to improve our understanding of the complex dynamics of measles epidemics.

Appendix

Proof of the lemma 2:

The lemma is a consequence of the definition of R_p and of the expressions $d_0\mathrm{,}{d}_1\mathrm{,}{d}_2$ :

$\begin{array}{l}{d}_0=\left(\frac{n\left(1-p\right)-\beta si-\delta s}{n\left(1-p\right)+\beta si+\delta s+\left(\mu +\gamma +\delta \right)i},\frac{\beta si-\left(\mu +\gamma +\delta \right)i}{n\left(1-p\right)+\beta si+\delta s+\left(\mu +\gamma +\delta \right)i}\right)\\ d_1=\left(\frac{n\left(1-p\right)-\delta s}{n\left(1-p\right)+\delta s},0\right)\\ d_2=\left(\frac{n\left(1-p\right)-\beta s_{0}i-\delta s_0}{n\left(1-p\right)+\beta s_{0}i+\delta s_0+\left(\mu +\gamma +\delta \right)i},\frac{\beta s_{0}i-\left(\mu +\gamma +\delta \right)i}{n\left(1-p\right)+\beta s_{0}i+\delta s_0+\left(\mu +\gamma +\delta \right)i}\right)\end{array}$ (10)

We can easily determine the signs of the abscissas and ordinates of $d_j$, indeed:

$s_0=\max \left(\frac{n\left(1-p\right)}{\delta },\frac{\mu +\gamma +\delta }{\beta }\right);R_p=\frac{\beta n\left(1-p\right)}{\delta \left(\mu +\gamma +\delta \right)}$

1) $d_{0x}<0<d_{0y}$ ; $d_{1x}=0$ ; $d_{2x}<0\le d_{2y}$

2) If $R_p\le 1$ we have: $0<d_{0y}<-d_{0x}$ ; $d_{1x}=0$ ; $d_{2x}<0=d_{2y}$

3) If $R_p>1$ we have: $0<d_{0y}<-d_{0x}$ ; $d_{1x}=0$ ; $d_{2x}<0<d_{2y}$ $\square$

Proof of the lemma 3:

In polar coordinates, we have $\nabla \varphi \left(r\mathrm{,}\theta \right)=\alpha r^{\alpha -1}\left(\cos \left(\alpha \theta -{\theta }_1\right)\mathrm{,}-\sin \left(\alpha \theta -{\theta }_1\right)\right)$. To establish the result, we distinguish the two cases $R_p\le 1$ and $R_p>1$.

· If $R_p\le 1$, the angle between $d_0$ and $\nabla \varphi \left(r\mathrm{,}\theta \right)$ is

$a\left(\theta \right)={\theta }_1-\alpha \theta -\left(\psi +\frac{\pi }{2}-\theta \right)$. From $-{\theta }_2-\psi -\frac{\pi }{2}\le a\left(\theta \right)\le \left({\theta }_1-\psi \right)-\frac{\pi }{2}$, the

angle between $d_1$ and $\nabla \varphi \left(r\mathrm{,0}\right)=\alpha r^{\alpha -1}\left(\cos {\theta }_1\mathrm{,}\sin {\theta }_1\right)$ is equal to $a_1={\theta }_1-\pi$. furthermore, we show that the angle between $d_2$ and $\nabla \varphi \left(r,\pi /2\right)=\alpha r^{\alpha -1}\left(\cos {\theta }_2\mathrm{,}-\sin {\theta }_2\right)$ is equal to $a_2=-{\theta }_2-\pi$.

The choice of ${\theta }_1$ and ${\theta }_2$ allows to have:

$-\frac{5\pi }{4}<{\theta }_1-\pi <-\frac{3\pi }{4},-\frac{3\pi }{2}<-{\theta }_2-\pi <{\theta }_1-\frac{3\pi }{2}<-\frac{\pi }{2}$

and for any $\theta ,0<\theta <\frac{\pi }{2}$, we obtain $-\frac{5\pi }{4}<a\left(\theta \right)<{\theta }_1-\frac{3\pi }{4}<-\frac{\pi }{2}$. As a result, we have inequalities $\cos a_1<\frac{-\sqrt{2}}{2}$, $\cos a_2<\cos \left({\theta }_1-\frac{3\pi }{2}\right)<0$ et $\cos \left(a\left(\theta \right)\right)<\cos \left({\theta }_1-\frac{3\pi }{4}\right)<0$.

· If $R_p>1$, we have good ${\theta }_1-\psi -\frac{\pi }{2}\le a\left(\theta \right)\le {\theta }_1-\psi$ et $a_1={\theta }_1-\pi$.

For $d_2$, we find $a_2=-{\theta }_2-{\psi }_2$. The choice of ${\theta }_1$ and ${\theta }_2$ leads to

$-\frac{3\pi }{2}<{\theta }_1-\pi <-\pi \text{et}-\frac{3\pi }{2}<-{\theta }_2-{\psi }_2<-{\theta }_2-\frac{\pi }{2}<-\frac{\pi }{2}$

and for any $\theta ,0<\theta <\frac{\pi }{2}$ et $-\frac{3\pi }{2}<{\theta }_1-\psi -\frac{\pi }{2}\le a\left(\theta \right)\le {\theta }_1-\psi <-\frac{\pi }{2}$. Thus

$\cos a_1<\frac{-\sqrt{2}}{2},\cos a_2<\cos \left(-{\theta }_2-\frac{\pi }{2}\right)<0\text{et}\cos \left(a\left(\theta \right)\right)<\cos \left({\theta }_1-\psi \right)<0.$

definitively, for any value of $R_p$ and for any $0\le \theta \le \pi /2$, we deduce that

$\nabla \varphi \left(r,\theta \right)\cdot d_0<0,\nabla \varphi \left(r,0\right)\cdot d_1<0\text{and}\nabla \varphi \left(r,\pi /2\right)\cdot d_2<0.$ (11)

So the assertions 1; 2. a) et 3. deduce.

To establish the assertion 2. (b), we consider the partial derivatives with respect to $x_1$ and $x_2$ of $\varphi$ :

$\begin{array}{l}{D}_1\varphi =\cos \theta {\varphi }_r-\frac{\sin \theta }{r}{\varphi }_{\theta }=\alpha r^{\alpha -1}\cos \left(\left(\alpha -1\right)\theta -{\theta }_1\right)\\ D_2\varphi =\sin \theta {\varphi }_r+\frac{\cos \theta }{r}{\varphi }_{\theta }=\alpha r^{\alpha -1}\sin \left(\left(\alpha -1\right)\theta -{\theta }_1\right)\end{array}$ (12)

where ${\varphi }_r$ and ${\varphi }_{\theta }$ are the partial derivatives with respect to r and $\theta$ of $\varphi$ of Jacobian matrix of $\varphi$ :

$D_{lj}\varphi =\alpha \left(\alpha -1\right)r^{\alpha -2}\left(\begin{array}{cc}\cos \left(\left(\alpha -2\right)\theta -{\theta }_1\right)& \sin \left(\left(\alpha -2\right)\theta -{\theta }_1\right)\\ -\sin \left(\left(\alpha -2\right)\theta -{\theta }_1\right)& \cos \left(\left(\alpha -2\right)\theta -{\theta }_1\right)\end{array}\right)$ (13)

The assertion 2. (c) follows from the definition of $\varphi$ and the fact that $\cos \left(\alpha \theta -{\theta }_1\right)>0$ ; indeed, for any $\theta \mathrm{,0}\le \theta \le \frac{\pi }{2}$ we have $-\frac{\pi }{2}<\alpha \theta -{\theta }_1<\frac{\pi }{2}$.

Hence the lemma 3. $\square$

Proof of the lemma 4:

The proof is analogous to that of the theorem 3 of [19]. The Taylor formula of the function $\varphi$ is:

$\varphi \left(x+h\right)-\varphi \left(x\right)=\nabla \varphi \left(x\right)\cdot h+R\left(x\mathrm{,}h\right)$

where $h=\left(h_1\mathrm{,}{h}_2\right)$ and $R\left(x,h\right)=\frac{1}{2}{\displaystyle {\sum }_{l,j=1,2}}{D}_{lj}\varphi \left(x+\eta h\right)h_l{h}_j$ is the remainder of Taylor with $0<\eta <1$.

For $l\in \left\{\mathrm{0,1,2}\right\}$, when we replace h by $A_l\left(x\right)$, we get:

$\mathbb{E}\left[\varphi \left(Y_{k+1}\right)-\varphi \left(Y_k\right)/Y_k=x\in D_l\right]=\nabla \varphi \left(x\right)\cdot d_l+\mathbb{E}\left[R\left(x,A_l\left(x\right)\right)\right]$

Applying the lemma 3 and the remark 4, we have

$\mathbb{E}\left[\varphi \left(Y_{k+1}\right)-\varphi \left(Y_k\right)/Y_k=x\in D_l\right]=\nabla \varphi \left(x\right)\cdot d_l+O\left({\Vert x\Vert }^{\alpha -2}\right)$

· If $R_p\le 1$, we have $1<\alpha <2$ and $\lim {\sup }_{\Vert x\Vert \to \infty }{\Vert x\Vert }^{-\alpha +1}\nabla \varphi \left(x\right)\cdot d_l\le C_0<0$ ; Therefore:

$\exists \epsilon >0\text{and}K>0,\mathbb{E}\left[\varphi \left(Y_{k+1}\right)-\varphi \left(Y_k\right)/Y_k=x\right]\le -\epsilon \forall \Vert x\Vert \ge K.$

· If $R_p>1$, it turns out that $0<\alpha <1$, we cannot conclude that

$\mathbb{E}\left[\varphi \left(Y_{k+1}\right)-\varphi \left(Y_k\right)/Y_k=x\right]\le 0\forall \Vert x\Vert \ge K.$

What completes the demonstration. $\square$

Proof of the lemma 5:

Let us show the recurrence in the case $R_p>1$.

We pose $B=\left\{x\mathrm{<mo>\; /\; </mo>}\Vert x\Vert \le K\right\}$, $T=\inf \left\{k\ge \mathrm{0,}{Y}_k\in B\right\}$ and $Z_k=\varphi \left(Y_k\right){\mathbb{I}}_{\left\{T>k\right\}}$ where ${\mathbb{I}}_A$ denotes the indicator map of A.

Let ${\left({\mathcal{F}}_k\right)}_{k\in \mathbb{N}}$ be the filtration associated to ${\left(Y_k\right)}_{k\in \mathbb{N}}$. Knowing that ${\mathbb{I}}_{\left\{T>k+1\right\}}\le {\mathbb{I}}_{\left\{T>k\right\}}$, we can write: $\mathbb{E}\left[Z_{k+1}/{\mathcal{F}}_k\right]\le \mathbb{E}\left[\varphi \left(Y_{k+1}\right){\mathbb{I}}_{\left\{T>k\right\}}/{\mathcal{F}}_k\right]$

$\mathbb{E}\left[\varphi \left(Y_{k+1}\right){\mathbb{I}}_{\left\{T>k\right\}}/{\mathcal{F}}_k\right]={\mathbb{I}}_{\left\{T>k\right\}}\mathbb{E}\left[\varphi \left(Y_{k+1}\right)/{\mathcal{F}}_k\right]\le {\mathbb{I}}_{\left\{T>k\right\}}\varphi \left(Y_k\right)=Z_k.$

In this last expression, the last inequality is obtained from the second assertion of the lemma 4. Thereafter ${\left(Z_k\right)}_{k\in \mathbb{N}}$ is a positive supermartingale and therefore $\mathbb{P}\left[{\lim }_{k\to +\infty }{Z}_k=0\right]=1$.

On the other hand, because the Markov chain ${\left(Y_k\right)}_{k\in \mathbb{N}}$ is irreducible, we have $\mathbb{P}\left[\lim {\sup }_{k\to +\infty }\Vert Y_k\Vert =\infty \right]=1$. In this case, on $\left\{T=+\infty \right\}$, it follows that

${\lim }_{k\to +\infty }{Z}_k={\lim }_{k\to +\infty }\varphi \left(Y_k\right)=+\infty$, thus $\mathbb{P}\left[T=+\infty \right]=0$. In other words, the finite set A is visited an infinite number of times by the Markov chain ${\left(Y_k\right)}_{k\in \mathbb{N}}$, which corresponds to recurrence. Finally, the last assertion of the lemma is a consequence of the first assertion of the lemma 4 and of Foster’s positive recurrence criterion [27]. $\square$

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	World Health Organization (2019) Fact about Measles.
[2]	Kitengeso, R.E., et al. (2015) A Mathematical Model for Control and Elimination of the Transmission Dynamics of Measles. Applied and Computational Mathematics, 4, 396-408. https://doi.org/10.11648/j.acm.20150406.12
[3]	Modèles compartimentaux en épidémiologie—Wikipédia, Oct 2019.
[4]	Alkassoum, S., et al. (2016) Surveillance épidémiologique de la rougeole au niger: Analyse de la base de données des maladies à déclaration obligatoire (mdo) de 2003 à 2015. International Journal of Innovation and Scientific Research, 17, 264-274.
[5]	Mitku, S.N. and Koya, P.R. (2017) Mathematical Modeling and Simulation Study for the Control and Transmission Dynamics of Measles. American Journal of Applied Mathematics, 5, 99-107. https://doi.org/10.11648/j.ajam.20170504.11
[6]	Moussa Tessa, O. (2006) Mathematical Model for Control of Measles by Vaccination. Proceedings of MSAS 06 Conference, Bamako, August 2006, 31-36.
[7]	Blake, A., DJibo, A., Guindo, O. and Bharti, N. (2020) Investigating Persistent Measles Dynamics in Niger and Associations with Rainfall. Journal of the Royal Society Interface, 17, 20200480. https://doi.org/10.1098/rsif.2020.0480
[8]	Moss, W.J. (2017) Measles. The Lancet, 390, 2490-2502. https://doi.org/10.1016/S0140-6736(17)31463-0
[9]	Moss, W.J. (2007) Measles Still Has a Devastating Impact in Unvaccinated Populations. PLOS Medicine, 4, e24. https://doi.org/10.1371/journal.pmed.0040024
[10]	Hethcote, H.W. (2000) The Mathematics of Infectious Diseases. SIAM Review, 42, 599-653. https://doi.org/10.1137/S0036144500371907
[11]	van der Driessche, P. and Watmough, J. (2002) Reproduction Numbers and Sub-Threshold Endemic Equilibria for Compartmental Models of Disease Transmission. Mathematical Biosciences, 180, 29-48. https://doi.org/10.1016/S0025-5564(02)00108-6
[12]	Finkenstädt, B.F. (2002) A Stochastic Model for Extinction and Recurrence of Epidemics Estimation and Inference for Measles Outbreaks. Biostatistics, 3, 493-510. https://doi.org/10.1093/biostatistics/3.4.493
[13]	Anderson, R.M. and Rhodes, C.J. (1996) Power Laws Governing Epidemics in Isolated Populations. Nature, 381, 600-602. https://doi.org/10.1038/381600a0
[14]	Anderson, R.M., Jensen, H.J. and Rhodes, C.J. (1997) On the Critical Behavior of Simple Epidemics. Proceedings of the Royal Society of London B, 264, 1639-1646. https://doi.org/10.1098/rspb.1997.0228
[15]	Cliff, A.D., Haggett, P. and Smallman-Raynor, M. (1993) Mesales: An Historical Geography of a Major Human Viral Disease from Global Expansion to Local Retreat. Blackwell, Oxford, 1840-1990.
[16]	Kitengeso, R.E. (2016) Stochastic Modelling of the Transmission Dynamics of Measles with Vaccination Control. International Journal of Theoretical and Applied Mathematics, 2, 60-73.
[17]	Graham, M., et al. (2019) Measles and the Canonical Path to Elimination. Science, 364, 584-587. https://doi.org/10.1126/science.aau6299
[18]	Lahrouz, A., Omari, L. and Kiouach, D. (2011) Global Analysis of a Deterministic and Stochastic Nonlinear SIRS Epidemic Model. Nonlinear Analysis: Modelling and Control, 16, 59-76. https://doi.org/10.15388/NA.16.1.14115
[19]	Rosenkrantz, W.A. (1989) Ergodicity Conditions for Two-Dimensional Markov Chains on the Positive Quadrant. Probability Theory and Related Fields, 83, 309-319. https://doi.org/10.1007/BF00964367
[20]	Meyn, S.P. and Tweddie, R.L. (1992) Stability of Markovian Processes I: Criteria for Discrete-Time Chains. Advances in Applied Probability, 24, 542-574. https://doi.org/10.1017/S000186780002440X
[21]	Meyn, S.P. and Tweddie, R.L. (1993) Markov Chains and Stochastic Stability. Springer-Verlag, London. https://doi.org/10.1007/978-1-4471-3267-7
[22]	Kemeny, J.G., Laurie Snell, J. and Knapp, A.W. (1982) Denumerable Markov Chains. Springer-Verlag, New York.
[23]	Allen, J.S. and Burgin, M. (2000) Comparison of Deterministic and Stochastic SIS and SIR Models in Discrete Time. Mathematical Biosciences, 163, 1-33. https://doi.org/10.1016/S0025-5564(99)00047-4
[24]	Darroch, J.N. and Seneta, E. (1967) On Quasi-Stationary Distributions in Absorbing Continuous-Time Finite Markov Chains. Journal of Applied Probability, 4, 192-196. https://doi.org/10.1017/S0021900200025341
[25]	Artalejo, J.R. and Lopez-Henero, M.J. (2010) Quasi-Stationary and Ratio of Expectations Distributions: A Comparative Study. Journal of Theoretical Biology, 266, 264-274. https://doi.org/10.1016/j.jtbi.2010.06.030
[26]	Bartlett, M.S. (1992) Measles Periodicity and Community Size. Journal of the Royal Statistical Society: Series A, 120, 48-70. https://doi.org/10.2307/2342553
[27]	Foster, F.G. (1982) On the Stochastic Matrices Associated with Certain Queueing Processes. Annals of Mathematical Statistics, 24, 295-308.