Asymptotic Analysis of a Stochastic Model of Mosquito-Borne Disease with the Use of Insecticides and Bet Nets

Boubacar Sidiki Kouyaté^1*, Modeste N’zi²

¹Department of Mathematics and Computer Science, University of Science, Technical and Technologies of Bamako, Bamako, Mali.
²Laboratory of Applied Mathematics and Computer Science, University Felix Houphouet-Boigny, Abidjan, Côte d’Ivoire.

DOI: 10.4236/jamp.2024.121024                

Abstract

Ross’ epidemic model describing the transmission of malaria uses two classes of infection, one for humans and one for mosquitoes. This paper presents a stochastic extension of a deterministic vector-borne epidemic model based only on the class of human infectious. The consistency of the model is established by proving that the stochastic delay differential equation describing the model has a unique positive global solution. The extinction of the disease is studied through the analysis of the stability of the disease-free equilibrium state and the persistence of the model. Finally, we introduce some numerical simulations to illustrate the obtained results.

Keywords

Vector-Borne Disease Epidemic Model, Stochastic Delay Differential Equations, Stochastic Stability, Lyapunov Functional Technique

1. Introduction

Stochastic differential equations (SDEs) in various modifications as descriptions of stochastic dynamical systems have been used in biological and medical problems [1] [2] , financial and economic problems [3] [4] , etc., in which significant uncertainty is present. In some real-world applications, many phenomena studied do not only depend on the present state, but also on the past ones (see [5] [6] [7] ). For instance, people infected with an infectious disease usually develop symptoms on average a few days or weeks (or even years for AIDS) after infection [8] . Stochastic delay differential equations (SDDEs) have been widely used to model such systems (see, e.g. [9] [10] ).

The SIRS epidemic model is used to describe the evolution of an endemic infection in which any infected person either dies or recovers from the disease and becomes temporarily immune [11] . This model relies on a linear incidence rate, which may not account for the saturation of effective contacts between infected and susceptible individuals due to the presence of already infected individuals (see, e.g. [12] [13] [14] ). Consequently, to capture some of the real characteristics of the infection process, non-linear incidence rates need to be taken into account.

Vector-borne diseases are one of the most dangerous diseases with more than 223 million cases and 700,000 deaths worldwide in 2020 [15] . However, many vector-borne diseases can be prevented through the application of protective measures. Mathematical modeling to study vector-borne diseases has a long history dating back to 1911 with the Ross model [16] , subsequently, important extensions were proposed by MacDonald [17] [18] . Since then, there have been numerous extensions and adaptations, such as the inclusion of acquired immunity proposed by Dietz, Molineaux and Thomas [19] . Some works have also included environmental effects [20] [21] , the spread of drug resistance [22] [23] , the treatment and impact of vaccination strategies [24] , as well as the timing of the period incubation [25] . Authors have also taken into account the effects of individual protection measures, such as the use of impregnated mosquito nets or repellents [26] [27] [28] , spatial dynamics [29] , the heterogeneity of hosts [30] , seasonality [31] , stochasticity [32] [33] and control [34] . The Ross epidemic model describes the dynamics of malaria transmission using two classes of infection, one for human hosts and the other for mosquito vectors. Another perspective for studying epidemic patterns for vector-borne diseases, such as malaria, is to approximately transform infectious mosquitoes in Ross-Mcdonald type model into infectious hosts by means of time-scale transformation. This way, we only need to study the transmission dynamics of the disease in the host population (see [35] [36] ). For example, Enatsu et al. [7] consider an endemic delayed epidemic model with nonlinear incidence rates in the form:

$\rho S\left(t\right){\displaystyle {\int }_0^{{\tau }_{\max }}}\Theta \left(\tau \right)\Psi \left(A\left(t-\tau \right)\right)\text{d}\tau$ (1)

where $S\left(t\right)$ , $A\left(t\right)$ and $R\left(t\right)$ denote respectively the fractions of susceptible, infective and recovered individuals at time t and Ψ, a nonlinear function satisfying some assumptions. This incidence rate is used to study the transmission of disease, which is caused by a pathogenic germ carried and inoculated by vectors that have an incubation time to become infectious. According to [36] , the vectors can be omitted from the equations by including a delay in the force of infection. Other authors (see [32] [37] ) have studied a stochastic model describing malaria transmission dynamics with compartments based solely on human host population. As mentioned in the annex of [37] , Wanduku considers a stochastic malaria epidemic model that is an extension of the model proposed in [35] [36] . But does not link the parameters of the proposed model to the parameters of the Ross-McDonald epidemic model from which it is derived nor the necessary hypotheses.

This paper deals with a delayed stochastic epidemic model with non-linear incidence rate describing the transmission dynamics of a vector-borne disease that is a stochastic extension of the model studied in [7] . Based on the results and assumptions in [36] , we linked the model parameters to the parameters of the corresponding Ross-MacDonald model. Furthermore, the practice of protective measures by part of the population is taken into account. The stochastic delay differential equation describing the model is obtained from a deterministic model by introducing noise into the contact rate due to environmental variations. We assume that the intensity of the stochastic disturbance is proportional to the number of infective individuals. Our study is broken down into three stages. First, we establish the consistency of the stochastic epidemic model by proving the existence of a unique global positive solution of the stochastic differential equations with delay describing the model. To study the extinction of disease, we secondly analyze the stability of the disease-free equilibrium point of the stochastic model under the condition $R_0<1$ . Where R₀ designates the reproduction number of the underlying deterministic model. In the third step, we study the persistence of the solution when the disease-free equilibrium is unstable. We conclude our study with numerical simulations to illustrate the importance of personal protective measures in disease control.

The remainder of the paper is structured as follows. In Section 2, we introduce the model and some preliminary definitions. In Section 3, we establish the existence and uniqueness of a global positive solution of our stochastic model. Section 4 deals with the stability of the disease-free equilibrium E₀, which is the unique equilibrium state of the stochastic model obtained. In the case where the disease-free equilibrium state is effectively unstable, the persistence of the stochastic model when $R_0>1$ is established in Section 5. In Section 6, some numerical simulations are given to illustrate the mathematical results. Finally, we conclude and propose some perspectives in Section 7.

2. Model Description

2.1. Deterministic Model Description

In this work, we propose to study an epidemic model of vector-borne disease transmission which makes it possible to take into account the effects of protective measures. The is We assume that the human population is subdivided into four compartments named SARB where S designates the class of individuals susceptible to infection, A the class of individuals affected by the pathogen after being bitten by infectious mosquitoes, R is the class of individuals treated and cleared of the pathogenic agent and who become susceptible again after a certain period of immunity, finally, B designates the class of individuals completely withdrawn from the epidemic process by using protective measures such as the use of insecticide-treated musketeers.

We recognize that the vector population is divided into two groups named susceptible and infectious. The size of the population of the vectors $N_v$ is constant so that the death and birth rates are the same ${\mu }_v$ . Susceptible mosquitoes are affected by the pathogen after biting an infectious human and will only be infectious after a certain latency period. An infectious mosquito stays infectious until it dies. Let N₀ be the average number of humans in the model such that $\frac{N_0}{N_v}\ll 1$ . $T_{\tau }$ the latency period for an infected mosquito to become infectious is distributed according to the probability distribution of density function Θ. ${\beta }_v$ denotes the contact rate between infected people and susceptible vectors such.

We make the following assumptions about the model:

· Bed net use hypothesis: Many efforts are made in the population to ensure that infection from infectious humans to susceptible vectors is low, so that $N_0{\beta }_v\approx 1$ .

· Insecticide use hypothesis: Many measures are implemented in the population to significantly reduce the vector population using insecticides or traps, so that ${\mu }_v$ is chosen large enough such that ${\epsilon }_v=\frac{{\text{e}}^{-{\mu }_v\mathbb{E}\left(T_{\tau }\right)}}{{\mu }_v}\ll 1$ .

Based on the works carried out in [36] , by using a time-scale transformation under the condition ${\epsilon }_v=\frac{{\text{e}}^{-{\mu }_v\mathbb{E}\left(T_{\tau }\right)}}{{\mu }_v}\ll 1$ , the vectors can be omitted from model equations by including a delay in the force in affected subgroup “A” such that the model can be described by the following flowchart (see Figure 1).

Using the non-linear incidence rate (1) proposed by Enatsu et al. [7] which is a generalization of such a class of model. Therefore, the differential system describing the model is given by:

$\{\begin{array}{l}\frac{\text{d}S\left(t\right)}{\text{d}t}=\varphi -\left({\mu }_1+b_1\right)S\left(t\right)-m_v{\epsilon }_v{\beta }_{h}S\left(t\right){\displaystyle {\int }_0^{{\tau }_{\max }}}\Theta \left(\tau \right)\Psi \left(A\left(t-\tau \right)\right)\text{d}\tau +\nu R\left(t\right),\hfill \\ \frac{\text{d}A\left(t\right)}{\text{d}t}=m_v{\epsilon }_v{\beta }_{h}S\left(t\right){\displaystyle {\int }_0^{{\tau }_{\max }}}\Theta \left(\tau \right)\Psi \left(A\left(t-\tau \right)\right)\text{d}\tau -\left({\mu }_2+\gamma +b_2\right)A\left(t\right),\hfill \\ \frac{\text{d}R\left(t\right)}{\text{d}t}=\gamma A\left(t\right)-\left({\mu }_3+\nu +b_3\right)R\left(t\right),\hfill \\ \frac{\text{d}B\left(t\right)}{\text{d}t}=b_{1}S\left(t\right)+b_{2}A\left(t\right)+b_{3}R\left(t\right)-{\mu }_{4}B\left(t\right),\hfill \end{array}$ (2)

the initial value of the model is given by:

$\{\begin{array}{l}S\left(\theta \right)={\eta }_1\left(\theta \right),A\left(\theta \right)={\eta }_2\left(\theta \right),R\left(\theta \right)={\eta }_3\left(\theta \right),B\left(\theta \right)={\eta }_4\left(\theta \right);\\ {\eta }_i\left(0\right)>0,{\eta }_i\in \mathcal{C}\left(\left[-{\tau }_{\max },0\right];{ℝ}^{+}\right)\text{for}i=1,2,3,4.\end{array}$

In this model, $\varphi$ is the rate at which new individuals appear in the susceptible compartment. The mortality rates of susceptible, affected, recovered and bed net use individuals are ${\mu }_i,i=1,2,3,4$ respectively, with the assumption ${\mu }_1\le \min \left\{{\mu }_2\mathrm{,}{\mu }_3{\mu }_4\right\}$ . The parameter ${\beta }_h$ denotes the disease contact rate between susceptible people and infected vectors, and $\gamma$ is the rate at which infected individuals recover from infection. $\nu$ is the rate at which recovered individuals lose their immunity and become susceptible again. The parameters

Figure 1. The chart flow of the model representing the different links between the compartments.

$b_1\le b_2\le b_3$ are the rates at which susceptible, infectious and recovered individuals, respectively, use disease protection practices. At last, ${\tau }_{\max }$ is the upper bound of the latency period in affected vectors. These parameters are assumed to be non-negative.

The functions Ψ and Θ satisfy the following conditions:

· (C₁) Θ is a probability density function with support $\left[0,{\tau }_{\max }\right]$ .

· (C₂) Ψ is Lipschitz and strictly increasing function on $\left[\mathrm{0,}\infty \right)$ with $\Psi \left(0\right)=0$ .

Therefore, $\Psi \left(x\right)\le x$ , $\forall x\ge 0$ .

· (C₃) Ψ is differentiable on $\left[\mathrm{0,}+\infty \right)$ such that ${\Psi }^{\prime }\left(0\right)=1$ .

In deterministic framework, the threshold that indicates if the disease persists in the population in large time or simple disappear is given by the reproduction number R₀ [38] . For this model,

$R_0=\frac{m_v{\epsilon }_v{\beta }_h\varphi {\Psi }^{\prime }\left(0\right)}{\left({\mu }_1+b_1\right)\left({\mu }_2+\gamma +b_2\right)}\mathrm{.}$

By a simple analysis we obtain that system (2) has a disease-free equilibrium $E_0=\left(s_0,i_0,r_0,b_0\right)=\left(\frac{\varphi }{{\mu }_1+b_1},0,0,\frac{b_1\varphi }{{\mu }_4\left({\mu }_1+b_1\right)}\right)$ . We will see in the following that E₀ is the unique equilibrium point of the stochastic model studied.

Recently many authors have studied a stochastic delayed epidemic model with perturbed parameter in various aspects (see, e.g. [37] [39] [40] ). The aim of this work is to study the stability analysis and the persistence of a stochastic version of model (2) by introducing noise in the contact rate.

2.2. Stochastic Model Derivation

Let’s consider a stochastic basis $\left(\Omega \mathrm{,}\mathcal{F}\mathrm{,}{\left\{{\mathcal{F}}_t\right\}}_{t\ge 0}\mathrm{,}P\right)$ with a filtration ${\left\{{\mathcal{F}}_t\right\}}_{t\ge 0}$ satisfying the usual conditions, on which are defined all random variables considered throughout this work. The contact process is inevitably affected by random perturbations due to stochastic environmental factors that can be modeled by a random variable ${\tilde{\beta }}_h$ with average value $\beta$ and variance ${\Xi }^2$ . Using the same ideas as in [41] [42] , the potentially infectious contacts ${\tilde{\beta }}_h\text{d}t$ made by each infected individual with each susceptible in the small time interval $\left[t\mathrm{,}t+\text{d}t\right]$ is approximately given by:

${\tilde{\beta }}_h\text{d}t={\beta }_h\text{d}t+\Xi \text{d}W\left(t\right)\mathrm{,}$

where W is a standard Brownian motion. On the other hand, we assume that, the increase in the number of infectious occurs with some spatial dispersions that increase the variability of contact processes. To take into account this situation, here we assume that the noise intensity at time t, depends on the infectious population size $A\left(t\right)$ . We obtain the stochastic model by replacing ${\beta }_h\text{d}t$ in (2) system by ${\tilde{\beta }}_h\text{d}t={\beta }_h\text{d}t+\sigma A\left(t\right)\text{d}W\left(t\right)$ , where $\sigma$ is a positive real. Then, the model is described by the stochastic delay differential equation:

$\{\begin{array}{l}\text{d}S\left(t\right)=\left[\varphi -\left({\mu }_1+b_1\right)S\left(t\right)-\rho S\left(t\right)O\left(A,t\right)+\nu R\left(t\right)\right]\text{d}t\\ -\sigma S\left(t\right)A\left(t\right)O\left(A,t\right)\text{d}W\left(t\right),\\ \text{d}A\left(t\right)=\left[\rho S\left(t\right)O\left(A,t\right)-\left({\mu }_2+\gamma +b_2\right)A\left(t\right)\right]\text{d}t+\sigma S\left(t\right)A\left(t\right)O\left(A,t\right)\text{d}W\left(t\right),\\ \text{d}R\left(t\right)=\left[\gamma A\left(t\right)-\left({\mu }_3+\nu +b_3\right)R\left(t\right)\right]\text{d}t,\\ \text{d}B\left(t\right)=\left[b_{1}A\left(t\right)+b_{2}S\left(t\right)+b_{3}R\left(t\right)-{\mu }_{4}B\left(t\right)\right]\text{d}t,\end{array}$ (3)

where, for all $t\ge 0$ ,

$O\left(A,t\right)={\displaystyle {\int }_0^{{\tau }_{\max }}}\Theta \left(\tau \right)\Psi \left(A\left(t-\tau \right)\right)\text{d}\tau \text{and}\rho =m_v{\epsilon }_v{\beta }_h.$

The description of the parameters is the same as in the deterministic model (2) with the same assumptions.

The initial condition is given by:

$\{\begin{array}{l}S\left(\theta \right)={\eta }_1\left(\theta \right),A\left(\theta \right)={\eta }_2\left(\theta \right),R\left(\theta \right)={\eta }_3\left(\theta \right),B\left(\theta \right)={\eta }_4\left(\theta \right),\theta \in \left[-{\tau }_{\max },0\right],\\ \eta =\left({\eta }_1\mathrm{,}{\eta }_2\mathrm{,}{\eta }_3\mathrm{,}{\eta }_4\right)\in {\mathcal{P}}^4\left(I\right)\mathrm{.}\end{array}$ (4)

Let $I$ be a compact and connected subset of ${ℝ}_{+}^n=\left\{\left(x_1\mathrm{,}\cdots \mathrm{,}{x}_n\right)\in {ℝ}^n\mathrm{<mo>\; :\; </mo>}{x}_1\ge \mathrm{0,}\cdots \mathrm{,}{x}_n\ge 0\right\}$ and $\left|x\right|={\left({\displaystyle {\sum }_{k=1}^n}{x}_i^2\right)}^{\frac{1}{2}}$ the Euclidean norm on ${ℝ}^n$ . Then, ${\mathcal{P}}^n\left(I\right)$ denotes the set of ${\mathcal{F}}_0$ -measurable $\mathcal{C}\left(\left[-{\tau }_{\max },0\right]\mathrm{<mo>\; ;\; </mo>}I\right)$ -valued random variables such that $\mathbb{E}\left({\Vert \eta \Vert }^2\right)<\infty$ where $\Vert \eta \Vert ={\sup }_{\theta \in \left[-{\tau }_{\max },0\right]}\left|\eta \left(\theta \right)\right|$ .

2.3. Definitions and Preliminary Results

Let $F\mathrm{<mo>\; :\; </mo>}\mathcal{C}\left(\left[-{\tau }_{\max },0\right]\mathrm{<mo>\; ;\; </mo>}I\right)\to {ℝ}^n$ be a n-dimensional functional and $G\mathrm{<mo>\; :\; </mo>}\mathcal{C}\left(\left[-{\tau }_{\max },0\right]\mathrm{<mo>\; ;\; </mo>}I\right)\to {\mathcal{M}}_{n\times m}\left(ℝ\right)$ be a $n\times m$ -matrix-valued functional. The functions F and G are Borel measurable.

Let $B={\left(B\left(t\right)\right)}_{t\ge 0}$ be a m-dimensional Brownian motion process. Consider the following n-dimensional stochastic system with time delay:

$\text{d}X\left(t\right)=F\left(X_t\right)\text{d}t+G\left(X_t\right)\text{d}B\left(t\right)\text{and}{X}_0=\eta \in {\mathcal{P}}^n\left(I\right)\mathrm{,}$ (5)

where $X_t=\left\{X\left(t+\theta \right),\theta \in \left[-{\tau }_{\max },0\right]\right\}$ is viewed as a $\mathcal{C}\left(\left[-{\tau }_{\max },0\right]\mathrm{<mo>\; ;\; </mo>}{ℝ}_{+}^n\right)$ -valued stochastic process with $X_t\left(0\right)=X\left(t\right)$ . Let $\left\{X\left(t\mathrm{,}\eta \right)\mathrm{<mo>\; :\; </mo>}t\ge 0\right\}$ the solution with initial value $X_0=\eta$ of the stochastic system (5).

Let $U\mathrm{<mo>\; :\; </mo>}\left[\mathrm{0,}\infty \right)\times {\mathcal{P}}^n\left(I\right)\to ℝ$ be a functional. The generating operator $\mathcal{L}$ of system (5) is defined (see, e.g. [10] [43] ) by the formula:

$\mathcal{L}U\left(t\mathrm{,}\eta \right)=\underset{\Delta \to 0}{lim}\frac{\mathbb{E}\left(U\left(t+\Delta \mathrm{,}{X}_{t+\Delta }\right)\mathrm{<mo>\; |\; </mo>}{X}_t=\eta \right)-U\left(t\mathrm{,}\eta \right)}{\Delta }\mathrm{.}$

Suppose that the functional U can be written in the form:

$U\left(t\mathrm{,}\eta \right)=U^0\left(t\mathrm{,}\eta \left(0\right)\mathrm{,}\eta \right),$

where $U^0$ is a $ℝ$ -valued functional defined on $\left[\mathrm{0,}\infty \right)\times {ℝ}^n\times {\mathcal{P}}^n\left(I\right)$ .

For any $\left(t\mathrm{,}x\right)\in \left[\mathrm{0,}\infty \right)\times {ℝ}^n$ and any $\eta \in {\mathcal{P}}^n\left(I\right)$ , we put:

$U_{\eta }\left(t\mathrm{,}x\right)=U^0\left(t\mathrm{,}x\mathrm{,}\eta \right),$

where $x=\eta \left(0\right)=X\left(t\right)$ and $\eta =X_t$ .

Let D be the class of all functional V for which functions $U_{\eta }\left(t\mathrm{,}x\right)$ are continuously twice differentiable in x and once in t. For functionals in D, the generating operator $\mathcal{L}$ of system (5) becomes:

$\begin{array}{c}\mathcal{L}U\left(t\mathrm{,}\eta \right)=\frac{\partial U_{\eta }\left(t\mathrm{,}\eta \left(0\right)\right)}{\partial t}+\nabla U_{\eta }^{\text{T}}\left(t\mathrm{,}\eta \left(0\right)\right)F\left(\eta \right)\\ +\frac{1}{2}trace\left[G^{\text{T}}\left(\eta \right){\nabla }^2{U}_{\eta }\left(t\mathrm{,}\eta \left(0\right)\right)G\left(\eta \right)\right]\mathrm{,}\end{array}$ (6)

where

$\nabla U_{\eta }\left(t\mathrm{,}x\right)=\left(\frac{\partial U_{\eta }\left(t\mathrm{,}x\right)}{\partial x_1}\mathrm{,}\cdots \mathrm{,}\frac{\partial U_{\eta }\left(t\mathrm{,}x\right)}{\partial x_n}\right)\mathrm{,}{\nabla }^2{U}_{\eta }\left(t\mathrm{,}x\right)={\left(\frac{{\partial }^2{U}_{\eta }\left(t\mathrm{,}x\right)}{\partial x_i\partial x_j}\right)}_{n\times n}\mathrm{.}$

The following theorem, which is a corollary of Theorem 3 in [44] , provides sufficient conditions for the stability of trivial solutions of system (5).

Theorem 1. Assume that both F and G satisfy the local Lipschitz condition and suppose that there exists a functional $U\left(t\mathrm{,}\eta \right)\in D$ such that:

$c_1{\left|\eta \left(0\right)\right|}^2\le U\left(t\mathrm{,}\eta \right)\le c_2{\Vert \eta \Vert }^2\text{and}\mathcal{L}V\left(\eta \mathrm{,}t\right)\le -\alpha {\left|\eta \left(0\right)\right|}^2,$

where $c_1\mathrm{,}{c}_2$ and $\alpha$ are positive constants. Then for all $\eta \in {\mathcal{P}}^n\left(I\right)$ , there exists a positive constant q such that the solution of system (5) satisfies:

$\underset{t\to \infty }{lim}\frac{1}{t}\ln \left(\left|X\left(t\mathrm{,}\eta \right)\right|\right)<-q\text{a}\text{.s}\mathrm{.}$

That is the trivial solution of (5) is almost surely exponentially stable.

Now, consider the following system formed by the first three equations of system (3):

$\{\begin{array}{l}\text{d}S\left(t\right)=\left[\varphi -\left({\mu }_1+b_1\right)S\left(t\right)-\rho S\left(t\right)O\left(A,t\right)+\nu R\left(t\right)\right]\text{d}t\\ -\sigma S\left(t\right)A\left(t\right)O\left(A,t\right)\text{d}W\left(t\right),\\ \text{d}A\left(t\right)=\left[\rho S\left(t\right)O\left(A,t\right)-\left({\mu }_2+\gamma +b_2\right)A\left(t\right)\right]\text{d}t+\sigma S\left(t\right)A\left(t\right)O\left(A,t\right)\text{d}W\left(t\right),\\ \text{d}R\left(t\right)=\left[\gamma A\left(t\right)-\left({\mu }_3+\nu +b_3\right)R\left(t\right)\right]\text{d}t,\end{array}$ (7)

with the initial condition

$\{\begin{array}{l}S\left(\theta \right)={\eta }_1\left(\theta \right),A\left(\theta \right)={\eta }_2\left(\theta \right),R\left(\theta \right)={\eta }_3\left(\theta \right),\theta \in \left[-{\tau }_{\max },0\right],\\ \eta =\left({\eta }_1\mathrm{,}{\eta }_2\mathrm{,}{\eta }_3\right)\in {\mathcal{P}}^3\left(I\right)\mathrm{.}\end{array}$ (8)

It’s straightforward to see that, $E_0=\left(\frac{\varphi }{{\mu }_1+b_1}\mathrm{,0,0,}\frac{b_1\varphi }{{\mu }_4\left({\mu }_1+b_1\right)}\right)$ and $E_0^{\mathrm{<mo>\; *\; </mo>}}=\left(\frac{\varphi }{{\mu }_1+b_1}\mathrm{,0,0}\right)$ are the disease-free equilibria of system (3) and (7), respectively.

Note that the first three equations of system (3) are independent of the fourth. In the following result, we prove that, the study of system (3) with initial condition (4) can be reduced to the study of system (7) with initial condition (8).

Lemma 2. Let assume that all component of the solution of the reduced system (7) with initial condition (8) are positive and ${\eta }_4\left(\theta \right)\ge 0$ , for all $\theta \in \left[-{\tau }_{\max },0\right]$ . Then,

1) Any solution of system (3) with initial condition (4) is positive.

2) The almost sure stability of the disease-free equilibrium $E_0^{\mathrm{<mo>\; *\; </mo>}}$ of the reduced system (7) with the initial condition (8) implies the almost sure stability of the disease-free equilibrium $E_0$ of system (3) with initial condition (4).

Proof.

1) Given the assumptions of the lemma, $S\left(t\right)\ge 0$ , $A\left(t\right)\ge 0$ , $R\left(t\right)\ge 0$ and ${\eta }_4\ge 0$ . By comparison theorem of ordinary differential equation, it follows that:

$\text{d}B\left(t\right)=\left(b_{1}S\left(t\right)+b_{2}A\left(t\right)+b_{3}R\left(t\right)-{\mu }_{4}B\left(t\right)\right)\text{d}t\ge -{\mu }_{4}B\left(t\right)\text{d}t,$

that is $B\left(t\right)\ge B\left(0\right){\text{e}}^{-{\mu }_{4}t}>0$ .

1) Now, let assume that the disease-free equilibrium $E_0^{\mathrm{<mo>\; *\; </mo>}}$ of the reduced system (7) with the initial condition (8) is stable, that is:

$\left(S\left(t\right)\mathrm{,}A\left(t\right)\mathrm{,}R\left(t\right)\right)\to \left(\frac{\varphi }{{\mu }_1+b_1}\mathrm{,0,0}\right)\text{a}\text{.s}\mathrm{.}$

Therefore, for all $\epsilon >0$ , there exists a real $t^{*}>0$ enough large, such that:

$\text{d}B\left(t\right)=\left(b_1\left(\frac{\varphi }{{\mu }_1+b_1}+\epsilon \right)-{\mu }_{4}B\left(t\right)\right)\text{d}t\text{for}\text{all}t\ge t^{*}.$

It follows that:

$B\left(t\right)=\frac{b_1}{{\mu }_4}\left(\frac{\varphi }{{\mu }_1+b_1}+\epsilon \right)+\left(B\left(0\right)-\frac{b_1\varphi }{{\mu }_4\left({\mu }_1+b_1\right)}\right){\text{e}}^{-{\mu }_{4}t}\text{a}\text{.s}\text{.}$

Consequently, by letting $\epsilon \to 0$ , we obtain that:

$\left(S\left(t\right)\mathrm{,}A\left(t\right)\mathrm{,}R\left(t\right)\mathrm{,}B\left(t\right)\right)\to \left(\frac{\varphi }{{\mu }_1+b_1}\mathrm{,0,0,}\frac{b_1\varphi }{{\mu }_4\left({\mu }_1+b_1\right)}\right)\text{a}\text{.s}\text{.}$ £

The stochastic model (7) with the initial condition (8) can be written in the form of (5) where $B\left(t\right)=\left(W\left(t\right)\mathrm{,}{W}^{\mathrm{<mo>\; *\; </mo>}}\left(t\right)\mathrm{,}{W}^{\mathrm{<mo>\; *\; </mo>\; <mo>\; *\; </mo>}}\left(t\right)\right)$ is a three-dimensional Brownian motion process. For any $\eta \in \mathcal{C}\left(\left[-{\tau }_{\max },0\right]\mathrm{<mo>\; ;\; </mo>}{ℝ}_{+}^3\right)$ , we have:

$F\left(\eta \right)=\left(\begin{array}{c}\varphi -\left({\mu }_1+b_1\right){\eta }_1\left(0\right)-\rho {\eta }_1\left(0\right){\displaystyle {\int }_{-{\tau }_{\max }}^0}\Theta \left(-\theta \right)\Psi \left({\eta }_2\left(\theta \right)\right)\text{d}\theta -\nu {\eta }_3\left(0\right)\\ \rho {\eta }_1\left(0\right){\displaystyle {\int }_{-{\tau }_{\max }}^0}\Theta \left(-\theta \right)\Psi \left({\eta }_2\left(\theta \right)\right)\text{d}\theta -\left({\mu }_2+\gamma +b_2\right){\eta }_2\left(0\right)\\ \gamma {\eta }_2\left(0\right)-\left({\mu }_3+\nu +b_3\right){\eta }_3(\; 0\; )\end{array}\right)$

$G\left(\eta \right)=\left(\begin{array}{ccc}-\sigma {\eta }_1\left(0\right){\eta }_2\left(0\right){\displaystyle {\int }_{-{\tau }_{\max }}^0}\Theta \left(-\theta \right)\Psi \left({\eta }_2\left(\theta \right)\right)\text{d}\theta & 0& 0\\ \sigma {\eta }_1\left(0\right){\eta }_2\left(0\right){\displaystyle {\int }_{-{\tau }_{\max }}^0}\Theta \left(-\theta \right)\Psi \left({\eta }_2\left(\theta \right)\right)\text{d}\theta & 0& 0\\ 0& 0& 0\end{array}\right).$

3. Consistency of the Model

In this section, we study the existence and uniqueness of global positive solution for model (3). To do this, we will first establish the existence and uniqueness of global positive solution for system (7) and conclude using Lemma 2. Under the assumptions (C₁ - C₃), the coefficients of system (7) are locally Lipschitz continuous. We therefore deduce that, for any initial condition $\eta \in {\mathcal{P}}^3\left(I\right)$ , system (7) has a unique local solution on $t\in \left[\mathrm{0,}{\vartheta }_e\right)$ , where ${\vartheta }_e$ denotes the explosion time (see, e.g. [9] Theorem 2.8 on page 154). In order to prove that the local solution is global and positive, we will first establish the existence and uniqueness of positive local solution. Therefore, we deduce that this solution does not explode towards infinity in a finite time, i.e. ${\vartheta }_e=\infty$ . Let us set:

$\Delta =\left\{\left(x,y,z\right)\in {ℝ}^3/x>0,y>0,z>0,x+y+z<\frac{\varphi }{{\mu }_1+b_1}\right\}$

and ${\mathcal{P}}^4\left(\Delta \right)$ be the class of ${\mathcal{F}}_0$ -measurable and $\mathcal{C}\left(\left[-{\tau }_{\max },0\right];\Delta \right)$ valued random variables.

Let $\mathcal{N}\left(t\right)=S\left(t\right)+A\left(t\right)+R\left(t\right)$ be the total population excluding the group of isolated people at time $t\in \left[-{\tau }_{\max },{\vartheta }_e\right[$ .

We consider the following stopping times:

${\vartheta }^{\mathrm{<mo>\; \{\; </mo>}-\mathrm{<mo>\; \}\; </mo>}}=\inf \left\{t\in \left[\mathrm{0,}\infty \right)\mathrm{,}\min \left\{S\left(t\right)\mathrm{,}A\left(t\right)\mathrm{,}R\left(t\right)\right\}\notin \left]\mathrm{0,}\infty \right[\right\}\mathrm{,}$

${\vartheta }_0=\inf \left\{t\in \left[\mathrm{0,}\infty \right)\mathrm{,}\min \left\{S\left(t\right)\mathrm{,}A\left(t\right)\mathrm{,}R\left(t\right)\right\}\notin \left]\mathrm{0,}\frac{\varphi }{{\mu }_1+b_1}\right[\right\}\mathrm{,}$

${\vartheta }_n=\inf \left\{t\in \left[\mathrm{0,}\infty \right)\mathrm{,}\min \left\{S\left(t\right)\mathrm{,}A\left(t\right)\mathrm{,}R\left(t\right)\right\}\notin \left]\frac{1}{n}\mathrm{,}\frac{\varphi }{{\mu }_1+b_1}\right[\right\}\mathrm{,}$

for all integers $n\ge n_0$ , where $n_0\in {\mathbb{N}}^{\ast }$ is such that $\min \left\{S\left(0\right)\mathrm{,}A\left(0\right)\mathrm{,}R\left(0\right)\right\}\in \left]\frac{1}{n_0}\mathrm{,}\frac{\varphi }{{\mu }_1+b_1}\right[$ , with the convention $\inf \varnothing =\infty$ .

In the following lemma, we establish that ${\vartheta }^{\mathrm{<mo>\; \{\; </mo>}-\mathrm{<mo>\; \}\; </mo>}}={\vartheta }_e$ , which means the positivity of the local solution of model (7) described above. We also establish the convergence of the sequence of stopping times $\left\{{\vartheta }_n\mathrm{<mo>\; :\; </mo>}n\in {\mathbb{N}}^{\mathrm{<mo>\; *\; </mo>}}\text{and}n\ge n_0\right\}$ towards the explosion time ${\vartheta }_e$ which will be used to establish that ${\vartheta }_e=\infty$ .

Lemma 3. Let us assume that the initial condition $\eta$ of (7) belongs to ${\mathcal{P}}^3\left(\Delta \right)$ . Then,

1) ${\vartheta }^{\mathrm{<mo>\; \{\; </mo>}-\mathrm{<mo>\; \}\; </mo>}}\ge {\vartheta }_e$ a.s.

2) ${\sup }_{t\in \left[-{\tau }_{\max },{\vartheta }_e\right)}\mathcal{N}\left(t\right)\le \frac{\varphi }{{\mu }_1+b_1}$ a.s.

3) Moreover ${\vartheta }_0={\vartheta }^{\{-\}}={\vartheta }_e$ a.s.

The sequence of stopping times ${\left({\vartheta }_n\right)}_{n\ge n_0}$ converges to ${\vartheta }_0$ a.s.

Proof. Obviously, for $t\in \left[-{\tau }_{\max },{\vartheta }_e\wedge {\vartheta }^{\{-\}}\right[$ , $\min \left\{S\left(t\right),A\left(t\right),R\left(t\right)\right\}>0$ . In view of Itô’s formula for all $t\in \left[-{\tau }_{\max },{\vartheta }_e\wedge {\vartheta }^{\{-\}}\right[$ , we have:

$\begin{array}{l}\ln \left[S\left(t\right)A\left(t\right)R\left(t\right)B\left(t\right)\right]-\ln \left[{\eta }_1\left(0\right){\eta }_2\left(0\right){\eta }_3\left(0\right)\right]\\ ={\displaystyle {\int }_0^t}\left[\frac{\varphi }{S\left(s\right)}-{\mu }_1-{\displaystyle \underset{i=1}{\overset{3}{\sum }}}{b}_i+\nu \frac{R\left(s\right)}{S\left(s\right)}-\rho O\left(A,s\right)+\rho \frac{S\left(s\right)}{A\left(s\right)}O\left(A,s\right)-\left({\mu }_2+\gamma \right)\right]\text{d}s\\ +{\displaystyle {\int }_0^t}\left[\gamma \frac{A\left(s\right)}{R\left(s\right)}-\left({\mu }_3+\nu \right)-\frac{{\sigma }^2}{2}{\left(A\left(s\right)O\left(A,s\right)\right)}^2-\frac{{\sigma }^2}{2}{\left(S\left(s\right)O\left(A,s\right)\right)}^2\right]\text{d}s\\ -\sigma {\displaystyle {\int }_0^t}A\left(s\right)O\left(A,s\right)\text{d}W\left(s\right)+\sigma {\displaystyle {\int }_0^t}S\left(s\right)O\left(A,s\right)\text{d}W\left(s\right).\end{array}$

It follows that:

$\begin{array}{l}\ln \left[S\left(t\right)A\left(t\right)R\left(t\right)\right]-\ln \left[{\eta }_1\left(0\right){\eta }_2\left(0\right){\eta }_3\left(0\right)\right]\\ \ge {\displaystyle {\int }_0^{t\wedge {\vartheta }_e\wedge {\vartheta }^{\{-\}}}}\left[-\left({\mu }_1+{\displaystyle \underset{i=1}{\overset{3}{\sum }}}{b}_i+{\mu }_2+\gamma +{\mu }_3+\nu \right)-\rho O\left(A,s\right)\right]\text{d}s\\ -\frac{{\sigma }^2}{2}{\displaystyle {\int }_0^{t\wedge {\vartheta }_e\wedge {\vartheta }^{\{-\}}}}\left(I^2\left(s\right)+S^2\left(s\right)\right)H^2\left(I,s\right)\text{d}s\\ +\sigma {\displaystyle {\int }_0^{t\wedge {\vartheta }_e\wedge {\vartheta }^{\{-\}}}}\left(S\left(s\right)+A\left(s\right)\right)O\left(A,s\right)\text{d}W\left(s\right)\\ =J\left(t\right).\end{array}$ (9)

Assume that $P\left(\left\{{\vartheta }^{\mathrm{<mo>\; \{\; </mo>}-\mathrm{<mo>\; \}\; </mo>}}<{\vartheta }_e\right\}\right)>0$ . By continuity of the solution of system (7), we have on the event $\left\{{\vartheta }^{\{-\}}<{\vartheta }_e\right\}$ :

$S\left({\vartheta }^{\mathrm{<mo>\; \{\; </mo>}-\mathrm{<mo>\; \}\; </mo>}}\right)A\left({\vartheta }^{\mathrm{<mo>\; \{\; </mo>}-\mathrm{<mo>\; \}\; </mo>}}\right)R\left({\vartheta }^{\mathrm{<mo>\; \{\; </mo>}-\mathrm{<mo>\; \}\; </mo>}}\right)=0.$

Hence,

$\underset{t\to {\vartheta }^{\mathrm{<mo>\; \{\; </mo>}-\mathrm{<mo>\; \}\; </mo>}}}{lim}\ln \left[S\left(t\right)A\left(t\right)R\left(t\right)\right]=-\infty \mathrm{.}$ (10)

Combining (9) and (10), we have on the event $\left\{{\vartheta }^{\{-\}}<{\vartheta }_e\right\}$ that $-\infty \ge J\left({\vartheta }^{\mathrm{<mo>\; \{\; </mo>}-\mathrm{<mo>\; \}\; </mo>}}\right)$ .

Therefore,

$\left\{{\vartheta }^{\{-\}}<{\vartheta }_e\right\}\subset \left\{-\infty \ge J\left({\vartheta }^{\{-\}}\right)\right\}.$

Since $J\left({\vartheta }^{\mathrm{<mo>\; \{\; </mo>}-\mathrm{<mo>\; \}\; </mo>}}\right)$ is finite on $\left\{{\vartheta }^{\{-\}}<{\vartheta }_e\right\}$ , we have a contradiction. So necessarily $P\left(\left\{{\vartheta }^{\mathrm{<mo>\; \{\; </mo>}-\mathrm{<mo>\; \}\; </mo>}}<{\vartheta }_e\right\}\right)=0$ and 1) is proved.

For any initial condition $\eta \in {\mathcal{P}}^3\left(\Delta \right)$ , the total size of the population $\mathcal{N}\left(t\right)$ at time $t\in \left[-{\tau }_{\max },{\vartheta }_e\right[$ is described by the equation:

$\{\begin{array}{l}\text{d}\mathcal{N}\left(t\right)=\left(\varphi -\left({\mu }_1+b_1\right)S\left(t\right)-\left({\mu }_2+b_2\right)A\left(t\right)-\left({\mu }_3+b_3\right)R\left(t\right)\right)\text{d}t,\\ \mathcal{N}\left(0\right)={\eta }_1\left(0\right)+{\eta }_2\left(0\right)+{\eta }_3\left(0\right).\end{array}$

In view of 1), for any $t\in \left[-{\tau }_{\max },{\vartheta }_e\right)$ , we see that $\min \left\{S\left(t\right),A\left(t\right),R\left(t\right)\right\}>0$ a.s.

Since ${\mu }_1+b_1\le \min \left\{{\mu }_2+b_2\mathrm{,}{\mu }_3+b_3\right\}$ , we get:

$\begin{array}{c}\text{d}\mathcal{N}\left(t\right)=[\varphi -\left({\mu }_2+b_2-{\mu }_1-b_1\right)A\left(t\right)-\left({\mu }_3+b_3-{\mu }_1-b_1\right)R\left(t\right)\\ -\left({\mu }_1+b_1\right)\mathcal{N}\left(t\right)]\text{d}t\\ =\left[-\Phi \left(t\right)+\left(\varphi -\left({\mu }_1+b_1\right)\mathcal{N}\left(t\right)\right)\right]\text{d}t\mathrm{,}\end{array}$

where $\Phi \left(t\right)=\left({\mu }_2+b_2-{\mu }_1-b_1\right)A\left(t\right)+\left({\mu }_3+b_3-{\mu }_1-b_1\right)R\left(t\right)>0$ .

Therefore, by virtue of a comparison theorem, we obtain that for any $t\in \left[\mathrm{0,}{\vartheta }_e\right)$ and $\mathcal{N}\left(0\right)\in \left(\mathrm{0,}\frac{\varphi }{{\mu }_1+b_1}\right)$ :

$\mathcal{N}\left(t\right)\le \left(\mathcal{N}\left(0\right)-\frac{\varphi }{{\mu }_1+b_1}\right){\text{e}}^{-\left({\mu }_1+b_1\right)t}+\frac{\varphi }{{\mu }_1+b_1}\text{a}\text{.s}\text{.}$

It follows that ${\sup }_{t\in \left[\mathrm{0,}{\vartheta }_e\right)}\mathcal{N}\left(t\right)\le \frac{\varphi }{{\mu }_1+b_1}$ which leads to ${\vartheta }^{\mathrm{<mo>\; \{\; </mo>}-\mathrm{<mo>\; \}\; </mo>}}\le {\vartheta }_e$ a.s. since $\max \left\{S\left(t\right),A\left(t\right),R\left(t\right)\right\}<\frac{\varphi }{{\mu }_1+b_1}$ a.s.

Hence, the solution might explode only toward $-\infty$ , which implies that ${\vartheta }_0={\vartheta }^{\{-\}}={\vartheta }_e$ a.s. Therefore, we have 2).

For $n\ge n_0$ , let set $F_n=\left]\frac{1}{n}\mathrm{,}\frac{\varphi }{{\mu }_1+b_1}\right[$ a sequence of real intervals. ${\left(F_n\right)}_{n\ge n_0}$ is increasing ( $F_n\subset F_{n+1}$ ) and converges to ${\displaystyle {\cup }_{n\ge n_0}}{F}_n=\left]0,\frac{\varphi }{{\mu }_1+b_1}\right[$ . Therefore, the sequence of stopping time ${\left({\vartheta }_n\right)}_{n\ge n_0}$ is increasing and there exists ${\vartheta }_{\infty }\in \left[\mathrm{0,}\infty \right[\cup \left\{\infty \right\}$ such that ${\lim }_{n\to \infty }{\vartheta }_n={\vartheta }_{\infty }$ . Since for any $n\ge n_0$ , $F_n\subset \left]\mathrm{0,}\frac{\varphi }{{\mu }_1+b_1}\right[$ , we have ${\vartheta }_{\infty }\le {\vartheta }_0$ . In particular, if ${\vartheta }_{\infty }=\infty$ then ${\vartheta }_{\infty }={\vartheta }_0=\infty$ .

Now, let us assume that ${\vartheta }_{\infty }<\infty$ and put $Y\left(t\right)=\min \left\{S\left(t\right),A\left(t\right),R\left(t\right)\right\}$ , for any $t\in \left[\mathrm{0,}{\vartheta }_e\right[$ . Since for all $n\ge n_0$ , $Y\left({\vartheta }_n\right)\notin F_n$ and $Y\left({\vartheta }_n\right)\to Y\left({\vartheta }_{\infty }\right)$ . It follows that $Y\left({\vartheta }_{\infty }\right)\notin \left]\mathrm{0,}\frac{\varphi }{{\mu }_1+b_1}\right[$ . Hence ${\vartheta }_{\infty }\ge {\vartheta }_0$ which gives 3).£

We complete the proof of the consistency of model (3) by establishing that the sequence of stopping time ${\vartheta }_n\mathrm{,}n\ge n_0$ converges to $+\infty$ , using stochastic calculus and absurdity reasoning.

Theorem 4. For any initial condition $\eta =\left({\eta }_1\mathrm{,}{\eta }_2\mathrm{,}{\eta }_3\right)\in {\mathcal{P}}^3\left(\Delta \right)$ , system (7) admits a unique solution $\left(S\left(t\right),A\left(t\right),R\left(t\right)\right)$ on $t\ge 0$ , and this solution remains in Δ with probability 1.

Proof. In view of Lemma 3, for any initial condition $\eta \in {\mathcal{P}}^3\left(\Delta \right)$ , system (7) has a unique local positive solution $\left(S\left(t\right),A\left(t\right),R\left(t\right)\right)$ on $t\in \left[\mathrm{0,}{\vartheta }_e\right)$ and ${\lim }_{n\to \infty }{\vartheta }_n={\vartheta }_e$ . In order to establish the existence and uniqueness of a global positive solution, it is enough to prove that ${\lim }_{n\to \infty }{\vartheta }_n=\infty$ .

Let $n_0\in {\mathbb{N}}^{\ast }$ such that $\min \left\{{\eta }_1\left(0\right)\mathrm{,}{\eta }_2\left(0\right)\mathrm{,}{\eta }_3\left(0\right)\right\}>\frac{1}{n_0}$ . Consider the function Q defined for any vector $x=\left(x_1\mathrm{,}{x}_2\mathrm{,}{x}_3\right)\in {ℝ}_{+}^3$ by:

$Q\left(x\right)=-\ln \left[\frac{\left({\mu }_1+b_1\right)x_1}{\varphi }\right]-\ln \left[\frac{\left({\mu }_1+b_1\right)x_2}{\varphi }\right]-\ln \left[\frac{\left({\mu }_1+b_1\right)x_3}{\varphi }\right]\mathrm{.}$

By virtue of Itô’s formula and $X\left(t\right)=\left(S\left(t\right),A\left(t\right),R\left(t\right)\right)$ , we get that for any $t\in \left[-{\tau }_{\max },{\vartheta }_e\right)$ :

$\begin{array}{c}\text{d}Q\left(X\left(t\right)\right)=[-\frac{\varphi }{S\left(t\right)}+{\mu }_1+b_1-\nu \frac{R\left(t\right)}{S\left(t\right)}+\rho O\left(A\mathrm{,}t\right)+\frac{{\sigma }^2}{2}{\left(A\left(t\right)O\left(A\mathrm{,}t\right)\right)}^2\\ -\gamma \frac{A\left(t\right)}{R\left(t\right)}+\left({\mu }_2+\gamma +b_2\right)+\left({\mu }_3+\nu +b_3\right)-\rho \frac{S\left(t\right)}{A\left(t\right)}O\left(A\mathrm{,}t\right)\\ +\frac{{\sigma }^2}{2}{\left(S\left(t\right)O\left(A\mathrm{,}t\right)\right)}^2]\text{d}t+\sigma O\left(A\mathrm{,}t\right)\left(S\left(t\right)-A\left(t\right)\right)\text{d}W\left(t\right)\mathrm{.}\end{array}$

Since for any $s\in \left[-{\tau }_{\max },t\wedge {\vartheta }_n\right)$ , $S\left(s\right)\mathrm{,}A\left(s\right)\mathrm{,}R\left(s\right)\in \left]\frac{1}{n}\mathrm{,}\frac{\varphi }{{\mu }_1+b_1}\right[$ a.s. Therefore, for all $n\ge n_0$ , we obtain that:

$\begin{array}{c}Q\left(X\left(t\wedge {\vartheta }_n\right)\right)\le Q\left(X\left(0\right)\right)+{\displaystyle {\int }_0^{t\wedge {\vartheta }_n}}\left({\mu }_1+{\mu }_2+\gamma +{\mu }_3+\nu +{\displaystyle \underset{i=1}{\overset{3}{\sum }}}{b}_i+\rho O\left(A,s\right)\right)\text{d}s\\ +\frac{{\sigma }^2}{2}{\displaystyle {\int }_0^{t\wedge {\vartheta }_n}}\left(A{\left(s\right)}^2+S{\left(s\right)}^2\right)O^2\left(I,s\right)\text{d}s\\ +\sigma {\displaystyle {\int }_0^{t\wedge {\vartheta }_n}}\left(S\left(s\right)-A\left(s\right)\right)O\left(A,s\right)\text{d}W\left(s\right).\end{array}$

Then in view of (C₁ - C₂), for any $s\in \left[\mathrm{0,}t\wedge {\vartheta }_k\right[$ , we have:

$O\left(A,s\right)={\displaystyle {\int }_{s-{\tau }_{\max }}^s}\Theta \left(s-z\right)\Psi \left(A\left(z\right)\right)\text{d}z\le \frac{\varphi }{{\mu }_1+b_1}{\displaystyle {\int }_0^{{\tau }_{\max }}}\Theta \left(z\right)\text{d}z=\frac{\varphi }{{\mu }_1+b_1}\text{a}\text{.s}.$ (11)

So, we get:

$Q\left(X\left(t\wedge {\vartheta }_n\right)\right)\le Q\left(X\left(0\right)\right)+C_{0}t\wedge {\vartheta }_n+\sigma {\displaystyle {\int }_0^{t\wedge {\vartheta }_n}}\left(S\left(s\right)-A\left(s\right)\right)O\left(A,s\right)\text{d}W\left(s\right)$ (12)

where $C_0=3\bar{\mu }+\gamma +\nu +{\displaystyle {\sum }_{i=1}^3}{b}_i+\rho \frac{\varphi }{{\mu }_1+b_1}+{\sigma }^2{\left(\frac{\varphi }{{\mu }_1+b_1}\right)}^4$ and $\bar{\mu }=\max \left\{{\mu }_1,{\mu }_2,{\mu }_3\right\}$ .

On the other hand, in view of (11), we have:

$\mathbb{E}\left[{\displaystyle {\int }_0^{t\wedge {\vartheta }_n}}\left(S\left(s\right)-A\left(s\right)\right)O\left(A\mathrm{,}s\right)\text{d}W\left(s\right)\right]=0.$

It follows from (12) that for any $t\ge 0$ :

$\begin{array}{c}\mathbb{E}\left[Q\left(X\left(t\wedge {\vartheta }_n\right)\right)\right]\le C_{0}t\wedge {\vartheta }_n+Q\left(X\left(0\right)\right)\\ \le C_{0}t+Q\left(X\left(0\right)\right)\mathrm{.}\end{array}$ (13)

Now, since $S\left(t\wedge {\vartheta }_n\right)\mathrm{,}A\left(t\wedge {\vartheta }_n\right)\mathrm{,}R\left(t\wedge {\vartheta }_n\right)$ are in $\left]\frac{1}{n}\mathrm{,}\frac{\varphi }{{\mu }_1+b_1}\right[$ , we have $Q\left(X\left(t\wedge {\vartheta }_n\right)\right)>0$ .

Therefore, we get:

$\begin{array}{c}\mathbb{E}\left[Q\left(X\left(t\wedge {\vartheta }_n\right)\right)\right]=\mathbb{E}\left[Q\left(X\left(t\wedge {\vartheta }_n\right)\right)1_{\left\{{\vartheta }_n\le t\right\}}\right]+\mathbb{E}\left[Q\left(X\left(t\wedge {\vartheta }_n\right)\right)1_{\left\{{\vartheta }_n>t\right\}}\right]\\ \ge \mathbb{E}\left[Q\left(X\left(t\wedge {\vartheta }_n\right)\right)1_{\left\{{\vartheta }_n\le t\right\}}\right]\mathrm{.}\end{array}$

In view of Lemma 3, for all $\eta \in {\mathcal{P}}^3\left(\Delta \right)$ , we see that:

$\min \left\{S\left(t\wedge {\vartheta }_n\right)\mathrm{,}A\left(t\wedge {\vartheta }_n\right)\mathrm{,}R\left(t\wedge {\vartheta }_n\right)\right\}<\frac{\varphi }{{\mu }_1+b_1}\mathrm{.}$

It follows that $\min \left\{S\left({\vartheta }_n\right)\mathrm{,}A\left({\vartheta }_n\right)\mathrm{,}R\left({\vartheta }_n\right)\right\}=\frac{1}{n}$ on $\left\{{\vartheta }_n\le t\right\}$ , which implies that $Q\left(X\left({\vartheta }_n\right)\right)\ge -\ln \left(\frac{{\mu }_1+b_1}{\varphi n}\right)$ .

Hence,

$\begin{array}{c}\mathbb{E}\left[Q\left(X\left(t\wedge {\vartheta }_n\right)\right)\right]\ge \mathbb{E}\left[Q\left(X\left(t\wedge {\vartheta }_n\right)\right)1_{\left\{{\vartheta }_n\le t\right\}}\right]\\ \ge -\ln \left(\frac{{\mu }_1+b_1}{\varphi n}\right)\mathbb{P}\left({\vartheta }_n\le t\right)\mathrm{.}\end{array}$ (14)

Combining (13) and (14), for any $t\ge 0$ , we get that:

$\mathbb{P}\left({\vartheta }_n\le t\right)\le \frac{C_{0}t+Q\left(X\left(0\right)\right)}{\ln \left(\frac{\varphi n}{{\mu }_1+b_1}\right)}\mathrm{.}$

By letting $n\to \infty$ , we obtain for any $t\ge 0$ , $P\left({\vartheta }_0\le t\right)=0$ . Consequently, $P\left({\vartheta }_0=\infty \right)=1$ . Now, since ${\vartheta }_e={\vartheta }_0$ a.s., we obtain that ${\vartheta }_e=\infty$ a.s.

Based on Lemma 2, we have the following result on the global positivity of model (3).

Corollary 1. For any initial condition $\eta =\left({\eta }_1\mathrm{,}{\eta }_2\mathrm{,}{\eta }_3\mathrm{,}{\eta }_4\right)\in {\mathcal{P}}^4\left(\Delta \times {ℝ}_{+}\right)$ , system (3) admits a unique solution $\left(S\left(t\right),A\left(t\right),R\left(t\right)\mathrm{,}B\left(t\right)\right)$ on $t\ge 0$ , and this solution remains in ${ℝ}_{+}^4$ with probability 1.

4. Extinction of the Disease

In mathematical modeling in epidemiology, one of the main questions is the determination of the conditions which ensure the disappearance of a disease within a population or to control its spread to a bearable level otherwise. Generally, the study of the extinction of a disease described by an epidemic model is carried out by analyzing the stability of the disease-free equilibrium point [10] [7] [39] . The basic reproduction number $R_0$ , defined as the average number of secondary cases produced by an infectious individual in a completely susceptible population, is a key indicator of whether a disease is becoming endemic in a population. The condition $R_0<1$ ensures the extinction of a disease in the deterministic framework [38] . The authors also studied the extinction of the disease described by a stochastic delay epidemic model under the condition $R_0<1$ in addition to a condition on the intensity of the noise $\sigma$ [33] [40] . But, considering only the equation describing the size of the infection in the model, Wanduku [37] establishes the extinction of the disease described by a delayed SIERS stochastic model under the sole condition $R_0<1$ as in the corresponding deterministic case.

In this section, we investigate the stability of the disease-free equilibrium

$E_0=\left(\frac{\varphi }{{\mu }_1+b_1}\mathrm{,0,0,}\frac{b_1\varphi }{{\mu }_4\left({\mu }_1+b_1\right)}\right)$ of model (3) through the stability analysis of the disease-free equilibrium $E_0^{\mathrm{<mo>\; *\; </mo>}}=\left(\frac{\varphi }{{\mu }_1+b_1}\mathrm{,0,0}\right)$ of model (7) based on Lemma 2.

We will first establish that the trivial solution of the stochastic delay differential equation describing the size of infectious individuals in model (7) is exponentially stable almost surely, i.e. ${\lim }_{t\to \infty }I\left(t\right)=0$ . From there, we deduce the almost surely stability of the disease-free equilibrium point $E_0^{\mathrm{<mo>\; *\; </mo>}}$ under the condition $R_0<1$ . This result is obtained by combined Lyapunov function technique and martingale convergence result (see, e.g. [45] [46] ). The stability of the disease-free equilibrium E₀ leads to the extinction of the disease described by model (3).

Theorem 5. Let $R_0<1$ , then the disease-free equilibrium $E_0^{\mathrm{<mo>\; *\; </mo>}}=\left(\frac{\varphi }{{\mu }_1+b_1}\mathrm{,0,0}\right)$ of model (7) is asymptotically almost surely stable for any initial condition $\eta =\left({\eta }_1\mathrm{,}{\eta }_2\mathrm{,}{\eta }_3\right)\in {\mathcal{P}}^3\left(\Delta \right)$ .

Proof of Theorem 5. We will first prove separately the almost sure asymptotic stability for every component of the solution $\left(S\left(t\right),A\left(t\right),R\left(t\right)\right)\mathrm{,}t\ge 0$ of system (7) and then conclude.

For any $\left(x_1\mathrm{,}{x}_2\mathrm{,}{x}_3\right)\in {ℝ}^3$ , let us put $P{r}_2\left(x_1\mathrm{,}{x}_2\mathrm{,}{x}_3\right)=x_2$ . It follows that the infectious size A(t) of model (7) is described by the following equation:

$\begin{array}{l}\text{d}A\left(t\right)=\left(\rho S\left(t\right)O\left(A\mathrm{,}t\right)-\left({\mu }_2+\gamma +b_2\right)A\left(t\right)\right)\text{d}t\\ +\sigma S\left(t\right)A\left(t\right)O\left(A\mathrm{,}t\right)\text{d}W\left(t\right)\mathrm{,}\end{array}$ (15)

with initial condition $P{r}_2\circ \eta ={\eta }_2$ where $\eta =\left({\eta }_1\mathrm{,}{\eta }_2\mathrm{,}{\eta }_3\right)\in {\mathcal{P}}^3\left(\Delta \right)$ .

In this equation, ${\left(S\left(t\right)\right)}_{t\ge 0}$ is considered as an adapted process and almost surely bounded by $s_0=\frac{\varphi }{{\mu }_1+b_1}$ . Let us consider the functional:

$U\left(t\mathrm{,}P{r}_2\circ \eta \right)=U_1\left(P{r}_2\circ \eta \right)+U_2\left(t\mathrm{,}P{r}_2\circ \eta \right)\mathrm{,}\text{for}\text{all}\eta \in {\mathcal{P}}^3\left(\Delta \right)\mathrm{,}$

where $U_1\left({\eta }_2\right)=\left|{\eta }_2\left(0\right)\right|$ , $U_2\left(t\mathrm{,}{\eta }_2\right)=\varpi {\displaystyle {\int }_{-{\tau }_{\max }}^0}\Theta \left(-\theta \right)\Psi \left({\eta }_2\left(\theta \right)\right)\text{d}\theta$ ,

and ${\eta }_2\left(\theta \right)=A\left(t+\theta \right)$ , $\theta \in \left[-{\tau }_{\max },0\right]$ .

In view of Theorem 4, for any ${\eta }_2$ such that $\eta =\left({\eta }_1\mathrm{,}{\eta }_2\mathrm{,}{\eta }_3\right)\in {\mathcal{P}}^3\left(\Delta \right)$ , we get that:

$\begin{array}{c}\mathcal{L}{U}_1\left({\eta }_2\right)=\rho S\left(t\right)O\left(A\mathrm{,}t\right)-\left({\mu }_2+\gamma +b_2\right)A\left(t\right)\\ \le \rho s_0{\displaystyle {\int }_0^{{\tau }_{\max }}}\Theta \left(\tau \right)\Psi \left(A\left(t-\tau \right)\right)\text{d}\tau -\left({\mu }_2+\gamma +b_2\right)A\left(t\right)\mathrm{.}\end{array}$

Let us put $\varpi =s_0\rho$ . In view of C₃, we have:

$\mathcal{L}{U}_2\left(t,{\eta }_2\right)\le s_0\rho \left[A\left(t\right)-{\displaystyle {\int }_0^{{\tau }_{\max }}}\Theta \left(\tau \right)\Psi \left(A\left(t-\tau \right)\right)\text{d}\tau \right].$

Finally, we get:

$\mathcal{L}U\left(t\mathrm{,}{\eta }_2\right)\le -\left({\mu }_2+\gamma +b_2-\rho s_0\right)\left|{\eta }_2\left(0\right)\right|\mathrm{.}$

Then, by virtue of Theorem 1, when $R_0=\frac{\rho s_0}{{\mu }_2+\gamma +b_2}<1$ , we have:

$\underset{t\to \infty }{lim}\frac{1}{t}\ln \left(A\left(t\right)\right)<-p\text{a}\text{.s}\mathrm{.,}$

where p is a positive constant. That is, there exists two positive constants $p_1$ and $p_2$ such that:

$A\left(t\right)<p_1\exp \left(-p_{2}t\right)\text{for}\text{any}t\ge 0\text{a}\text{.s}.$ (16)

Now, consider the third equation of model (7). From the well-know variation of constants approach, we obtain:

$R\left(t\right)=R\left(0\right){\text{e}}^{-\left({\mu }_3+\nu +b_3\right)t}+\gamma {\displaystyle {\int }_0^t}A\left(s\right){\text{e}}^{\left({\mu }_3+\nu \right)\left(s-t\right)}\text{d}s.$

In view of (16), for all $\epsilon >0$ , there exists $T\left(\epsilon \right)>0$ such that for any $t>T\left(\epsilon \right)$ , $A\left(t\right)<\epsilon \frac{{\mu }_3+\nu +b_3}{\gamma }$ :

$\begin{array}{c}R\left(t\right)=R\left(0\right){\text{e}}^{-\left({\mu }_3+\nu \right)t}+\gamma {\displaystyle {\int }_0^t}A\left(s\right){\text{e}}^{\left({\mu }_3+\nu +b_3\right)\left(s-t\right)}\text{d}s\\ \le \left(R\left(0\right)-\epsilon \right){\text{e}}^{-\left({\mu }_3+\nu +b_3\right)t}+\epsilon ,\end{array}$ (17)

By letting $\epsilon \to 0$ , we have:

$\underset{t\to \infty }{lim}R\left(t\right)\le \underset{t\to \infty }{lim}R\left(0\right){\text{e}}^{-\left({\mu }_3+\nu +b_3\right)t}\mathrm{.}$

Hence,

$\underset{t\to \infty }{lim}R\left(t\right)=0\text{a}\text{.s}$

Let us now prove that ${\lim }_{t\to \infty }\left(\frac{\varphi }{{\mu }_1+b_1}-S\left(t\right)\right)=0$ . From the first equation of model (7), we get:

$\begin{array}{c}\frac{\varphi }{{\mu }_1+b_1}-S\left(t\right)=\frac{\varphi }{{\mu }_1+b_1}-S\left(0\right)+\rho {\displaystyle {\int }_0^t}S\left(s\right)O\left(A,t\right)\text{d}s\\ -{\displaystyle {\int }_0^t}\left[\left({\mu }_1+b_1\right)\left(\frac{\varphi }{{\mu }_1+b_1}-S\left(s\right)\right)+\nu R\left(s\right)\right]\text{d}s\\ +\sigma {\displaystyle {\int }_0^t}S\left(s\right)A\left(s\right)O\left(A,t\right)\text{d}W\left(s\right).\end{array}$ (18)

In view of Theorem 4, Hölder inequality and (16), we obtain:

$\begin{array}{l}\underset{t\to \infty }{lim}\rho {\displaystyle {\int }_0^t}S\left(s\right){\displaystyle {\int }_0^{{\tau }_{\max }}}\Theta \left(\tau \right)\Psi \left(A\left(s-\tau \right)\right)\text{d}\tau \text{d}s\\ \le \underset{t\to \infty }{lim}\rho {\displaystyle {\int }_0^t}S\left(s\right){\displaystyle {\int }_{s-{\tau }_{\max }}^s}\Theta \left(s-u\right)A\left(u\right)\text{d}u\text{d}s\\ \le h{p}_1\rho s_0\exp \left(p_2{\tau }_{\max }\right)\left(\underset{t\to \infty }{lim}{\displaystyle {\int }_0^t}\exp \left(-C_{2}s\right)\text{d}s\right)<\infty ,\end{array}$

and

$\frac{\varphi }{{\mu }_1+b_1}-S\left(t\right)\ge 0\text{a}\text{.s}\text{.}$

Therefore, by virtue of the non-negative semimartingale convergence result established in Liptser and Shiryayev ( [47] , Theorem 7, p. 139), we get from (18):

$\underset{t\to \infty }{lim}\left(\frac{\varphi }{{\mu }_1+b_1}-S\left(t\right)\right)<\infty$ and

$\underset{t\to \infty }{lim}{\displaystyle {\int }_0^t}\left[\left({\mu }_1+b_1\right)\left(\frac{\varphi }{{\mu }_1+b_1}-S\left(s\right)\right)+\nu R\left(s\right)\right]\text{d}s<\infty \text{a}\text{.s}\mathrm{.}$

Since $R\left(t\right)$ and $\frac{\varphi }{{\mu }_1+b_1}-S\left(s\right)$ are positives for all $t\ge 0$ , we get:

$\underset{t\to \infty }{lim}{\displaystyle {\int }_0^t}\left(\frac{\varphi }{{\mu }_1+b_1}-S\left(s\right)\right)\text{d}s={\displaystyle {\int }_0^{\infty }}\left(\frac{\varphi }{{\mu }_1+b_1}-S\left(s\right)\right)\text{d}s<\infty \mathrm{.}$ (19)

Assume that $\frac{\varphi }{{\mu }_1+b_1}-S\left(s\right)$ does not converge almost surely to 0. Then there is a set ${\Omega }_1\subset \Omega$ with $P\left({\Omega }_1\right)>0$ such that for all $\omega \in {\Omega }_1$ :

$\underset{t\to \infty }{\lim \inf }\left(\frac{\varphi }{{\mu }_1+b_1}-S\left(t,\omega \right)\right)=\tau \left(\omega \right)>0.$

Then, there exists a $T>0$ such that $\frac{\varphi }{{\mu }_1+b_1}-S\left(t,\omega \right)>\frac{1}{2}\tau \left(\omega \right)$ for all $t\ge T$ . It follows that:

$\begin{array}{l}\underset{t\to \infty }{\lim }{\displaystyle {\int }_0^t}\left(\frac{\varphi }{{\mu }_1+b_1}-S\left(s,\omega \right)\right)\text{d}s\\ ={\displaystyle {\int }_0^T}\left(\frac{\varphi }{{\mu }_1+b_1}-S\left(s,\omega \right)\right)\text{d}s+{\displaystyle {\int }_T^{\infty }}\left(\frac{\varphi }{{\mu }_1+b_1}-S\left(s,\omega \right)\right)\text{d}s\\ \ge {\displaystyle {\int }_T^{\infty }}\left(\frac{\varphi }{{\mu }_1+b_1}-S\left(s,\omega \right)\right)\text{d}s=\infty .\end{array}$

Therefore, ${\Omega }_1\subset {\Omega }_2$ , where ${\Omega }_2=\left\{\omega ,{\displaystyle {\int }_T^{\infty }}\left(\frac{\varphi }{{\mu }_1+b_1}-S\left(s,\omega \right)\right)\text{d}s=\infty \right\}$ . Hence, $P\left({\Omega }_2\right)>0$ , which contradicts (19). So, we have:

$\underset{t\to \infty }{lim}\left(\frac{\varphi }{{\mu }_1+b_1}-S\left(t\right)\right)=0\text{a}\text{.s}\mathrm{.}$

Finally, we obtain that, when $t\to \infty$ , then

$\left(S\left(t\right),A\left(t\right),R\left(t\right)\right)\to \left(\frac{\varphi }{{\mu }_1+b_1}\mathrm{,0,0}\right)\text{a}\text{.s}\mathrm{.}$ £

Considering the previous result in Theorem 5 and Lemma 2, we get the following result.

Corollary 2. Under the condition $R_0<1$ , the disease-free equilibrium $E_0=\left(\frac{\varphi }{{\mu }_1+b_1}\mathrm{,0,0,}\frac{b_1\varphi }{{\mu }_4\left({\mu }_1+b_1\right)}\right)$ of model (3) is globally asymptotically almost surely stable for any initial condition $\eta =\left({\eta }_1\mathrm{,}{\eta }_2\mathrm{,}{\eta }_3\mathrm{,}{\eta }_4\right)\in {\mathcal{P}}^4\left(\Delta \times {ℝ}_{+}\right)$ .

5. Persistence When $R_0>1$

The stochastic model (3) obtained from the deterministic system (2) has a single equilibrium position which is the disease-free equilibrium E₀. However, even if it does not admit endemic equilibrium, it is interesting to understand the asymptotic behavior of the solution of the stochastic model when $R_0>1$ . In the following results, we establish the persistence of the solution of the stochastic model (7) when $R_0>1$ since the persistence of model (3) follows by Lemma 2. Recall that the solution of the stochastic model (3) is said to be persistent with probability one if, for each initial value $\eta =\left({\eta }_1\mathrm{,}{\eta }_2\mathrm{,}{\eta }_3\mathrm{,}{\eta }_4\right)\in {\mathcal{P}}^4\left(\Delta \times {ℝ}_{+}\right)$ . We have the property:

$\underset{t\to \infty }{\lim \inf }S\left(t\right)>0,\underset{t\to \infty }{\lim \inf }I\left(t\right)>0,\underset{t\to \infty }{\lim \inf }R\left(t\right)>0,\underset{t\to \infty }{\lim \inf }B\left(t\right)>0\text{a}\text{.s}.$

This persistence property translates into the endemicity of the disease described by the model in the population.

Theorem 6. Assume that $R_0>1$ . If the disease-free equilibrium $E_0^{\mathrm{<mo>\; *\; </mo>}}$ of model (7) is unstable in ${\mathcal{P}}^3\left(\Delta \right)$ , then for any initial condition

$\left({\eta }_1\mathrm{,}{\eta }_2\mathrm{,}{\eta }_3\right)\in {\mathcal{P}}^3\left(\Delta \right)$ , the solution of system (7) is persistent with probability 1, that is there exists a constant $\upsilon \in \left(\mathrm{0,}\frac{\varphi }{{\mu }_1+b_1}\right)$ such that,

$\underset{t\to \infty }{\lim \inf }S\left(t\right)\ge \upsilon ,\underset{t\to \infty }{\lim \inf }I\left(t\right)\ge \upsilon ,\underset{t\to \infty }{\lim \inf }R\left(t\right)\ge \upsilon .$

Proof. Let us assume that for any initial condition $\eta \in {\mathcal{P}}^3\left(\Delta \right)$ the disease-free equilibrium E₀ of system (7) is unstable and the trivial solution of Equation (15) describing the infectious size A(t) with initial condition ${\eta }_2=P{r}_2\left(\eta \right)$ is stable.

It follows that:

$\begin{array}{l}\forall \epsilon >\mathrm{0,}\exists \eta \left(\epsilon \right)>0\text{and}\exists T\left(\epsilon \right)>0\text{such}\text{that}\\ \Vert A\left(t\right)\Vert <\epsilon \frac{{\mu }_3+\nu }{\gamma }\mathrm{,}\forall t\ge T\left(\epsilon \right)\text{and}\Vert P{r}_2\left(\eta \right)\Vert <\eta \left(\epsilon \right)\mathrm{.}\end{array}$

From (17), we obtain:

$R\left(t\right)\le \left(R\left(0\right)-\epsilon \right){\text{e}}^{-\left({\mu }_3+\nu \right)t}+\epsilon \mathrm{,}$

By letting $\epsilon \to 0$ , for any $t\ge 0$ , we have:

${\lim }_{t\to \infty }R\left(t\right)\le {\lim }_{t\to \infty }R\left(0\right){\text{e}}^{-\left({\mu }_3+\nu \right)t}$ .

Hence,

$\underset{t\to \infty }{\lim }R\left(t\right)=\underset{t\to \infty }{\lim }A\left(t\right)=0\text{a}\text{.s}.$ (20)

Since ${\mu }_1\le \min \left\{{\mu }_2\mathrm{,}{\mu }_3\right\}$ , the size of the whole population $\mathcal{N}\left(t\right)$ in model (7) verifies:

$\mathcal{N}\left(t\right)={\text{e}}^{-\left({\mu }_1+b_1\right)t}\left(\mathcal{N}\left(0\right)+{\displaystyle {\int }_0^t}{\text{e}}^{\left({\mu }_1+b_1\right)s}\left(\varphi -{\alpha }_{1}A\left(s\right)+{\alpha }_{2}R\left(s\right)\right)\text{d}s\right)\mathrm{,}$

where ${\alpha }_1={\mu }_2-{\mu }_1$ and ${\alpha }_2={\mu }_3-{\mu }_1$ .

In view of (20) for any $\epsilon >0$ , there exists $T\left(\epsilon \right)>0$ such that for any $t>T\left(\epsilon \right)$ , $A\left(t\right)<\epsilon$ and $R\left(t\right)<\epsilon$ , we have:

$\mathcal{N}\left(t\right)\ge \mathcal{N}\left(0\right){\text{e}}^{-\left({\mu }_1+b_1\right)t}+\frac{\varphi -\left({\alpha }_1+{\alpha }_2\right)\epsilon }{{\mu }_1+b_1}\left(1-{\text{e}}^{-\left({\mu }_1+b_1\right)t}\right)\mathrm{.}$

By letting $\epsilon \to 0$ , we obtain:

$\underset{t\to \infty }{\lim \inf }\mathcal{N}\left(t\right)\ge \frac{\varphi }{{\mu }_1+b_1}\text{a}\text{.s}\text{.}\forall t\mathrm{<mo>\; >\; </mo>}T(\; \epsilon \; )$

Since $\mathcal{N}\left(t\right)=S\left(t\right)+A\left(t\right)+R\left(t\right)$ , by virtue of Lemma 3 and (20), we have:

$\underset{t\to \infty }{\lim }\mathcal{N}\left(t\right)=\frac{\varphi }{{\mu }_1+b_1}\text{and}\underset{t\to \infty }{\lim }S\left(t\right)=\frac{\varphi }{{\mu }_1+b_1}\text{a}\text{.s}\text{.}$

So, the disease-free equilibrium E₀ is stable, which is a contradiction since by hypothesis the disease-free equilibrium E₀ is assumed to be unstable. Therefore, the trivial solution of Equation (15) describing the infectious size A(t) is unstable. Finally, there exists a constant $\upsilon >0$ such that:

$\underset{t\to \infty }{\lim \inf }S\left(t\right)>\upsilon ,\underset{t\to \infty }{\lim \inf }A\left(t\right)>\upsilon ,\underset{t\to \infty }{\lim \inf }R\left(t\right)>\upsilon .$ £

6. Numerical Simulation and Discussion

Her, we propose some numerical simulations to understand the results on the extinction and persistence of the disease described by model (3). The main objective being to understand the effect of individual protective parameters b₁, b₂ and b₃ against the disease on the extinction of the disease in the population as well as the effect of the magnitude of the basic reproduction number R₀ on the level of disease endemicity. In a first scenario, we give a simulation of a sample of paths of model (3) in the case $R_0<1$ and with high values of the isolation parameters. In a second scenario, we lower the values of the isolation protection parameters and we actually observe an increase in the basic reproduction number, which generates endemicity of the disease with $R_0>1$ . Finally, from the previous simulation with $R_0>1$ , we increase the value of the random noise intensity and conclude. We use the Euler-Maruyama method (see, e.g. [48] ) to simulate the path of model (3). The parameters used in the simulation are all identical except the for the isolation protection parameters. The following parameters are used in all simulations presented below:

· The nonlinear function in the incidence rate is $\Psi \left(x\right)=x/\left(1+x\right)$ for all $x\in \left[\mathrm{0,}\infty \right)$ .

· The latency period $T_{\tau }$ is distributed according to a uniform law $\mathcal{U}\left(\left[0,{\tau }_{\max }\right]\right)$ where ${\tau }_{\max }=15$ and $\mathbb{E}\left(T_{\tau }\right)=7.5$ . That is $\Theta \left(s\right)=1/{\tau }_{\max }$ for all $s\in \left[0,{\tau }_{\max }\right]$ , null otherwise.

· The initial values are: $S\left(\theta \right)=50$ , $A\left(\theta \right)=8$ , $R\left(\theta \right)=0$ , $B\left(\theta \right)=0$ for $\theta \in \left[-\mathrm{15,0}\right]$ .

· The birth and death rate in mosquito population is ${\mu }_v=0.765$ , therefore ${\epsilon }_v=\frac{{\text{e}}^{-{\mu }_v\mathbb{E}\left(T_{\tau }\right)}}{{\mu }_v}=0.00421$ .

· The density of mosquitoes per human is $mv=350$ and the contact rate from infectious humans to susceptible vectors is ${\beta }_h=0.02$ . Therefore, $\rho =m_v{\epsilon }_v{\beta }_h=0.0295$ .

· The recruitment rate of susceptible humans $\varphi =20$ , the mortality rates of humans are ${\mu }_1={\mu }_3={\mu }_4=0.0095$ , ${\mu }_2=0.0115$ .

· The recovery rate and the immunity loss rate are respectively $\gamma =0.1$ and $\nu =0.01$ . The noise intensity coefficient is $\sigma =0.01$ .

Scenario 1

In Figure 2, we give a sample path of the stochastic epidemic model (3) under the conditions $R_0=0.8891<1$ with high level protection rates $b_1=0.72$ , $b_2=0.78$ , $b_3=0.78$ . We see that this numerical simulation agrees with the analytical results of Theorem 5, that is the condition $R_0=0.8891<1$ is sufficient to ensure the asymptotic stability of the disease free equilibrium E₀.

Scenario 2

In Figure 3, we give a sample path of the stochastic model (3) under the condition $R_0=2.7604>1$ due to a reduction in protection rates against disease to $b_1=0.4$ , $b_2=0.4$ , $b_3=0.4$ . In this case, we see that the solution $\left(S\left(t\right),A\left(t\right),R\left(t\right)\right)$ of the model is persistent, that is $S\left(t\right)>0$ , $A\left(t\right)>0$ , $R\left(t\right)>0$ , $B>0$ .

Scenario 3

In this case, we give an example trajectory of the stochastic epidemic model (3) with a high basic reproduction number. Figure 4 in which $R_0=24.9663$ presents a higher endemicity level of the disease than the case of Figure 3 where $R_0=2.7604$ .

Figure 2. Sample paths of the stochastic SIRS epidemic model (7) with initial values are: $S\left(\theta \right)=50$ , $A\left(\theta \right)=8$ , $R\left(\theta \right)=0$ , $B\left(\theta \right)=0$ for $\theta \in \left[-\mathrm{15,0}\right]$ . ${\epsilon }_v=\frac{{\text{e}}^{-{\mu }_v\mathbb{E}\left(T_{\tau }\right)}}{{\mu }_v}=0.00421$ and $\rho =0.0064$ . The remaining parameters are given by: $\varphi =20$ , ${\mu }_1={\mu }_3=0.0095$ , ${\mu }_2=0.0115$ , $b_1=0.72$ , $b_2=0.78$ , $b_3=0.78$ , $\gamma =0.1$ , $\sigma =0.01$ , $\nu =0.01$ . The condition $R_0=0.8891<1$ of Theorem 5 is checked.

Figure 3. Sample paths of the stochastic SIRS epidemic model (3) under the condition of $R_0=1.8444>1$ . The protection rates are $b_1=0.4$ , $b_2=0.4$ , $b_3=0.4$ . The rest of the parameters are as in Figure 2 and $R_0=2.7604$ .

Figure 4. Sample paths of the stochastic SIRS epidemic model (3) under the condition of $R_0=24.9663$ . Here, the protection rates are $b_1=0.1$ , $b_2=0.1$ , $b_3=0.1$ . The rest of the parameters are as in Figure 2.

7. Conclusion and Perspective

In this work, we consider a stochastic delay differential equation representing a stochastic model (3) describing a mosquito-borne disease in a randomly varying environment where insecticides and mosquito nets are used. First, we proved the global positivity of the solution (see Corollary 1). Using a Lyapunov functional technique, we established the almost sure stability of the disease-free equilibrium E₀ of the stochastic model (3) under the condition $R_0<1$ (see Corollary 2). In Theorem 6, under the condition $R_0>1$ , we proved the persistence of the stochastic model solution. As shown in the numerical simulations in Figure 2 and Figure 3, the disease can be fully controlled by only acting on the protection rate parameters. In terms of perspective, model (3) can be improved by allowing variable parameters, which make it possible to take into account the effects of seasonal variations on the model. On the other hand, the increase in certain parameters such as the recovery rate $\gamma$ or the decrease in certain parameters such as the disease contact rate $\beta$ is accompanied by a certain cost due to financial efforts necessary for such action. It would be more realistic to take these costs into account by writing for example $\beta =\beta \left(cost\right)$ . In this way, we can carry out optimization to select the parameter values that will enable us to control the disease and minimize the cost of control.

Acknowledgements

Sincere thanks to the members of JAMP for their professional performance, and special thanks to managing editor Hellen XU for a rare attitude of high quality.

Conflicts of Interest

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

References

[1]	Ludwing, D. (1974) Stochastic Population Theories. Springer, Heidelberg.
[2]	Capasso, V. and Bakstein, D. (2005) An Introduction to Continuous-Time Stochastic Processes: Theory, Models, and Applications to Finance, Biology, and Medicine. Birkhäuser, New York.
[3]	Black, F. and Scholes, M. (1976) The Pricing of Options and Corporate Liabilities. Journal of Political Economy, 81, 637-654. https://doi.org/10.1086/260062
[4]	Arriojas, M., Hu, Y., Salah-Eldin Mohammed, S.E. and Pap, G. (2007) A Delayed Black and Scholes Formula. Stochastic Analysis and Applications, 25, 471-492. https://doi.org/10.1080/07362990601139669
[5]	Kolmanovskii, V. and Myshkis, A. (1992) Applied Theory of Functional Differential Equations. Kluwer Academic Publishers Group, Dordrecht. https://doi.org/10.1007/978-94-015-8084-7
[6]	Kolmanovskii, V. and Shaikhet, L. (1996) Control of Systems with Aftereffect. American Mathematical Society, Providence.
[7]	Enatsu, Y., Muroya, Y. and Nakata, Y. (2012) Lyapunov Functional Techniques for the Global Stability Analysis of a Delayed SIRS Epidemic Model. Nonlinear Analysis: Real World Applications, 13, 2120-2133. https://doi.org/10.1016/j.nonrwa.2012.01.007
[8]	Centers for Disease Control and Prevention (2012) Lesson 1, Section 9: Natural History and Spectrum of Disease. https://archive.cdc.gov/www_cdc_gov/csels/dsepd/ss1978/lesson1/section9.html
[9]	Mao, X. (2008) Stochastic Differential Equations and Applications. 2nd Edition, Horwood Publications, Chichester. https://doi.org/10.1533/9780857099402
[10]	Beretta, E., Kolmanovski, V. and Shaikhet, L. (1998) Stability of Epidemic Model with Time Delays Influenced by Stochastic Perturbations. Mathematic and Computers in Simulation, 45, 269-277. https://doi.org/10.1016/S0378-4754(97)00106-7
[11]	Kermack, W.O. and McKendrick, A.G. (1993) Contributions to the Matuhematical Theory of Epidemics (Part III), Further Studies of the Problem of Endemicity. Proceedings of the Royal Society London. Series A, 141, 94-122. https://doi.org/10.1098/rspa.1933.0106
[12]	Capasso, V. and Serio, G. (1968) A Generalization of the Kermack-Mckendrick Deterministic Epidemic Model. Mathematical Biosciences, 42, 43-61. https://doi.org/10.1016/0025-5564(78)90006-8
[13]	Xu, R. and Ma, Z. (2009) Stability of a Delayed SIRS Epidemic Model with a Nonlinear Incidence Rate. Chaos Solitons and Fractals, 41, 2319-2325. https://doi.org/10.1016/j.chaos.2008.09.007
[14]	Hethcote, H.W. and van den Driessche, P. (1991) Some Epidemiological Models with Nonlinear Incidence. Journal of Mathematical Biology, 29, 271-287. https://doi.org/10.1007/BF00160539
[15]	WHO (2020) Vector-Borne Diseases. https://www.who.int/news-room/fact-sheets/detail/vector-borne-diseases
[16]	Ross, R. (1911) Some Quantitative Studies in Epidemiology. Nature, 87, 466-467. https://doi.org/10.1038/087466a0
[17]	MacDonald, G. (1957) The Epidemiology and Control of Malaria. Oxford University Press, Oxford.
[18]	Smith, D.L., Barlle, K.E., Barker, C.M., et al. (2010) Ross, Macdonald, and a Theory for the Dynamics and Control of Mosquito—Transmitted Pathogens. PLOS Pathogens, 8, e1002588. https://doi.org/10.1371/journal.ppat.1002588
[19]	Dietz, K., Molineaux, L. and Thomas, A. (1974) A Malaria Model Tested in the African Savannah. Bulletin of the World Health Organization, 50, 347-357. https://www.ncbi.nlm.nih.gov/pmc/articles/PMC2481186
[20]	Li, J., Welch, R.M., Nair, U.S., et al. (2002) Dynamic Malaria Models with Environmental Changes. Proceedings Thirty-Fourth Southeastern Symposium on System Theory, Huntsville, 19 March 2002, 396-400.
[21]	Yang, H.M. (2000) Malaria Transmission Model for Different Levels of Acquired Immunity and Temperature-Dependent Parameters (Vector). Revista de Saúde Pública, 34, 223-231. https://doi.org/10.1590/S0034-89102000000300003
[22]	Bacäer, N. and Sokhna, C. (2005) A Reaction-Diffusion System Modeling the Spread of Resistance to an Antimalarial Drug. Mathematical Biosciences and Engineering, 2, 227-238. https://doi.org/10.3934/mbe.2005.2.227
[23]	Koella, J.C. and Antia, R. (2003) Epidemiological Models for the Spread of Anti-Malarial Resistance. Malaria Journal, 2, Article No. 3. https://doi.org/10.1186/1475-2875-2-3
[24]	Chiyaka, C., Tchuenche, J.M., Garira, W. and Dube, S. (2008) A Mathematical Analysis of the Effects of Control Strategies on the Transmission Dynamics of Malaria. Applied Mathematics and Computation, 195, 641-662. https://doi.org/10.1016/j.amc.2007.05.016
[25]	Ruan, S., Xiao, D. and Beier, J.C. (2008) On the Delayed Ross-Macdonald Model for Malaria Transmission. Bulletin of Mathematical Biology, 70, 1098-1114. https://doi.org/10.1007/s11538-007-9292-z
[26]	Demers, J., Bewick, S. and Fagan, W.F. (2018) Dynamic Modelling of Personal Protection Control Strategies for Vector-Borne Disease Limits the Role of Diversity Amplification. Journal of the Royal Society Interface, 15, Article ID: 20180116. https://doi.org/10.1098/rsif.2018.0166
[27]	Dordević, J.A. (2022) A Stochastic Model for Malaria and Its Behavior under Insecticide-Treated Nets. Studies in Applied Mathematics, 149, 631-656. https://doi.org/10.1111/sapm.12515
[28]	Tchoumi, S.Y., Kouakep-Tchaptchie, Y., Fotsa-Mbogne, D.J., et al. (2021) Optimal Control of a Malaria Model with Long-Lasting Insecticide-Treated Nets. Mathematical Modelling and Control, 1, 188-207. https://doi.org/10.3934/mmc.2021018
[29]	Favier, C., Schmit, D., Muller-Graf, C.D.M., et al. (2005) Influence of Spatial Heterogeneity on an Emerging Infectious Disease: The Case of Dengue Epidemics. Proceedings of the Royal Society B: Biological Sciences, 272, 1171-1177. https://doi.org/10.1098/rspb.2004.3020
[30]	Dushoff, J.G. (1997) Modeling the Effects of Host Heterogeneity on the Spread of Human Diseases. Ph.D. Thesis, Princeton University, Princeton.
[31]	Alitzer, S., Dobson, A. and Hosseini, P. (2006) Seasonality and the Dynamics of Infectious Diseases. Ecology Letters, 9, 467-484. https://doi.org/10.1111/j.1461-0248.2005.00879.x
[32]	Tornatore, E. and Buccellato, S.M. (2009) Parasite Population Delay Model of Malaria Type with Stochastic Perturbation and Environmental Criterion for Limitation of Disease. Journal of Mathematical Analysis and Applications, 360, Article ID: 624630. https://doi.org/10.1016/j.jmaa.2009.06.078
[33]	Krstic, M. (2011) The Effect of Stochastic Perturbation on a Nonlinear Delay Malaria Epidemic Model. Mathematics and Computers in Simulation, 82, 558-569. https://doi.org/10.1016/j.matcom.2011.09.003
[34]	Agusto, F.B., Marcus, N. and Okosun, K.O. (2012) Application of Optimal Control to the Epidemiology of Malaria. Electronic Journal of Differential Equations, 2012, 1-22.
[35]	Takeuchi, Y., Ma, W. and Beretta, E. (2000) Global Asymptotic Properties of a Delay SIR Epidemic Model with Finite Incubation Times. Nonlinear Analysis: Theory, Methods & Applications, 42, 931-947. https://doi.org/10.1016/S0362-546X(99)00138-8
[36]	Zhen, J., Ma, Z. and Han, M. (2006) Global Stability of an SIRS Epidemic Model with Delays. Acta Mathematica Scientia, 26, 291-306. https://doi.org/10.1016/S0252-9602(06)60051-9
[37]	Wanduku, D. (2019) The Stochastic Extinction and Stability Conditions for Nonlinear Malaria Epidemics. Mathematical Biosciences and Engineering, 16, 3771-3806. https://doi.org/10.3934/mbe.2019187
[38]	Kamgang, J.C. and Gauthier, S. (2008) Computation of Threshold Conditions for Epidemiological Models and Global Stability of the Disease-Free Equilibrium (DFE). Mathematical Biosciences, 213, 1-12. https://doi.org/10.1016/j.mbs.2008.02.005
[39]	Lui, Q., Jiang, D., Hayat, T. and Ahmad, B. (2018) Analysis of a Delayed Vaccinated SIR Epidemic Model with Temporary Immunity and Lévy Jumps. Nonlinear Analysis: Hybrid Systems, 27, 29-43. https://doi.org/10.1016/j.nahs.2017.08.002
[40]	Kiouach, D. and Sabbar, Y. (2022) The Long-Time Behavior of a Stochastic SIR Epidemic Model with Distributed Delay and Multidimensional Levy Jumps. International Journal of Biomathematics, 15, Article ID: 2250004. https://doi.org/10.1142/S1793524522500048
[41]	Gray, A., Greenhalgh, D., Hu, L., Mao, X. and Pan, J. (2011) A Stochastic Differential Equation SIS Epidemic Model. SIAM Journal on Applied Mathematics, 71, 876-902. https://doi.org/10.1137/10081856X
[42]	Mao, X., Yuan, C. and Zou, J. (2005) Stochastic Differential Delay Equations in Population Dynamics. Journal of Mathematical Analysis and Applications, 304, 296-320. https://doi.org/10.1016/j.jmaa.2004.09.027
[43]	Shaikhet, L. (2013) Lyapunov Functionals and Stability of Stochastic Functional Differential Equations. Springer-Verlag, Berlin. https://doi.org/10.1007/978-3-319-00101-2
[44]	Luo, Q., Mao, X. and Shen, Y. (2011) Generalised Theory on Asymptotic Stability and Boundedness of Stochastic Functional Differential Equations. Automatica, 47, 2075-2081. https://doi.org/10.1016/j.automatica.2011.06.014
[45]	Lahrouz, A., Omari, L. and Kiouach, D. (2018) 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
[46]	Dalal, N., Greenhalgh, D. and Mao, X. (2007) A Stochastic Model of AIDS and Condom Use. Journal of Mathematical Analysis and Applications, 325, 36-53. https://doi.org/10.1016/j.jmaa.2006.01.055
[47]	Liptser, R.S. and Shiryayev, A.N. (1989) Theory of Martingales. Kluwer Academic Publishers, Dordrecht. https://doi.org/10.1007/978-94-009-2438-3
[48]	Mao, X. (2003) Numerical Solution of Stochastic Functional Differential Equations. LMS Journal of Computation and Mathematics, 6, 141-161. https://doi.org/10.1112/S1461157000000425