Loading [MathJax]/jax/output/SVG/jax.js
Research article

Global stability of multi-group SEIQR epidemic models with stochastic perturbation in computer network

  • Received: 06 February 2023 Revised: 03 April 2023 Accepted: 24 April 2023 Published: 29 May 2023
  • In this paper, a class of multi-group SEIQR models with random perturbation in computer network is investigated. The existence and uniqueness of global positive solution with any positive initial value are obtained. The sufficient conditions on the asymptotic behavior of solutions around the disease-free equilibrium and endemic equilibrium of the corresponding deterministic model are established. Furthermore, the existence and uniqueness of stationary distribution are also obtained. Lastly, the analytical results are illustrated by the numerical simulations.

    Citation: Ramziya Rifhat, Kai Wang, Lei Wang, Ting Zeng, Zhidong Teng. Global stability of multi-group SEIQR epidemic models with stochastic perturbation in computer network[J]. Electronic Research Archive, 2023, 31(7): 4155-4184. doi: 10.3934/era.2023212

    Related Papers:

    [1] Xueyong Zhou, Xiangyun Shi . Stability analysis and backward bifurcation on an SEIQR epidemic model with nonlinear innate immunity. Electronic Research Archive, 2022, 30(9): 3481-3508. doi: 10.3934/era.2022178
    [2] Yanjiao Li, Yue Zhang . Dynamic behavior on a multi-time scale eco-epidemic model with stochastic disturbances. Electronic Research Archive, 2025, 33(3): 1667-1692. doi: 10.3934/era.2025078
    [3] Meng Gao, Xiaohui Ai . A stochastic Gilpin-Ayala nonautonomous competition model driven by mean-reverting OU process with finite Markov chain and Lévy jumps. Electronic Research Archive, 2024, 32(3): 1873-1900. doi: 10.3934/era.2024086
    [4] Pinglan Wan . Dynamic behavior of stochastic predator-prey system. Electronic Research Archive, 2023, 31(5): 2925-2939. doi: 10.3934/era.2023147
    [5] Miaomiao Gao, Yanhui Jiang, Daqing Jiang . Threshold dynamics of a stochastic SIRS epidemic model with transfer from infected individuals to susceptible individuals and log-normal Ornstein-Uhlenbeck process. Electronic Research Archive, 2025, 33(5): 3037-3064. doi: 10.3934/era.2025133
    [6] Wenhui Niu, Xinhong Zhang, Daqing Jiang . Dynamics and numerical simulations of a generalized mosquito-borne epidemic model using the Ornstein-Uhlenbeck process: Stability, stationary distribution, and probability density function. Electronic Research Archive, 2024, 32(6): 3777-3818. doi: 10.3934/era.2024172
    [7] Erhui Li, Qingshan Zhang . Global dynamics of an endemic disease model with vaccination: Analysis of the asymptomatic and symptomatic groups in complex networks. Electronic Research Archive, 2023, 31(10): 6481-6504. doi: 10.3934/era.2023328
    [8] Wenxuan Li, Suli Liu . Dynamic analysis of a stochastic epidemic model incorporating the double epidemic hypothesis and Crowley-Martin incidence term. Electronic Research Archive, 2023, 31(10): 6134-6159. doi: 10.3934/era.2023312
    [9] Jing Yang, Shaojuan Ma, Dongmei Wei . Dynamical analysis of SIR model with Gamma distribution delay driven by Lévy noise and switching. Electronic Research Archive, 2025, 33(5): 3158-3176. doi: 10.3934/era.2025138
    [10] Tao Zhang, Mengjuan Wu, Chunjie Gao, Yingdan Wang, Lei Wang . Probability of disease extinction and outbreak in a stochastic tuberculosis model with fast-slow progression and relapse. Electronic Research Archive, 2023, 31(11): 7104-7124. doi: 10.3934/era.2023360
  • In this paper, a class of multi-group SEIQR models with random perturbation in computer network is investigated. The existence and uniqueness of global positive solution with any positive initial value are obtained. The sufficient conditions on the asymptotic behavior of solutions around the disease-free equilibrium and endemic equilibrium of the corresponding deterministic model are established. Furthermore, the existence and uniqueness of stationary distribution are also obtained. Lastly, the analytical results are illustrated by the numerical simulations.



    As is well known, the internet world has brought great changes in the society. In reality, we know that cyber world is being threatened by the attack of malicious objects. Malicious object is a code that infects computer systems. There are different kinds of malicious objects such as: Worm, Virus, Trojan horse, etc., which differ according to the way they attack computer systems and the malicious actions they perform (see [1,2,3]). With the development of the computer network, malicious objects be widely spread through a network, through an online service, through shared computer software or through a mobile storage tool, and so on. Because of the similarity between the transmission of human infectious diseases and transmission of malicious objects in the computer network, some authors employ the epidemic models to describe the transmission of malicious objects in the cyber world (see [2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17]).

    Considering different contact patterns, different anti-virus software, or distinct number of contacts etc., it is more appropriate to divide individual hosts into groups in modeling epidemic disease. Therefore, it is reasonable to propose multi-group models to describe the transmission dynamics of malicious objects in heterogeneous host populations on computer network. At present, many scholars have focused their study on various forms of multi-group epidemic models (see [18,19,20,21,22,23]). They have also proved the global stability of the unique endemic equilibrium through Lyapunov function, which is one of the main mathematical challenges in analyzing multi-group models. Particularly, Wang et al. [23] proposed the following multi-group SEIQR epidemic model for describing the transmission of malicious objects in computer network

    {dSk(t)=[Λknj=1βkjSk(t)Ij(t)dSkSk]dt,dEk(t)=[nj=1βkjSk(t)Ij(t)(dEk+ϵk)Ek]dt,dIk(t)=[ϵkEk(dIk+αk+δk+γk)Ik]dt,dQk(t)=[δkIk(dQk+αk+μk)Qk]dt,dRk(t)=[γkIk+μkQkdRkRk]dt,1kn, (1.1)

    where the total network nodes are divided into n groups of nodes, n2 is an integer. Sk(t), Ek(t), Ik(t), Qk(t) and Rk(t) express the numbers of susceptible nodes, exposed (infected but not yet infectious) nodes, infectious nodes, quarantined nodes and recovered nodes at time t in the k-th group (1kn), respectively. The definitions of all parameters in model (1.1) are listed in Table 1. We assume that the parameters dSk, dEk, dIk, dQk, dRk and Λk are positive and the rest of parameters in model (1.1) ia nonnegative for all k. In particular, βkj=0 if there is no transmission of the disease between compartments Sk and Ij. In model (1.1), the basic reproduction number R0=ρ(M0), the spectral radius of matrix M0=(βkjϵkΛkdk(dEk+ϵk)(dIk+αk+δk+γk))n×n, is a threshold which completely determines the persistence or extinction of the disease. It is shown that, if R01, the disease-free equilibrium E0 is globally stable in the feasible region and the disease always dies out, and if R0>1, a unique endemic equilibrium E exists and is globally stable in the interior of the feasible region, and once the disease appears, it eventually persists at the unique endemic equilibrium level.

    Table 1.  Description of parameters in model (1.1).
    Symbol Description
    Λk influx of individuals into the kth group
    βkj transmission coefficient between compartments Sk and Ij
    dSk,dEk,dIk,dQk,dRk natural death rates of Sk,Ek,Ik,Qk,Rk compartments in the kth group
    ϵk the rate constant for nodes leaving the exposed class Ek for infective
    compartment in the kth group
    δk the rate constant for nodes leaving the infective compartment Ik
    for quarantine compartment in the kth group
    αk the disease related death rate(crashing of nodes due to the attack
    of malicious objects)constant in the compartments
    γk and μk the rates at which nodes recover temporarily after the runof anti-malicious
    software and return to recovered class R from compartments Ik and Qk
    in the kth group

     | Show Table
    DownLoad: CSV

    On the other hand, there exist uncertainties and random phenomena everywhere in nature [23,24,25,26,27]. Environmental noises are usually considered to be harmful, which will lead to the disorder of the dynamics [20,21]. Nevertheless, the noises also play a positive role in the dynamics of complex nonlinear systems, especially in interdisciplinary physical models and biomathematics models, such as noise induced resonances, noise enhanced stability (NES) and so on [22,23,24,28,29,30]. According to the noise source, the noises can be divided into the additive noise and the multiplicative noise. The former is not controlled by the system and can be directly introduced to the system, while the latter is related to system parameters and variables. The multiplicative noises can always ensure the nonnegativity of the solution. The two main peculiarities of the presence of the multiplicative noise are the presence of the absorbing barrier in zero population density and the phenomenon of the anomalous fluctuations [25,31]. The noise existing in biological systems is caused by environmental fluctuations, which is usually considered as the multiplicative white noise. For example, Caruso et al. [26] described the dynamic behavior of an ecosystem of two competing species by a stochastic Lotka-Volterra model with the multiplicative white noise. The multiplicative noise models the interaction between the environment and the species.

    For human disease related epidemics, the nature of epidemic growth and spread is random due to the unpredictability in person to person contacts. Because of environmental noises, the deterministic approach has some limitations in the mathematical modeling transmission of an infectious disease, several authors began to consider the effect of white noise on the computer network systems (see [23,24,25,26,27]).

    There are different approaches used in the literature to introduce random perturbations into population models, both from a mathematical and biological perspective (see [23,24,25,26,27,28,29,31]). One is to perturb the positive equilibria in order for making robust the equilibria of deterministic models. In this situation, the essence of the investigation using the approach is to check if the asymptotic stability of the positive equilibria of deterministic models can be preserved. For example, Wang et al. [23] investigated a multi-group SEIQR model with random perturbation around the positive equilibrium of corresponding deterministic model, which revealed that the stochastic stability of endemic equilibrium depends on the magnitude of the intensity of noise as well as the parameters involved within the model. The other important approach is with parameters perturbation. We find that there are many literatures on this approach, see [25,26,27] and the references cited therein. In epidemic models, the natural death rate and the disease transmission rate are two of the key parameters to disease transmission. And in the real situation, the natural death rate and the disease transmission rate always fluctuate around some average value due to continuous fluctuation in the environment. For example, El Ansari et al. [25] considered a stochastic version of model (1.1) with noises introduced in the rate at which nodes are crashed due to reasons other than the attacks of viruses and the transmission rate, and they proved the various conditions that control the extinction and stability of a nonlinear mathematical spread model with stochastic perturbations.

    We now turn to a continuous time SEIQRS model which takes random effects into account. In SEIQRS model (1.1), the natural death rate dXik, where 1kn and (X1,X2,X3,X4,X5)=(S,E,I,Q,R), is one of the key parameters to disease transmission. May [30] pointed out that all the parameters involved in the population model exhibit random fluctuation as the factors controlling them are not constant. And in the real situation, the natural death rate d always fluctuate around some average value due to continuous fluctuation in the environment. In this sense, dXik can seem as a random variable ˜dXik. More precisely, in [t,t+dt),

    ˜dXikdt=dXikdt+σikdBik(t),1kn,i=1,2,3,4,5,

    where Bik(t)(1kn,i=1,2,3,4,5) are the independent standard Brownian motion defined on the complete probability space (Ω,{Ft}t0,P) with a filtration {Ft}t0 satisfying the usual conditions, and σ2ik is the intensity of Bik(t). The reason of adopting σ2ik(1kn,i=1,2,3,4,5) as the intensity of the noise for the group Sk, Ek, Ik, Qk and Rk, respectively, is considering the difference between the group mobility response to infection risks. And then, in [t,t+dt), ˜dXikdt is normally distributed with mean E(˜dXikdt)=dXikdt and variance Var(˜dXikdt)=σ2idt. Due to Var(˜dXikdt)=σ2idt0 as dt0, this is a biologically reasonable assumption. Indeed this is a well-established way of introducing stochastic environmental noise into biologically realistic population dynamic models.

    Therefore, replace dXikdt in model (1.1) with ˜dXikdt=dXikdt+σikdBik(t)(1kn,i=1,2,3,4,5), and for simplicity, we replace ˜dXik with dXik again, then we can obtain the same SDE epidemic model as the following model (1.2) that is analog to its deterministic version model (1.1) by introducing stochastic perturbation terms to the growth equations of susceptible, infectious, recovered individuals to incorporate the effect of randomly fluctuating environments:

    {dSk=[Λknj=1βkjSk(t)Ij(t)dSkSk]dt+σ1kSkdB1k,dEk=[nj=1βkjSk(t)Ij(t)(dEk+ϵk)Ek]dt+σ2kEkdB2k,dIk=[ϵkEk(dIk+αk+δk+γk)Ik]dt+σ3kIkdB3k,dQk=[δkIk(dQk+αk+μk)Qk]dt+σ4kQkdB4k,dRk=[γkIk+μkQkdRkRk]dt+σ5kRkdB5k,1kn. (1.2)

    Throughout this paper, we always assume that model (1.2) is defined on a complete probability space (Ω,{Ft}t0,P) with a filtration {Ft}t0 satisfying the usual conditions (i.e., it is right continuous and F0 contain all P-null sets). Furthermore, we also always assume that the infection rate matrix B=(βkj)n×n in model (1.2) is irreducible.

    In this paper, we will study the asymptotic behavior of positive solutions of model (1.2) around the disease-free and endemic equilibria of corresponding deterministic model (1.1) in probability meaning by using the theory of graphs, Lyapunov functions method, Itˆo's formula and the theory of stochastic analysis. Then by using the theory of stationary distributions of stochastic process we will study the existence of stationary distribution of model (1.2).

    The paper is organized as follows. In Section 2, the criterion on the asymptotic behavior of positive solutions of model (1.2) around the disease-free equilibrium of the corresponding deterministic model is stated and proved. In Section 3, the sufficient condition the asymptotic behavior of positive solutions of model (1.2) around the endemic equilibrium of corresponding deterministic model and the existence of stationary distribution are stated and proved. In Section 4, we make some numerical simulations to illustrate our analytical results. Finally, in Section 5, we give a brief conclusion.

    We first give a lemma to show that for any positive initial value model (1.2) has a unique positive solution defined on [0,).

    Lemma 1. For any initial value in R5n+ model (1.2) has a unique positive solution defined for all t0 and the solution remain in R5n+ with probability one.

    This lemma can be easily proved by using the standard arguments as in [14,18] and with the help of Lyapunov function

    V(Sk,Ek,Ik,Qk,Rk,1kn)=nk=1[(SkaalogSkack)+(Ek1logEk)+(Ik1logIk)+(Qk1logQk)+(Rk1logRk)],

    where positive constant a satisfies amin{dIk+αknj=1βjk,k=1,2,,n}.

    For deterministic model (1.1), in [23] the authors have obtained that there is a disease-free equilibrium E0=(S01,0,0,0,0,S02,0,0,0,0,,S0n,0,0,0,0), where S0k=ΛkdSk, and if R01, then E0 is globally asymptotically stable, which means the disease will die out. Therefore, it is interesting to study the stability of disease-free equilibrium for controlling the spread of infectious disease. However, for stochastic model (1.2) there is not any disease-free equilibrium. Therefore, it is natural to ask how we can consider the disease will be extinct. In this section we mainly through estimating the asymptotic oscillation around equilibrium E0 of any positive solutions of stochastic model (1.2) to reflect whether the disease in stochastic model (1.2) will die out. We have the following result.

    Theorem 1. Assume that R01 and the following conditions hold

    dSk>σ21k,dIk+αk+δk+γk>12σ23k,dEk+ϵk>12σ22k,dQk+αk+μk>12σ24k,dRk>12σ25k,1kn. (2.1)

    Then for any positive solution (Sk(t),Ek(t),Ik(t),Qk(t),Rk(t),1kn) of model (1.2) one has

    lim supt1tEt0nk=1{Ak(Sk(r)ΛkdSk)2+BkE2k(r)+CkI2k(r)+DkQ2k(r)+FkR2k(r)}drnk=1(bak+1)(σ1kΛkdSk)2,

    where Ak=(dSkσ21k), Bk=14(dEk+ϵk12σ22k) and

    Ck=ck[2(dIk+αk+δk+γk)σ23k4ckϵ2k[(dEk+ϵk)12σ22k]]dkδ2kdQk+αk+μk12σ24kekγ2kdRk12σ25k,Dk=dk(dQk+αk+μk12σ24k),Fk=ek(dRk12σ25kekμ2kdk(dQk+αk+μk12σ24k)),

    and positive constants ak,dk,ck,ek(1kn) and b will be confirmed in the proof of the theorem.

    Proof. Let uk=SkΛkdSk,vk=Ek,wk=Ik,yk=Qk,zk=Rk(1kn), then model (1.2) becomes into

    {duk=[nj=1βkjuk(t)wj(t)nj=1βkjwj(t)ΛkdSkdSkuk]dt+σ1k(uk+ΛkdSk)dB1k,dvk=[nj=1βkjuk(t)wj(t)+nj=1βkjwj(t)ΛkdSk(dEk+ϵk)vk]dt+σ2kvkdB2k,dwk=[ϵkvk(dIk+αk+δk+γk)wk]dt+σ3kwkdB3k,dyk=[δkwk(dQk+αk+μk)yk]dt+σ4kykdB4k,dzk=[γkwk+μkykdRkzk]dt+σ5kzkdB5k.

    Since B=(βkj)n×n is irreducible, then M0 is also nonnegative and irreducible. Hence, by Lemma A.1 in [3], M0 has a positive left eigenvector η=(η1,η2,,ηn) such that

    (η1,η2,,ηn)ρ(M0)=(η1,η2,,ηn)M0. (2.2)

    Define a Lyapunov function as follows.

    V=V1+b(V2+V3)+V4+V5+V6

    with V1=12nk=1(uk+vk)2, V2=12nk=1aku2k, V3=nk=1ϵkηk(dEk+ϵk)(dIk+αk+δk+γk)(vk+dEk+ϵkϵkwk), V4=nk=1ckw2k, V5=nk=1dky2k and V6=nk=1ekz2k, where positive constants ak,ck,dk,ek(1kn) and b will be determined later. By Itˆo's formula, we get

    dV=LVdt+nk=1σ1k(uk+ΛkdSk)[(1+bak)uk+vk]dB1k+nk=1σ2kvk[uk+vk+bωkϵk(dEk+ϵk)(dIk+αk+δk+γk)]dB2k+nk=1σ3kwk[bωkdIk+αk+δk+γk+ckwk]dB3k+nk=1dkσ4ky2kdB4k+nk=1ekσ5kz2kdB5k (2.3)

    with LV=LV1+b(LV2+LV3)+LV4+LV5+LV6, where

    LV1=nk=1(uk+vk)[nj=1βkjuk(t)wj(t)nj=1βkjwj(t)ΛkdSkdSkuk+nj=1βkjuk(t)wj(t)+nj=1βkjwj(t)ΛkdSk(dEk+ϵk)vk]+nk=1[σ21k(uk+ΛkdSk)2+σ22kv2k]nk=1{(dSkσ21k)u2k+[dEk+ϵk12σ22k]v2k+(dSk+dEk+ϵk)ukvk(ΛkdSk)2σ21k}, (2.4)
    LV2=nk=1akuk[nj=1βkjuk(t)wj(t)nj=1βkjwj(t)ΛkdSkdSkuk]+nk=1akσ21ku2k+nk=1akσ21k(ΛkdSk)2nk=1ak[(dSkσ21k)u2k+nj=1βkjΛkdSkuk(t)wj(t)(σ1kΛkdSk)2] (2.5)

    and

    LV3=nk=1ωkϵk(dEk+ϵk)(dIk+αk+δk+γk)[nj=1βkjuk(t)wj(t)+nj=1βkjwj(t)ΛkdSk(dEk+ϵk)vk+ϵkvkdEk+ϵkϵk(dIk+αk+δk+γk)wk]nk=1nj=1βkjωkϵk(dEk+ϵk)(dIk+αk+δk+γk)uk(t)wj(t)nk=1ωkwk+nk=1nj=1βkjωkϵk(dEk+ϵk)(dIk+αk+δk+γk)ΛkdSkwj(t).

    Note from (2.2) that

    ηkwk+nk=1nj=1βkjωkϵk(dEk+ϵk)(dIk+αk+δk+γk)ΛkdSkwj(t)=(R01)ηw,

    where w=(w1,w2,,wn)T. If R01, then

    LV3nk=1nj=1βkjηkϵk(dEk+ϵk)(dIk+αk+δk+γk)uk(t)wj(t). (2.6)

    Furthermore, we also have

    LV4=nk=1ck[2(dIk+αk+δk+γk)σ23k]w2k+2nk=1ckϵkwkvk,LV5=nk=1dk[2(dQk+αk+μk)σ24k]y2k+2nk=1dkδkwkyk,LV6=nk=1ek[2dRkσ25k]z2k+2nk=1ekγkwkzk+2nk=1ekμkykzk. (2.7)

    and

    2ckϵkwkvk14[(dEk+ϵk)12σ22k]v2k+4c2kϵ2k[(dEk+ϵk)12σ22k]w2k,2dkδkwkykdk[(dQk+αk+μk)12σ24k]y2k+dkδ2kdQk+αk+μk12σ24kw2k,2ekγkwkzkek(dRk12σ25k)z2k+ekγ2kdRk12σ25kw2k,2ekμkykzkdk(dQk+αk+μk12σ24k)y2k+e2kμ2kdk(dQk+αk+μk12σ24k)z2k. (2.8)

    Choosing ak=dSkηkϵk(dEk+ϵk)(dIk+αk+δk+γk)Λk(1kn) and b=max1kn{(dSk+dEk+ϵk)22ak(dEk+ϵk12σ22k)}, then from (2.4)–(2.8) we finally obtain

    LVnk=1{Aku2k+Bkv2k+Ckw2k+Dky2k+Fkz2k}+nk=1(bak+1)(σ1kΛkdSk)2, (2.9)

    where Ak,Bk,Ck,Dk and Fk are given in the above.

    If (2.1) holds, then Ak>0, Bk>0 and Dk>0. Further, we can choose ck, dk and ek such that

    0<ck<[(dEk+ϵk)12σ22k]4ϵ2k[2(dIk+αk+δk+γk)σ23k],0<dk<ckηk[2(dIk+αk+δk+γk)σ23k4ckϵ2k[(dEk+ϵk)12σ22k]],0<ek<dk(dQk+αk+μk12σ24k)(dRk12σ25k)μ2k.

    Particularly, we can take

    ck=[(dEk+ϵk)12σ22k]8ϵ2k[2(dIk+αk+δk+γk)σ23k],dk=ck2ηk[2(dIk+αk+δk+γk)σ23k4ckϵ2k[(dEk+ϵk)12σ22k]],ek=dk(dQk+αk+μk12σ24k)(dRk12σ25k)2μ2k,

    where ηk=δ2kdQk+αk+μk12σ24k+γ2k(dQk+αk+μk12σ24k)μ2k>0. Thus, we have

    Ck=ck[2(dIk+αk+δk+γk)σ23k4ckϵ2k[(dEk+ϵk)12σ22k]]dkδ2kdQk+αk+μk12σ24kekγ2kdRk12σ25k>ck[2(dIk+αk+δk+γk)σ23k4ckϵ2k[(dEk+ϵk)12σ22k]]dk[δ2kdQk+αk+μk12σ24k+γ2k(dQk+αk+μk12σ24k)μ2k]>0

    and Fk=ek(dRk12σ25kekμ2kdk(dQk+αk+μk12σ24k))>0. By integration and taking expectation of both sides of (2.3), from (2.9) we obtain

    E(V(t))E(V(0))=E[t0LV(r)dr]Et0nk=1{Aku2k(r)+Bkv2k(r)+Ckw2k(r)+Dky2k(r)+Fkz2k(r)}dr+nk=1(bak+1)(σ1kΛkdSk)2.

    Therefore,

    lim supt1tEt0nk=1{Aku2k(r)+Bkv2k(r)+Ckw2k(r)+Dky2k(r)+Fkz2k(r)}drnk=1(bak+1)(σ1kΛkdSk)2.

    Consequently,

    lim supt1tEt0nk=1{Ak(Sk(r)ΛkdSk)2+BkE2k(r)+CkI2k(r)+DkQ2k(r)+FkR2k(r)}drnk=1(bak+1)(σ1kΛkdSk)2.

    This completes the proof.

    Remark 1. From Theorem 1, we see that under some conditions the solution of model (1.2) will oscillates around the disease-free equilibrium of deterministic model (1.1), and the intensity of fluctuation is only relation to the intensity of the white noise B1k(t), but do not relation to the intensities of the other white noises. In a biological interpretation, as the intensity of stochastic perturbations is small, the solution of model (1.2) will be close to the disease-free equilibrium of model (1.1) most of the time.

    As a special case of model (1.2), when σ1k=0, then model (1.2) becomes into

    {dSk=[Λknj=1βkjSk(t)Ij(t)dSkSk]dt,dEk=[nj=1βkjSk(t)Ij(t)(dEk+ϵk)Ek]dt+σ2kEkdB2k,dIk=[ϵkEk(dIk+αk+δk+γk)Ik]dt+σ3kIkdB3k,dQk=[δkIk(dQk+αk+μk)Qk]dt+σ4kQkdB4k,dRk=[γkIk+μkQkdRkRk]dt+σ5kRkdB5k. (2.10)

    Obviously, E0 is also the disease-free equilibrium of model (2.10). From the proof of Theorem 2, we get

    LVnk=1{2akdSk(Sk(r)ΛkdSk)2+BkE2k(r)+CkI2k(r)+DkQ2k(r)+FkR2k(r)},

    which is negative definite if for each 1kn

    dIk+αk+δk+γk>12σ23k,dRk>12σ25k,dEk+ϵk>12σ22k,dQk+αk+μk>12σ24k. (2.11)

    Therefore, as a consequence of Theorem 1 we have the following result.

    Corollary 1. Assume that R01 and condition (2.11) holds. Then disease-free equilibrium E0 of model (2.9) is globally stochastically asymptotically stable.

    Firstly, we introduce some concepts and conclusions of graph theory (see [10]). A directed graph g=(V,E) contains a set V={1,2,,n} of vertices and a set E of arcs (k,j) leading from initial vertex k to terminal vertex j. A subgraph H of g is said to be spanning if H and g have the same vertex set. A directed digraph g is weighted if each arc (k,j) is assigned a positive weight akj. Given a weighted digraph g with n vertices, define the weight matrix A=(akj)n×n whose entry akj equals the weight of arc (k,j) if it exists, and 0 otherwise. A weighted digraph is denoted by (g,A). A digraph g is strongly connected if for any pair of distinct vertices, there exists a directed path from one to the other and it is well known that a weighted digraph (g,A) is stronly connected if and only if the weight matrix A is irreducible (see [32]).

    The Laplacian matrix of graph (g,A) is defined by

    LA=(k1a1ka12a1na21k2a2ka2n............an1an2knank).

    Let ck(1kn) denote the cofactor of the k-th diagonal element of LA. The following lemmas are the classical results of graph theory (see [21,33]) which will be used in this paper.

    Lemma 2. Assume that A is a irreducible matrix and n2. Then ck>0 for all 1kn.

    Lemma 3. Assume that A is a irreducible matrix and n2. Then the following equality holds

    nk=1nj=1ckakjGk(xk)=nk=1nj=1ckakjGj(xj),

    where Gk(xk)(1kn) are arbitrary functions.

    For model (1.2), we see that there is not any endemic equilibrium. Therefore, in order to study the persistence of disease in model (1.2), we need to study the asymptotic behavior of the endemic equilibrium of model (1.2) which is surrounding the deterministic model (1.1), we obtain the following result.

    Theorem 2. Assume that R0>1 and the following conditions hold

    σ21k<dSk,σ22k<12dEk,σ23k<12(dIk+αk+δk+γk),σ24k<12(dQk+αk+μk),σ25k<12dRk,1kn. (3.1)

    Then for any positive solution (Sk(t),Ek(t),Ik(t),Qk(t),Rk(t),1kn) of model (1.2) one has

    limt1tt0{nk=1{ckrdSkσ21kSk2akDk}(Sk(s)Sk)2+2nk=1(dEk2σ22k)(Ek(s)Ek)2+nk=1{ak(dIk+αk+δk+γk2σ23k)bkδ2kdQk+αk+μkdkγ2kdRk}(Ik(s)Ik)2+nk=1{bk(dQk+αk+μk2σ24k)μ2k}(Qk(s)Qk)2+nk=1dk{(dRk2σ25k)dk}(Rk(s)Rk)2}dsnk=1ρk,

    where E=(Sk,Ek,Ik,Qk,Rk,1kn) be the endemic equilibrium of model (1.1), and

    ρk=2nk=1ak{σ21k(Sk)2+σ22k(Ek)2+(1+dEk+dIk+αk+δk+γkϵk)σ23k(Ik)2}+2nk=1bkσ24k(Qk)2+2nk=1dkσ25k(Rk)2+12nk=1ck[(K+2)σ21kSk+(K+1)σ22kEk+(K+1)dEk+ϵkϵkσ23kIk],

    and positive constants r, ak,bk,ck and Dk(1kn) will be confirmed in the proof of the theorem.

    Proof. When R0>1, from [23] there exits an endemic equilibrium E of model (1.1), then

    Λk=nj=1βkjSkIj+dSkSk,nj=1βkjSkIj=(dEk+ϵk)Ek,ϵkEk=(dIk+αk+δk+γk)Ik,δkIk=(dQk+αk+μk)Qk,γkIk+μkQk=dRkRk,1kn.

    Let matrix A=(akj)n×n with akj=βkjSkIj,k,j=1,2,,n. Since B=(βkj)n×n is irreducible, then A also is irreducible.

    Firstly, define the C2-function V1:R3n+R+ by

    V1(Sk,Ek,Ik,1kn)=nk=1ck[(SkSkSklogSkSk)+(EkEkEklogEkEk)+dEk+ϵkϵk(IkIkIklogIkIk)],

    where ck(1kn) are the cofactor of the k-th diagonal element of LA. V1 is positive definite. From Itˆo's formula, by calculating we can get

    LV1=nk=1ck[3nj=1βkjSkIj+2dSkSkdSkSk(Sk)2dSkSknj=1βkj(Sk)2IjSk+nj=1βkjSkIjnj=1βkjSkIjSkIjEkSkEkIjnj=1βkjSkIjIkEkEkIknj=1βkjSkIjIkIk]+12nk=1ck(σ21kSk+σ22kEk+dEk+ϵkϵkσ23kIk)=nk=1ckdSkSk(2SkSkSkSk)+nk=1ck[3nj=1βkjSkIjnj=1βkjSkIjSkSknj=1βkjSkIjSkIjEkSkEkIjnj=1βkjSkIjIkEkEkIk]+nk=1ck[nj=1βkjSkIjnj=1βkjSkIjIkIk]+12nk=1ck(σ21kSk+σ22kEk+dEk+ϵkϵkσ23kIk). (3.2)

    By Lemma 2, we obtain

    nk=1ck[nj=1βkjSkIjnj=1βkjSkIjIkIk]=nk=1nj=1ckβkjSkIjIjIjnk=1nj=1ckβkjSkIjIkIk=nk=1nj=1ckβkjSkIjIkIknk=1nj=1ckβkjSkIjIkIk=0. (3.3)

    Similarly, we also get

    nk=1nj=1ckβkjSkIjIkEkEkIk=nk=1nj=1ckβkjSkIjIjEjEjIj.

    Hence

    nk=1ck[3nj=1βkjSkIjnj=1βkjSkIjSkSknj=1βkjSkIjSkIjEkSkEkIjnj=1βkjSkIjIkEkEkIk]=nk=1cknj=1βkjSkIj[3SkSkSkIjEkSkEkIjIjEjEjIj]nk=1cknj=1βkjSkIj[33lnEkEklnEjEj]=nk=1cknj=1βkjSkIjlnEkEknk=1cknj=1βkjSkIjlnEjEj=0, (3.4)

    where the last equality is derived from Lemma 3. Substituting (3.3) and (3.4) into (3.2), we have

    LV1nk=1ckdSkSk(2SkSkSkSk)+12nk=1ck(σ21kSk+σ22kEk+dEk+ϵkϵkσ23kIk). (3.5)

    Secondly, define the C2-function V2:R2n+R+ as follows.

    V2(Ek,Ik,1kn)=nk=1ck[(EkEkEklogEkEk)+dEk+ϵkϵk(IkIkIklogIkIk)],

    where ck(1kn) are given as in V1. V2 is positive definite. It follows from Itˆo's formula that

    LV2=nk=1ck[nj=1βkjSkIj(dEk+ϵk)(dIk+αk+δk+γk)ϵkIknj=1βkjSkIjEkEk+(dEk+ϵk)Ek(dEk+ϵk)EkIkIk+(dEk+αk+δk+ϵk)(dIk+αk+δk+γk)ϵkIk]+12nk=1ck(σ22kEk+dEk+ϵkϵkσ23kIk)]=nk=1nj=1ckβkj(SkSk)(IjIj)+nk=1nj=1ckβkjSkIj[1+SkSkSkEkIjSkEkIjEkIkEkIk]+nk=1nj=1ckβkjSkIjIjIjnk=1nj=1ckβkjSkIjIkIk+12nk=1ck(σ22kEk+dEk+ϵkϵkσ23kIk)]. (3.6)

    We have

    nk=1nj=1ckβkjSkIj[1+SkSkSkEkIjSkEkIjEkIkEkIk]nk=1nj=1ckβkjSkIj[SkSk1logSkEkIjSkEkIjlogEkIkEkIk]=nk=1nj=1ckβkjSkIj[SkSk1logSkSklogIjIjlogIkIk]nk=1nj=1ckβkjSkIj[SkSk+SkSk2]nk=1nj=1ckβkjSkIj[logIjIj+logIkIk]=nk=1nj=1ckβkjSkIj[SkSk+SkSk2], (3.7)

    where the last equality is derived from Lemma 3 such that

    nk=1nj=1ckβkjSkIklogIjIjnk=1nj=1ckβkjSkIklogIkIk=0.

    We further get

    nk=1nj=1ckβkjSkIjIjIjnk=1nj=1ckβkjSkIjIkIk=0. (3.8)

    Substituting (3.7) and (3.8) into (3.6), we have

    LV2nk=1nj=1ckβkj(SkSk)(IjIj)+nk=1nj=1ckβkjSkIk[SkSk+SkSk2]+12nk=1ck(σ22kEk+dEk+ϵkϵkσ23kIk)]. (3.9)

    Thirdly, define the C2-function V3:Rn+R+ by

    V3(Sk,1kn)=nk=1ck(SkSk)22Sk,

    where ck(1kn) are given as in V1. We obtain

    LV3=nk=1ckdSk(SkSk)2Sknk=1nj=1ckβkj(SkSk)2IjSk+12nk=1ckS2kσ21knk=1nj=1ckβkj(SkSk)(IjIj)nk=1ck(dSkσ21k)(SkSk)2Sknk=1nj=1ckβkj(SkSk)(IjIj)+nk=1ckSkσ21k. (3.10)

    Choose K=nj=1βkjIkdSk, then (3.5) together with (3.9) and (3.10) implies

    L(KV1+V2+V3)nk=1KckdSkSk(2SkSkSkSk)+12nk=1Kck(σ21kSk+σ22kEk+dEk+ϵkϵkσ23kIk)+nk=1nj=1ckβkj(SkSk)(IjIj)+nk=1nj=1ckβkjSkIk[SkSk+SkSk2]+12nk=1ck(σ22kEk+dEk+ϵkϵkσ23kIk)]nk=1ck(dSkσ21k)(SkSk)2Sknk=1nj=1ckβkj(SkSk)(IjIj)+nk=1ckSkσ21knk=1ck(dSkσ21k)(SkSk)2Sk+Ak, (3.11)

    where Ak=12nk=1ck[(K+2)σ21kSk+(K+1)σ22kEk+(K+1)dEk+ϵkϵkσ23kIk].

    Next, define the C2-function V4:R3n+R+ by

    V4(Sk,Ek,Ik,1kn)=nk=1ak(SkSk+EkEk+IkIk)2,

    where ak(1kn) are positive constants to be determined later. By calculating, we can get

    LV4=2nk=1ak[dSk(SkSk)2+dEk(EkEk)2+(dIk+αk+δk+γk)(IkIk)2]2nk=1{ak(dSk+dEk)(SkSk)(EkEk)+(dSk+dIk+αk+δk+γk)×(SkSk)(IkIk)+(dEk+dIk+αk+δk+γk)(EkEk)(IkIk)}+nk=1ak(σ21kS2k+σ22kEk+σ23kIk).

    Since 2(dSk+dEk)(SkSk)(EkEk)(dSk+dEk)2dEk(SkSk)2+dEk(EkEk)2 and

    2(dSk+dIk+αk+δk+γk)(SkSk)(IkIk)(dSk+dIk+αk+δk+γk)2(dIk+αk+δk+γk)(SkSk)2+(dIk+αk+δk+γk)(IkIk)2,

    we further obtain

    LV42nk=1ak[Dk(SkSk)2(dEk2σ22k)(EkEk)2(dIk+αk+δk+γk2σ23k)(IkIk)2]2nk=1ak(dEk+dIk+αk+δk+γk)(EkEk)(IkIk)+2nk=1ak(σ21k(Sk)2+σ22k(Ek)2+σ23k(Ik)2), (3.12)

    where Dk=dSk+dEk+(dSk)2dEk+(dSk+dIk+αk+δk+γk)2dIk+αk+δk+γk+σ21k.

    Further, define the C2-function V5:Rn+R+ by

    V5(Ik,1kn)=nk=1ak(dEk+dIk+αk+δk+γk)ϵk(IkIk)2.

    We obtain

    LV5=2nk=1ak[(dEk+dIk+αk+δk+γk)ϵk(dIk+αk+δk+γkσ23k)(IkIk)2(dEk+dIk+αk+δk+γk)(EkEk)(IkIk)]+2nk=1ak(dEk+dIk+αk+δk+γk)ϵkσ23k(Ik)2. (3.13)

    Finally, define the C2 functions V6 and V7:Rn+R+ as follows.

    V6(Qk,1kn)=nk=1bk(QkQk)2,V7(Rk,1kn)=nk=1dk(RkRk)2,

    where bk, dk(1kn) are positive constants to be determined later. We get

    LV6=2nk=1bk(dQk+αk+μkσ24k)(QkQk)2+2nk=1bkδk(QkQk)(IkIk)+2nk=1bkσ24k(Qk)2nk=1bk(dQk+αk+μk2σ24k)(QkQk)2+nk=1bkδ2kdQk+αk+μk(IkIk)2+2nk=1bkσ24k(Qk)2 (3.14)

    and

    LV7=2nk=1dk(dRkσ25k)(RkRk)2+2nk=1dkγk(RkRk)(IkIk)+2nk=1dkμk(QkQk)(RkRk)+2nk=1dkσ25k(Rk)2nk=1dk(dRk2σ25k)(RkRk)2+nk=1dkγ2kdRk(IkIk)2+nk=1μ2k(QkQk)2+nk=1d2k(RkRk)2+2nk=1dkσ25k(Rk)2, (3.15)

    where the last equality is derived by the inequality 2aba2+b2.

    From (3.12)–(3.15) we obtain

    L(V4+V5+V6+V7)2nk=1akDk(SkSk)22nk=1(dEk2σ22k)(EkEk)2nk=1{ak(dIk+αk+δk+γk2σ23k)bkδ2kdQk+αk+μkdkγ2kdRk}(IkIk)2nk=1{bk(dQk+αk+μk2σ24k)μ2k}(QkQk)2nk=1dk{(dRk2σ25k)dk}(RkRk)2+nk=1Ck, (3.16)

    where

    Ck=2nk=1ak{σ21k(Sk)2+σ22k(Ek)2+(1+dEk+dIk+αk+δk+γkϵk)σ23k(Ik)2)}+2nk=1bkσ24k(Qk)2+2nk=1dkσ25k(Rk)2.

    From condition (3.1), we can choose positive numbers r, ak, bk and dk for k=1,2,,n satisfying dk<dRk2σ25k and

    r>2SkDkak(dSkσ21k)ck,ak>[bkδ2kdQk+αk+μk+dkγ2kdRk](dIk+αk+δk+γk2σ23k),bk>μ2k(dQk+αk+μk2σ24k)

    such that for each 1kn

    ak(dIk+αk+δk+γk2σ23k)[bkδ2kdQk+αk+μk+dkγ2kdRk]>0,dRk2σ25kdk>0,bk(dQk+αk+μk2σ24k)μ2k>0,ckrSkdSkσ21k2akDk>0.

    Lastly, define a Lyapunov function as follows

    V=r(KV1+V2+V3)+V4+V5+V6+V7.

    By Itˆo's formula, we obtain

    dV=LVdt+nk=1σ1k[ckr(K+SkSk)(SkSk)+2ak(SkSk+EkEk+IkIk)Sk]dB1k+2nk=1σ2k[ckrK(EkEk)+ak(SkSk+EkEk+IkIk)Ek]dB2k+nk=1σ3k{[r(K+1)ckdEk+ϵkϵk+akdEk+dIk+αk+δk+γkϵkIk](IkIk)+2ak×(SkSk+IkIk+EkEk)Ik}dB3k+nk=1σ4kbk(QkQk)QkdB4k+nk=1σ5kdk(RkRk)RkdB5k, (3.17)

    where (3.11) together with (3.16) implies

    LVnk=1[{ckr(dSkσ21k)Sk2akDk}(SkSk)2+2(dEk2σ22k)(EkEk)2+{ak(dIk+αk+δk+γk2σ23k)bkδ2kdQk+αk+μkdkγ2kdRk}(IkIk)2+{bk(dQk+αk+μk2σ24k)μ2k}(QkQk)2+{(dRk2σ25k)dk}(RkRk)2]+nk=1ρk. (3.18)

    By integration and taking expectation of both sides of (3.17), we obtain

    E(V(t))E(V(0))=E[t0LV(r)dr]Et0nk=1[{ckr(dSkσ21k)Sk2akDk}(SkSk)2+2(dEk2σ22k)(EkEk)2+{ak(dIk+αk+δk+γk2σ23k)bkδ2kdQk+αk+μkdkγ2kdRk}(IkIk)2+{bk(dQk+αk+μk2σ24k)μ2k}(QkQk)2+{(dRk2σ25k)dk}(RkRk)2]dr+tnk=1ρk.

    Therefore,

    lim supt1tEt0nk=1[{ckr(dSkσ21k)Sk2akDk}(SkSk)2+2(dEk2σ22k)(EkEk)2+{ak(dIk+αk+δk+γk2σ23k)bkδ2kdQk+αk+μkdkγ2kdRk}(IkIk)2+{bk(dQk+αk+μk2σ24k)μ2k}(QkQk)2+{(dRk2σ25k)dk}(RkRk)2]drnk=1ρk.

    This completes the proof.

    As a consequence of Theorem 2, we have the following result on the existence and uniqueness of stationary distribution for model (1.2).

    Theorem 3. Assume that all conditions in Theorem 2 hold. Then model (1.2) has a unique stationary distribution μ() in R5n+.

    Proof. Choose region Ω in ([34], Lemma 2.5) by Ω=R5n+. Consider the following inequality

    nk=1{ck(dSkσ21k)Sk2akBk}(SkSk)2+2nk=1(dEk2σ22k)(EkEk)2+nk=1{ak(dIk+αk+δk+γk2σ23k)bkδ2kdQk+αk+μkdkγ2kdRk}(IkIk)2+nk=1{bk(dQk+αk+μk2σ24k)μ2k}(QkQk)2+nk=1dk{(dRk2σ25k)dk}(RkRk)2H.

    Let region U1 denote all points (Sk,Ek,Ik,Qk,Rk,1kn) which satisfy the above inequality with H=2nk=1ρk and region U2 denote all points (Sk,Ek,Ik,Qk,Rk,1kn) which satisfy the above inequality with H=3nk=1ρk. Obviously, U2 is a neighborhood of U1 and the closure ˉU2Ω. Then from (3.18), for any xΩU1,

    LVnk=1[{ckr(dSkσ21k)Sk2akDk}(SkSk)2+2(dEk2σ22k)(EkEk)2+{ak(dIk+αk+δk+γk2σ23k)bkδ2kdQk+αk+μkdkγ2kdRk}(IkIk)2+{bk(dQk+αk+μk2σ24k)μ2k}(QkQk)2+{(dRk2σ25k)dk}(RkRk)2]+nk=1ρknk=1ρk,

    which implies condition (ⅱ) in ([35], Lemma 2.5) is satisfied.

    For model (1.2), the diffusion matrix is

    A(x)=diag(σ21kS2k,σ22kE2k,σ23kI2k,σ24kQ2k,σ25kR2k,1kn).

    Choose a positive constant MinfˉU2{σ21iS2i,σ22iE2i,σ23iI2i,σ24iQ2i,σ25iR2i,1in}. Then,

    5ni,j=1aijξiξj=ni=1σ21iS2iξ25i4+ni=1σ22iE2iξ25i3+ni=1σ23iI2iξ25i2+ni=1σ24iQ2iξ25i1+ni=1σ25iR2iξ25iMξ2,

    for all (Si,Ei,Ii,Qi,Ri,1in)ˉU2 and ξR5n. This implies condition (ⅰ) in ([34], Lemma 2.5) is also satisfied. Therefore, by ([34], Lemma 2.5), model (1.2) has a unique stationary distribution μ in R5n+. This completes the proof.

    In this section, we analyse the stochastic behaviour of model (1.2) by means of the numerical simulations in order to make readers understand our results more better. The numerical simulation method can be found in [36]. The corresponding discretization system of

    {Sk,i+1=Sk,i+[Λkβk1Sk,iI1,iβk2Sk,iI2,idSkSk,i]Δt+σ1kSk,iΔtε1k,i+σ21kSk,i2(ε21k,iΔtΔt),Ek,i+1=Ek,i+[βk1Sk,iI1,i+βk2Sk,iI2,i(dEk+ϵk)Ek,i]Δt+σ2kEk,iΔtε2k,i+σ22kEk,i2(ε22k,iΔtΔt),Ik,i+1=Ik,i+[ϵkEk,i(dIk+αk+δk+γk)Ik,i]Δt+σ3kIk,iΔtε3k,i+σ23kIk,i2(ε23k,iΔtΔt),Qk,i+1=Qk,i+[δkIk,i(dQk+αk+μk)Qk,i]Δt+σ4kQk,iΔtε4k,i+σ24kQk,i2(ε24k,iΔtΔt),Rk,i+1=Rk,i+[γkIk,i+μkQk,idRkRk,i]Δt+σ5kRk,iΔtε5k,i+σ25kRk,i2(ε25k,iΔtΔt),

    where time increment Δt>0, and ε1k,i, ε2k,i, ε3k,i, ε4k,i, ε5k,i for 1kn are N(0,1)-distributed independent random variables which be generated numerically by pseudo-random number generators.

    Example 1. In model (1.2), we choose n=2 and the parameters Λ1=3.2, ϵ1=0.1, α1=0.1, β11=0.409, dS1=0.9, dE1=0.7, dI1=0.81, dQ1=0.2, dR1=0.65, μ1=0.3, γ1=0.04, β12=0.02, δ1=0.1, σ11=0.15, σ21=0.1, σ31=0.41, σ41=0.2, σ51=0.3, Λ2=7.5, ϵ2=2.4, α2=0.2, β21=0.05, dS2=0.49, dE2=0.25, dI2=0.15, dQ2=0.25, dR2=0.39, μ2=0.5, γ2=0.15, β22=0.0014, δ2=0.43, σ12=0.2, σ22=0.6, σ32=0.5, σ42=0.8 and σ52=0.8.

    By computing, we have R00.8675<1 and disease-free equilibrium E0=(3.56,0,0,0,0,15.31,0,0,0,0) for corresponding deterministic model (1.1), and the conditions in Theorem 1 are satisfied. Therefore, according to the conclusion in Theorem 1 by numerical calculation we can obtain that for the solution (Sk(t),Ek(t),Ik(t),Qk(t),Rk(t),k=1,2) satisfying the initial values (S1(0),E1(0),I1(0),Q1(0),R1(0))=(0.75,0.8,0.8,0.2,0.2) and (S2(0),E2(0),I2(0),Q2(0),R2(0))=(1.7,4.5,2.7,4.3,5) one has

    lim supt1tEt02k=1{Ak(Sk(r)S0k)2+BkE2k(r)+CkI2k(r)+DkQ2k(r)+FkR2k(r)}dt10.49, (4.1)

    where S01=3.56, S02=15.31, A1=0.8775, A2=0.48, B1=0.1988, B2=0.6175, C1=18.21, C2=0.04295, D1=0.1482, D2=0.0077, F1=0.1507 and F2=3.7551×105.

    From the numerical simulations given in Figure 1 we easily see that the above formula (4.1) holds. That is, the solution of stochastic model (1.2) asymptotically oscillates in probability around disease-free equilibrium E0.

    Figure 1.  The numerical simulations of asymptotic oscillation in probability around disease-free equilibrium E0 for the solution (Sk(t),Ek(t),Ik(t),Qk(t),Rk(t),k=1,2) of stochastic model with initial values (S1(0),E1(0),I1(0),Q1(0),R1(0))=(0.75,0.8,0.8,0.2,0.2) and (S2(0),E2(0),I2(0),Q2(0),R2(0))=(1.7,4.5,2.7,4.3,5).

    In addition, from Figure 1 we also easily see that the mean of susceptible Sk(t)(k=1,2) tend to S0k and all exposed Ek, infectious Ik, quarantined Qk and recovered Rk for k=1,2 tend to zero in probability as t.

    Example 2. In model (1.2), we choose n=2 and the parameters Λ1=0.8, ϵ1=0.1, α1=0.1, β11=0.109, dS1=0.19, dE1=1.107, dI1=0.081, dQ1=0.2, dR1=0.65, μ1=0.3, γ1=0.04, β12=0.02, δ1=0.01, σ11=1.15, σ21=1.1, σ31=1.41, σ41=01.2, σ51=1.3, Λ2=1.5, ϵ2=2.4, α2=0.2, β21=0.05, dS2=0.49, dE2=0.25, dI2=0.15, dQ2=0.25, dR2=0.39, μ2=0.5, γ2=0.15, β22=0.0014, δ2=0.043, σ12=1.2, σ22=1.6, σ32=0.5, σ42=0.8 and σ52=0.8.

    By computing, we have R00.51741. Since dS1σ211=1.13<0, dS2σ212=0.33<0, dR112σ251=0.2<0 and dR212σ252=0.46<0, the condition (2.1) in Theorem 1 does not hold. However, from the numerical simulations given in Figure 2, we can see that the solution (Sk(t),Ek(t),Ik(t),Qk(t),Rk(t),k=1,2) of stochastic model (1.2) with initial values (S1(0),E1(0),I1(0),Q1(0),R1(0))=(0.75,0.8,0.8,0.2,0.2) and (S2(0),E2(0),I2(0),Q2(0),R2(0))=(1.7,4.5,2.7,4.3,5) asymptotically oscillates in probability around the disease-free equilibrium E0=(4.21,0,0,0,0,3.06,0,0,0,0) of corresponding deterministic model (1.1). This example seems to indicate that the condition (2.1) in Theorem 1 can be weakened or taken out.

    Figure 2.  The numerical simulations of asymptotic oscillation in probability around disease-free equilibrium E0 for the solution (Sk(t),Ek(t),Ik(t),Qk(t),Rk(t),k=1,2) of stochastic model (1.2) with initial values (S1(0),E1(0),I1(0),Q1(0),R1(0))=(0.75,0.8,0.8,0.2,0.2) and (S2(0),E2(0),I2(0),Q2(0),R2(0))=(1.7,4.5,2.7,4.3,5).

    Example 3. In model (1.2), we choose n=2 and the parameters Λ1=4.5, ϵ1=1, α1=0.1, β11=1.55, dS1=0.5, dE1=0.15, dI1=0.1, dQ1=0.2, dR1=0.65, μ1=0.3, γ1=0.4, β12=1.35, δ1=0.6, σ11=0.3, σ21=0.5, σ31=0.4, σ41=0.2, σ51=0.4, Λ2=7.5, ϵ2=2.4, α2=0.2, β21=1.5, dS2=0.49, dE2=0.25, dI2=0.15, dQ2=0.25, dR2=0.39, μ2=0.5, γ2=0.15, β22=1.24, δ2=0.43, σ12=0.2, σ22=0.6, σ32=0.5, σ42=0.8 and σ52=0.3.

    By computing, we have R01.1032>1 and the conditions in Theorem 2 are satisfied. The numerical simulations are given in Figures 3 and 4. Figure 3 shows that the solution (Sk(t),Ek(t),Ik(t),Qk(t),Rk(t),k=1,2) of stochastic model (1.2) satisfying the initial values (S1(0),E1(0),I1(0),Q1(0),R1(0))=(0.75,0.8,0.8,0.2,0.2) and (S2(0),E2(0),I2(0),Q2(0),R2(0))=(1.7,4.5,2.7,4.3,5) asymptotically oscillates in probability around the endemic equilibrium E=(0.37,3.35,0.27,0.38,0.19,0.79,2.68,6.93,3.14,6.68) of corresponding deterministic model (1.1). Figure 4 shows that the solution has a unique stationary distribution. Therefore, the conclusions of Theorem 3 are validated by the numerical example.

    Figure 3.  The numerical simulations of asymptotic oscillation in probability around endemic equilibrium E for the solution (Sk(t),Ek(t),Ik(t),Qk(t),Rk(t),k=1,2) of stochastic model (1.2) with initial values (S1(0),E1(0),I1(0),Q1(0),R1(0))=(0.75,0.8,0.8,0.2,0.2) and (S2(0),E2(0),I2(0),Q2(0),R2(0))=(1.7,4.5,2.7,4.3,5).
    Figure 4.  The stationary distribution of the solution (Sk(t),Ek(t),Ik(t),Qk(t),Rk(t),k=1,2) for the stochastic model (1.2).

    In addition, from Figure 3 we also easily see that the mean value of the solution for stochastic model (1.2) asymptotically oscillates in probability around the endemic equilibrium E of corresponding deterministic model (1.1). From Figure 5 we can find the relationship between variances of the solution (Sk(t),Ek(t),Ik(t),Qk(t),Rk(t),k=1,2) and the intensities of noises (σ2ik,σ22k,σ23k,σ24k,σ25k,k=1,2) as time t is enough large.

    Figure 5.  The numerical simulations of variances for the solution (Sk(t),Ek(t),Ik(t),Qk(t),Rk(t),k=1,2) of stochastic model (1.2) with initial values (S1(0),E1(0),I1(0),Q1(0),R1(0))=(0.75,0.8,0.8,0.2,0.2) and (S2(0),E2(0),I2(0),Q2(0),R2(0))=(1.7,4.5,2.7,4.3,5).

    Example 4. In model (1.2), we choose n=2 and the parameters Λ1=4.5, ϵ1=0.1, α1=0.1, β11=1.55, dS1=2.05, dE1=1.015, dI1=0.51, dQ1=0.02, dR1=0.65, μ1=0.3, γ1=0.04, β12=1.35, δ1=0.6, σ11=2.3, σ21=1.5, σ31=0.5, σ41=0.4, σ51=0.4, Λ2=7.5, ϵ2=2.4, α2=0.2, β21=1.5, dS2=0.49, dE2=0.25, dI2=0.15, dQ2=0.25, dR2=0.39, μ2=0.5, γ2=0.15, β22=0.24, δ2=0.43, σ12=1.2, σ22=0.6, σ32=0.5, σ42=0.8 and σ52=0.3.

    By computing, we have R01.09013>1. Since dS1σ211=10.12<0, dS2σ221=0.43<0 and dE112σ221=0.11<0, the condition (3.1) in Theorem 2 does not hold. However, from the numerical simulations are given in Figures 6 we can see that the solution (Sk(t),Ek(t),Ik(t),Qk(t),Rk(t),k=1,2) of stochastic model (1.2) with initial values (S1(0),E1(0),I1(0),Q1(0),R1(0))=(0.75,0.8,0.8,0.2,0.2) and (S2(0),E2(0),I2(0),Q2(0),R2(0))=(1.7,4.5,2.7,4.3,5) asymptotically oscillates in probability around the endemic equilibrium E=(0.43,3.24,0.26,0.37,0.19,3.34,2.22,5.7,2.59,5.51) of corresponding deterministic model (1.1). This example seems to indicate that the condition (3.1) in Theorem 2 can be weakened or taken out.

    Figure 6.  The numerical simulations of asymptotic oscillation in probability around endemic equilibrium E for the solution (Sk(t),Ek(t),Ik(t),Qk(t),Rk(t),k=1,2) of stochastic model (1.2) with initial values (S1(0),E1(0),I1(0),Q1(0),R1(0))=(0.75,0.8,0.8,0.2,0.2) and (S2(0),E2(0),I2(0),Q2(0),R2(0))=(1.7,4.5,2.7,4.3,5).

    In this research we consider a class of stochastic multi-group SEIQR (susceptible, exposed, infectious, quarantined and recovered) models in computer network. For the deterministic system, if the reproduction number R0>1, the system has unique endemic equilibrium which is globally stable, this means that the disease will persist at the endemic equilibrium level if it is initially present. It is clear that when the disease is endemic, the recovery nodes increases with the increasing quarantine nodes, and finally both reach the steady state values. Thus, it will be of great importance for one to run anti-malicious software to quarantine infected nodes. In order to study the asymptotic behavior of model (1.2), we first introduce the global existence of a positive solution. Then by using the theory of graphs, stochastic Lyapunov functions method, Itˆo's formula and the theory of stochastic analysis, we carry out a detailed analysis on the asymptotic behavior of model (1.2). If R01, the solution of model (1.2) oscillates around the disease-free equilibrium, while if R0>1, the solution of model (1.2) fluctuates around the endemic equilibrium. The investigation of this stochastic model revealed that the stochastic stability of E depends on the magnitude of the intensity of noise as well as the parameters involved within the model system. finally, numerical methods are employed to illustrate the dynamic behavior of the model. The effect of quarantine on recovered nodes is also analyzed in the stochastic model.

    Some interesting topics deserve further consideration. On the one hand, we can solve the corresponding probability density function of various stochastic epidemic models. On the other hand, we need to establish a more complete and systematic theory to obtain more accurate conditions and density function. The reader is referred to [37,38,39,40,41,42,43,44,45]. These problems are expected to be studied and solved as planned future work.

    This research is supported by the Natural Science Foundation of Xinjiang of China (Grant Nos. 2020D01C178) and the National Natural Science Foundation of China (Grant Nos. 12101529, 12061079, 72163033, 72174175, 11961071).

    The authors declare there is no conflicts of interest.



    [1] R. M. May, A. L. Lloyd, Infection dynamics on scale-free networks, Phys. Rev. E, 64 (2001), 1–3. https://doi.org/10.1103/PhysRevE.64.066112 doi: 10.1103/PhysRevE.64.066112
    [2] B. K. Mishra, D. K. Saini, SEIRS epidemic model with delay for transmission of malicious objects in computer network, Appl. Math. Comput., 188 (2007), 1476–1482. https://doi.org/10.1016/j.amc.2006.11.012 doi: 10.1016/j.amc.2006.11.012
    [3] B. K. Mishra, D. Saini, Mathematical models on computer virus, Appl. Math. Comput., 187 (2007), 929–936. https://doi.org/10.1016/j.amc.2006.09.062 doi: 10.1016/j.amc.2006.09.062
    [4] B. K. Mishra, N. Jha, Fixed period of temporary immunity after run of anti-malicious software on computer nodes, Appl. Math. Comput., 190 (2007), 1207–1212. https://doi.org/10.1016/j.amc.2007.02.004 doi: 10.1016/j.amc.2007.02.004
    [5] H. Yuan, G. Chen, Network virus-epidemic model with the point-to-group information propagation, Appl. Math. Comput., 206 (2008), 357–367. https://doi.org/10.1016/j.amc.2008.09.025 doi: 10.1016/j.amc.2008.09.025
    [6] J. B. Shukla, G. Singh, P. Shukla, A. Tripathi, Modeling and analysis of the effects of antivirus software on an infected computer network, Appl. Math. Comput., 227 (2014), 11–18. https://doi.org/10.1016/j.amc.2013.10.091 doi: 10.1016/j.amc.2013.10.091
    [7] M. Sun, Q. Liu, An SIS epidemic model with time delay and stochastic perturbation on heterogeneous networks, Math. Biosci. Eng., 18 (2021), 6790–6805. https://doi.org/10.3934/mbe.2021337 doi: 10.3934/mbe.2021337
    [8] R. Zhao, Q. Liu, M. Sun, Dynamical behavior of a stochastic SIQS epidemic model on scale-free networks, J. Appl. Math. Comput., 68 (2022), 813–838. https://doi.org/10.1007/s12190-021-01550-9 doi: 10.1007/s12190-021-01550-9
    [9] Y. E. Ansari, A. E. Myr, L. Omari, Deterministic and stochastic study for an infected computer network model powered by a system of antivirus programs, Discrete Dyn. Nat. Soc., (2017), 3540278. https://doi.org/10.1155/2017/3540278 doi: 10.1155/2017/3540278
    [10] Y. Xu, X. Sun, H. Hu, Extinction and stationary distribution of a stochastic SIQR epidemic model with demographics and non-monotone incidence rate on scale-free networks, J. Appl. Math. Comput., 68 (2022), 3367–3395. https://doi.org/10.1007/s12190-021-01645-3 doi: 10.1007/s12190-021-01645-3
    [11] R. Ikram, A. Khan, M. Zahri, A. Saeed, M. Yavuz, P. Kumam, Extinction and stationary distribution of a stochastic COVID-19 epidemic model with time-delay, Comput. Biol. Med., 141 (2022), 105115. https://doi.org/10.1016/j.compbiomed.2021.105115 doi: 10.1016/j.compbiomed.2021.105115
    [12] M. Naim, Y. Sabbar, A. Zeb, Stability characterization of a fractional-order viral system with the non-cytolytic immune assumption, Math. Modell. Numer. Simul. Appl., 2 (2022), 164–176. https://doi.org/10.53391/mmnsa.2022.013 doi: 10.53391/mmnsa.2022.013
    [13] H. Joshi, M. Yavuz, S. Townley, B. K. Jha, Stability analysis of a non-singular fractional-order covid-19 model with nonlinear incidence and treatment rate, Phys. Scr., 98 (2023), 045216. https://doi.org/10.1088/1402-4896/acbe7a doi: 10.1088/1402-4896/acbe7a
    [14] M. Yavuz, F. Özköse, M. Susam, M. Kalidass, A new modeling of fractional-order and sensitivity analysis for Hepatitis-B disease with real data, Fractal Fract., 7 (2023), 165. https://doi.org/10.3390/fractalfract7020165 doi: 10.3390/fractalfract7020165
    [15] A. Din, M. Z. Abidin, Analysis of fractional-order vaccinated Hepatitis-B epidemic model with Mittag-Leffler kernels, Math. Modell. Numer. Simul. Appl., 2 (2022), 59–72. https://doi.org/10.53391/mmnsa.2022.006 doi: 10.53391/mmnsa.2022.006
    [16] H. Joshi, B. K. Jha, M. Yavuz, Modelling and analysis of fractional-order vaccination model for control of COVID-19 outbreak using real data, Math. Biosci. Eng., 20 (2023), 213–240. https://doi.org/10.3934/mbe.2023010 doi: 10.3934/mbe.2023010
    [17] Y. Sabbar, M. Yavuz, F. Özköse, Infection eradication criterion in a general epidemic model with logistic growth, quarantine strategy, media intrusion, and quadratic perturbation, Mathematics, 10 (2022), 4213. https://doi.org/10.3390/math10224213 doi: 10.3390/math10224213
    [18] B. K. Mishra, N. Jha, SEIQRS model for the transmission of malicious objects in computer network, Appl. Math. Modell., 34 (2010), 710–715. https://doi.org/10.1016/j.apm.2009.06.011 doi: 10.1016/j.apm.2009.06.011
    [19] B. K. Mishra, A. K. Singh, Two quarantine models on the attack of malicious objects in computer network, Math. Probl. Eng., (2012), 407064. https://doi.org/10.1155/2012/407064 doi: 10.1155/2012/407064
    [20] H. Guo, M. Y. Li, Z. Shuai, A graph-theoretic approach to the method of global Lyapunov functions, Proc. Am. Math. Soc., 136 (2008), 2793–2802. https://doi.org/10.1090/S0002-9939-08-09341-6 doi: 10.1090/S0002-9939-08-09341-6
    [21] M. Y. Li, Z. Shuai, Global-stability problem for coupled systems of differential equation on networks, J. Differ. Equations, 248 (2010), 1–20. https://doi.org/10.1016/j.jde.2009.09.003 doi: 10.1016/j.jde.2009.09.003
    [22] H. Guo, M. Y. Li, Z. Shuai, Global stability of the endemic equilibrium of multigroup SIR epidemic models, Can. Appl. Math. Q., 14 (2006), 259–284.
    [23] Z. Wang, X. Fan, Q. Han, Global stability of deterministic and stochastic multigroup SEIQR models in computer network, Appl. Math. Modell., 37 (2013), 8673–8686. https://doi.org/10.1016/j.apm.2013.07.037 doi: 10.1016/j.apm.2013.07.037
    [24] A. L. Krause, L. Kurowski, K. Yawar, R. A. Van Gorder, Stochastic epidemic metapopulation models on networks: SIS dynamics and control strategies, J. Theor. Biol., 449 (2018), 35–52. https://doi.org/10.1016/j.jtbi.2018.04.023 doi: 10.1016/j.jtbi.2018.04.023
    [25] Y. El Ansari, A. El Myr, L. Omari, Deterministic and stochastic study for an infected computer network model powered by a system of antivirus programs, Discrete Dyn. Nat. Soc., 2017, 3540278. https://doi.org/10.1155/2017/3540278 doi: 10.1155/2017/3540278
    [26] A. Caruso, M. E. Gargano, D. Valenti, A. Fiasconaro, B. Spagnolo, Cyclic fluctuations, climatic changes and role of noise in planktonic foraminifera in the mediterranean sea, Fluctuation Noise Lett., 5 (2005), 349–355. https://doi.org/10.1142/S0219477505002768 doi: 10.1142/S0219477505002768
    [27] C. Zhang, Y. Zhao, Y. Wu, S. Deng, A stochastic dynamic model of computer viruses, Discrete Dyn. Nat. Soc., (2012), 264874. https://doi.org/10.1155/2012/264874 doi: 10.1155/2012/264874
    [28] C. Yuan, D. Jiang, D. O'Regan, R. P. Agarwal, Stochastically asymptotically stability of the multi-group SEIR and SIR models with random perturbation, Commun. Nonlinear Sci. Numer. Simul., 17 (2012), 2501–2516. https://doi.org/10.1016/j.cnsns.2011.07.025 doi: 10.1016/j.cnsns.2011.07.025
    [29] J. Yu, D. Jiang, N. Shi, Global stability of two-group SIR model with random perturbation, J. Math. Anal. Appl., 360 (2009), 235–244. https://doi.org/10.1016/j.jmaa.2009.06.050 doi: 10.1016/j.jmaa.2009.06.050
    [30] R. M. May, Stability and Complexity in Model Ecosystems, Princeton University Press, Princeton, 1973. https://doi.org/10.1515/9780691206912
    [31] C. Ji, D. Jiang, N. Shi, Multigroup SIR epidemic model with stochastic perturbation, Phys. A, 390 (2011), 1747–1762. https://doi.org/10.1016/j.physa.2010.12.042 doi: 10.1016/j.physa.2010.12.042
    [32] A. Berman, R. J. Plemmons, Nonnegative Matrices in the Mathematical Sciences, Academic Press, New York, 1979.
    [33] J. D. Brown, Microprobe analysis, X-Ray Spectrom., 3 (1974), 28. https://doi.org/10.1002/xrs.1300030415 doi: 10.1002/xrs.1300030415
    [34] R. Z. Khas'Miniskii, Stochastic Stability of Differential Equation, Springer Dordrecht, The Netherlands, 1980.
    [35] Q. Liu, D. Jiang, N. Shi, T. Hayat, A. Alsaedi, Asymptotic behavior of multi-group epidemic models with distributed delays, Phys. A, 467 (2017), 527–541. https://doi.org/10.1016/j.physa.2016.10.034 doi: 10.1016/j.physa.2016.10.034
    [36] D. Higham, An algorithmic introduction to numerical simulation of stochastic differential equations, SIAM Rev., 43 (2001), 525–546. https://doi.org/10.1137/S0036144500378302 doi: 10.1137/S0036144500378302
    [37] B. Zhou, X. Zhang, D. Jiang, Dynamics and density function analysis of a stochastic SVI epidemic model with half saturated incidence rate, Chaos Solitons Fractals, 137 (2020), 109865. https://doi.org/10.1016/j.chaos.2020.109865 doi: 10.1016/j.chaos.2020.109865
    [38] B. Zhou, D. Jiang, Y. Dai, T. Hayat, A. Alsaedi, Stationary distribution and probability density function of a stochastic SVIS epidemic model with standard incidence and vaccination strategies, Chaos Solitons Fractals, 143 (2021), 110601. https://doi.org/10.1016/j.chaos.2020.110601 doi: 10.1016/j.chaos.2020.110601
    [39] J. Ge, W. Zuo, D. Jiang, Stationary distribution and density function analysis of a stochastic epidemic HBV model, Math. Comput. Simul., 191 (2022), 232–255. https://doi.org/10.1016/j.matcom.2021.08.003 doi: 10.1016/j.matcom.2021.08.003
    [40] B. Zhou, B. Han, D. Jiang, T. Hayat, A. Alsaedi, Stationary distribution, extinction and probability density function of a stochastic vegetation-water model in arid ecosystems, J. Nonlinear Sci., 32 (2022), 30. https://doi.org/10.1007/s00332-022-09789-7 doi: 10.1007/s00332-022-09789-7
    [41] D. H. Nguyen, G. Yin, C. Zhu, Long-term analysis of a stochastic SIRS model with general incidence rates, SIAM J. Appl. Math., 80 (2020), 814–838. https://doi.org/10.1137/19M1246973 doi: 10.1137/19M1246973
    [42] S. Pan, Q. Zhang, M. Anke, Near-optimal control of a stochastic vegetation-water system with reaction diffusion, Math. Meth. Appl. Sci., 43 (2020), 6043–6061. https://doi.org/10.1002/mma.6346 doi: 10.1002/mma.6346
    [43] S. Pan, Q. Zhang, A. Meyer-Baese, Stationary distribution of a stochastic vegetation-water system with reaction-diffusion, Appl. Math. Lett., 123 (2022), 107589. https://doi.org/10.1016/j.aml.2021.107589 doi: 10.1016/j.aml.2021.107589
    [44] B. Han, D. Jiang, T. Hayat, A. Alsaedi, B. Ahmad, Stationary distribution and extinction of a stochastic staged progression AIDSdel with staged treatment and second-order perturbation, Chaos Solitons Fractals, 140 (2020), 110238. https://doi.org/10.1016/j.chaos.2020.110238 doi: 10.1016/j.chaos.2020.110238
    [45] B. Han, D. Jiang, B. Zhou, T. Hayat, A. Alsaedi, Stationary distribution and probability density function of a stochastic SIRSI epidemic model with saturation incidence rate and logistic growth, Chaos Solitons Fractals, 142 (2021), 110519. https://doi.org/10.1016/j.chaos.2020.110519 doi: 10.1016/j.chaos.2020.110519
  • This article has been cited by:

    1. Shirali Kadyrov, Farkhod Haydarov, Khudoyor Mamayusupov, Komil Mustayev, Endemic coexistence and competition of virus variants under partial cross-immunity, 2025, 33, 2688-1594, 1120, 10.3934/era.2025050
  • Reader Comments
  • © 2023 the Author(s), licensee AIMS Press. This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0)
通讯作者: 陈斌, bchen63@163.com
  • 1. 

    沈阳化工大学材料科学与工程学院 沈阳 110142

  1. 本站搜索
  2. 百度学术搜索
  3. 万方数据库搜索
  4. CNKI搜索

Metrics

Article views(1487) PDF downloads(126) Cited by(1)

Figures and Tables

Figures(6)  /  Tables(1)

Other Articles By Authors

/

DownLoad:  Full-Size Img  PowerPoint
Return
Return

Catalog