Approximations of Quasi-Stationary Distributions of the Stochastic SVIR Model for the Measles

Moussa Seydou, Moussa Tessa

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

DOI: 10.4236/jamp.2021.99145                

Abstract

In this paper, we analyze the quasi-stationary distribution of the stochastic SVIR (Susceptible, Vaccinated, Infected, Recovered) model for the measles. The quasi-stationary distributions, as discussed by Danoch and Seneta, have been used in biology to describe the steady state behaviour of population models which exhibit discernible stationarity before to become extinct. The stochastic SVIR model is a stochastic SIR (Susceptible, Infected, Recovered) model with vaccination and recruitment where the disease-free equilibrium is reached, regardless of the magnitude of the basic reproduction number. But the mean time until the absorption (the disease-free) can be very long. If we assume the effective reproduction number R_p < 1 or , the quasi-stationary distribution can be closely approximated by geometric distribution. β and δ stands respectively, for the disease transmission coefficient and the natural rate.

Keywords

Compartment Models, SIR, Markov Chains, Stochastic Simulation, Basic Reproduction Number, Quasi-Stationary Distribution, Measles

1. Introduction

Measles is a highly contagious viral infection that manifests as a rash associated with signs of respiratory infections. It is caused by a virus of the paramyxovirus family whose reservoir is exclusively human [1] and is transmitted by direct contact with secretions from the nose, throat and through the air [2] [3]. The virus primarily infects the respiratory tract. Upon infection, the patient passes through a latent period of 6 to 9 days, followed by 6 to 7 day infective period [4]. 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, hence its essential role in any measles control initiative.

The children under 5 years remain the most affected. 90% who die have less than 5 years. In developing countries, like Niger where children under one year old represent 4.32%, those under 5 years old 19.73% and those under 15 years 51.18%. The measles remains one of the main causes of infant mortality [2] [5] [6] [7].

Our stochastic model is a stochastic SVIR (Susceptible, Vaccinated, Infected, Recovered) model for the measles [8], where the process $X_t={\left(S_t\mathrm{,}{I}_t\right)}_{t\ge 0}$ is a continuous-time Markov chain resulting from a set of transient states $E_0$ which evolves until it escapes to a set of absorbing states corresponding to disease-free equilibrium. $S\left(t\right)\mathrm{,}I\left(t\right)$ denote respectively the number of susceptible and infected. When the process reaches the set of absorbing states, it remains there permanently. However, before the instant of absorption (which is relatively long), the process passes through a quasi-stationary state.

To understand this phenomenon, we study the long time behavior of the process conditioned on non extinction, which leads us to consider the quasi-stationary distribution introduced by Danoch and Seneta in biology. It allows to describe the steady state behaviour of population models which exhibit discernible stationarity before to become extinct [9]. The term quasi-stationarity refers to the distribution of the Markov chain by conditioning on the event that absorption has not occurred yet [10]. It gives a good measure of the behavior before absorption when the absorption time is very long. But this measure has a number of flaws. Indeed if the set of transient states is finite and irreducible, it is well known that the quasi-stationary distribution exists [11]. 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.

For the continuous time SVIR model under some conditions on the effective reproduction number [12], the quasi-stationary distribution of the number of infected exists and can be closely approximated by geometric distribution.

The main results are stated in theorem 4 and theorem 5. Precisely, let $R_p$ be the effective reproduction number. In theorem 4, we prove that, if $R_p<1$, the quasi-stationary distribution of the number of infected can be closely approximated by geometric distribution with parameter $1-R_p$. However in theorem 5, we note that if $R_p>1+\frac{\beta }{\delta }$ this latter distribution is approximatively geometric with parameter $1/I_e^{\mathrm{<mo>\; *\; </mo>}}$, where $\beta$, $\delta$ and $I_e^{\mathrm{<mo>\; *\; </mo>}}$ stands respectively, for the disease transmission coefficient, the natural death rate and the endemic equilibrium point, for the number of infected.

