Stability of a Delayed Stochastic Epidemic COVID-19 Model with Vaccination and with Differential Susceptibility

Modeste N’zi¹, Boubacar Sidiki Kouyaté^2*, Ilimidi Yattara², Modibo Diarra²

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

DOI: 10.4236/jamp.2024.122034                

Abstract

In this paper, we treat the spread of COVID-19 using a delayed stochastic SVIRS (Susceptible, Infected, Recovered, Susceptible) epidemic model with a general incidence rate and differential susceptibility. We start with a deterministic model, then add random perturbations on the contact rate using white noise to obtain a stochastic model. We first show that the delayed stochastic differential equation that describes the model has a unique global positive solution for any positive initial value. Under the condition R₀ ≤ 1, we prove the almost sure asymptotic stability of the disease-free equilibrium of the model.

Keywords

SIRS Delayed Epidemic Model, Nonlinear Incidence rate, Lyapunov Function, Asymptotic Stability in Probability

1. Introduction

Coronavirus disease 2019 (COVID-19) is a contagious disease caused by a virus, SARS-CoV-2. The first known case was identified in the city of Wuhan, China in December 2019. The disease has spread around the world, leading to the COVID-19 pandemic. The World Health Organization (WHO) declared the outbreak a public health emergency of international concern on January 30, 2020, and a pandemic on March 11, 2020. As of June 25, 2022, the pandemic has caused over 543 million cases and 6.32 million confirmed deaths, making it one of the deadliest in history.

In an attempt to control the COVID-19 pandemic and its consequences, many countries resorted to public health measures such as isolation, containment, and barrier measures that had disastrous long-term economic consequences [1] [2] [3] . Thus, vaccination has remained the most appropriate means of controlling a pandemic, particularly COVID-19 [4] [5] [6] . Since the first COVID-19 vaccines became available in 2021, several countries have implemented vaccination campaigns focusing on high-risk groups such as the elderly at high risk of exposure [7] [8] . Since, studies have shown that the severity of the disease increases with age and that the risk of developing symptoms increases by 4% per year in adults between 30 and 60 years of age (see [9] ). Studies on the efficacy of COVID-19 vaccination suggest primarily increased protection against severe cases rather than protection against infection. Thus, a vaccinated person may well be infected with COVID-19 but will only develop a severe form with a very low probability (see e.g. [7] ). What’s more, in difficult economic environments such as those in developing countries, vaccines may not be accessible to everyone. It is therefore highly appropriate to implement vaccination strategies targeting high-risk groups, such as the elderly.

Mathematical modeling has played an important role in controlling an epidemic in a population [10] [11] [12] . At the beginning of the COVID-19 pandemic, authors proposed mathematical models to study the spread of SARS-CoV-2 in some regions of the world [13] [14] [15] [16] . In [17] , it is reported that individuals infected with Sars-cov2 usually develop symptoms on average 5 - 6 days after infection. The median time to symptom onset for SARS-CoV-2 is estimated to be 3 days, the shortest 1 day and the longest 24 days [18] . In view of [19] , a person infected with Sars-cov2 has an average of 5.5 (95% CI: 5.1 - 5.9) days of latency period, which corresponds to the time between infection and when it becomes infectious. Thus, the average incubation period is a good approximation of the average latency period, with a difference of about one to two days. In the literature, many authors have used delay stochastic differential equations to model the dynamics of the spread of an infectious disease with a latency or incubation period (see e.g. [20] [21] [22] ). In [23] , the authors used a stochastic delayed model to study a deconfinement strategy in Morocco. In [24] , a stochastic epidemic SIRC (Susceptible-Infectious-Recovered-Cross immune) model with delay was proposed to analyze the effect of cross-immunity in the spread of COVID-19. In [25] , the authors studied the effectiveness of quarantine using a stochastic delayed SIAQR (Susceptible-Infectiouse-Asymptomatic-Quarantine-Recovered) model. In [26] [27] , the authors used stochastic SEIRV and SIRV models that account for vaccination in the transmission of COVID-19 to study its spread.

In deterministic framework, epidemics are commonly studied by using deterministic compartmental models where the population is divided into several classes, namely susceptible, infected, and recovered groups. For example in [28] , a delayed SIRS epidemic model is proposed, the authors use a nonlinear and functional incidence rate $\beta S\left(t\right){\displaystyle {\int }_0^{h}g\left(\tau \right)\mathcal{K}\left(I\left(t-\tau \right)\right)\text{d}\tau }$ where $S\left(t\right)$ is the number of susceptible at time t, $I\left(t-\tau \right)$ the number of infected individual at time $t-\tau$ , g a probability density on $\left[\mathrm{0,}h\right]$ and $\mathcal{K}$ a function from ${ℝ}^{+}$ to ${ℝ}^{+}$ and satisfying certain assumptions. This class of incidence rate is used to model the spread of a disease in which the transmission of infection occurs through vectors that have an incubation time $\tau$ to become infectious [29] [30] . This incidence rate stipulates that the susceptible at time t can only be infected by infectious people who have been infected at time $t-\tau$ , $\tau \in \left[\mathrm{0,}h\right]$ . Then, other authors have considered this type of incidence rate as a generalization of the standard bilinear incidence rate and which can be used to model diseases that are not vector-borne (see e.g. [31] [32] [33] ). In the stochastic framework, some authors have used this type of incidence rate with a random perturbation modeled by Brownian motion. In particular, in [26] [27] the authors use a delayed stochastic model with this class of incidence rate to study the spread of COVID-19.

In the transmission process of some infectious diseases, the susceptibility to infection differs between groups of people, for example, age groups, immunological status, or fragile subgroups such as pregnant women for malaria. In the deterministic framework, many researchers have considered an epidemic model in which the susceptibility varies from individual to individual [34] [35] [36] [37] . In the reality of the spread of the COVID-19 disease, the susceptibility to infection of individuals and the burden depends on the age group. In [38] , the authors compiled several serological studies from around the world and deduced the following observations: young adults under the age of 35 had the highest seroprevalence of almost any age group; infection rates in people over the age of 55 were significantly lower than in people aged 18 - 54; the highest infection rates in New York State were in people aged 45 - 54. They also note that the age group for which the seroprevalence estimate is highest varies according to location. An analysis of the number of COVID-19 cases in Mali as of June 5, 2022 [39] confirms the same trend. In Mali, the age group under 34 years represented up to 44.7% of confirmed cases, the age group 35 - 54 years represented 36.1% of confirmed cases while the age group Age 55 and over represented 19.2% of confirmed cases. Likewise in Mali, the risk of hospitalization or death following infection by SARS-CoV-2 increases significantly with age. Recently, many studies on COVID-19 transmission have emphasized the heterogeneity in the number of cases and the severity across age groups [9] [40] [41] .

In the current literature, to our knowledge, there is no stochastic delayed SIRS model with vaccination that takes into account this differential susceptibility aspect for COVID-19. This work therefore allows the development of a new stochastic delayed SIRS with vaccination subgroup epidemic model taking into account not only the differential susceptibility but also a general transmission rate for the spread of COVID-19 and which can be used by public health authorities to adopt control strategies for COVID-19 and any other similar epidemic. In this work, we first propose a deterministic model describing the spread of COVID-19 under the hypothesis of differential susceptibility according to age groups and recourse to vaccination of the oldest subgroup. More precisely, we assume that the population is subdivided into age classes (1 - 35, 36 - 54 and >54). In order to take into account the effect of random variations in the environment on the contact process, we add a random disturbance in the contact rate of the deterministic model. We thus obtain a stochastic epidemic model described by a stochastic differential equation with delay. To ensure that the model is well-posed, we first prove the existence and uniqueness of a positive global solution. Based on the Lyapunov technique combined with stochastic analysis, we establish disease extinction below the threshold $R_0<1$ , where $R_0$ is the basic reproduction number of the deterministic model. Finally, numerical simulations are carried out to illustrate the theoretical results in a practical context.

The remainder of this work is structured as follows. In Section 2 we describe the model and Section 3 presents some definitions and notation. In Section 4 we study the consistency of the model, i.e. the well-posedness. In Section 5 we analyze the almost stability of the disease-free-equilibrium state of the model and in Section 6 we illustrate our theoretical results with numerical simulations. Finally, in section 7 we conclude and present some perspectives.

2. Model Formulation

In this work, we propose a stochastic SVIRS (susceptible-vaccinated-infectious-retired-susceptible) epidemic model with delay and different types of susceptible individuals in which the incidence rate is

$\beta S\left(t\right){\displaystyle {\int }_0^{h}g\left(\tau \right)\mathcal{K}\left(I\left(t-\tau \right)\right)\text{d}\tau }$ .

2.1. Deterministic Model

⁴As mentioned above, the severity and the transmission of COVID-19 to a susceptible individual by an infectious individual depends on several factors, for instance, the behaviour of susceptible individuals and the age groups. In the previous section, we explained that vaccination does not protect against SARS-CoV-2 infection but does protect against severe forms of the disease [7] [8] . It is therefore entirely appropriate to implement vaccination strategies targeting high-risk groups such as the elderly. On the other hand, studies carried out in [9] [38] show us that the severity and probability of catching COVID-19 depends on the age group. In some cases, such as Mali, the variability in susceptibility according to age group (1 - 35, 36 - 54 and >54) was relevant. To take into account some of these specificities, we then assume that the entire population $N\left(t\right)$ at time t is divided into six compartments, namely susceptible individuals $S_1\left(t\right)$ , $S_2\left(t\right)$ , $S_3\left(t\right)$ , vaccinated individuals $V\left(t\right)$ , infected individuals $I\left(t\right)$ and recovered individuals $R\left(t\right)$ .

2.1.1. Model Assumption

1) The population is subdivided into three age groups represented by the compartments S₁, S₂, S₃ and defined as follows:

S₁: the sub-population of susceptibles who are 35 years old or younger;

S₂: the sub-population of susceptibles over 35 and under 55 years old;

S₃: the sub-population of susceptibles 55 years or older.

2) Only the susceptible sub-population aged 55 or older are vaccinated and denoted by V. Moreover, the vaccinated individuals develop immunity related to the vaccination and move into the R compartment of recovered individuals.

3) Births occur only in the susceptible class S₁ at the rate Λ, since the newborns are less than one year old.

4) The functions $\mathcal{K}$ and g satisfy the following assumptions A₁-A₃:

(A₁) $\mathcal{K}$ is Lipschitz continuous on $\left[\mathrm{0,}+\infty \right)$ and satisfies $0<\mathcal{K}\left(x\right)\le x$ , $\forall x\ge 0$ .

(A₂) $\mathcal{K}$ is monotone increasing on $\left[\mathrm{0,}+\infty \right)$ , with $\mathcal{K}\left(0\right)=0$

(A₃) g is a probability density function with support $\left[\mathrm{0,}h\right]$ .

2.2.2. Parameter Description and the Model Chart Flow

In this model, the births occur only in the susceptible class S₁ at the rate Λ. The parameters ${\mu }_1$ , ${\mu }_2$ , ${\mu }_3$ , ${\mu }_4$ , ${\mu }_I$ and ${\mu }_R$ are respectively the mortality rates in the sub-populations S₁, S₂ and S₃ of susceptible individuals and those in the compartments V, I and R. ${\beta }_1$ , ${\beta }_2$ , ${\beta }_3$ , ${\beta }_4$ are the contact coefficients (disease transmission rate). $\gamma$ is the rate of recovery of infectious and the incubation period $\tau$ of the disease is assumed to be distributed in $\left[\mathrm{0,}h\right]$ . The transfer rate from compartment S₁ to compartment S₂ and from compartment S₂ to compartment S₃ are ${\theta }_1$ and ${\theta }_2$ , respectively. ${\theta }_3$ is the vaccination rate of those in compartment S₃, i.e., those who are older than 55 years. The parameter ${\theta }_4$ is the rate at which the vaccinated individuals develop immunity related to the vaccination. $\eta ={\displaystyle {\sum }_{i=1}^3{p}_i\eta }$ is the rate of loss of immunity of the recovered individuals while $p_i\eta ,i=1,2,3$ are the transfer rates from the compartment of recovered individuals to the compartments of susceptible individuals S₁, S₂, and S₃ respectively. The model is described by the following flowchart.

2.2. Stochastic Model