The rest of the paper is organized as follows: The Section 2 describes the deterministic SVIR model by a system of differential equations. The equilibrium points of the system of differential equations are also given. In Section 3, we use the continuous time Markov chains model to form our stochastic SVIR model [8]. Section 4 is devoted to the study of the quasi-stationary distribution of the stochastic SVIR 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 susceptible (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. p is the probability that a newborn will acquire immunity after being vaccinated. In addition, we assume that:

● the natural death rate is $\delta$ for each compartment.

● infectious patients recover at the rate of $\gamma$.

● infectious patients have an additional $\mu$ death rate from measles.

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

The compartment diagram of the transitions in the SVIR model is in Figure 1.

The dynamics of a well-mixed population can be described by the following system of 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.

Figure 1. Compartment diagram of model SVIR.

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

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

$R_0=\frac{\beta n}{\delta \left(\delta +\mu +\gamma \right)}$ is the basic reproduction number [13] and the effective

reproduction number is defined by $R_p=\left(1-p\right)R_0$ [12]. Recall that $R_0$ is defined as average number of secondary infections produced when one infected individual is introduced into a host population where everyone is susceptible [14] [15]. If $R_p<1$ this equilibrium point is asymptotically stable [13]. In addition, we have $R_p<R_0$ and $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 equilibrium point appears ( ${\stackrel{\mathrm{<mo>\; *\; </mo>}}{S}}_e\mathrm{,}{\stackrel{\mathrm{<mo>\; *\; </mo>}}{I}}_e\mathrm{,}{\stackrel{\mathrm{<mo>\; *\; </mo>}}{R}}_e$ ) asymptotically stable [13], where

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

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\mathrm{,}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\mathrm{,}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\mathrm{,}l\right)\mathrm{<mo>\; /\; </mo>}{X}_t=\left(s\mathrm{,}i\right)\right\}$

Figure 2. States transition.

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)

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 Equations (5), we obtain the Kolmogorov Forward equations, for all $s\ge 0$ and $i\ge 0$

$\begin{array}{c}\frac{\text{d}{P}_{s\mathrm{,}i}}{\text{d}t}=n\left(1-p\right)\left[P_{s-\mathrm{1,}i}-P_{s\mathrm{,}i}\right]+\beta \left[\left(s+1\right)\left(i-1\right)P_{s+\mathrm{1,}i-1}-si{P}_{s\mathrm{,}i}\right]\\ +\left(\mu +\gamma +\delta \right)\left[\left(i+1\right)P_{s\mathrm{,}i+1}-i{P}_{s\mathrm{,}i}\right]+\delta \left[\left(s+1\right)P_{s+\mathrm{1,}i}-s{P}_{s\mathrm{,}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,i}\left(t\right),\bar{I}\left(t\right)={\displaystyle \underset{s=0}{\overset{+\infty }{\sum }}}{\displaystyle \underset{i=0}{\overset{+\infty }{\sum }}}i{P}_{s,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,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}$

In the previous article, see [8], we established the following result:

Theorem 1. 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$.

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

If $R_p\le 1$, then (1) $\mathbb{E}\left[T_0\right]<+\infty$ and (2) $\underset{t\to +\infty }{\lim }\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)$.

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

If $R_p>1$, then (1) $\mathbb{E}\left[T_0\right]=+\infty$ and (2) $\underset{t\to +\infty }{\lim }\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)$

Unlike the deterministic approach, we note 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. However, before the instant of absorption (which is relatively long) the process passes through a quasi-stationary state. To understand this phenomenon, we study the long time behavior of the process conditioned on non extinction, which leads us to consider the quasi-stationary distribution introduced by Danoch and Seneta in biology.

4. Quasi-Stationary Distribution

The term quasi-stationarity refers to the distribution of the Markov chain by conditioning on the event that absorption has not occurred yet [10]. It gives a good measure of the behavior before absorption when the absorption time is very long. In all the following ${\mathbb{P}}_x$ refer to the probability measure conditional on $X_0=x$ and ${\mathbb{P}}_{\alpha }={\displaystyle {\sum }_x}\alpha \left(x\right){\mathbb{P}}_x$ for any probability measure $\alpha$. ${\mathbb{E}}_x$ and ${\mathbb{E}}_{\alpha }$ are the corresponding expectations

Definition 4.1. A probability distribution $\pi$ on the set of transient states $E_T$ is called a quasi-stationary distribution for the process ${\left(X_t\right)}_{t\ge 0}$ if for all $t\ge 0$ and any measurable set $A\subset E$ we have $\pi \left(A\right)={\mathbb{P}}_{\pi }\left(X_t\in A/T_0>t\right)$ [10].

Remark

1) Equivalently, $\pi$ is the unique limiting conditional probability distribution such that $\forall t\ge 0$, $\forall A\subset E_T$ ${\lim }_{t\to +\infty }{\mathbb{P}}_{\alpha }\left\{X_t\in A/T>t\right\}=\pi \left(A\right)$ independently of initial distribution $\alpha$ [10].

2) If the set of transient states is finite and irreducible, it is well known that the quasi-stationary distribution exist. But if this set is infinite the existence of quasi-stationary distribution is not guaranteed, furthermore even if it exist, it is typically impossible to evaluate it explicitly. One is therefore lead to consider iterative methods [16] or asymptotic solutions by diffusion processes [17] [18] for the quasi-stationary distribution.

Theorem 4. Let $Q_{s\mathrm{,}i}^{\mathrm{<mo>\; *\; </mo>}}$ be the quasi-stationary distribution of the process ${\left(X_t\right)}_{t\ge 1}$ and $Q_{.i}^{*}={\displaystyle \underset{s\ge 0}{\sum }{Q}_{s,i}^{*}}$ the marginal distribution of the number $I^{\mathrm{<mo>\; *\; </mo>}}$ of infected in a quasi-stationary regime. If $R_p<1$, for all $i\ge 1$, $Q_{\mathrm{.}i}^{\mathrm{<mo>\; *\; </mo>}}\approx \left(1-R_p\right)R_p^{i-1}$

Proof. For all $i\mathrm{,}j\ge 1$, setting