Throughout this work, we let $\left(\Omega ,\mathcal{F},{\left\{{\mathcal{F}}_t\right\}}_{t>0},P\right)$ be a complete probability space with a filtration ${\left\{{\mathcal{F}}_t\right\}}_{t>0}$ satisfying the usual conditions (i.e., ${\mathcal{F}}_0$ contains all P-null sets of $\mathcal{F}$ and ${\mathcal{F}}_{t_{+}}:={\displaystyle {\cap }_{s>t}{\mathcal{F}}_s}={\mathcal{F}}_t$ ) and we let ${\left\{W\left(t\right)\right\}}_{t\ge 0}$ be a scalar Brownian motion defined on the probability space. To obtain the stochasticity of the model, we integrate random fluctuations in the form of white noise to reveal the effect of random environmental perturbations on the parameters. Previously, we assume that $i=1,2,3,4$ each contact rate ${\tilde{\beta }}_i$ is given as a random variable with average value ${\beta }_i$ plus some random fluctuation term ${\epsilon }_i$ with a mean zero. So, in the small time interval $\left[t\mathrm{,}t+\text{d}t\right]$ , each infected individual makes ${\tilde{\beta }}_i\text{d}t={\beta }_i\text{d}t+{\epsilon }_i\text{d}t$ potentially infectious contact with susceptible individuals. Based on the same argument as in [42] , we fund that ${\epsilon }_i\text{d}t~N\left(\mathrm{0,}{\mathcal{D}}_i^2\text{d}t\right)$ or ${\epsilon }_i\text{d}t~{\mathcal{D}}_i\text{d}W\left(t\right)$ , where $\text{d}W\left(t\right)=W\left(t+\text{d}t\right)-W\left(t\right)$ is the increment of a standard scalar Brownian motion ${\left(W\left(t\right)\right)}_{t\ge 0}$ that follows $N\left(\mathrm{0,}\text{d}t\right)$ . It follows that, ${\tilde{\beta }}_i\left(t\right)\text{d}t\mathrm{<mo>\; :\; </mo>}={\beta }_i\text{d}t+{\mathcal{D}}_i\text{d}W\left(t\right)$ where ${\mathcal{D}}_i$ is the noise intensity with ${\mathcal{D}}_i\ge 0$ .

Naturally, the increase in the number of infectious individuals occurs with certain spatial dispersion of these infectious individuals that increase the level of the variability (variance) of the contact process due to the change of environment. To include this in the model, we assume that the noise intensity ${\mathcal{D}}_i$ at time t, depends on infectious population size $I\left(t\right)$ . Therefore, by replacing ${\beta }_i\text{d}t$ in the deterministic model described by the flow chart (see Figure 1) with ${\tilde{\beta }}_i\text{d}t={\beta }_i\text{d}t+{\sigma }_{i}I\left(t\right)\text{d}W\left(t\right)$ and by setting $V=S_4$ , we obtain the following delayed stochastic differential equation describing the model,