$\{\begin{array}{l}{P}_s\left(i,j,\Delta t\right)=\mathbb{P}\left(I\left(t+\Delta t\right)=j/S\left(t\right)=s,I\left(t\right)=i\right)\\ P_I\left(i,j,t,\Delta t\right)=\mathbb{P}\left(I\left(t+\Delta t\right)=j/I\left(t\right)=i\right)\end{array}$ (8)

we have

$P_I\left(i\mathrm{,}j\mathrm{,}t\mathrm{,}\Delta t\right)={\displaystyle \underset{s\ge 0}{\sum }}\mathbb{P}\left(S\left(t\right)=s\right)P_s\left(i\mathrm{,}j\mathrm{,}\Delta t\right)$ (9)

and according to the process definition ${\left(X_t\right)}_{t>0}$,

$P_s\left(i,j,\Delta t\right)=\{\begin{array}{ll}\beta is\Delta t+o\left(\Delta t\right)\hfill & \text{if}j=i+1\hfill \\ \left(\mu +\gamma +\delta \right)i\Delta t+o\left(\Delta t\right)\hfill & \text{if}j=i-1\hfill \\ 1-\left[\beta is+\left(\mu +\gamma +\delta \right)i\right]\Delta t+o\left(\Delta t\right)\hfill & \text{if}j=i\hfill \end{array}$ (10)

we deduce that

$P_I\left(i,j,t,\Delta t\right)=\{\begin{array}{ll}\left(\mu +\gamma +\delta \right)i\Delta t+o\left(\Delta t\right)\hfill & \text{if}j=i+1\hfill \\ 1-\left[\beta i\bar{S}\left(t\right)+\left(\mu +\gamma +\delta \right)i\right]\Delta t+o\left(\Delta t\right)\hfill & \text{if}j=i-1\hfill \\ 1-\left[\beta i\bar{S}\left(t\right)+\left(\mu +\gamma +\delta \right)i\right]\Delta t+o\left(\Delta t\right)\hfill & \text{if}j=i\hfill \end{array}$ (11)

As in the case of disease-free equilibrium, $\underset{t\to +\infty }{\lim }\bar{S}\left(t\right)={\stackrel{*}{S}}_0=\frac{n\left(1-p\right)}{\delta }$, thanks to the Equation (11), we have

$P_I\left(i,j,t,\Delta t\right)=\{\begin{array}{ll}\beta i{\stackrel{*}{S}}_0\Delta t+o\left(\Delta t\right)\hfill & \text{if}j=i+1\hfill \\ \left(\mu +\gamma +\delta \right)i\Delta t+o\left(\Delta t\right)\hfill & \text{if}j=i-1\hfill \\ 1-\left[\beta i{\stackrel{*}{S}}_0+\left(\mu +\gamma +\delta \right)i\right]\Delta t+o\left(\Delta t\right)\hfill & \text{if}j=i\hfill \end{array}$ (12)

thus asymptotically the process $I\left(t\right)$ is a linear birth-death process with infinitesimal generator:

$q_{ij}=\{\begin{array}{l}\lambda i\text{if}j=i+1\\ \nu i\text{if}j=i-1\end{array}\text{where}\lambda =\frac{\beta n\left(1-p\right)}{\delta }\text{and}\nu =\left(\mu +\gamma +\delta \right)$ (13)

In this case, under the condition $\lambda <\nu$, it is well known [16] there is a unique quasi-stationary distribution for the process which follows the geometric law with parameter $1-\frac{\lambda }{\nu }=1-R_p$. Hence if $R_p<1$ for all $i\ge 1$, we obtain:

$Q_{\mathrm{.}i}^{\mathrm{<mo>\; *\; </mo>}}\approx \left(1-R_p\right)R_p^{i-1}$ or $Q_{\mathrm{.}i}^{\mathrm{<mo>\; *\; </mo>}}={\displaystyle \underset{s\ge 0}{\sum }}{Q}_{s\mathrm{,}i}^{\mathrm{<mo>\; *\; </mo>}}$ (14)

The proof is completed for the theorem 4. □

Remark. Under the condition $R_p\le 1$, the irreducible Markov chain ${\left(X_t\right)}_{t>0}$ is positive recurrent. Then an unique invariant probability measure $\pi$ exists and

$\pi \left(s,i\right)=\frac{1}{{\mathbb{E}}_{\left(s,i\right)}\left({\tau }_{s,i}\right)}\text{where}{\tau }_{s,i}=\inf \left\{t>0/X_t=\left(s,i\right)\right\}.$

Thus, the theorem 4 simply states that for all $i\ge 1$, ${\pi }_i={\displaystyle {\sum }_{s\ge 0}}\pi \left(s,i\right)\approx \left(1-R_p\right)R_p^{i-1}$ if $R_p<1$.

Theorem 5. Let $Q_{s\mathrm{,}i}^{\mathrm{<mo>\; *\; </mo>}}$ be the quasi-stationary distribution of the process ${\left(X_t\right)}_{t\ge 1}$ and $Q_{\mathrm{.}i}^{\mathrm{<mo>\; *\; </mo>}}={\displaystyle \underset{s\ge 0}{\sum }{Q}_{s\mathrm{,}i}^{\mathrm{<mo>\; *\; </mo>}}}$ the marginal distribution of the number $I^{\mathrm{<mo>\; *\; </mo>}}$ of infected in a quasi-stationary regime. If $R_p>1+\frac{\beta }{\delta }$, for all $i\ge 1$, $Q_{\mathrm{.}i}^{\mathrm{<mo>\; *\; </mo>}}\approx \frac{1}{I_e^{\mathrm{<mo>\; *\; </mo>}}}{\left(1-\frac{1}{I_e^{\mathrm{<mo>\; *\; </mo>}}}\right)}^{i-1}$

Proof. If $R_p>1$, we have $\underset{t\to +\infty }{\lim }\bar{S}\left(t\right)={\stackrel{*}{S}}_e=\frac{\mu +\gamma +\delta }{\beta }=\frac{\nu }{\beta }$.

The approximation of the process $I\left(t\right)$ by a birth-death process does not lead to a satisfactory result. In fact we will use the recursive method of Nåssell [17] and estimate the characteristics of quasi-stationary distribution of the process by those of a diffusion in stationary regime.

Let $Q_{s\mathrm{,}i}\left(t\right)$ the conditional distribution of $X_t=\left(S\left(t\right)\mathrm{,}I\left(t\right)\right)$ given that the epidemic has not extinguished:

$Q_{s,i}\left(t\right)=\frac{P_{s,i}\left(t\right)}{1-P_{.0}\left(t\right)}\text{et}{Q}_{.i}\left(t\right)={\displaystyle \underset{s\ge 0}{\sum }}{Q}_{s,i}\left(t\right)$ (15)

$Q_{\mathrm{.}i}\left(t\right)$ denotes the marginal distribution of $I\left(t\right)$, conditional on non-extinction. From Kolmogorov forward equations (Equation (6)) we obtain the following system:

$\begin{array}{c}{Q^{\prime }}_{s,i}\left(t\right)=n\left(1-p\right)Q_{s-1,i}\left(t\right)+\beta \left(s+1\right)\left(i-1\right)Q_{s+1,i-1}\left(t\right)\\ +\left(\mu +\gamma +\delta \right)\left(i+1\right)Q_{s,i+1}\left(t\right)+\delta \left(s+1\right)Q_{s+1,i}\left(t\right)\\ -\tau \left(s,i\right)Q_{s,i}\left(t\right)+\left(\mu +\gamma +\delta \right)Q_{.1}\left(t\right)Q_{s,i}\left(t\right)\end{array}$ (16)

The quasi-stationary distribution is the stationary probability distribution ${\left\{Q_{s\mathrm{,}i}^{\mathrm{<mo>\; *\; </mo>}}\right\}}_{s\ge \mathrm{0,}i\ge 0}$ satisfying:

$\begin{array}{l}0=n\left(1-p\right)Q_{s-1,i}^{*}+\beta \left(s+1\right)\left(i-1\right)Q_{s+1,i-1}^{*}+\left(\mu +\gamma +\delta \right)\left(i+1\right)Q_{s,i+1}^{*}\\ +\delta \left(s+1\right)Q_{s+1,i}^{*}-\tau \left(s,i\right)Q_{s,i}^{*}+\left(\mu +\gamma +\delta \right)Q_{.1}^{*}{Q}_{s,i}^{*}\end{array}$ (17)

the recursive method of Nåsell [17] gives the following recursion relationship for $Q_{\mathrm{.}i}^{\mathrm{<mo>\; *\; </mo>}}$ :

$\left(i+1\right)Q_{.i+1}^{*}=\frac{\delta R_0}{n}i{Q}_{.i}^{*}{e}_{S^{*}}\left(i\right)+Q_{.1}^{*}\left(1-{\displaystyle \underset{k=1}{\overset{i}{\sum }}}{Q}_{.k}^{*}\right)$ (18)

where $e_{S^{\mathrm{<mo>\; *\; </mo>}}}\left(i\right)={\displaystyle \underset{s\ge 0}{\sum }s\frac{Q_{s\mathrm{,}i}^{\mathrm{<mo>\; *\; </mo>}}}{Q_{\mathrm{.}i}^{\mathrm{<mo>\; *\; </mo>}}}}$ is the conditional expectation of $S^{\mathrm{<mo>\; *\; </mo>}}$ given that $I^{*}=i$. Let denote $\overline{S^{\mathrm{<mo>\; *\; </mo>}}}$ the expectation of $S^{\mathrm{<mo>\; *\; </mo>}}$, by summing the Equation (18) over $i\ge 1$ we obtain:

$Q_{.1}^{\mathrm{<mo>\; *\; </mo>}}=1-\frac{\delta R_0}{n}\left(\overline{S^{\mathrm{<mo>\; *\; </mo>}}}+\frac{cov\left(S^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}{I}^{\mathrm{<mo>\; *\; </mo>}}\right)}{\overline{I^{\mathrm{<mo>\; *\; </mo>}}}}\right)$ (19)

It follows that the marginal distribution $Q_{\mathrm{.}i}^{\mathrm{<mo>\; *\; </mo>}}$ satisfies the recursion relationship

$\{\begin{array}{l}{Q}_{.i+1}^{*}=\frac{\delta R_0}{n}\frac{i}{i+1}{Q}_{.i}^{*}{e}_{S^{*}}\left(i\right)+\frac{Q_{.1}^{*}}{i+1}\left(1-{\displaystyle \underset{k=1}{\overset{i}{\sum }}}{Q}_{.k}^{*}\right)\\ Q_{.1}^{*}=1-\frac{\delta R_0}{n}\left(\overline{S^{*}}+\frac{cov\left(S^{*},I^{*}\right)}{\overline{I^{*}}}\right)\end{array}$ (20)

Let denote ${\hat{X}}_e=\left(S_e^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}{I}_e^{\mathrm{<mo>\; *\; </mo>}}\right)$ the endemic equilibrium point in Equation (3) and $\alpha =\frac{\mu +\gamma +\delta }{\delta }$. We approximate the process $X_t-{\hat{X}}_e$ by a diffusion process with drift matrix $A_e$ and covariance matrix $co{v}_e$

$A_e=\delta \left(\begin{array}{cc}-R_p& -\alpha \\ R_p-1& 0\end{array}\right)\text{et}co{v}_e=\frac{n\left(1-p\right)}{R_p}\left(\begin{array}{cc}2R_p& -\left(R_p-1\right)\\ -\left(R_p-1\right)& 2\left(R_p-1\right)\end{array}\right)$ (21)

its stationary distribution is approximately bivariate normal with mean 0 and covariance matrix [17]:

$\Sigma =\frac{n\left(1-p\right)}{\delta {\left(R_p\right)}^2}\left(\begin{array}{cc}\alpha +R_p& -R_p\\ -R_p& R_p-1+\frac{R_p^2}{\alpha }\end{array}\right)$ (22)

from Equation (18) we approximate $cov\left(S^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}{I}^{\mathrm{<mo>\; *\; </mo>}}\right)$ by ${\Sigma }_{12}={\Sigma }_{21}=-\frac{n\left(1-p\right)}{\delta R_p}$ and $\left(\overline{S^{\mathrm{<mo>\; *\; </mo>}}}\mathrm{,}\overline{I^{\mathrm{<mo>\; *\; </mo>}}}\right)$ by $\left(S_e^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}{I}_e^{\mathrm{<mo>\; *\; </mo>}}\right)$, thus we obtain $Q_{.1}^{*}=\frac{1}{I_e^{*}}$.

In other hand since the stochastic means are close to equilibrium points in stationary regime we have set $e_{S^{*}}\left(i\right)=c{S}_e^{*}$ where c is a constant, and from the Equation (18) necessarily $c=\left(1-Q_{.1}^{\mathrm{<mo>\; *\; </mo>}}\right)$. It follows that the marginal distribution $Q_{\mathrm{.}i}^{\mathrm{<mo>\; *\; </mo>}}$ satisfies the recursion relationship

$\{\begin{array}{l}{Q}_{.i+1}^{*}=\frac{i}{i+1}{Q}_{.i}^{*}\left(1-Q_{.1}^{*}\right)+\frac{Q_{.1}^{*}}{i+1}\left(1-{\displaystyle \underset{k=1}{\overset{i}{\sum }}}{Q}_{.k}^{*}\right)\\ Q_{.1}^{*}=\frac{1}{I_e^{*}}\end{array}$ (23)

The only solution of which is $Q_{.i}^{*}=\frac{1}{I_e^{*}}{\left(1-\frac{1}{I_e^{*}}\right)}^{i-1}$, pour tout $i\ge 1$. The condition

of the theorem ensures that $I_e^{*}>1$ so the probability distribution ${\left(Q_{\mathrm{.}i}^{\mathrm{<mo>\; *\; </mo>}}\right)}_{i\ge 1}$ is well defined. The proof of theorem 4 is complete. □

5. Simulation

Two sample paths of $I\left(t\right)$ in Figure 3 are simulated over the interval $\left[0,25\right]$,

Figure 3. Two sample paths of $I\left(t\right)$ in quasi-stationary regime for the parameter values $S_0=100$ ; $\beta =0.69$ ; $\delta =0.25$ ; $\mu =0.02$ ; $\gamma =0.5$ ; $n=3.5$ ; $p=0.51$ ; $t\in \left[0,25\right]$. $R_p=6.1473$. All simulations started from the quasi-stationary distribution (Equation (23)).

Figure 4. approximations $Q_{\mathrm{.}i}^{\mathrm{<mo>\; *\; </mo>}}$ and $O_{\mathrm{.}i}^{\mathrm{<mo>\; *\; </mo>}}$ (dashed curve), for the parameter values $\beta =0.69$ ; $\delta =0.25$ ; $\mu =0.02$ ; $\gamma =0.5$ ; $n=3.5$ ; $p=0.51$. ${\mathcal{R}}_p=6.1473$.

with initial distributions the quasi-stationary distribution estimated at Equation (23). In the case $R_p>1$, the approximation of the quasi-stationary distribution of I by a diffusion process gives a normal distribution whit mean ${\mu }_I$ and variance ${\sigma }_I^2$ [17] [18].

${\mu }_I=I_e\text{and}{\sigma }_I^2=\frac{n\left(1-p\right)}{\delta {\left(R_p\right)}^2}\left(R_p-1+{\left(R_p\right)}^2\right).$

since $I>0$, the approximation of this distribution is

$O_{.i}^{*}=\frac{\varphi \left(\frac{i-{\mu }_I}{{\sigma }_I}\right)}{{\sigma }_I\Phi \left(\frac{{\mu }_I-0.5}{{\sigma }_I}\right)}$.

where $\varphi$ and $\Phi$ denote the normal density function and the normal cumulative distribution function respectively. Figure 4 gives $Q_{\mathrm{.}i}^{\mathrm{<mo>\; *\; </mo>}}$ and $O_{\mathrm{.}i}^{\mathrm{<mo>\; *\; </mo>}}$, for the values of $i\in \left\{\mathrm{1,2},\cdots ,100\right\}$.

6. Conclusion and Discussion

The asymptotic analysis of our model identifies three regions in parameter space with qualitatively different behaviors of the quasi-stationary distribution. $R_p$ is significantly greater or less than the deterministic threshold value 1 in the first region or the third region respectively and that $R_p$ is in a second region can be the transition region close to the deterministic threshold value 1 namely

$1\le R_p\le 1+\frac{\beta }{\delta }$ or $\frac{1}{1-p}\le R_0\le \frac{\delta +\beta }{\delta \left(1-p\right)}$. This result is analogous to that of Nåsell [17] for models SIS and SIR with finite population size N where the transition region is determined by a parameter $-3<\rho \left(N,R_0\right)<3$. The fact that a given population belongs to one of these three regions is governed by the size N of the population and the number $R_0$. In our study, this is a function of $\rho \left(\beta ,\delta \right)=1+\frac{\beta }{\delta }$. For example, for parameters $\beta =0.69$ and $\delta =0.25$ the transition region corresponds to $1\le R_p\le 3.76$, and if in addition $p=0.51$ we have $2.04\le R_0\le 7.67$. As Nåsell [17] noted, the approximation of the quasi-stationary distribution in the transition region remains a relatively complicated problem.

Note that $\frac{1}{\delta }$ is average life length in the target population. For measles the target group remains the children of 0 to 5 years. The mortality rate of children before five years of age in Niger is 280% [19]. This means that one in three children does not reach the age of five years and therefore we have $\frac{1}{\delta }<5$ years. The parameter $\beta$ in our model hardly exceeds 1. We see by this estimate that the transition region for $R_p$ is in the interval $\left[1,6\right[$ regardless of the size of the population.

Another approach for understanding the dynamics of the system before absorption is the ratio of means approach. Given the initial state, the ratio of expectations distribution (RE) is defined as a ration between the time that the process spends at each transient state and the expected time to absorption (provided that the expected time to absorption is finite). Precisely let $T_i$ be the time that the process ${\left(X_t\right)}_{t\ge 0}$, starting from state j, spending in state i before absorption, the

RE is defined by [10]: ${\hat{Q}}_{ji}\stackrel{def}{=}\frac{{\mathbb{E}}_j\left(T_i\right)}{{\mathbb{E}}_j\left(T\right)}$, where T denotes the absorption time. If ${\left(X_t\right)}_{t\ge 0}$ is a birth and death process with state space $\left\{\mathrm{0,1,}\cdots ,N\right\}$ and infinitesimal generator

$q_{ik}=\{\begin{array}{cc}\lambda \left(i\right)& \text{if}k=i+1\\ \mu \left(i\right)& \text{if}k=i-1\end{array}$ (24)

we have [10]:

${\hat{Q}}_{ji}=\frac{\frac{1}{\mu \left(i\right)}{\displaystyle \underset{k=1}{\overset{\min \left(j,i\right)}{\sum }}}{\displaystyle \underset{n=1}{\overset{i-1}{\prod }}}\frac{\lambda \left(n\right)}{\mu \left(n\right)}}{{\displaystyle \underset{i=1}{\overset{N}{\sum }}}\frac{1}{\mu \left(i\right)}{\displaystyle \underset{k=1}{\overset{\min \left(j,i\right)}{\sum }}}{\displaystyle \underset{n=1}{\overset{i-1}{\prod }}}\frac{\lambda \left(n\right)}{\mu \left(n\right)}},1\le i,j\le N$ (25)

as $N\to \infty$ and ${\left(X_t\right)}_{t\ge 0}$ is a linear birth-death process with infinitesimal generator given in Equation (13), It is easy to see that

${\hat{Q}}_{ji}=\left(1-\frac{\lambda }{\nu }\right){\left(\frac{\lambda }{\nu }\right)}^{i-1}=Q_{.i}^{*}$

Indeed, the RE is a another natural measure of the behavior of absorbing Markov chains before absorption, but the approximation is good only if the convergence to quasi-stationarity is relatively fast [10]. This measure is useful in studying the steady state of outbreaks measles epidemic.

Conflicts of Interest

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

References

[1]	Bartlett, T. (1987) The Molecular Biology of the Morbillivirus (Measles) Group. Biochemical Society Symposium, 53, 25-37.
[2]	World Health Organization (2019) Fact about Measles.
[3]	Edward, S., Kitengeso, R.E., Kiria, G.T., Felician, N., Mwema, G.G., 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
[4]	Wikipédia (2019) Modèles compartimentaux en épidémiologie. https://fr.wikipedia.org/wiki/Mod/%60eles_compartimentaux_en_/'epid/'emiologie
[5]	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.
[6]	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
[7]	Moussa Tessa, O. (2006) Mathematical Model for Control of Measles by Vaccination. Proceedings of Mali Symposium on Applied Sciences, Bamako, August 2006, 31-36.
[8]	Seydou, M. and Moussa Tessa, O. (2021) A Stochastic SVIR Model for Measles. Applied Mathematics, 12, 209-223. https://doi.org/10.4236/am.2021.123013
[9]	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.2307/3212311
[10]	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
[11]	Jacquez, J.A. and O’Neill, P. (1991) Reproduction Numbers and Thresholds in Stochastic Epidemic Models I. Homogeneous Populations. Mathematical Biosciences, 107, 161-186. https://doi.org/10.1016/0025-5564(91)90003-2
[12]	Anderson, R.M. and May, R.M. (1991) Infectious Diseases of Humans: Dynamics and Control. Oxford University Press, Oxford.
[13]	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
[14]	Hethcote, H.W. (2000) The Mathematics of Infectious Diseases. SIAM Review, 42, 599-653. https://doi.org/10.1137/S0036144500371907
[15]	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
[16]	Nasell, I. (2001) Extinction and Quasi-Stationarity in the Verhulst Logistic Model. Journal of Theoretical Biology, 211, 11-27. https://doi.org/10.1006/jtbi.2001.2328
[17]	Nasell, I. (1999) On the Time to Extinction in Recurrent Epidemics. Journal of the Royal Statistical Society: Series B (Statistical Methodology), 61, 309-330. https://doi.org/10.1111/1467-9868.00178
[18]	Gillespie, D.T. (1996) Exact Numerical Simulation of the Ornstein-Uhlenbeck Process and Its Integral. Physical Review E, 54, 2084-2091. https://doi.org/10.1103/PhysRevE.54.2084
[19]	C. Une Mortalité Infantile Très élevée, Liée Notamment à une Malnutrition Endémique (2021) https://www.senat.fr/rap/r04-512/r04-5126.html