$\{\begin{array}{l}\text{d}{S}_1\left(t\right)=\left[\Lambda -{\beta }_1{S}_1\left(t\right)H\left(I\left(t-\cdot \right)\right)-\left({\theta }_1+{\mu }_1\right)S_1\left(t\right)+\eta p_{1}R\left(t\right)\right]\text{d}t\\ -{\sigma }_{1}I\left(t\right)S_1\left(t\right)H\left(I\left(t-\cdot \right)\right)\text{d}W\left(t\right),\\ \text{d}{S}_2\left(t\right)=\left[{\theta }_1{S}_1\left(t\right)-{\beta }_2{S}_2\left(t\right)H\left(I\left(t-\cdot \right)\right)-\left({\theta }_2+{\mu }_2\right)S_2\left(t\right)+\eta p_{2}R\left(t\right)\right]\text{d}t\\ -{\sigma }_{2}I\left(t\right)S_2\left(t\right)H\left(I\left(t-\cdot \right)\right)\text{d}W\left(t\right),\\ \text{d}{S}_3\left(t\right)=\left[{\theta }_2{S}_2\left(t\right)-{\beta }_3{S}_3\left(t\right)H\left(I\left(t-\cdot \right)\right)-\left({\mu }_3+{\theta }_3\right)S_3\left(t\right)+\eta p_{3}R\left(t\right)\right]\text{d}t\\ -{\sigma }_{3}I\left(t\right)S_3\left(t\right)H\left(I\left(t-\cdot \right)\right)\text{d}W\left(t\right),\\ \text{d}V\left(t\right)=\left[{\theta }_3{S}_3\left(t\right)-{\beta }_{4}V\left(t\right)H\left(I\left(t-\cdot \right)\right)-\left({\mu }_4+{\theta }_4\right)V\left(t\right)\right]\text{d}t\\ -{\sigma }_{4}I\left(t\right)V\left(t\right)H\left(I\left(t-\cdot \right)\right)\text{d}W\left(t\right),\\ \text{d}I\left(t\right)=\left[{\displaystyle \underset{k=1}{\overset{4}{\sum }}}{\beta }_k{S}_k\left(t\right)H\left(I\left(t-\cdot \right)\right)-\left({\mu }_I+\gamma \right)I\left(t\right)\right]\text{d}t\\ +I\left(t\right){\displaystyle \underset{k=1}{\overset{4}{\sum }}}{\sigma }_k{S}_k\left(t\right)H\left(I\left(t-\cdot \right)\right)\text{d}W\left(t\right),\\ \text{d}R\left(t\right)=\left[\gamma I\left(t\right)+{\theta }_{4}V\left(t\right)-\left({\mu }_R+\eta \right)R\left(t\right)\right]\text{d}t\end{array}$ (1)

with initial condition

$\{\begin{array}{l}{S}_1\left(\theta \right)={\phi }_1\left(\theta \right),S_2\left(\theta \right)={\phi }_2\left(\theta \right),S_3\left(\theta \right)={\phi }_3\left(\theta \right),V\left(\theta \right)={\phi }_4\left(\theta \right),I\left(\theta \right)={\phi }_5\left(\theta \right),\\ R\left(\theta \right)={\phi }_6\left(\theta \right),\phi =\left({\phi }_1,{\phi }_2,{\phi }_3,{\phi }_4,{\phi }_5,{\phi }_6\right)\text{and}\phi \in {\mathcal{C}}_{\left[-h,0\right]}\left({ℝ}_{+}^6\right)\cap {\mathcal{F}}_0,\end{array}$

Figure 1. The chart flow of model describing the transfer rule between the model compartments.

where the set ${\mathcal{C}}_{\left[-h\mathrm{,0}\right]}\left({ℝ}_{+}^6\right)\cap {\mathcal{F}}_0$ will be defined in the next section and the functional $H(\cdot )$ is given by

$H\left(I\left(t-\cdot \right)\right)={\displaystyle {\int }_0^{h}g}\left(\tau \right)\mathcal{K}\left(I\left(t-\tau \right)\right)\text{d}\tau \mathrm{.}$

3. Definition and Notation

Let ${ℝ}_{+}^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 a real $h>0$ , we define ${\mathcal{C}}_{\left[-h\mathrm{,0}\right]}\left({ℝ}^n\right)$ , $x\in {\mathcal{C}}_{\left[-h\mathrm{,}+\infty \right)}\left({ℝ}_{+}^n\right)$ and ${\mathcal{C}}_{\left[-h\mathrm{,}+\infty \right)}\left({ℝ}^n\right)$ respectively as the spaces of continuous functions from $\left[-h\mathrm{,0}\right]$ to ${ℝ}^n$ , from $\left[-h\mathrm{,0}\right]$ to ${ℝ}_{+}^n$ and from $\left[-h\mathrm{,}+\infty \right)$ to ${ℝ}^n$ , endowed with the supremum norm. For any $x\in {\mathcal{C}}_{\left[-h\mathrm{,}+\infty \right)}\left({ℝ}^n\right)$ , $x_t$ denotes the segment process of $x$ given by $x_t\left(\theta \right)=x\left(t+\theta \right)$ , $\theta \in \left[-h\mathrm{,0}\right]$ , $t\ge 0$ . For any vector $v\in {ℝ}^n$ , we denote by $\left|v\right|:={\left({\displaystyle {\sum }_{i=1}^n{v}_i^2}\right)}^{1/2}$ the Euclidean norm and for any $x\in {\mathcal{C}}_{\left[-h\mathrm{,0}\right]}\left({ℝ}^n\right)$ , we denote by $\Vert x\Vert \mathrm{<mo>\; :\; </mo>}={\sup }_{\theta \in \left[-h\mathrm{,0}\right]}\left|x\left(\theta \right)\right|$ the supremum norm.

Consider the general n-dimensional stochastic functional differential equation

$\text{d}X\left(t\right)=D\left(X_t\mathrm{,}t\right)\text{d}t+F\left(X_t\mathrm{,}t\right)\text{d}W\left(t\right)\mathrm{,}{X}_0=\phi \in {\mathcal{C}}_{\left[-h\mathrm{,0}\right]}\left({ℝ}^n\right)\cap {\mathcal{F}}_0$ (2)

where $D\mathrm{<mo>\; :\; </mo>}{\mathcal{C}}_{\left[-h\mathrm{,0}\right]}\left({ℝ}^n\right)\times \left[\mathrm{0,}\infty \right)\to {ℝ}^n$ , $F\mathrm{<mo>\; :\; </mo>}{\mathcal{C}}_{\left[-h\mathrm{,0}\right]}\left({ℝ}^n\right)\times \left[\mathrm{0,}\infty \right)\to {ℝ}^{n\times m}$ and ${\left\{W\left(t\right)\right\}}_{t\ge 0}$ is an m-dimensional Brownian motion $\left(\Omega \mathrm{,}\mathcal{F}\mathrm{,}{\left\{{\mathcal{F}}_t\right\}}_{t\ge 0}\mathrm{,}P\right)$ . Moreover, ${\mathcal{C}}_{\left[-h\mathrm{,0}\right]}\left({ℝ}^n\right)\cap {\mathcal{F}}_0$ denote the set of ${\mathcal{C}}_{\left[-h\mathrm{,0}\right]}\left({ℝ}^n\right)$ -valued random variables that are ${\mathcal{F}}_0$ -measurable.

Let $C^{\mathrm{2,1}}\left({ℝ}^n\times \left[\mathrm{0,}\infty \right)\mathrm{<mo>\; ;\; </mo>}{ℝ}_{+}\right)$ denote the family of all continuous non-negative functions $V\left(x\mathrm{,}t\right)$ defined on ${ℝ}^n\times \left[\mathrm{0,}\infty \right)$ such that they are continuously twice differentiable in x and once in t. For any $V\left(x\mathrm{,}t\right)\in C^{\mathrm{2,1}}\left({ℝ}^n\times \left[\mathrm{0,}\infty \right)\mathrm{<mo>\; ;\; </mo>}{ℝ}_{+}\right)$ , we define the function $\mathcal{L}V\mathrm{<mo>\; :\; </mo>}{\mathcal{C}}_{\left[-h\mathrm{,0}\right]}\left({ℝ}^n\right)\times {ℝ}_{+}\to ℝ$ by

$\begin{array}{c}\mathcal{L}V\left(X_t\mathrm{,}t\right)=V_t\left(X\left(t\right)\mathrm{,}t\right)+V_x\left(X\left(t\right)\mathrm{,}t\right)b\left(X_t\mathrm{,}t\right)\\ +\frac{1}{2}trace\left[{\sigma }^{\text{T}}\left(X_t\mathrm{,}t\right)V_{xx}\left(X\left(t\right)\right)\sigma \left(X_t\mathrm{,}t\right)\right]\end{array}$ (3)

where $V_t=\frac{\partial V}{\partial t}$ , $V_x=\left(\begin{array}{c}\frac{\partial V}{\partial x_1},\cdots ,\frac{\partial V}{\partial x_n}\end{array}\right)$ , $V_{xx}\left(x\right)={\left(\begin{array}{c}\frac{{\partial }^{2}V\left(x\right)}{\partial x_i\partial x_j}\end{array}\right)}_{1\le i,j\le n}$ .

In what follows, we consider the stochastic system (1) which is of the form (2) with dimension $n=6$ . We always assume that the initial value $\phi =\left({\phi }_1\mathrm{,}{\phi }_2\mathrm{,}{\phi }_3\mathrm{,}{\phi }_4\mathrm{,}{\phi }_5\mathrm{,}{\phi }_6\right)\in {\mathcal{C}}_{\left[-h\mathrm{,0}\right]}\left({ℝ}_{+}^6\right)\cap {\mathcal{F}}_0$ which is the set of ${\mathcal{C}}_{\left[-h\mathrm{,0}\right]}\left({ℝ}_{+}^6\right)$ -valued random variables and is ${\mathcal{F}}_0$ -measurable.

4. Existence of Unique Global Positive Solution

In general, we are interested in positive solutions since our study concerns processes describing the size of population compartments. The drift and diffusion coefficients of the system (1) are locally Lipschitz continuous under the assumptions A₁-A₃, for any given initial value. Then system (1) has a unique local solution $X\left(t\right)$ $\left[-h\mathrm{,}{\tau }_e\right)$ , where the explosion time ${\tau }_e$ (see Mao [43] ) is defined by

${\tau }_e=\sup \left\{t\ge 0;\underset{\left[0,t\right]}{\sup }\left|X\left(s\right)\right|<\infty \right\}.$

In order to guarantee that the unique local solution is global, it is necessary to establish its non-explosion in a finite time. The following result assures us of the existence and uniqueness of the global positive solution.

Theorem 1. Let’s assume that A₁-A₃ is valid. Then, for any initial value $\phi =\left({\phi }_1\mathrm{,}{\phi }_2\mathrm{,}{\phi }_3\mathrm{,}{\phi }_4\mathrm{,}{\phi }_5\mathrm{,}{\phi }_6\right)\in {\mathcal{C}}_{\left[-h\mathrm{,0}\right]}\left({ℝ}_{+}^6\right)\cap {\mathcal{F}}_0$ , the model (1) has a unique positive solution on $t\in \left[\mathrm{0,}{\tau }_e\right)$ ,

$X\left(t\right)=\left(S_1\left(t\right)\mathrm{,}{S}_2\left(t\right)\mathrm{,}{S}_3\left(t\right)\mathrm{,}V\left(t\right)\mathrm{,}I\left(t\right)\mathrm{,}R\left(t\right)\right)\in {ℝ}_{+}^6\text{a}\text{.s}\mathrm{.}$

Proof. As mentioned above, under assumptions A₁-A₃, model (1) has a unique local solution $X\left(t\right)$ on $\left[-h\mathrm{,}{\tau }_e\right)$ , for any initial value $\phi =\left({\phi }_1\mathrm{,}{\phi }_2\mathrm{,}{\phi }_3\mathrm{,}{\phi }_4\mathrm{,}{\phi }_5\right)\in {\mathcal{C}}_{\left[-h\mathrm{,0}\right]}\left({ℝ}_{+}^6\right)\cap {\mathcal{F}}_0$ . Let us define the following stopping time as in [44]

${\tau }_{S_1}=\inf \left\{t\in \left[-h,{\tau }_e\right):S\left(t\right)\le 0\right\}.$

In a similar way, we define ${\tau }_{S_2}\mathrm{,}{\tau }_{S_3}\mathrm{,}{\tau }_V\mathrm{,}{\tau }_I\mathrm{,}{\tau }_R$ groups S₂, S₃, V, I and R respectively.

We can see from (1) that $S_1\left(t\right)$ satisfies the linear stochastic differential equation. Thus, by (Mao [43] , Chap.3.) one has

$\begin{array}{l}\text{d}{S}_1\left(t\right)=\left[\Lambda -{\beta }_1{S}_1\left(t\right)H(I_t)-\left({\theta }_1+{\mu }_1\right)S_1\left(t\right)+{\eta }_{1}R\left(t\right)\right]\text{d}t\\ -{\sigma }_{1}I\left(t\right)S_1\left(t\right)H\left(I_t\right)\text{d}W\left(t\right),t\in \left[0,{\tau }_{S_1}\right),\end{array}$

where ${\left(R\left(t\right)\right)}_{t\ge 0}$ is considered as ${\left\{{\mathcal{F}}_t\right\}}_{t\ge 0}$ -adapted and almost surely locally bounded process.

We have

$S_1\left(t\right)=Z_{S_1}\left(t\right)\left(S_1\left(0\right)+{\displaystyle {\int }_0^t}\frac{\Lambda +{\eta }_{1}R\left(u\right)}{Z_{S_1}\left(u\right)}\text{d}u\right)\mathrm{,}t\in \left[\mathrm{0,}{\tau }_{S_1}\right)\mathrm{,}$

where

$\begin{array}{c}{Z}_{S_1}\left(t\right)=\exp [-\left({\theta }_1+{\mu }_1\right)t-{\displaystyle {\int }_0^t\left({\beta }_{1}H\left(I_s\right)-\frac{1}{2}{\left({\sigma }_{1}I\left(t\right)H\left(I_s\right)\right)}^2\right)\text{d}s}\\ \stackrel{}{}-{\sigma }_1{\displaystyle {\int }_0^{t}I\left(s\right)H\left(I_s\right)\text{d}W\left(s\right)}]\\ >0\text{a}\text{.s}\text{.}\end{array}$

Therefore, we deduce that ${\tau }_{S_1}\ge {\tau }_R$ almost surely.

For the second group of susceptible, we have

$S_2\left(t\right)=Z_{S_2}\left(t\right)\left(S_2\left(0\right)+{\displaystyle {\int }_0^t}\frac{{\theta }_1{S}_2\left(u\right)+{\eta }_{2}R\left(u\right)}{Z_{S_2}\left(u\right)}\text{d}u\right),t\in \left[0,{\tau }_{S_2}\right),$

where

$\begin{array}{c}{Z}_{S_2}\left(t\right)=\exp [-\left({\theta }_2+{\mu }_2\right)t-{\displaystyle {\int }_0^t\left({\beta }_{2}H\left(I_s\right)-\frac{1}{2}{\left({\sigma }_{2}I\left(t\right)H\left(I_s\right)\right)}^2\right)\text{d}s}\\ \stackrel{}{}-{\sigma }_2{\displaystyle {\int }_0^{t}I\left(s\right)H\left(I_s\right)\text{d}W\left(s\right)}]\\ >0\text{a}\text{.s}\text{.}\end{array}$

It follows that ${\tau }_{S_2}\ge {\tau }_R$ a.s.

For the third group of susceptible, we have

$S_3\left(t\right)=Z_{S_2}\left(t\right)\left(S_3\left(0\right)+{\displaystyle {\int }_0^t}\frac{{\theta }_2{S}_2\left(t\right)+{\eta }_{3}R\left(t\right)}{Z_{S_3}\left(u\right)}\text{d}u\right),t\in \left[0,{\tau }_{S_3}\right),$

where

$\begin{array}{c}{Z}_{S_3}\left(t\right)=\exp [-\left({\mu }_3+{\theta }_3\right)t-{\displaystyle {\int }_0^t\left({\beta }_{3}H\left(I_s\right)-\frac{1}{2}{\left({\sigma }_{3}I\left(s\right)H\left(I_s\right)\right)}^2\right)\text{d}s}\\ \stackrel{}{}-{\sigma }_3{\displaystyle {\int }_0^{t}I\left(s\right)H\left(I_s\right)\text{d}W\left(s\right)}]\\ >0\text{a}\text{.s}\text{.}\end{array}$

So, we get ${\tau }_{S_3}\ge {\tau }_R$ a.s. In the same way as before, we have ${\tau }_V\ge {\tau }_R$ a.s.

The recovered population $R\left(t\right)$ satisfies the linear stochastic differential equation

$\text{d}R\left(t\right)=\left[\gamma I\left(t\right)+{\theta }_{4}V\left(t\right)-\left({\mu }_R+\eta \right)R\left(t\right)\right]\text{d}t\mathrm{,}$

where ${\left(I\left(t\right)\right)}_{t\ge 0}$ is an ${\left\{\mathcal{F}\left(t\right)\right\}}_{t\ge 0}$ -adapted and almost surely locally bounded process.

$R\left(t\right)=Z_R\left(t\right)\left(R\left(0\right)+{\displaystyle {\int }_0^t}\frac{\gamma I\left(u\right)+{\theta }_{4}V\left(u\right)}{Z_R\left(u\right)}\text{d}u\right)\mathrm{,}t\in \left[\mathrm{0,}{\tau }_R\right)$

where

$Z_R\left(t\right)=\exp \left[-\left({\mu }_R+\eta \right)t\right]>0a.s.$

Since ${\tau }_V\ge {\tau }_R$ , thus $R\left(t\right)$ become negative after only $I\left(t\right)$ is negative. That is ${\tau }_I\le \min \left\{{\tau }_{S_1}\mathrm{,}{\tau }_{S_2}\mathrm{,}{\tau }_{S_3}\mathrm{,}{\tau }_V\mathrm{,}{\tau }_R\right\}$ almost surely.

For the infected group, we have

$I\left(t\right)=Z_I\left(t\right)\left(I\left(0\right)+{\displaystyle {\int }_0^t}\frac{H\left(I_t\right){\displaystyle {\sum }_{k=1}^4{\beta }_k{S}_k\left(t\right)}}{Z_I\left(u\right)}\text{d}u\right),t\in \left[0,{\tau }_e\right)$

and

$\begin{array}{c}{Z}_I\left(t\right)=\exp \left[-\left({\mu }_I+\gamma \right)t+\frac{1}{2}{\sigma }^2{\displaystyle {\int }_0^t{K}^2}\left(I_s,S\left(s\right)\right)\text{d}s+\sigma {\displaystyle {\int }_0^{t}K}\left(I_s,S\left(s\right)\right)\text{d}W\left(s\right)\right]\\ >0\text{a}\text{.s}.\end{array}$

with $K\left(I_t,S\left(t\right)\right)=H\left(I_t\right){\displaystyle {\sum }_{k=1}^4{\sigma }_k{S}_k\left(t\right)}$ .

We propose to show in the following step that ${\tau }_e\le \min \left\{{\tau }_{S_1}\mathrm{,}{\tau }_R\mathrm{,}{\tau }_{S_2}\mathrm{,}{\tau }_{S_3}\mathrm{,}{\tau }_I\right\}$ almost surely by establishing that ${\tau }_e\le {\tau }_I$ almost surely. To do this, we proceed by contradiction.

Suppose that there exists a Borel set B of Ω with $P\left(B\right)>0$ and for all $\omega \in B$ we have ${\tau }_I\left(\omega \right)\le {\tau }_e\left(\omega \right)$ . By definition of ${\tau }_I$

$I\left(u\right)>0,\forall u\in \left[0,{\tau }_I\right)\text{and}I\left({\tau }_I\right)=0.$

In view of assumption (A₂) and the fact that ${\tau }_I\le \min \left\{{\tau }_{S_1}\mathrm{,}{\tau }_{S_2}\mathrm{,}{\tau }_{S_3}\mathrm{,}{\tau }_V\right\}$ a.s., for all $\omega \in B$ and for all $u\in \left[\mathrm{0,}{\tau }_I\left(\omega \right)\right)$ , it follows that

$H\left(I_u\left(\cdot ,\omega \right)\right)>0\text{and}{S}_k\left(u,\omega \right)>0,k=1,2,3,4$

Therefore

$\begin{array}{c}0=\underset{t\to {\tau }_I\left(\omega \right)}{\lim }I\left(t,\omega \right)\\ =Z_I\left({\tau }_I\left(\omega \right),\omega \right)\left(I\left(0,\omega \right)+{\displaystyle {\int }_0^t}\frac{H\left(I_u\left(\cdot ,\omega \right)\right){\displaystyle {\sum }_{k=1}^4{\beta }_k{S}_k\left(u,\omega \right)}}{Z_I\left(u,\omega \right)}\text{d}u\right)>0,\end{array}$

which leads to a contradiction. Necessary, we must have ${\tau }_e\le {\tau }_I$ almost surely.

□

Let us put

$\begin{array}{l}{\Gamma }_{\epsilon }=\left\{\left(u_1,u_2,u_3,u_4,u_5,u_6\right)\in {ℝ}_{+}^6|u_1+u_2+u_3+u_4+u_5+u_6<\frac{\Lambda }{{\mu }_1}+\epsilon \right\}\\ \text{and}\Gamma ={\displaystyle {\cap }_{\epsilon >0}{\Gamma }_{\epsilon }}\end{array}$

The following result gives us the boundness of any local solution of the model (1) and achieves the proof of existence and uniqueness of global positive solution.

Corollary 1. Assume that A₁-A₃ hold. Then, the system (1) has a unique global positive solution $X\left(t\right)=\left(S_1\left(t\right)\mathrm{,}{S}_2\left(t\right)\mathrm{,}{S}_3\left(t\right)\mathrm{,}V\left(t\right)\mathrm{,}I\left(t\right)\mathrm{,}R\left(t\right)\right)$ for any initial value $\phi =\left({\phi }_1\mathrm{,}{\phi }_2\mathrm{,}{\phi }_3\mathrm{,}{\phi }_4\mathrm{,}{\phi }_5\mathrm{,}{\phi }_6\right)\in {\mathcal{C}}_{\left[-h\mathrm{,0}\right]}\left({ℝ}_{+}^6\right)\cap {\mathcal{F}}_0$ . Moreover, there exists a $T=T\left(\epsilon \right)>0$ for any sufficiently small $\epsilon >0$ such that for all $t>T$ this solution remains in ${\Gamma }_{\epsilon }$ with probability 1. In particular, if $\phi \left(\left[-h\mathrm{,0}\right]\right)\subset \Gamma$ this solution lies in Γ almost surly.

Proof. By Theorem 1, for any initial value $\phi =\left({\phi }_1\mathrm{,}{\phi }_2\mathrm{,}{\phi }_3\mathrm{,}{\phi }_4\mathrm{,}{\phi }_5\mathrm{,}{\phi }_6\right)\in {\mathcal{C}}_{\left[-h\mathrm{,0}\right]}\left({ℝ}_{+}^6\right)\cap {\mathcal{F}}_0$ , the model (1) has a unique positive solution $t\in \left[\mathrm{0,}{\tau }_e\right)$ . To prove that this solution is global, it suffices to prove that it is bounded. We have $\max \left\{S_1\left(t\right)\mathrm{,}{S}_2\left(t\right)\mathrm{,}{S}_3\left(t\right)\mathrm{,}V\left(t\right)\mathrm{,}I\left(t\right)\mathrm{,}R\left(t\right)\right\}\le N\left(t\right)$ a.s. on $\left[\mathrm{0,}{\tau }_e\right)$ , where $N\left(t\right)=S_1\left(t\right)+S_2\left(t\right)+S_3\left(t\right)+V\left(t\right)+I\left(t\right)+R\left(t\right)$ . Therefore

$\begin{array}{c}\text{d}N\left(t\right)=\left[\Lambda -{\mu }_1{S}_1\left(t\right)-{\mu }_2{S}_2\left(t\right)-{\mu }_3{S}_3\left(t\right)-{\mu }_{4}V\left(t\right)-{\mu }_{I}I\left(t\right)-{\mu }_{R}R\left(t\right)\right]\text{d}t\\ =[\Lambda -\left({\mu }_2-{\mu }_1\right)S_1\left(t\right)-\left({\mu }_3-{\mu }_1\right)S_3\left(t\right)-\left({\mu }_4-{\mu }_1\right)S_3\left(t\right)\\ -\left({\mu }_4-{\mu }_1\right)V\left(t\right)-\left({\mu }_I-{\mu }_1\right)I\left(t\right)-\left({\mu }_R-{\mu }_1\right)R\left(t\right)-{\mu }_{1}N\left(t\right)]\text{d}t,\\ \le \left[\Lambda -{\mu }_{1}N\left(t\right)\right]\text{d}t,t\in \left[0,{\tau }_e\right)\end{array}$

because naturally ${\mu }_1\le \min \left\{{\mu }_2\mathrm{,}{\mu }_3\mathrm{,}{\mu }_4\mathrm{,}{\mu }_I\mathrm{,}{\mu }_R\right\}$ .

We denote by $\underset{\_}{N}\left(t\right)$ the solution of the following differential equation

$\text{d}\underset{\_}{N}\left(t\right)=\left[\Lambda -{\mu }_{1}N\left(t\right)\right]\text{d}t$ , with the same initial value $\underset{\_}{N}\left(0\right)=N\left(0\right)$ .

Therefore, by the comparison theorem [45] , we have

$N\left(t\right)\le \underset{\_}{N}\left(t\right)=\left(N\left(0\right)-\frac{\Lambda }{{\mu }_1}\right){\text{e}}^{-{\mu }_{1}t}+\frac{\Lambda }{{\mu }_1}<\infty \text{on}\left[0,\infty \right)\text{a}\text{.s}.$

For any $\epsilon >0$ , there exists $T=T\left(\epsilon \right)>0$ such that for all $t>T$

$N\left(t\right)<\frac{\Lambda }{{\mu }_1}$

□

5. Stability Analysis and Asymptotic Behavior

This section deals with the stability analysis of the stochastic epidemic model (1). A simple analysis established that this model has a unique disease free-equilibrium $E^0=\left(s_1^0\mathrm{,}{s}_2^0\mathrm{,}{s}_3^0\mathrm{,}{s}_4^0\mathrm{,0,}{r}^0\right)$ where $v^0=s_4^0$ , which is given by:

$\begin{array}{l}{s}_1^0=\frac{\Lambda }{{\mu }_1+{\theta }_1-{\eta }_1{\eta }_2{\alpha }_1{\alpha }_2{\alpha }_3{\alpha }_4},s_2^0={\alpha }_4{s}_1^0,s_3^0={\alpha }_3{s}_2^0,v^0={\alpha }_3{s}_3^0,\\ i^0=0\text{and}{r}^0={\alpha }_1{v}^0\end{array}$

where

$\begin{array}{l}{\alpha }_1=\frac{{\theta }_4}{{\mu }_R+{\eta }_1+{\eta }_2+{\eta }_3},{\alpha }_2=\frac{{\theta }_3}{{\mu }_4+{\theta }_4},{\alpha }_3=\frac{{\theta }_2}{{\mu }_3+{\theta }_3-{\eta }_3{\alpha }_1{\alpha }_2}\\ \text{and}{\alpha }_4=\frac{{\theta }_1}{{\mu }_2+{\theta }_2-{\eta }_2{\alpha }_1{\alpha }_2{\alpha }_3}\mathrm{.}\end{array}$

By using a change of variable, we first reduce the analysis of the stability of the equilibrium point $E^0$ to the study of the stability of the trivial equilibrium point zero of another system. We first establish that the solution of the system obtained by change of variable remains in suitable subset ${ℝ}^6$ . Then, by using a Lyapunov functional technique and a local martingale convergence result, we deduce the almost sure stability of the disease-free equilibrium $E^0$ of the model (1) under the condition $R_0\le 1$ . In the following, we consider the class of initial conditions $\phi =\left({\phi }_1\mathrm{,}{\phi }_2\mathrm{,}{\phi }_3\mathrm{,}{\phi }_4\mathrm{,}{\phi }_5\mathrm{,}{\phi }_6\right)\in {\mathcal{C}}_{\left[-h\mathrm{,0}\right]}\left({ℝ}_{+}^6\right)\cap {\mathcal{F}}_0$ such that $\phi \left(\left[-h\mathrm{,0}\right]\right)\subset \Gamma$ .

Let’s put

$y_1=S_1-s_1^0,y_2=S_2-s_2^0,y_3=S_3-s_3^0,y_4=V-v^0,y_5=I,y_6=R-r^0.$ (4)

Then, by virtue of Itô’s formula, we get the following system

$\{\begin{array}{l}\text{d}{y}_1\left(t\right)=\left[-\left({\mu }_1+{\theta }_1\right)y_1\left(t\right)-{\beta }_1\left(y_1\left(t\right)+s_1^0\right)H\left(y_5\left(t-\cdot \right)\right)+p_1\eta y_6\left(t\right)\right]\text{d}t\\ -\sigma y_5\left(t\right)\left(y_1\left(t\right)+s_1^0\right)H\left(y_5\left(t-\cdot \right)\right)\text{d}W\left(t\right),\\ \text{d}{y}_2\left(t\right)=\left[{\theta }_1{y}_1\left(t\right)-\left({\mu }_1+{\theta }_2\right)y_2\left(t\right)-{\beta }_2\left(y_2\left(t\right)+s_2^0\right)H\left(y_5\left(t-\cdot \right)\right)+p_2\eta y_6\left(t\right)\right]\text{d}t\\ -\sigma y_5\left(t\right)\left(y_2\left(t\right)+s_2^0\right)H\left(y_5\left(t-\cdot \right)\right)\text{d}W\left(t\right),\\ \text{d}{y}_3\left(t\right)=\left[{\theta }_2{y}_2\left(t\right)-\left({\mu }_3+{\theta }_3\right)y_3\left(t\right)-{\beta }_3\left(y_3\left(t\right)+s_3^0\right)H\left(y_5\left(t-\cdot \right)\right)+p_3\eta y_6\left(t\right)\right]\text{d}t\\ -\sigma y_4\left(t\right)\left(y_3\left(t\right)+s_3^0\right)H\left(y_5\left(t-\cdot \right)\right)\text{d}W\left(t\right),\\ \text{d}{y}_4\left(t\right)=\left[{\theta }_3{y}_3\left(t\right)-\left({\mu }_4+{\theta }_4\right)y_4\left(t\right)-{\beta }_4\left(y_4\left(t\right)+s_4^0\right)H\left(y_5\left(t-\cdot \right)\right)+p_3\eta y_3\left(t\right)\right]\text{d}t\\ -\sigma y_4\left(t\right)\left(y_3\left(t\right)+s_3^0\right)H\left(y_5\left(t-\cdot \right)\right)\text{d}W\left(t\right),\\ \text{d}{y}_5\left(t\right)=\left[{\displaystyle \underset{k=1}{\overset{4}{\sum }}}{\beta }_k\left(y_k\left(t\right)+s_k^0\right)H\left(y_5\left(t-\cdot \right)\right)-\left({\mu }_I+\gamma \right)y_5\left(t\right)\right]\text{d}t\\ +\sigma y_5\left(t\right){\displaystyle \underset{k=1}{\overset{4}{\sum }}}\left(y_k\left(t\right)+s_k^0\right)H\left(y_5\left(t-\cdot \right)\right)\text{d}W\left(t\right),\\ \text{d}{y}_6\left(t\right)=\left[\gamma y_5\left(t\right)+{\mu }_4{y}_4\left(t\right)-\left({\mu }_3+\eta \right)y_6\left(t\right)\right]\text{d}t,\end{array}$ (5)

with initial condition

$\{\begin{array}{l}{y}_1\left(\theta \right)={\psi }_1\left(\theta \right),y_2\left(\theta \right)={\psi }_2\left(\theta \right),y_3\left(\theta \right)={\psi }_3\left(\theta \right),y_4\left(\theta \right)={\psi }_4\left(\theta \right),y_5\left(\theta \right)={\psi }_5\left(\theta \right),\\ \psi =\left({\psi }_1,{\psi }_2,{\psi }_3,{\psi }_4,{\psi }_5\right)\text{and}\psi +E^0\in {\mathcal{C}}_{\left[-h,0\right]}\left({ℝ}_{+}^6\right)\cap {\mathcal{F}}_0,\end{array}$

where

$H\left(y_5\left(t-\cdot \right)\right)={\displaystyle {\int }_0^{h}g}\left(\tau \right)\mathcal{K}\left(y_5\left(t-\tau \right)\right)\text{d}\tau \mathrm{.}$

It is easy to see that the stability analysis of the disease-free equilibrium $E^0$ of the model (1) can be obtained from the stability analysis of the trivial solution $E_y^0=\left(\mathrm{0,0,0,0,0,0}\right)$ of the model (5). Before studying the stability analysis of the trivial solution of the model (5), we need some information about the sign of the components of its solution.

Theorem 2. Either $\left(S_1\left(t\right)\mathrm{,}{S}_2\left(t\right)\mathrm{,}{S}_3\left(t\right)\mathrm{,}V\left(t\right)\mathrm{,}I\left(t\right)\mathrm{,}R\left(t\right)\right)$ the solution of the system (1) with initial condition $\phi =\left({\phi }_1\mathrm{,}{\phi }_2\mathrm{,}{\phi }_3\mathrm{,}{\phi }_4\mathrm{,}{\phi }_5\right)\in {\mathcal{C}}_{\left[-h\mathrm{,0}\right]}\left({ℝ}_{+}^6\right)\cap {\mathcal{F}}_0$ . Suppose that

$\begin{array}{l}{\phi }_1\left(\theta \right)+{\phi }_2\left(\theta \right)+{\phi }_3\left(\theta \right)+{\phi }_4\left(\theta \right)+{\phi }_6\left(\theta \right)<s_1^0+s_2^0+s_3^0+s_4^0+r^0,\\ \text{forall}\theta \in \left[-h\mathrm{,0}\right]\end{array}$ (6)

and

$\mathcal{K}\left(s_0\right)\ge \underset{1\le k\le 4}{\min }\frac{2{\beta }_k}{{\left({\sigma }_k{s}^0\right)}^2}\mathrm{,}\text{where}{s}^0=\frac{\Lambda }{\mu }\mathrm{.}$ (7)

Then, we have

$S_1\left(t\right)<s_1^0,S_2\left(t\right)<s_2^0,S_3\left(t\right)<s_3^0,V\left(t\right)<v^0\text{and}R\left(t\right)<r^0\text{a}\text{.s}.,\text{for}\text{all}t\ge 0.$

Proof. To arrive at the result, reformulate the equilibrium states $s_1^0$ , $s_2^0$ , $s_3^0$ , $v^0$ . A simple analysis gives us

$s_1^0=\frac{\Lambda +{\eta }_1{r}^0}{{\mu }_1+{\theta }_1},s_2^0=\frac{{\theta }_1{s}_1^0+{\eta }_1{r}^0}{{\mu }_2+{\theta }_2},s_3^0=\frac{{\theta }_2{s}_2^0+{\eta }_2{r}^0}{{\mu }_3+{\theta }_3},v^0=\frac{{\theta }_3{s}_3^0}{{\mu }_4+{\theta }_4}.$

Based on the proof of the Theorem 1 and the Corollary 1, we get

$S_1\left(t\right)=Z_{S_2}\left(t\right)\left(S_1\left(0\right)+{\displaystyle {\int }_0^t}\frac{\Lambda +{\eta }_{1}R\left(u\right)}{Z_{S_1}\left(u\right)}\text{d}u\right),t\ge 0$

where

$\begin{array}{c}{Z}_{S_1}\left(t\right)=\exp [-\left({\theta }_1+{\mu }_1\right)t-{\displaystyle {\int }_0^t\left({\beta }_{2}H\left(I_s\right)-\frac{1}{2}{\left({\sigma }_{2}I\left(t\right)H\left(I_s\right)\right)}^2\right)\text{d}s}\\ \stackrel{}{}-{\displaystyle {\int }_0^t{\sigma }_{2}I\left(s\right)H\left(I_s\right)\text{d}W\left(s\right)}],\end{array}$

Let $Y\left(s\right)={\sigma }_2{\displaystyle {\int }_u^{t}I\left(s\right)H\left(I_s\right)\text{d}W\left(s\right)}$ , the quadratic variation

$\langle Y\left(s\right)\rangle ={\sigma }_2^2{\displaystyle {\int }_u^t}{\left(I\left(s\right)H\left(I_s\right)\right)}^2\text{d}s$ of $Y\left(s\right)$ is locally bounded by Corollary 1. So, by the strong law of large numbers for local martingales (see e.g. [43] ) we have

$\frac{Y\left(t\right)}{t}\to 0$ when $t\to \infty$ . Therefore, there exists $T^0>0$ large enough, such that for all $t>T^0$ such that

$\frac{Y\left(s\right)}{t}\approx 0\text{and}\frac{1}{t}{\displaystyle {\int }_{T_0}^t}\left({\beta }_{2}H\left(I_s\right)-\frac{1}{2}{\left({\sigma }_{k}I\left(t\right)H\left(I_s\right)\right)}^2\right)\text{d}s\approx 0.$

It follows for all $t>T^0$ that

$\begin{array}{c}{Z}_{S_k}\left(u\right)=\exp \left[-\left({\theta }_1+{\mu }_1\right)t-{\displaystyle {\int }_0^t}\left({\beta }_{1}H\left(I_s\right)-\frac{1}{2}{\left({\sigma }_{1}I\left(t\right)H\left(I_s\right)\right)}^2\right)\text{d}s-Y\left(t\right)\right]\\ =\exp [t(-\left({\theta }_1+{\mu }_1\right)-\frac{1}{t}{\displaystyle {\int }_0^{T_0}}\left({\beta }_{1}H\left(I_s\right)-\frac{1}{2}{\left({\sigma }_{k}I\left(t\right)H\left(I_s\right)\right)}^2\right)\text{d}s\\ -\frac{1}{t}{\displaystyle {\int }_{T_0}^u}\left({\beta }_{1}H\left(I_s\right)-\frac{1}{2}{\left({\sigma }_{k}I\left(t\right)H\left(I_s\right)\right)}^2\right)\text{d}s-\frac{Y\left(t\right)}{t})]\\ \approx \exp \left[-\left({\theta }_1+{\mu }_1\right)t-{\displaystyle {\int }_{T_0}^t}\left({\beta }_{1}H\left(I_s\right)-\frac{1}{2}{\left({\sigma }_{k}I\left(t\right)H\left(I_s\right)\right)}^2\right)\text{d}s\right],\end{array}$

On the other hand, it is easy to see that

${\beta }_{k}x-\frac{1}{2}{\left(\frac{{\sigma }_1\Lambda }{\mu }\right)}^2{x}^2\ge 0\text{implies}x\ge \underset{1\le k\le 4}{\min }\frac{2{\beta }_1}{{\left({\sigma }_1{s}^0\right)}^2}\mathrm{.}$ (8)

In view of Corollary (1) and assumptions A₁-A₂, we have

$I\left(t\right)\le s_0\text{and}\mathcal{K}\left(I\left(t\right)\right)\le s_0\mathrm{,}\forall t\ge 0.$ (9)

Therefore

$H\left(I_t\right)={\displaystyle {\int }_0^{h}g}\left(\tau \right)\mathcal{K}\left(I\left(t-\tau \right)\right)\text{d}\tau \le K\left(s_0\right){\displaystyle {\int }_0^{h}g}\left(\tau \right)\text{d}\tau \le \mathcal{K}\left(s_0\right)$ (10)

By combining (8), (9) and (10), we deduce that

${\beta }_{2}H\left(I_t\right)-\frac{1}{2}{\left({\sigma }_{2}I\left(t\right)H\left(I_t\right)\right)}^2\ge \mathrm{0,}\forall t\ge 0\text{implies}\mathcal{K}\left(s_0\right)\ge \underset{1\le k\le 4}{\min }\frac{2{\beta }_1}{{\left({\sigma }_1{s}^0\right)}^2}$

Consequently, under the condition (7), for all $t,u>T^0$ , we get

$\begin{array}{c}{Z}_{S_1}\left(t\right)\approx \exp \left[-\left({\theta }_1+{\mu }_1\right)t-{\displaystyle {\int }_{T_0}^t}\left({\beta }_{1}H\left(I_s\right)-\frac{1}{2}{\left({\sigma }_{1}I\left(t\right)H\left(I_s\right)\right)}^2\right)\text{d}s\right]\\ \le \exp \left[-\left({\theta }_1+{\mu }_1\right)t\right],\end{array}$ (11)

And

$\begin{array}{c}\frac{Z_{S_1}\left(t\right)}{Z_{S_1}\left(u\right)}\approx \exp \left[-\left({\theta }_k+{\mu }_k\right)\left(t-u\right)-{\displaystyle {\int }_u^t}\left({\beta }_{2}H\left(I_s\right)-\frac{1}{2}{\left({\sigma }_{k}I\left(t\right)H\left(I_s\right)\right)}^2\right)\text{d}s\right]\\ \le \exp \left[-\left({\theta }_k+{\mu }_k\right)\left(t-u\right)\right],\end{array}$ (12)

Now, let put $\lim {\sup }_{t\to \infty }R\left(t\right)=r^{*}$ . So, there exists two positive reals $\epsilon$ and $T^{\epsilon }$ such ${\sup }_{t>T^{\epsilon }}R\left(t\right)=\epsilon +r^{*}$ . Therefore, for all $t>T^{\epsilon }$ , we have

$\begin{array}{c}{S}_1\left(t\right)=Z_{S_1}\left(t\right)\left(S_1\left(0\right)+{\displaystyle {\int }_0^t}\frac{\Lambda +{\eta }_{1}R\left(t\right)}{Z_{S_1}\left(u\right)}\text{d}u\right)\\ \le Z_{S_1}\left(t\right)S_1\left(0\right)+\left(\Lambda +{\eta }_1\left(\epsilon +r^{*}\right)\right)Z_{S_1}\left(t\right){\displaystyle {\int }_{T^{\epsilon }}^t}{Z}_{S_1}^{-1}\left(u\right)\text{d}u\\ +Z_{S_1}\left(t\right){\displaystyle {\int }_0^{T^{\epsilon }}}\frac{\Lambda +{\eta }_{1}R\left(t\right)}{Z_{S_1}\left(u\right)}\text{d}u.\end{array}$

Let us put $u_0=\max \left\{T^0,T^{\epsilon }\right\}$ , it is straightforward to see that ${\displaystyle {\int }_0^{u_0}}\frac{\Lambda +{\eta }_{1}R\left(t\right)}{Z_{S_1}\left(u\right)}\text{d}u=x_0<\infty$ . Based on (11) and (12), we get for all $t>u_0$

$\begin{array}{c}{S}_1\left(t\right)\le Z_{S_1}\left(t\right)\left(S_1\left(0\right)+x_0\right)+\left(\Lambda +{\eta }_1\left(r^{*}+\epsilon \right)\right)Z_{S_1}\left(t\right){\displaystyle {\int }_{u_0}^t}{Z}_{S_1}^{-1}\left(u\right)\text{d}u\\ \approx \left(S_1\left(0\right)+x_0\right){\text{e}}^{-\left({\theta }_1+{\mu }_1\right)t}+\left(\Lambda +{\eta }_1\left(r^{*}+\epsilon \right)\right){\displaystyle {\int }_0^t}{\text{e}}^{-\left({\theta }_1+{\mu }_1\right)\left(t-u\right)}\text{d}u\\ =\left(S_1\left(0\right)+x_0\right){\text{e}}^{-\left({\theta }_1+{\mu }_1\right)t}+\frac{\Lambda +{\eta }_1\left(r^{*}+\epsilon \right)}{{\theta }_1+{\mu }_1}-\frac{\left(\Lambda +{\eta }_1\left(r^{*}+\epsilon \right)\right){\text{e}}^{-\left({\theta }_1+{\mu }_1\right)t}}{{\theta }_1+{\mu }_1}.\end{array}$

As

${\displaystyle {\int }_0^t{\text{e}}^{-\left({\theta }_1+{\mu }_1\right)\left(t-u\right)}\text{d}u}={\text{e}}^{-\left({\theta }_1+{\mu }_1\right)t}{\displaystyle {\int }_0^t{\text{e}}^{\left({\theta }_1+{\mu }_1\right)u}\text{d}u}=\frac{1}{{\theta }_1+{\mu }_1}-\frac{{\text{e}}^{-\left({\theta }_1+{\mu }_1\right)t}}{{\theta }_1+{\mu }_1}.$

Letting $\epsilon \to 0$ , given the fact that the population may be without infectious ( $I=0$ ), we deduce that

$\underset{t\to \infty }{\lim \sup }{S}_1\left(t\right)=\frac{\Lambda +{\eta }_1{r}^{*}}{{\theta }_1+{\mu }_1}=s_1^{*}\text{a}\text{.s}\mathrm{..}$ (13)

Taking this result into account and repeating the same reasoning on the following expressions

$S_2\left(t\right)=Z_{S_2}\left(t\right)\left(S_2\left(0\right)+{\displaystyle {\int }_0^t}\frac{{\theta }_1{S}_2\left(u\right)+{\eta }_{2}R\left(u\right)}{Z_{S_2}\left(u\right)}\text{d}u\right),t\ge 0$

where

$\begin{array}{c}{Z}_{S_2}\left(t\right)=\exp [-\left({\theta }_2+{\mu }_2\right)t-{\displaystyle {\int }_0^t}\left({\beta }_{2}H\left(I_s\right)-\frac{1}{2}{\left({\sigma }_{2}I\left(t\right)H\left(I_s\right)\right)}^2\right)\text{d}s\\ \stackrel{}{}-{\displaystyle {\int }_0^t}{\sigma }_{2}H\left(I_s\right)\text{d}W\left(s\right)],\end{array}$

$S_3\left(t\right)=Z_{S_2}\left(t\right)\left(S_3\left(0\right)+{\displaystyle {\int }_0^t}\frac{{\theta }_2{S}_2\left(t\right)+{\eta }_{3}R\left(t\right)}{Z_{S_3}\left(u\right)}\text{d}u\right),t\ge 0$

where

$\begin{array}{c}{Z}_{S_3}\left(t\right)=\exp [-\left({\mu }_3+{\theta }_3\right)t-{\displaystyle {\int }_0^t}\left({\beta }_{3}H\left(I_s\right)-\frac{1}{2}{\left({\sigma }_{3}I\left(s\right)H\left(I_s\right)\right)}^2\right)\text{d}s\\ \stackrel{}{}-{\displaystyle {\int }_0^t}\sigma H\left(I_s\right)\text{d}W\left(s\right)],\end{array}$

and

$V\left(t\right)=Z_V\left(t\right)\left(V\left(0\right)+{\displaystyle {\int }_0^t}\frac{{\theta }_2{S}_3\left(t\right)}{Z_V\mathrm{<mo\; stretchy="false">\; (\; </mo>}u\mathrm{<mo\; stretchy="false">\; )\; </mo>}}\text{d}u\right)\mathrm{,}t\ge 0$

where

$\begin{array}{c}{Z}_V\left(t\right)=\exp [-\left({\mu }_4+{\theta }_4\right)t-{\displaystyle {\int }_0^t}\left({\beta }_{4}H\left(I_s\right)-\frac{1}{2}{\left({\sigma }_{4}I\left(s\right)H\left(I_s\right)\right)}^2\right)\text{d}s\\ \stackrel{}{}-{\displaystyle {\int }_0^t}\sigma H\left(I_s\right)\text{d}W\left(s\right)],\end{array}$

respectively, we deduce almost surly, that

$\begin{array}{l}\underset{t\to \infty }{\lim \sup }{S}_2\left(t\right)=\frac{{\theta }_1{s}_1^{*}+{\eta }_1{r}^{*}}{{\mu }_2+{\theta }_2}=s_2^{*},\underset{t\to \infty }{\lim \sup }{S}_3\left(t\right)=\frac{{\theta }_2{s}_2^{*}+{\eta }_2{r}^{*}}{{\mu }_3+{\theta }_3}=s_3^{*},\\ \underset{t\to \infty }{\lim \sup }V\left(t\right)=\frac{{\theta }_3{s}_3^{*}}{{\mu }_4+{\theta }_4}=v^{*}.\end{array}$ (14)

Let us put $y_1=S_1-s_1^0$ , $y_2=S_2-s_2^0$ , $y_3=S_3-s_3^0$ , $y_4=V-v^0$ , $y_5=I$ et $y_6=R-r^0$ , therefore $Y\left(t\right)=y_1\left(t\right)+y_2\left(t\right)+y_3\left(t\right)+y_4\left(t\right)+y_5\left(t\right)+y_6\left(t\right)=N\left(t\right)-s_1^0-s_2^0-s_3^0-v^0-r^0$ . In view of (5), we get

$\begin{array}{l}dY\left(t\right)=-{\mu }_1{y}_1\left(t\right)-{\mu }_2{y}_2\left(t\right)-{\mu }_3{y}_3\left(t\right)-{\mu }_4{y}_4\left(t\right)-{\mu }_I{y}_5\left(t\right)-{\mu }_R{y}_6\left(t\right)\\ \le -\min \left\{{\mu }_1\mathrm{,}{\mu }_2\mathrm{,}-{\mu }_3\mathrm{,}{\mu }_4\mathrm{,}{\mu }_I\mathrm{,}{\mu }_R\right\}Y\left(t\right)\mathrm{,}\text{where}\mu =\min \left\{{\mu }_1\mathrm{,}{\mu }_2\mathrm{,}-{\mu }_3\mathrm{,}{\mu }_4\mathrm{,}{\mu }_I\mathrm{,}{\mu }_R\right\}\mathrm{.}\end{array}$

Under the condition (6), it follows that $Y\left(0\right)=N\left(0\right)-s_1^0-s_2^0-s_3^0-v^0-r^0\le 0$

Therefore

$Y\left(t\right)\le Y\left(0\right){\text{e}}^{-\max \left\{{\mu }_1,{\mu }_2,-{\mu }_3,{\mu }_4,{\mu }_I,{\mu }_R\right\}t}\le 0,$

so we get

$N\left(t\right)\le s_1^0+s_2^0+s_3^0+v^0+r^0\mathrm{.}$

Consequently, for all $t\ge 0$ ,

$\left(s_1^{\mathrm{<mo>\; *\; </mo>}}-s_1^0\right)+\left(s_2^{\mathrm{<mo>\; *\; </mo>}}-s_2^0\right)+\left(s_3^{\mathrm{<mo>\; *\; </mo>}}-s_3^0\right)+\left(v^{\mathrm{<mo>\; *\; </mo>}}-v^0\right)+R\left(t\right)\le r^0\text{a}\text{.s}\mathrm{.}$ (15)

On the other hand,

$s_1^{\mathrm{<mo>\; *\; </mo>}}-s_1^0=K_1\left(r^{\mathrm{<mo>\; *\; </mo>}}-r^0\right)$ , $s_2^{\mathrm{<mo>\; *\; </mo>}}-s_2^0=K_2\left(r^{\mathrm{<mo>\; *\; </mo>}}-r^0\right)$ , $s_3^{\mathrm{<mo>\; *\; </mo>}}-s_3^0=K_3\left(r^{\mathrm{<mo>\; *\; </mo>}}-r^0\right)$ , $v^{\mathrm{<mo>\; *\; </mo>}}-v^0=K_4\left(r^{\mathrm{<mo>\; *\; </mo>}}-r^0\right)$ ,

where

$K_1=k_1$ , $K_2=k_1{k}_2+k_3$ , $K_3=k_1{k}_2{k}_4+k_3{k}_4+k_5$ and $K_4=k_6\left(k_1{k}_2{k}_4+k_3{k}_4+k_5\right)$ ,

with

$k_1=\frac{{\eta }_1}{{\mu }_1+{\theta }_1}$ , $k_2=\frac{{\theta }_1}{{\mu }_2+{\theta }_2}$ , $k_3=\frac{{\eta }_2}{{\mu }_2+{\theta }_2}$ , $k_4=\frac{{\theta }_2}{{\mu }_3+{\theta }_3}$ , $k_5=\frac{{\eta }_3}{{\mu }_3+{\theta }_3}$ , $k_6=\frac{{\theta }_3}{{\mu }_4+{\theta }_4}$ .

Therefore, based on (15), we have almost surly $\left(K_1+K_2+K_3+K_4\right)\left(r^{\mathrm{<mo>\; *\; </mo>}}-r^0\right)+R\left(t\right)\le r^0$ for all $t\ge 0$ .

Let us assume $r^{*}>r^0$ , that is $\left(K_1+K_2+K_3+K_4\right)\left(r^{*}-r^0\right)=e>0$ . It follows that $e+R\left(t\right)\le r^0$ , for all $t\ge 0$ .

In particular $e+r^{\mathrm{<mo>\; *\; </mo>}}\le r^0$ that leads to a contradiction, necessarily we must have $r^{*}<r^0$ . By (13) and (14), finally obtain almost surly

$\underset{t\to \infty }{\lim \sup }{S}_2\left(t\right)=\frac{\Lambda +{\eta }_1{r}^0}{{\mu }_1+{\theta }_1}=s_1^0,\underset{t\to \infty }{\lim \sup }{S}_2\left(t\right)=\frac{{\theta }_1{s}_1^0+{\eta }_2{r}^0}{{\mu }_2+{\theta }_2}=s_2^0,$

$\underset{t\to \infty }{\lim \sup }{S}_3\left(t\right)=\frac{{\theta }_2{s}_2^0+{\eta }_2{r}^0}{{\mu }_3+{\theta }_3}=s_3^0,\underset{t\to \infty }{\lim \sup }V\left(t\right)=\frac{{\theta }_3{s}_3^0}{{\mu }_4+{\theta }_4}=v^0.$

□

The following corollary which is necessary to establish our stability result is a consequence of the previous result of Theorem (2).

Corollary 2. Assume that the assumptions A₁-A₂ and the condition (6) in Theorem (2) are satisfied. Then, the system (5) has a unique global solution $y\left(t\right)=\left(y_1\left(t\right)\mathrm{,}{y}_2\left(t\right)\mathrm{,}{y}_3\left(t\right)\mathrm{,}{y}_4\left(t\right)\mathrm{,}{y}_5\left(t\right)\mathrm{,}{y}_6\left(t\right)\right)$ for any initial value $\psi =\left({\psi }_1\mathrm{,}{\psi }_2\mathrm{,}{\psi }_3\mathrm{,}{\psi }_4\mathrm{,}{\psi }_5\mathrm{,}{\psi }_6\right)$ such that $\psi +E^0\in {\mathcal{C}}_{\left[-h\mathrm{,0}\right]}\left({ℝ}_{+}^6\right)\cap {\mathcal{F}}_0$ . Moreover, for any $\theta \in \left[\mathrm{0,}h\right]$ such that for $i=1,2,3,4$ ,

$-s_i^0\le {\psi }_i\left(\theta \right)\le \mathrm{0,}{\psi }_5\left(\theta \right)\ge 0\text{and}-r^0\le {\psi }_6\left(\theta \right)\le 0\text{a}\text{.s}\mathrm{.,}$

we have

for all $t\ge 0$ , $y_i\left(t\right)<0$ a.s. $i=1,2,3,4,6$ .

Proof. Let us put ${\phi }_1={\psi }_1+s_1^0$ , ${\phi }_2={\psi }_1+s_1^0$ , ${\phi }_3={\psi }_1+s_1^0$ , ${\phi }_4={\psi }_1+s_1^0$ , ${\phi }_5={\psi }_1+s_1^0$ , ${\phi }_6={\psi }_1+s_1^0$ . It’s easy to see that the system (1) with initial condition $\phi =\left({\phi }_1\mathrm{,}{\phi }_2\mathrm{,}{\phi }_3\mathrm{,}{\phi }_4\mathrm{,}{\phi }_5\mathrm{,}{\phi }_6\right)\in {\mathcal{C}}_{\left[-h\mathrm{,0}\right]}\left({ℝ}_{+}^6\right)\cap {\mathcal{F}}_0$ is equivalent to the system (5) with the initial condition $\psi =\left({\psi }_1\mathrm{,}{\psi }_2\mathrm{,}{\psi }_3\mathrm{,}{\psi }_4\mathrm{,}{\psi }_5\mathrm{,}{\psi }_6\right)$ such that $\psi +E^0\in {\mathcal{C}}_{\left[-h\mathrm{,0}\right]}\left({ℝ}_{+}^6\right)\cap {\mathcal{F}}_0$ . Therefore, in view of the Corollary 1, the system (5) has a unique solution global solution. Moreover, the condition, $\theta \in \left[\mathrm{0,}h\right]$ , for $i=1,2,3,4$ ,

$-s_i^0\le {\psi }_i\left(\theta \right)\le \mathrm{0,}{\psi }_5\left(\theta \right)\ge 0\text{and}-r^0\le {\psi }_6\left(\theta \right)\le 0\text{a}\text{.s}\mathrm{.,}$

imply that $0\le {\phi }_i\left(\theta \right)\le s_i^0$ , ${\phi }_5\left(\theta \right)\ge 0$ and $0\le {\phi }_6\left(\theta \right)\le r^0$ a.s. That is

$\begin{array}{l}{\phi }_1\left(\theta \right)+{\phi }_2\left(\theta \right)+{\phi }_3\left(\theta \right)+{\phi }_4\left(\theta \right)+{\phi }_5\left(\theta \right)+{\phi }_6\left(\theta \right)<s_1^0+s_2^0+s_3^0+s_4^0+r^0\text{a}\text{.s}\text{.},\\ \text{forall}\theta \in \left[-h\mathrm{,0}\right]\end{array}$

In view of the Theorem 2, the solution

$y\left(t\right)+E^0=\left(y_1\left(t\right)+s_1^0\mathrm{,}{y}_2\left(t\right)+s_2^0\mathrm{,}{y}_3\left(t\right)+s_3^0\mathrm{,}{y}_4\left(t\right)+s_4^0\mathrm{,}{y}_5\left(t\right)\mathrm{,}{y}_6\left(t\right)+r^0\right)$

of the system (1) is such that

$\begin{array}{l}{y}_1\left(t\right)+s_1^0<s_1^0,y_2\left(t\right)+s_2^0<s_2^0,y_3\left(t\right)+s_3^0<s_3^0,y_4\left(t\right)+s_4^0<v^0\text{and}\\ y_6\left(t\right)+s_6^0<r^0\text{a}\text{.s}.,\text{for}\text{all}t\ge 0.\end{array}$

Therefore for all $t\ge 0$ , $y_i\left(t\right)<0$ a.s. $i=1,2,3,4,6$ . □

Now, we establish a stability result for the trivial solution $E_y^0=\left(\mathrm{0,0,0,0,0,0}\right)$ of the model (5) by combining a stochastic Lyapunov technique and martingale convergence theory (see [42] [46] ).

Theorem 3. Let’s assume that $R_0=\frac{{\displaystyle {\sum }_{k=1}^4}{\beta }_k{s}_k^0}{{\mu }_I+\gamma }<1$ , then the disease-free equilibrium $E_y^0=\left(\mathrm{0,0,0,0,0,0}\right)$ of model (5) is globally asymptotically stable almost surely for any initial condition $\psi =\left({\psi }_1\mathrm{,}{\psi }_2\mathrm{,}{\psi }_3\mathrm{,}{\psi }_4\mathrm{,}{\psi }_5\mathrm{,}{\psi }_6\right)$ such that $\psi +E^0\in {\mathcal{C}}_{\left[-h\mathrm{,0}\right]}\left({ℝ}_{+}^6\right)\cap {\mathcal{F}}_0$ .

The proof of this Theorem requires the useful non-negative semimartingale convergence result ( [43] Theorem 1.3.9, p.14).

Lemma 4. Let A₁ and A₂ be two continuous adapted increasing processes on $t\ge 0$ with $A_1\left(0\right)=A_2\left(0\right)=0$ a.s. Let M be a real-valued continuous local martingale with $M\left(0\right)=0$ a.s. Let Z be a non-negative measurable random variable such that $\mathbb{E}\left(Z\right)<\infty$ . Define

$X\left(t\right)=Z+A_1\left(t\right)-A_2\left(t\right)+M\left(t\right)\text{for}t\ge 0.$

If X is non-negative, then

$\left\{\underset{t\to \infty }{\lim }{A}_1\left(t\right)<\infty \right\}\subset \left\{\underset{t\to \infty }{\lim }X\left(t\right)<\infty \right\}\cap \left\{\underset{t\to \infty }{\lim }{A}_2\left(t\right)<\infty \right\}\text{a}\text{.s}.,$

where $E\subset F$ a.s., means $P\left(E\cap F^c\right)=0$ .

In particular, if ${\lim }_{t\to \infty }{A}_1\left(t\right)<\infty$ a.s., then

$\underset{t\to \infty }{\lim }X(t)<\infty ,\underset{t\to \infty }{\lim }{A}_2(t)<\infty \text{and}\underset{t\to \infty }{\lim }\left|M\left(t\right)\right|<\infty \text{a}\text{.s}$

That is, all of the processes X, A₂, and M converge to finite random variables.

Proof of Theorem 3. We will first establish separately the almost sure asymptotic stability of the trivial solution of the component $y_5$ of the system (5), then we deduce that the trivial solution $E_y^0$ of the system (5) is asymptotically stable almost surly.

For any $\left(x_1\mathrm{,}{x}_2\mathrm{,}{x}_3\mathrm{,}{x}_4\mathrm{,}{x}_5\mathrm{,}{x}_6\right)\in {ℝ}^5$ , let us define $P{r}_5\left(x_1\mathrm{,}{x}_2\mathrm{,}{x}_3\mathrm{,}{x}_4\mathrm{,}{x}_5\mathrm{,}{x}_6\right)=x_5$ . So the component $y_5$ of the system (5) is described by the following equation

$\begin{array}{c}\text{d}{y}_5\left(t\right)=\left[{\displaystyle \underset{k=1}{\overset{4}{\sum }}}{\beta }_k\left(y_k+s_k^0\right)H\left(y_5,t\right)-\left({\mu }_I+\gamma \right)y_5\left(t\right)\right]\text{d}t\\ +y_5\left(t\right){\displaystyle \underset{k=1}{\overset{4}{\sum }}}{\sigma }_k\left(y_k+s_k^0\right)H\left(y_5,t\right)\text{d}W\left(t\right),\end{array}$ (16)

with initial condition $P{r}_5\circ \psi ={\psi }_5$ for any $\psi =\left({\psi }_1\mathrm{,}{\psi }_2\mathrm{,}{\psi }_3\mathrm{,}{\psi }_4\mathrm{,}{\psi }_5\mathrm{,}{\psi }_6\right)$ such that $\psi +E^0\in {\mathcal{C}}_{\left[-h\mathrm{,0}\right]}\left({ℝ}_{+}^6\right)\cap {\mathcal{F}}_0$ .

Let us consider the Lyapunov functional

$\bar{V}\left(\phi \right)=P{r}_5\circ \phi (\; 0\; )$

where $P{r}_5\circ \psi ={\psi }_5\left(\theta \right)=y_4\left(t+\theta \right)$ , $\theta \in \left[-h\mathrm{,0}\right]$ .

In view of (3) and Corollary (2), we get that

$\begin{array}{c}\mathcal{L}\bar{V}({\phi }_2)={\displaystyle \underset{k=1}{\overset{4}{\sum }}}{\beta }_k\left(y_k+s_k^0\right){\displaystyle {\int }_0^h}g\left(\tau \right)\mathcal{K}\left(y_5\left(t-\tau \right)\right)\text{d}\tau -\left({\mu }_I+\gamma \right)y_4\left(t\right)\\ \le {\displaystyle \underset{k=1}{\overset{4}{\sum }}}{\beta }_k{s}_k^0{\displaystyle {\int }_0^h}g\left(\tau \right)\mathcal{K}\left(y_5\left(t-\tau \right)\right)\text{d}\tau -\left({\mu }_I+\gamma \right)y_4\left(t\right).\end{array}$

In view of Theorems 3 in [47] , if ${\displaystyle {\sum }_{k=1}^4}{\beta }_k{s}_k^0<{\mu }_I+\gamma$ , we obtain that

${\lim }_{t\to \infty }\frac{1}{t}\ln \left(y_4\left(t\right)\right)<-p$ , where p is a positive constant. That is, if $R_0\le 1$ there exists two positive constants p₁ and p₂ such that

$y_4\left(t\right)<p_1\exp \left(-p_{2}t\right)\text{for}\text{any}t\ge 0.$ (17)

In this step we will prove that ${\lim }_{t\to \infty }Y\left(t\right)=0$ , where $Y\left(t\right)=-y_1\left(t\right)-y_2\left(t\right)-y_3\left(t\right)-y_4\left(t\right)-y_6\left(t\right)$ . From Corollary 2, it is clear that $Y\left(t\right)\ge 0$ almost surly.

Let put $\mathcal{M}\left(t\right)=\sigma {\displaystyle {\int }_0^t}{y}_5\left(v\right){\displaystyle {\sum }_{k=1}^4}\left(y_k\left(v\right)+s_k^0\right)H\left(y_5\left(v-\cdot \right)\right)\text{d}W\left(v\right)$ . According to Corollary 1, the quantity $y_5\left(v\right){\displaystyle {\sum }_{k=1}^4}\left(y_k\left(v\right)+s_k^0\right)H\left(y_5\left(v-\cdot \right)\right)$ is bounded for all $v\ge 0$ , consequently $\mathcal{M}\left(t\right)$ is a local martingale.

From the four first equations of the model (5) and the fact that ${\mu }_1\le \min \left\{{\mu }_2\mathrm{,}{\mu }_3\mathrm{,}{\mu }_4\mathrm{,}{\mu }_I\mathrm{,}{\mu }_R\right\}$ , we obtain

$Y\left(t\right)\le Y\left(0\right)+{\displaystyle \underset{k=1}{\overset{4}{\sum }}}{\beta }_k{s}_k^0{\displaystyle {\int }_0^t}H\left(y_5,v\right)\text{d}v-{\mu }_1{\displaystyle {\int }_0^t}Y\left(v\right)\text{d}v+\mathcal{M}\left(t\right).$

In view of Corollary 1, Hölder inequality and (17), we obtain

$\begin{array}{c}\underset{t\to \infty }{\lim }{\displaystyle \underset{k=1}{\overset{4}{\sum }}}{\beta }_k{s}_k^0{\displaystyle {\int }_0^t}H\left(y_4,s\right)\text{d}s\le \underset{t\to \infty }{\lim }{\displaystyle \underset{k=1}{\overset{4}{\sum }}}{\beta }_k{s}_k^0{\displaystyle {\int }_0^t}{\displaystyle {\int }_{s-h}^s}g\left(s-u\right)y_5\left(u\right)\text{d}u\text{d}s\\ \le \underset{t\to \infty }{\lim }h{\displaystyle \underset{k=1}{\overset{4}{\sum }}}{\beta }_k{s}_k^0{\displaystyle {\int }_0^t}\underset{u\in \left[s-h,s\right]}{\sup }{y}_5\left(u\right)\text{d}s\\ \le h{p}_1{\displaystyle \underset{k=1}{\overset{4}{\sum }}}{\beta }_k{s}_k^0{\text{e}}^{p_{2}h}\left(\underset{t\to \infty }{\lim }{\displaystyle {\int }_0^t}{\text{e}}^{-p_{2}s}\text{d}s\right)<\infty .\end{array}$

Therefore, by virtue of Lemma 4, we get

$\underset{t\to \infty }{\lim }Y\left(t\right)<\infty \text{and}\underset{t\to \infty }{\lim }{\displaystyle {\int }_0^t}Y\left(s\right)\text{d}s<\infty a.s.$

In accordance with Theorems 2, $Y\left(t\right)$ is positive for all $t\ge 0$ . Therefore, we get

$\underset{t\to \infty }{\lim }{\displaystyle {\int }_0^t}Y\left(s\right)\text{d}s={\displaystyle {\int }_0^{\infty }}Y\left(s\right)\text{d}s<\infty .$ (18)

Assume that $Y\left(t\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$ ,

$\lim {\inf }_{t\to \infty }Y\left(t,\omega \right)=\tau \left(\omega \right)>0.$

Then, there exists a $T>0$ such that $Y\left(t,\omega \right)>\frac{1}{2}\tau \left(\omega \right)$ for all $t\ge T$ . It follows that

$\begin{array}{c}\underset{t\to \infty }{\lim }{\displaystyle {\int }_0^t}Y\left(s,\omega \right)\text{d}s={\displaystyle {\int }_0^T}Y\left(s,\omega \right)\text{d}s+{\displaystyle {\int }_T^{\infty }}Y\left(s,\omega \right)\text{d}s\\ \ge {\displaystyle {\int }_T^{\infty }}Y\left(s,\omega \right)\text{d}s=\infty .\end{array}$

Therefore, ${\Omega }_1\subset {\Omega }_2$ , where ${\Omega }_2=\left\{\omega ,{\displaystyle {\int }_T^{\infty }}Y\left(s,\omega \right)\text{d}s=\infty \right\}$ . Hence $P\left({\Omega }_2\right)>0$ , which contradicts (18). So, we have

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

Finally, we have proved that, when $t\to \infty$ ,

$\left(y_1\left(t\right)\mathrm{,}{y}_2\left(t\right)\mathrm{,}{y}_3\left(t\right)\mathrm{,}{y}_4\left(t\right)\mathrm{,}{y}_5\left(t\right)\right)\to \left(\mathrm{0,0,0,0,0}\right)$ a.s. □

Following the result of Theorem 3 and the change of variable (4), we deduce the following corollary which gives us the almost sure stability of the disease-free equilibrium $E^0$ of the model (1) under the condition $R_0<1$ .

Corollary 3. If $R_0=\frac{{\displaystyle {\sum }_{k=1}^4}{\beta }_k{s}_k^0}{{\mu }_I+\gamma }<1$ , then the disease-free equilibrium

$E^0=\left(s_1^0\mathrm{,}{s}_2^0\mathrm{,}{s}_3^0\mathrm{,}{s}_3^0\mathrm{,0,}{r}^0\right)$ of model (1) is stable almost surly for any initial condition $\phi =\left({\phi }_1\mathrm{,}{\phi }_2\mathrm{,}{\phi }_3\mathrm{,}{\phi }_4\mathrm{,}{\phi }_5\mathrm{,}{\phi }_6\right)\in {\mathcal{H}}_b^5\left(\left[-h\mathrm{,0}\right]\right)$ .

6. Numerical Simulation and Commentary

In this section, we give an illustration of the stability result and the effect of noise intensity in the model by numerical simulation. We use the Euler-Maruyama method (see e.g [48] ) to simulate the sample paths of the model (1) with $G\left(x\right)=x/\left(1+x\right)$ (i.e. $G\left(x\right)\le 1$ ) for all $x\in \left[\mathrm{0,}\infty \right)$ and $f\left(s\right)=1/h$ for all $s\in \left[\mathrm{0,}h\right]$ , null otherwise. Ten sample paths of the stochastic model (1) under the condition $R_0<1$ given in Figure 2, we effectively observe the stability of the disease-free equilibrium $E_0$ . In Figure 3, we represent a sample path of model (1) with $R_0>1$ , in this case, the disease persists in the population ( $I\left(t\right)>0$ , $\forall t\ge 0$ ). We therefore see that these numerical simulations (Figure 2 and Figure 3) agree well with the analytical results of theorem 3. Finally, in Figure 4, we illustrate model (1) under the condition $R_0>1$ and with higher noise intensity compared to the case of Figure 2. Thus, we observe that the increase in noise intensity increases the intensity of fluctuations in the model with larger extreme values.

Figure 2. Ten (3) sample paths of the stochastic SVIRS epidemic models (1) with $G\left(x\right)=x/\left(1+x\right)$ . The initial values are: $S_1\left(\theta \right)=80$ , $S_2\left(\theta \right)=15$ , $S_3\left(\theta \right)=1$ , $V\left(\theta \right)=1$ , $I\left(\theta \right)=4$ , $R\left(\theta \right)=0$ for $\theta \in \left[-\mathrm{5,0}\right]$ . The values of the parameters are given by: $h=5$ , $\Lambda =30$ , ${\mu }_1={\mu }_3={\mu }_2=0.2$ , ${\theta }_1=0.08$ , ${\theta }_2=0.18$ , ${\theta }_3=0.6$ , ${\theta }_4=0.85$ , ${\beta }_1=0.006$ , ${\beta }_2=0.005$ , ${\beta }_3=0.0045$ , ${\beta }_4=0.035$ , $\gamma =0.8$ , $\sigma =0.1$ , ${\eta }_1={\eta }_2={\eta }_3=0.01$ . The conditions $R_0=0.7940<1$ of the theorem 3 is checked. We observe then that the disease-free equilibrium E₀ is asymptotically stable almost surly.

Figure 3. Sample paths of the stochastic SVIRS epidemic models (1) with $G\left(x\right)=x/\left(1+x\right)$ . The initial values are: $S_1\left(\theta \right)=101$ , $S_2\left(\theta \right)=20$ , $S_3\left(\theta \right)=3$ , $V\left(\theta \right)=2$ , $I\left(\theta \right)=1$ , $R\left(\theta \right)=8$ for $\theta \in \left[-\mathrm{5,0}\right]$ . The values of the parameters are given by: $h=5$ , $\Lambda =30$ , ${\mu }_1={\mu }_3={\mu }_2=0.2$ , ${\theta }_1=0.08$ , ${\theta }_2=0.18$ , ${\theta }_3=0.6$ , ${\theta }_4=0.85$ , ${\beta }_1=0.04$ , ${\beta }_2=0.037$ , ${\beta }_3=0.035$ , ${\beta }_4=0.035$ , $\gamma =0.0035$ , $\sigma =0.0035$ , ${\eta }_1={\eta }_2={\eta }_3=0.01$ . In this case $R_0=18.5002$ .

Figure 4. Sample paths of the stochastic SVIRS epidemic models (1) with $G\left(x\right)=x/\left(1+x\right)$ . The initial values are: $S_1\left(\theta \right)=101$ , $S_2\left(\theta \right)=20$ , $S_3\left(\theta \right)=3$ , $V\left(\theta \right)=2$ , $I\left(\theta \right)=1$ , $R\left(\theta \right)=8$ for $\theta \in \left[-\mathrm{5,0}\right]$ . The values of the parameters are given by: $h=5$ , $\Lambda =30$ , ${\mu }_1={\mu }_3={\mu }_2=0.2$ , ${\theta }_1=0.16$ , ${\theta }_2=0.1$ , ${\theta }_3=0.2$ , ${\theta }_4=0.3$ , ${\beta }_1=0.04$ , ${\beta }_2=0.037$ , ${\beta }_3=0.035$ , ${\beta }_4=0.035$ , $\gamma =0.1$ , $\sigma =0.0098$ , ${\eta }_1={\eta }_2={\eta }_3=0.01$ . In this case $R_0=18.1255$ .

7. Conclusion and Perspective

We considered a stochastic epidemic SIRS model represented by a delayed stochastic differential equation to describe the spread of COVID-19 in a population where susceptibility to the disease varies across age groups and where a fraction of people over 55 years of age are vaccinated. First, we established the consistency of the model, i.e. the existence and uniqueness of the global and positive solution (see Corollary 1). Then, we have established the almost sure asymptotic stability of the $E^0$ disease-free equilibrium of the model when $R_0<1$ (see Theorem 3). The work performed in this paper could be improved by taking into account time-related parameters, which would allow taking into account seasonal effects or times of the year favoring large gatherings, where contact rates can increase. We can also imagine the possibility of using delay-dependent contact rates. This is relevant in certain situations where supporting measures are required, such as Contact tracing of an infected person can reduce the number of infectious contacts over time, so an increase in the length of the latency period can reduce the contact rate. Finally, given the large amount of data on COVID-19 globally, we can improve this work by adding parameter estimation methods to adapt the model to reality.

Acknowledgement

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]	Choi, Y., Kim, H.J. and Lee, Y. (2022) Economic Consequences of the COVID-19 Pandemic: Will It Be a Barrier to Achieving Sustainability? Sustainability, 14, Article 1629. https://doi.org/10.3390/su14031629
[2]	Pak, A., Adegboye, O.A., Adekunle, A.I., et al. (2020) Economic Consequences of the COVID-19 Outbreak: The Need for Epidemic Preparedness. Frontiers in Public Health, 8, Article 546036. https://doi.org/10.3389/fpubh.2020.00241
[3]	Aragie, E., Taffesse, A.S. and Thurlow, J. (2020) The Short-Term Economywide Impacts of COVID-19 in Africa: Insights from Ethiopia. African Development Review, 33, S152-S164. https://doi.org/10.1111/1467-8268.12519
[4]	Anderson, R.M. and May, R.M. (1982) Directly Transmitted Infections Diseases: Control by Vaccination. Science, 215, 1053-1060. https://doi.org/10.1126/science.7063839
[5]	Anderson, R.M. and May, R.M. (1985) Vaccination and Herd Immunity to Infectious Diseases. Nature, 318, 323-329. https://doi.org/10.1038/318323a0
[6]	Hethcote, H.W. (2000) The Mathematics of Infectious Diseases. SIAM Review, 42, 599-653. https://doi.org/10.1137/S0036144500371907
[7]	Matrajt, L., Eaton, J., Leung, T. and Brown, E.R. (2021) Vaccine Optimization for COVID-19: Who to Vaccinate First? Science Advances, 7, eabf1374. https://doi.org/10.1126/sciadv.abf1374
[8]	Zhao, Z.Y., Niu, Y., Luo, L., et al. (2021) The Optimal Vaccination Strategy to Control COVID-19: A Modeling Study in Wuhan City, China. Infectious Diseases of Poverty, 10, Article No. 140. https://doi.org/10.1186/s40249-021-00922-4
[9]	Griette, Q., Magal, P. and Seydi, O. (2020) Unreported Cases for Age Dependent COVID-19 Outbreak in Japan. Biology, 9, Article 132. https://doi.org/10.1101/2020.05.07.20093807
[10]	Anderson, R. and May, R. (1991) Infectious Disease of Humans, Dynamics and Control. Oxford University Press, Oxford. https://doi.org/10.1093/oso/9780198545996.001.0001
[11]	Issanov, A., Amanbek, Y., Abbay, A., et al. (2021) COVID-19 Outbreak in Post-Soviet States: Modeling the Best and Worst Possible Scenarios. Electronic Journal of General Medicine, 17, em256. https://doi.org/10.1101/2020.04.19.20071704
[12]	Beretta, E. and Takeuchi, Y. (1995) Global Stability of an SIR Epidemic Model with Time Delays. Journal of Mathematical Biology, 33, 250-260. https://doi.org/10.1007/BF00169563
[13]	Coulibaly, M. and N’zi, M. (2021) A Stochastic Model with Jumps for the COVID-19 Epidemic in the Greater Abidjan Region during Public Health Measures. Journal of Infectious Diseases and Epidemiology, 7, Article 196. https://doi.org/10.23937/2474-3658/1510196
[14]	Liu, Z., Magal, P., Seydi, O. and Webb, G. (2020) Understanding Unreported Cases in the COVID-19 Epidemic Outbreak in Wuhan, China, and the Importance of Major Public Health Interventions. Biology, 7, Article 50. https://doi.org/10.3390/biology9030050
[15]	Otunuga, O.M. (2021) Estimation of Epidemiological Parameters for COVID-19 Cases Using a Stochastic SEIRS Epidemic Model with Vital Dynamics. Results in Physics, 28, Article ID: 104664. https://doi.org/10.1016/j.rinp.2021.104664
[16]	Ndaïrou, F., Area, I., Nieto, J.J. and Torres, D.F.M. (2020) Mathematical Modeling of COVID-19 Transmission Dynamics with a Case Study of Wuhan. Chaos, Solitons & Fractals, 135, Article ID: 109846. https://doi.org/10.1016/j.chaos.2020.109846
[17]	WHO-China (2020) Report of the WHO-China Joint Mission on Coronavirus Disease 2019 (COVID-19). https://www.who.int/docs/default-source/coronaviruse/who-china-joint-mission-on-covid-19-final-report.pdf
[18]	Yang, Z., Zeng, Z., et al. (2020) Modified SEIR and AI Prediction of the Epidemics Trend of COVID-19 in China under Public Health Interventions. Journal of Thoracic Disease, 12, 165-174. https://doi.org/10.21037/jtd.2020.02.64
[19]	Xin, H., Li, Y., Wu, P., Li, Z., et al. (2022) Estimating the Latent Period of Coronavirus Disease 2019 (COVID-19). Clinical Infectious Diseases, 74, 1678-1681. https://doi.org/10.1093/cid/ciab746
[20]	Beretta, E., Kolmanovskii, V. and Shaikhet, L. (1998) Modified SEIR and AI Prediction of the Epidemics Trend of COVID-19 in China under Public Health Interventions. Mathematics and Computers in Simulation, 45, 269-277. https://doi.org/10.1016/S0378-4754(97)00106-7
[21]	Lui, Q., Jiang, D., Shi, N., Hayat, T. and Alsaedi, A. (2017) Asymptotic Behavior of Stochastic Multi-Group Epidemic Models with Distributed Delays. Physica A: Statistical Mechanics and its Applications, 467, 527-541. https://doi.org/10.1016/j.physa.2016.10.034
[22]	N’zi, M. and Yattara, I. (2018) The Dynamics of a Delayed Deterministic and Stochastic Epidemic SIR Models with a Saturated Incidence Rate. Random Operators and Stochastic Equation, 26, 235-245. https://doi.org/10.1515/rose-2018-0021
[23]	Zine, H., Boukhouima, A., et al. (2020) A Stochastic Time-Delayed Model for the Effectiveness of Moroccan COVID-19 Deconfinement Strategy. Mathematical Modelling of Natural Phenomena, 15, Article No. 50.
[24]	Rihan, F., Alsakaji, H. and Rajivganthi, C. (2020) Stochastic SIRC Epidemic Model with Time-Delay for COVID-19. Advances in Difference Equations, 502, 235-245. https://doi.org/10.1186/s13662-020-02964-8
[25]	Rihan, F. and Alsakaji, H. (2021) Dynamics of a Stochastic Delay Differential Model for COVID-19 Infection with Asymptomatic Infected and Interacting People: Case Study in the UAE. Results in Physics, 28, Article ID: 104658. https://doi.org/10.1016/j.rinp.2021.104658
[26]	Alsakaji, H., Rihan, F. and Hashish, A. (2022) Dynamics of a Stochastic Epidemic Model with Vaccination and Multiple Time-Delays for COVID-19 in the UAE. Complexity, 2022, Article ID: 4247800. https://doi.org/10.1155/2022/4247800
[27]	Ikram, R., Khan, A., et al. (2012) Extinction and Stationary Distribution of a Stochastic COVID-19 Epidemic Model with Time-Delay. Computers in Biology and Medicine, 141, Article ID: 105115. https://doi.org/10.1016/j.compbiomed.2021.105115
[28]	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
[29]	Cooke, K.L. (1979) Stability Analysis for a Vector Disease Model. Rocky Mountain Journal of Mathematics, 9, 31-42. https://doi.org/10.1216/RMJ-1979-9-1-31
[30]	Han, M., Ma, Z. and Zhen, J. (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
[31]	Zhang, T. and Teng, Z. (2008) Global Behavior and Permanence of SIRS Epidemic Model with Time Delay. Nonlinear Analysis: Real World Applications, 9, 1409-1424. https://doi.org/10.1016/j.nonrwa.2007.03.010
[32]	Wang, J., Takeuchi, Y. and Liu, S. (2012) A Multi-Group SVEIR Epidemic Model Distributed Delay and Vaccination. International Journal of Biomathematics, 5, Article ID: 1260001. https://doi.org/10.1142/S1793524512600017
[33]	Sekiguchia, M., Ishiwata, E. and Ishiwata, E. (2010) Global Dynamics of a Discretized SIRS Epidemic Model with Time Delay. Journal of Mathematical Analysis and Applications, 371, 195-202. https://doi.org/10.1016/j.jmaa.2010.05.007
[34]	Hyman, J. and Li, J. (2010) Differential Susceptibility Epidemic Models. Journal of Mathematical Biology, 50, 626-644. https://doi.org/10.1007/s00285-004-0301-7
[35]	Hyman, J. and Li, J. (2009) Epidemic Models with Differential Susceptibility and Staged Progression and Their Dynamics. Mathematical Biosciences and Engineering, 6, 321-332. https://doi.org/10.3934/mbe.2009.6.321
[36]	Bonzi, B., Fall, A.A., Iggidr, A. et al. (2011) Global Dynamics of a Discretized SIRS Epidemic Model with Time Delay. Journal of Mathematical Biology, 62, 39-64. https://doi.org/10.1007/s00285-010-0327-y
[37]	Zhang, X.B., Huo, H.F., Sun, X.K., et al. (2010) The Differential Susceptibility SIR Epidemic Model with Time Delay and Pulse Vaccination. Journal of Applied Mathematics and Computing, 34, 287-298. https://doi.org/10.1007/s12190-009-0321-y
[38]	Goldstein, E., Lipsitch, M. and Cevik, M. (2021) On the Effect of Age on the Transmission of SARS-CoV-2 in Households, Schools, and the Community. The Journal of Infectious Diseases, 223, 362-369. https://doi.org/10.1101/2020.07.19.20157362
[39]	WHO (2022) Rapport de Situation COVID-19 au Mali, 30 mai au 05 juin 2022/No 202. https://reliefweb.int/node/3854229
[40]	Fields, R., Humphrey, L., et al. (2021) Age-Stratified Transmission Model of COVID-19 in Ontario with Human Mobility during Pandemic’s First Wave. Heliyon, 7, E07905. https://doi.org/10.1016/j.heliyon.2021.e07905
[41]	Zhao, Z.Y., Zhu, Y.Z. and Xu, J.W., et al. (2020) A Five-Compartment Model of Age-Specific Transmissibility of SARS-CoV-2. Infectious Diseases of Poverty, 9, Article No. 117. https://doi.org/10.1186/s40249-020-00735-x
[42]	Dalal, N., Greenhalgh, D. and Mao, X. (2020) 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
[43]	Mao, X. (2007) Stochastic Differential Equations and Applications. Horwood Publishing Limited. https://doi.org/10.1533/9780857099402
[44]	Fu, X.M. (2019) On Invariant Measures and the Asymptotic Behavior of a Stochastic Delayed SIRS Epidemic Model. Physica A: Statistical Mechanics and Its Applications, 523, 1008-1023. https://doi.org/10.1016/j.physa.2019.04.181
[45]	Ikeda, N. and Watanabe, S. (1977) A Comparison Theorem for Solutions of Stochastic Differential Equations and Its Applications. Osaka Journal of Mathematics, 14, 619-633.
[46]	Lahrouz, A., Omari, L. and Omari, L. (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
[47]	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
[48]	Mao, X. (2003) A Stochastic Model of AIDS and Condom Use. LMS Journal of Computation and Mathematics, 6, 141-161. https://doi.org/10.1112/S1461157000000425