Research article Special Issues

Dynamics and numerical simulations of a generalized mosquito-borne epidemic model using the Ornstein-Uhlenbeck process: Stability, stationary distribution, and probability density function

  • In this paper, we proposed a generalized mosquito-borne epidemic model with a general nonlinear incidence rate, which was studied from both deterministic and stochastic insights. In the deterministic model, we proved that the endemic equilibrium was globally asymptotically stable when the basic reproduction number R0 was greater than unity and the disease free equilibrium was globally asymptotically stable when R0 was lower than unity. In addition, considering the effect of environmental noise on the spread of infectious diseases, we developed a stochastic model in which the infection rates were assumed to satisfy the mean-reverting log-normal Ornstein-Uhlenbeck process. For this stochastic model, two critical values, known as Rs0 and RE0, were introduced to determine whether the disease will persist or die out. Additionally, the exact probability density function of the stationary distribution near the quasi-equilibrium point was obtained. Numerical simulations were conducted to validate the results obtained and to examine the impact of stochastic perturbations on the model.

    Citation: 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[J]. Electronic Research Archive, 2024, 32(6): 3777-3818. doi: 10.3934/era.2024172

    Related Papers:

    [1] Yifan Wu, Xiaohui Ai . Analysis of a stochastic Leslie-Gower predator-prey system with Beddington-DeAngelis and Ornstein–Uhlenbeck process. Electronic Research Archive, 2024, 32(1): 370-385. doi: 10.3934/era.2024018
    [2] 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
    [3] Wenjie Zuo, Mingguang Shao . Stationary distribution, extinction and density function for a stochastic HIV model with a Hill-type infection rate and distributed delay. Electronic Research Archive, 2022, 30(11): 4066-4085. doi: 10.3934/era.2022206
    [4] Jitendra Singh, Maninder Singh Arora, Sunil Sharma, Jang B. Shukla . Modeling the variable transmission rate and various discharges on the spread of Malaria. Electronic Research Archive, 2023, 31(1): 319-341. doi: 10.3934/era.2023016
    [5] Bing Zhao, Shuting Lyu, Qimin Zhang . Dynamics and density function for a stochastic anthrax epidemic model. Electronic Research Archive, 2024, 32(3): 1574-1617. doi: 10.3934/era.2024072
    [6] 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
    [7] Bo Zheng, Lijie Chang, Jianshe Yu . A mosquito population replacement model consisting of two differential equations. Electronic Research Archive, 2022, 30(3): 978-994. doi: 10.3934/era.2022051
    [8] Qin Guo, Binlei Cai . Learning capability of the rescaled pure greedy algorithm with non-iid sampling. Electronic Research Archive, 2023, 31(3): 1387-1404. doi: 10.3934/era.2023071
    [9] Mingzhan Huang, Xiaohuan Yu, Shouzong Liu . Modeling and analysis of release strategies of sterile mosquitoes incorporating stage and sex structure of wild ones. Electronic Research Archive, 2023, 31(7): 3895-3914. doi: 10.3934/era.2023198
    [10] Feng Wang, Taotao Li . Dynamics of a stochastic epidemic model with information intervention and vertical transmission. Electronic Research Archive, 2024, 32(6): 3700-3727. doi: 10.3934/era.2024168
  • In this paper, we proposed a generalized mosquito-borne epidemic model with a general nonlinear incidence rate, which was studied from both deterministic and stochastic insights. In the deterministic model, we proved that the endemic equilibrium was globally asymptotically stable when the basic reproduction number R0 was greater than unity and the disease free equilibrium was globally asymptotically stable when R0 was lower than unity. In addition, considering the effect of environmental noise on the spread of infectious diseases, we developed a stochastic model in which the infection rates were assumed to satisfy the mean-reverting log-normal Ornstein-Uhlenbeck process. For this stochastic model, two critical values, known as Rs0 and RE0, were introduced to determine whether the disease will persist or die out. Additionally, the exact probability density function of the stationary distribution near the quasi-equilibrium point was obtained. Numerical simulations were conducted to validate the results obtained and to examine the impact of stochastic perturbations on the model.



    Mosquito-borne infectious diseases, a common type of vector-borne infections diseases, are primarily caused by pathogens, which are transmitted from mosquitoes to humans or other animals [1]. Mosquitoes are one of the most widely distributed vectors in the world, transmitting a variety of parasites and viruses. While the vector organisms may not develop the disease themselves, they serve as a means for the pathogen to spread among hosts. With mosquitoes found in various locations from tropical to temperate zones, mosquito-borne infections have a wide and diverse geographic distribution [2]. Diseases such as malaria, lymphatic filariasis, West Nile Virus, Zika virus, and dengue fever are commonly transmitted by mosquitoes, with no specific vaccine or medication available for treatment. These diseases pose a significant threat to human health and socio-economic development worldwide.

    Malaria is an infectious disease caused by a protozoan parasite that is mainly transmitted to humans through the bite of a mosquito [3]. The disease is widespread in tropical and subtropical regions, particularly in poorer areas of Africa, Asia, and Latin America. Globally, malaria remains a major public health problem, resulting in millions of infections and hundreds of thousands of deaths annually. Apart from malaria, as the classic mosquito-borne infectious disease, dengue fever also attracts considerable attention. Dengue fever is a serious infectious disease caused by the dengue virus, primarily transmitted to humans by the Aedes aegypti mosquito, but also through blood transfusion, organ transplantation, and vertical transmission [4]. In addition, both yellow fever and West Nile disease are also caused by mosquitoes carrying the corresponding viral infections [5,6]. Yellow fever is a rare disease among U.S. travelers. Conversely, West Nile virus is the primary mosquito-borne disease in the continental U.S., commonly transmitted through the bite of an infected mosquito.

    Research on mosquito-borne diseases has proliferated [7,8,9,10]. Mathematical models have been increasingly used for experimental and observational studies of different biological phenomena, and a wide range of techniques and applications have been developed to study epidemic diseases. For example, Newton and Reiter [8] developed a deterministic Susceptibility, Exposure, Infection, Resistance, and Removal (SEIR) model of dengue transmission to explore the behavior of epidemics and realistically reproduce epidemic transmission in immunologically unmade populations. Moreover, Pandey et al. [9] proposed a Caputo fractional order derivative mathematical model of dengue disease to study the transmission dynamics of the disease and to make reliable conclusions about the behavior of dengue epidemics. In addition, in order to investigate the effect of the vector on the dynamics of the disease, Shi and Zhao et al. [10] proposed a differential system with saturated incidence to model a vector-borne disease

    {˙SH=μK+dIHμSH(¯β1IV1+η1IV+¯β2IH1+η2IH)SH,˙IH=(¯β1IV1+η1IV+¯β2IH1+η2IH)SH(d+μ+γ)IH,˙RH=γIHμRH,˙SV=Λβ3IHSV1+η3IHmSV,˙IV=β3IHSV1+η3IHmIV, (1.1)

    the biological significance of the model (1.1) parameters are shown in Table 1. In this model, there are two populations, namely mosquito vector population and human host population. The mosquito vector population is divided into two categories, SV and IV, and N=SV+IV, while human host population is divided into three categories, SH,IH and RH. According to [10], it is reasonable to assume that the total number of human K=SH+IH+RH is a positive constant.

    Table 1.  Variables and parameters in model (1.1).
    Variables and Parameters Description
    SH number of the susceptible human host
    IH number of the infected human host
    RH number of the recovered human host
    K sum of the total human host
    SV density of the susceptible mosquito vectors
    IV density of the infected mosquito vectors
    N sum of the total mosquito vectors density
    ¯β1 biting rate of an infected vector on the susceptible human
    ¯β2 infection incidence between infected and susceptible hosts
    β3 infection ratio between infected hosts and susceptible vectors
    η1 determines the level at which the force of infection saturates
    η2 determines the level at which the force of infection saturates
    η3 determines the level at which the force of infection saturates
    γ the conversion rate of infected hosts to recovered hosts
    μ natural death rate of human
    Λ birth or immigration of human
    m natural death rate of mosquito vectors
    d disease-induced mortality of infected hosts

     | Show Table
    DownLoad: CSV

    Obviously, we can get that there always exists a compact positively invariant set for model (1.1) as follows

    Γ0={(SH,IH,RH,SV,IV)R5+:SH+IH+RHK,SV+IVΛm}. (1.2)

    The incidence rate has various forms and plays an important role in the study of epidemic dynamics. In addition to the saturated incidence used in model (1.1) of this paper, other forms of the incidence rate have been widely used. For example, Chong, Tchuenche, and Smith [11] studied a mathematical model of avian influenza with half-saturated incidence rate SbIbHb+Ib, SHIaHa+Ia and ShImHm+Im. In addition, Li et al. [12] also carried out numerical analysis of friteral order pine wilt disease model with bilinear incidence rate ShIv and IhSv. In order to make the model (1.1) have wider research significance and apply to more infectious diseases, we consider replacing the saturated incidence rate with the general incidence rate ϕ1(IV), ϕ2(IH) and ϕ3(IH), and then give the epidemic model with the general incidence as follows

    {˙SH=μ(KSH)¯β1ϕ1(IV)SH¯β2ϕ2(IH)SH+dIH,˙IH=¯β1ϕ1(IV)SH+¯β2ϕ2(IH)SHωIH,˙RH=γIHμRH,˙SV=Λβ3Φ3(IH)SVmSV˙IV=β3ϕ3(IH)SVmIV. (1.3)

    Furthermore, the incidence rate in model (1.3) are assumed to meet the following conditions

    (A1) ϕ1(0)=ϕ2(0)=ϕ3(0)=0,

    (A2) ϕ1(IV)0,ϕ2(IH)0,ϕ3(IH)0, IV,IH0,

    (A3) 0(IVϕ1(IV))m0,0(IHϕ2(IH))m0,0(IHϕ3(IH))m0, where m0 is a positive constant.

    By looking at the model (1.3), the following equation is valid, dNdt=ΛmN, that indicates, N(t)Λm as t. Note that SH+IH+RH=K, SV+IV=Λm, this means that RH=KSHIH, SV=ΛmIV, let ω=d+μ+γ, thus the host population and pathogen population system are equivalent to the following system

    {˙SH=μ(KSH)¯β1ϕ1(IV)SH¯β2ϕ2(IH)SH+dIH,˙IH=¯β1ϕ1(IV)SH+¯β2ϕ2(IH)SHωIH,˙IV=β3ϕ3(IH)(ΛmIV)mIV. (1.4)

    Noise is ubiquitous in real life, and the spread of infectious diseases will inevitably be affected by the environment and other external factors. With the unpredictable environment, some key parameters in the infectious disease model are inevitably affected by external environmental factors. Therefore, in order to more accurately describe the transmission process, we use a stochastic model to describe and predict the epidemic trend of diseases. For the perturbation term of the parameter, two methods are commonly used, including linear function of Gaussian noise and mean-reverting stochastic process [13,14,15,16,17]. For Gaussian noise, when the time interval is very small, the variance of the parameter will become infinite, indicating that the parameter has changed greatly in a short time, which is unreasonable [18]. So we mainly consider the mean-reverting Ornstein–Uhlenbeck process. Let βi(i=1,2) take the following form

    dβi(t)=αi(¯βiβi(t))dt+σidBi(t), (1.5)

    where αi represent the speed of reversal and σi represent the intensity of fluctuation. Solve the Eq (1.5), we can get

    βi(t)=eαitβi(0)+¯βi(1eαit)+σit0eαi(ts)dBi(s),

    where βi(0) is the initial value of βi(t). For arbitrary initial value βi(0), βi(t) follows a Gaussian distribution βi(t)N(¯βi,σ2i2αi)(t). Furthermore, setting βi(0)=¯βi, then the average value of βi(t) satisfies

    ¯βi(t)=1tt0βi(s)ds=¯βi+1tt0σiαi(1eαi(st))dBi(s),

    it is known that the mathematical expectation and variance of βi(t) are ¯βi and σ2it3+O(t2), respectively, where O(t2) is the higher order infinitesimal of t2. Obviously, the variance becomes zero instead of infinity as t0. This shows the universality of the Ornstein–Uhlenbeck process. Moreover, in order to ensure the positivity of the parameter values after adding the perturbation, the log-normal Ornstein-Uhlenbeck process for the noise perturbation to the transmission rates β1 and β2 of the system (1.4) is used, then following stochastic model is obtained

    {dSH(t)=[μ(KSH(t))β1ϕ1(IV(t))SH(t)β2(t)ϕ2(IH(t))SH(t)+dIH(t)]dt,dIH(t)=[β1(t)ϕ1(IV(t))SH(t)+β2(t)ϕ2(IH(t))SH(t)ωIH(t)]dt,dIV(t)=[β3ϕ3(IH(t))(ΛmIV(t))mIV(t)]dt,dlogβ1(t)=α1(log¯β1logβ1(t))dt+σ1dB1(t),dlogβ2(t)=α2(log¯β2logβ2(t))dt+σ2dB2(t). (1.6)

    In this paper, we extend the saturated incidence rate of model (1.1) to the general incidence ϕ1(IV), ϕ2(IH) and ϕ3(IH) to obtain models (1.3) and (1.4), and investigate the global asymptotic stability of the equilibrium point of model (1.3). Furthermore, we choose to modify the parameter β1 and β2 to satisfy the log-normal Ornstein-Uhlenbeck process to obtain the stochastic model (1.6), and study its stationary distribution, exponential extinction, probability density function near the quasi-equilibrium point and other dynamic properties.

    The rest of this article is organized as follows. In Section 2, some necessary mathematical symbols and lemmas are introduced. In Section 3, some conclusions of deterministic model (1.3) are obtained and the global stability of equilibrium point in this model is proved. In Section 4, we obtain some theoretical results for the stochastic system (1.6) where we prove the existence of a unique global positive solution for the stochastic system (1.6). In addition, through the ergodic properties of parameters βi(t),i=1,2 and the construction of a series of suitable Lyapunov functions, sufficient criterion for the existence of stationary distribution is obtained, which indicates that the disease in the system will persist. Next, we have sufficient conditions for the disease to go extinct. Further, we solve the corresponding matrix equation to obtain an expression for the probability density function near the quasi-local equilibrium point of the stationary distribution. Next, in Section 5, some theoretical results are verified by several numerical simulations. Finally, several conclusions are given in Section 6.

    To make it easier to understand, denote Rn+={(y1,y2,...yn)Rn|yj>0,1jn}. In represents the n-dimensional unit matrix. IA denotes the indicator function of set A, and it means that when xA,IA=1, otherwise, IA=0. If A is a matrix or vector, then AT stands for its inverse matrix, and A1 stands for its inverse matrix.

    Lemma 2.1. (Itˆo's formula [19]) Consider the n-dimensional stochastic differential equation

    dx(t)=f(t)dt+g(t)dB(t), (2.1)

    where B(t)=(B1(t),B2(t),...,Bn(t)) and it represents n-dimensional Brownian motion defined on a complete probability space, let L act on a function VC2,1(Rn×R+;R), then we have

    dV(x(t),t)=LV(x(t),t)dt+Vx(x(t),t)g(t)dB(t),a.s.,

    where

    LV(x(t),t)=Vt(x(t),t)+Vx(x(t),t)f(t)+12trace(gT(t)Vxx(x(t),t)g(t)),

    it represents the differential operator, and

    Vt=Vt,Vx=(Vx1,...,Vxn),Vxx=(2Vxixj)n×n.

    Lemma 2.2. (Ma et al. [20]) Letting ϕ(λ)=λn+a1λn1+a2λn2++an1λ+an is the characteristic polynomial of the square matrix A, the matrix A is called a Hurwitz matrix if and only if all characteristic roots of A are negative real parts, that is equivalent to the following conditions being true

    Hk=|a1a3a5a2k11a2a4a2k20a1a3a2k301a2a2k4000ak|>0,

    k=1,2,,n, among them j>n, replenishing definition aj=0.

    Lemma 2.3. ([21]) For five-dimension algebraic equation G20+LΘ+ΘLT=0, where Θ is a symmetric matrix, G0=diag(1,0,0,0,0) and

    L=(l1l2l3l4l51000001000001000000l6).

    If l1>0,l2>0,l3>0,l4>0 and l1l2l3l23l21l4>0, then the symmetric matrix Θ is a positive semi-definite matrix. Thus, we have

    Θ=(l1l4l2l3l0l3l000l3l0l1l0l3l0l1l000l1l0l3l1l2ll4000000),

    where l=2(l4l21l1l2l3+l23).

    Lemma 2.4. ([22,23]) For n-dimension stochastic process (1.6), X(t)Rn and its initial value X(0)Rn, if there is a bounded closed domain U in Rn with a regular boundary and

    lim inft+1tt0P(τ,X(0),U)dτ>0,a.s.,

    in which P(τ,X(0),U) represents the transition probability of X(t), then X(t) has an invariant probability measure on Rn, then it admits at least one stationary distribution.

    In this section, we focus on the local stability of the equilibrium point of the deterministic model (1.3). Initially, we verify the existence and uniqueness of equilibria model (1.3). We can calculate the basic reproduction number of the deterministic model (1.3) by the next generation method [24], define

    F=(¯β2ϕ2(0)K¯β1ϕ1(0)K00),V=(ω0β3Λϕ3(0)mm),

    therefore, the next generation matrix is

    FV1=(¯β2Kϕ2(0)ω+¯β1β3ΛKϕ1(0)ϕ3(0)m2ω¯β1ϕ1(0)Km00),

    then, the basic reproduction number for system (1.3) is obtained

    R0=ρ(FV1)=¯β2Kϕ2(0)ω+¯β1β3ΛKϕ1(0)ϕ3(0)m2ω.

    Based on the key value of the basic reproduction number R0, the conditions for the existence of local equilibrium point for the model (1.3) can be found.

    Theorem 3.1. The disease-free equilibrium E0 of model (1.3) is E0=(SH0, 0, 0, SV0, 0)=(K, 0, 0, Λm, 0) which always exists. If R0>1, there is a unique local equilibrium E=(SH, IH, RH, SV, IV).

    After finding the conditions for the existence of the equilibrium point of model (1.3), next we verify that the global stability of equilibria.

    Theorem 3.2. (i) If R0<1, the disease-free equilibrium point E0 is globally asymptotically stable. If R0>1, E0 is unstable. (ii) If R0>1, the endemic equilibrium point E is globally asymptotically stable.

    Remark 3.2 The global stability of E0 in this theorem can be referred to the method in the literature [10]. By constructing a series of suitable Lyapunov functions, we prove the global stability of E.

    Initially, we verify that the stochastic model (1.6) has a unique global positive solution. This provides preparation for the dynamic behavior of the model. For model (1.6), it is easy to see that

    Γ={(SH,IH,IV,β1,β2)R5+:SH+IH<K,IV<Λm}

    is the positive invariant set, and the subsequent research will be discussed in Γ.

    Theorem 4.1. For any initial value (SH(0),IH(0),IV(0),β1(0),β2(0))Γ, there exists a unique solution (SH(t),IH(t),IV(t),β1(t),β2(t)) of system (1.6) and the solution will remain in Γ with probability one (a.s.).

    The stationary distribution of stochastic model (1.6) plays a key role in regulating the dynamics of disease and analyzing the sustainable development of disease. Next, sufficient conditions for the existence of stationary distribution will be obtained. We define

    Rs0=~β2Kϕ2(0)ω+~β1β3ΛKϕ1(0)ϕ3(0)m2ω,

    where

    ~β1=¯β1eσ2120α1,~β2=¯β2eσ2212α2.

    Theorem 4.2. If Rs0>1, then stochastic system (1.6) has a stationary distribution.

    Remark 4.2 The Theorem 4.2 is proved by constructing the Lyapunov functions. It is observed that when Rs0>1, a stationary distribution exists and the disease will be endemic for a long period of time.

    Furthermore, disease propagation and extinction are two major areas of research in stochastic system dynamics. After establishing the conditions under which a disease reaches a stable state, it is also essential to understand the conditions under which the disease becomes extinct. In addition, we discuss the sufficient condition for disease extinction in the model (1.6), define

    RE0=R0+mR0(eσ21α12eσ214α1+1)12+ϕ2(0)K¯β2(eσ22α22eσ224α2+1)12min{m,¯β2ϕ2(0)KR0}.

    Theorem 4.3. If RE0<1, then the disease of the system (1.6) will become exponentially extinct with probability 1.

    Remark 4.3 Theorem 4.3 gives the sufficient condition for the exponential extinction of diseases IH, IV. From the expressions of Rs0 and RE0, the relationships Rs0R0 and RE0R0 are deduced, and the equal sign holds if and only if σ1=σ2=0.

    Additionally, the density function of a continuous distribution is essential to understanding a stochastic system, making the precise determination of its expression a crucial challenge. Through the matrix analysis method, we have successfully derived the expression for the probability density function in the vicinity of the equilibrium point for model (1.6). Linearize the model (1.6) before calculate the probability density function. Define a quasi-endemic equilibrium point P=(SH,IH,IV,logβ1,logβ2), it satisfies

    {μKμSβ1ϕ1(IV)SHβ2ϕ2(IH)S+dIH=0,β1ϕ1(IV)SH+β2ϕ2(IH)SωIH=0,β3ϕ3(IH)(ΛmIV)mIV=0,α1(log¯β1logβ1)=0,α2(log¯β2logβ2)=0. (4.1)

    Let (z1,z2,z3,x1,x2)T=(SHSH,IHIH,IVIV,logβ1logβ1,logβ2logβ2)T, then system (4.1) can be linearized around P as follows

    {dz1=[a11z1a12z2a13z3a14x1a15x2]dt,dz2=[a21z1a22z2+a13z3+a14x1+a15x2]dt,dz3=[a32z2a33z3]dt,dx1=a44x1dt+σ1dB1(t),dx2=a55x2dt+σ2dB2(t), (4.2)

    where a11=μ+¯β1ϕ1(IV)+¯β2ϕ2(IH), a12=¯β2ϕ2(IH)SHd, a13=¯β1ϕ1(IV)SH, a14=¯β1ϕ1(IV)SH, a15=¯β2ϕ2(IH)SH, a21=¯β1ϕ1(IV)+¯β2ϕ2(IH), a22=ω¯β2ϕ2(IH)SH, a32=β3ϕ3(IH)(ΛmIV), a33=β3ϕ3(IH)+m, a44=α1, a55=α2. Apparently, aij>0.

    By denoting

    A=(a11a12a13a14a15a21a22a13a14a150a32a3300000a4400000a55), G=(000σ1σ2),
    B(t)=(0,0,0,B1(t),B2(t))T, Z(t)=(z1,z2,z3,x1,x2)T.

    Then system (4.1) can be expressed as

    dZ(t)=AZ(t)dt+GdB(t), (4.3)

    then, the solution of (4.3) can be calculated as

    X(t)=eAtX(0)+t0eA(ts)GdB(t),

    since t0eA(ts)GdB(t) obeys a normal distribution N(0,ˆΣ(t)) at time t, where ˆΣ(t)=t0eA(ts)GGTeAT(ts)ds, then, we can get X(t)N(eAtX(0),ˆΣ(t)).

    First, we need to verify that matrix A is Hurwitz matrix [25], the characteristic polynomial of A can be obtained as follows

    φA(λ)=|λ+a11a12a13a14a15a21λ+a22a13a14a150a32λ+a3300000λ+a4400000λ+a55|=(λ+a44)(λ+a55)|λ+a11a12a13a21λ+a22a230a32λ+a33|.

    Obviously there are λ1=a44,λ2=a55, other characteristic roots can also be found with negative real part according to Section 3. Therefore, matrix A is a Hurwitz matrix by Lemma 2.2. According to the stability theory of zero solution to the general linear equation [20], we have

    limt+eAtX(0)=0,
    Σlimt+ˆΣ=limt+t0eA(ts)GGTeAT(ts)ds=+0eATtG2eAtdt.

    Based on the solution of Gardiner [26], it follows

    G2+AΣ+ΣAT=0. (4.4)

    Second, we will solve the Eq (4.4) to get the exact expression of probability density function near the quasi-endemic equilibrium.

    Theorem 4.4. For any initial value (SH(0),IH(0),IV(0),β1(0),β2(0))Γ, if R0>1 and mμ+β3ϕ3(IH)>0, then the stationary distribution of stochastic system (1.6) near the P approximately admits a normal probability density function as follows

    Φ(SH,IH,IV,logβ1,logβ2)=(2π)52|Σ|12exp[12(SHSH,IHIH,IVIV,logβ1log¯β1,logβ2log¯β2)Σ1(SHSH,IHIH,IVIV,logβ1log¯β1,logβ2log¯β2)T].

    The exact expression of covariance matrix Σ is shown in the proof.

    Remark 4.4 In Theorem 4.4, by defining the quasi-endemic equilibrium P, we derive an exact expression of probability density function of the stationary distribution around a quasi-positive equilibrium P.

    In order to illustrate the above theoretical results, we perform several numerical simulations in this section. Consider bilinear incidence rate ϕ1(IV)=IV, ϕ2(IH)=IH, ϕ3(IH)=IH, letting xi(t)=logβi(t)logˉβi,i=1,2, then according to the method in [27], the following is the corresponding discretized equation for the system (1.6)

    {Sj+1H=SjH+[μ(KSjH)¯β1exj1IjVSjH+¯β2exj2IjHSjH+dIjH]Δt,Ij+iH=IjH+[¯β1exj1IjVSjH+¯β2exj2IjHSjHωIjH]Δt,Ij+1V=IjV+[β3IjH(ΛmIjV)mIjV]Δt,xj+11=xj1α1xj1Δt+σ1Δtξ1,j+σ212(ξ21,j1)Δt,xj+12=xj2α2xj2Δt+σ2Δtξ2,j+σ222(ξ22,j1)Δt,

    let ξi,j be random variables that follow a Gaussian distribution N(0,1) for i=1,2 and j=1,2,...,n. The time interval is denoted by Δt>0. The values (SjH,IjH,IjV,xj1,xj2) correspond to the j-th iteration of the discretization equation.

    First and foremost, taking into consideration the importance of parameter selection, rationality, and the visual effectiveness of theoretical results, we choose the following appropriate parameters, referring to [10,28,29], and denoted them as Number 1

    μ=0.05,K=100,¯β1=0.15,¯β2=0.1,β3=0.1,ω=0.8,d=0.5,
    γ=0.25,Λ=5,m=0.5,α1=0.8,α2=0.8,σ1=0.1,σ2=0.1,

    after calculation, we can get mμ+β3IH=2.0391>0 and the indexes of deterministic system and stochastic system can be obtained respectively, as shown below

    R0=¯β2Kω+¯β1β3Λkm2ω=50>1, Rs0=~β2Kω+~β1β3Λkm2ω=50.0365>1,

    which satisfies the condition in Theorem 4.2. More importantly, we can calculate the quasi-endemic equilibrium and the covariance matrix Σ, which have the following forms

    (SH, IH, IV, logβ1, logβ2)=(4.6565, 15.8906, 7.6066, log0.15, log0.1),
    Σ=(0.05290.03770.00380.00930.01300.03770.03520.00310.00720.01000.00380.00310.00040.00060.00080.00930.00720.00060.006200.01300.01000.000800.0062).

    Therefore, we can get that the solution (SH(t),IH(t),IV(t),β1(t),β2(t)) obeys the normal density function

    Φ(SH,IH,IV,β1,β2)N((4.6565, 15.8906, 7.6066,0.15,0.1)T,Σ).

    The marginal density functions are as follows

    ΦSH=1.7345e9.4518(SH4.6565)2, ΦIH=2.1264e14.2045(IH15.8906)2, ΦIV=19.9471e1250(IV7.6066)2.

    Using the above parameters, we can get trajectories of SH(t), IH(t) and IV(t) respectively, which are present in Figure 1. It is used to represent the variation of the solution (SH(t),IH(t),IV(t)) in the deterministic model (1.4) and the stochastic model (1.6). Frequency histograms and marginal density function curves for SH(t),IH(t) and IV(t) are also given in the right column of the Figure 1.

    Figure 1.  The left and right columns show the trajectories of the solutions (SH(t),IH(t),IV(t)) of the stochastic and deterministic systems under perturbations σ1=0.1,σ2=0.1, as well as histograms of the solutions and the marginal density functions, respectively.

    In addition, the frequency fitted density functions and the marginal density functions for SH(t),IH(t) and IV(t) are given in Figure 2, respectively, which are highly consistent. Therefore, we deduce that the solutions (SH(t),IH(t),IV(t)) have a smooth distribution and their density functions follow a normal distribution. As we can see, the disease eventually spreads, which is consistent with Theorem 4.2.

    Figure 2.  The frequency fitting density functions and marginal density functions of SH(t),IH(t) and IV(t), respectively.

    On the other hand, we select that a part of parameters are shown below, and the remaining parameters are consistent with Number 1, ¯β1=0.015,¯β2=0.001,β3=0.015. These are denoted as Number 2. The crucial value RE0 takes the form

    RE0=R0+(eσ21α12eσ214α1+1)12+K¯β2m(eσ22α22eσ224α2+1)12=0.7986<1

    which satisfies the condition in Theorem 4.3. Figure 3 represents the trajectory of the solution (SH(t),IH(t),IV(t)), and it is clearly visible that the eventual trend of the disease is towards extinction.

    Figure 3.  The trajectory of solution (SH(t),IH(t),IV(t)) under the condition RE0<1.

    Further, by choosing the parameters in Numbers 1 and 2, respectively, the left panel in Figure 4 satisfies the condition Rs0>1 and the right panel satisfies RE0<1. It can be seen that the disease exhibits a trend towards stabilization and extinction as conditions Theorems 4.2 and 4.3 are satisfied, respectively.

    Figure 4.  The left panel shows the trajectories of the solution (SH(t),IH(t),IV(t)) in the stochastic model (1.6) for Rs0>1, and the right panel shows the trajectories of the solutions of the stochastic model (1.6) for RE0<1.

    Now, we study the effect of perturbations for a mosquito-borne epidemic model. Assuming that all parameters take the values in Number 1, we choose different reversion speed α and volatility intensity σ to plot the graphs, respectively. Taking α1=α2=0.8 and different volatility intensity as shown in the Figure 5, the icons are red line σi=0.05, blue line σi=0.1 and green line σi=0.15, i=1,2, the trends of the solution (SH(t),IH(t),IV(t)) of the stochastic model (1.6) are represented by the figure. It shows that the fluctuation decrease as the volatility intensity decreases. Then, we set the volatility intensity σ1=σ2=0.1, the reversion speed αi=0.1 as shown by the red line in the figure, the blue line shows αi=1.0, and similarly, the green line indicates αi=1.5, i=1,2, then the same insightful changes in Figure 6 indicate that the fluctuation decreases with the increase of the reversion speed.

    Figure 5.  Trajectory plots of the solution (SH(t),IH(t),IV(t)) of the stochastic model (1.6) at the reversion speed αi=0.8, i=1,2 and different volatility intensity is shown in the icon with the red line σ1=σ2=0.05, the blue line σ1=σ2=0.1 and the green line σ1=σ2=0.15.
    Figure 6.  Trajectory plots of the solution (SH(t),IH(t),IV(t)) of the stochastic model (1.6) at the volatility intensity σi=0.1, i=1,2 and different reversion is shown in the icon with the red line α1=α2=0.1, the blue line α1=α2=1 and the green line α1=α2=1.5.

    Further, we make the rest of the parameter assumptions consistent with Numbers 1 and 2, respectively, except for the the volatility intensity and reversion speed. Figures 7 and 8 depict the trends of R0,Rs0, and RE0 under different volatility intensity and different reversion speed, respectively, and the range of the two variables we choose is [0,1]. Combining the information in the two figures, it can be concluded that higher reversion speed and lower volatility intensity can make RE0 and Rs0 smaller.

    Figure 7.  Trend plots of R0,Rs0, and RE0 at fixed reversion speed α1=α2=0.8 and volatility intensity σ1[0.01,1], σ2=0.1. The rest of the parameter values in the left figure are consistent with those in Number 1, and the rest of the parameter values in the right figure are consistent with those in Number 2.
    Figure 8.  Trend plots of R0,Rs0, and RE0 at fixed volatility intensity σ1=σ2=0.1 and reversion speed α1[0.01,1], α2=0.8. The rest of the parameter values in the left figure are consistent with those in Number 1, and the rest of the parameter values in the right figure are consistent with those in Number 2.

    Next, we continue to discuss the effects of reversion speed α and volatility intensity σ on Rs0 and RE0 and their magnitude relationships under different conditions.

    (i) Assuming that α1, σ1 are the variables, and the other parameters are consistent with Number 1, The Figure 9 shows the three-dimensional chromatograms of Rs0 and RE0, which are consistent with the results of Figures 7 and 8, where Rs0 increases with increasing σ1, and decreases with increasing α1. This indicates that the disease stabilizes as the reversion speed decreases or volatility intensity increases. In addition, it can be seen that both Rs0 and RE0 are greater than 1 and RE0 is greater than Rs0 in the range where α1 belongs to [0.5,0.6] and σ1 belongs to [0.005,0.09];

    Figure 9.  Color plot of the trend of RE0 and Rs0 with variables (α1,σ1)[0.5,0.6]×[0.005,0.09], with the rest of the parameters consistent with those in Number 1.

    (ii) Conditionally the same as in (i), we set the other parameters are consistent with Number 2, it is noted that RE0 increases with the increase of σ1 and decreases with the increase of α1 in Figure 10, implying that the diseases in the stochastic model (1.6) tend to become extinct when the volatility intensity decreases. Moreover, in the parameter range of the plot, both Rs0 and RE0 are less than 1, and RE0 is greater than Rs0.

    Figure 10.  Color plot of the trend of RE0 and Rs0 with variables (α1,σ1)[0.5,0.6]×[0.005,0.09], with the rest of the parameters consistent with those in Number 2.

    Next, we will discuss the mean first passage time, the moment a stochastic process first transitions from one state to another is termed the first passage time (FPT) [30]. The mean first passage time (MFPT) is then defined as the average of these first passage times [31]. Starting with an initial value of (SH(0),IH(0),IV(0)), we aim to examine the time it takes for the system to evolve from this initial state to either a stationary state (MFPT1) or to an extinction state (MFPT2).

    Then we define τ1 as the FPT from the initial state to the persistent state, and τ2 as the FPT from the initial state to the extinct state

    τ1=inf{t:SH<SH,IH>IH,IV>IV},
    τ2=inf{t:IH0<0.0001,IV0<0.0001}.

    Then we have

    MFPT1=E(τ1), MFPT2=E(τ2).

    Using Monte Carlo numerical simulation method, if SH(nΔt)<SH, IH(nΔt)>IH, IV(nΔt)>IV, then τ1=nΔt, assuming that the number of simulations is N, then

    MFPT1=Ni=1niΔtN.

    Similarly, if IH0<0.0001 and IV0<0.0001, then τ2=mΔt and

    MFPT2=Ni=1niΔtN.

    Here, we set N = 2000, σi and αi, i=1,2 are random variables. Figures 11 and 12 depict the relationship between MFPT1 and MFPT2 and the speed of reversion αi and the volatility intensity σi, i=1,2 in stochastic system (1.6) with the bilinear incidence rate, respectively. Figure 11 reveals the values of MFPT1 with N=2000,σi[0.01,0.1] and αi=4,5,6, i=1,2, respectively. It shows that MFPT1 decreases with decreasing reversion speed αi or increasing voltility intensity σi, implying that the disease is much easier to arrive the stable state. Similarly, Figure 12 shows the trend of MFPT2 at N=2000, σi[0.01,0.1] and αi=0.1,0.15,0.2. Through the figure it can be noted that MFPT2 increases with αi decrease and σi increase.

    Figure 11.  The mean first passage time for transitioning from the initial state values (SH(0),IH(0),IV(0)) = (3, 1, 1) to the state of stationary with σi[0.01,0.1], αi=4,5,6, i=1,2. The other fixed parameter values are consistent with those in Number 1.
    Figure 12.  The mean first passage time for transitioning from the initial state values (SH(0),IH(0),IV(0)) = (5,200, 10) to the state of extinction with σi[0.01,0.1], αi=0.1,0.15,0.2, i=1,2. The other fixed parameter values are consistent with those in Number 2.

    We mainly develop a stochastic model, coupled with the general incidence rate and Ornstein-Uhlenbeck process, to study the dynamic of infectious disease spread, which includes the stationary distribution and probability density function. In view of our analysis, we can draw the following conclusions

    (i) For the deterministic epidemic model, two equilibria and the basic reproduction number R0 are obtained, and their global asymptotic stability are deduced. Specificaly, the endemic equilibrium is globally asymptotically stable if R0>1, and the disease free equilibrium is globally asymptotically stable if R0<1.

    (ii) Considering that the spread of infectious disease is inevitably affected by environmental perturbations, we propose a stochastic model with general incidence and the Ornstein-Uhlenbeck process. By constructing a series of appropriate Lyapunov functions, the stationary distribution of model (1.6) is derived and we establish the sufficient criterion for the existence of the extinction. Specifically, the innovation of this paper is that we obtain a precise expression of the distribution around its quasi-positive equilibrium P by solving a difficult five-dimensional matrix equation, which is quite challenging.

    (iii) We also verify some conclusions of this paper by several numerical simulations.

    When Rs0>1, i.e., the parameters satisfy the condition of Theorem 4.2, we obtain trajectory plots of the solutions of the deterministic model (1.4) and the stochastic model (1.6), as well as the corresponding frequency histograms and edge density functions, and as shown in Figures 1 and 2, the disease eventually persists. This provides some verification of Theorem 4.2 of the theoretical results.

    Also, from Figure 3, we can further find that when RE0<1, the population strengths decrease with time and eventually converge to zero, which implies extinction of the disease. Figure 4 also further illustrates these points. The theoretical result of Theorem 4.3 is visualized through Figures 3 and 4.

    In addition, we also investigate the effect of perturbations and give Figures 5 and 6 to depict the effect of trends with different reversion speed and different volatility intensity. It can be seen that the fluctuation decreases as the volatility intensity decrease and the reversion speed increase.

    Then, correlation plots of Rs0, RE0 with reversion speed and volatility intensity are obtained in Figures 710. We conculde that higher reversion speed and lower volatility intensity can make RE0 and Rs0 more smaller.

    Finally, Figures 11 and 12 visually demonstrate the relationship between MFPT and the αi and σi, i=1,2 in a stochastic system (1.6) with bilinear incidence. We can see that if voltility intensity σi is much bigger (or the reversion speed αi is much smaller), then the disease is more easier to arrive the stable state. This is consistent with the results of the above conclusions.

    The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

    The authors would like to thank the anonymous referee for his/her useful suggestions. This work is supported by the Fundamental Research Funds for the Central Universities (No. 24CX03002A).

    The authors declare that they have no conflict of interest.

    Proof. The equilibria of system (1.3) satisfy

    {μ(KSH)¯β1ϕ1(IV)SH¯β2ϕ2(IH)SH+dIH=0,¯β1ϕ1(IV)SH+¯β2ϕ2(IH)SHωIH=0,γIHμRH=0,Λβ3Φ3(IH)SVmSV=0,β3ϕ3(IH)SVmIV=0. (A.1)

    Notice that

    RH=KSHIH, SV=ΛmIV,

    and we have

    {μ(KSH)¯β1ϕ1(IV)SH¯β2ϕ2(IH)SH+dIH=0,¯β1ϕ1(IV)SH+¯β2ϕ2(IH)SHωIH=0,β3ϕ3(IH)(ΛmIV)mIV=0. (A.2)

    Obviously E0=(SH0, 0, 0, SV0, 0)=(K, 0, 0, Λm, 0) always exists.

    On the other hand, we have

    {SH=K(1+γμ)IH,IV=β3ϕ3(IH)Λm2+mβ3ϕ3(IH)IH=¯β1ϕ1(IV)SHω+¯β2ϕ2(IH)SHω.

    Therefore, IH[0,μKμ+γ], let

    H(IH)=¯β1ϕ1(IV)SHω+¯β2ϕ2(IH)SHωIH,

    then by calculation, H(0)=0,H(μKμ+γ)=μKμ+γ<0, and

    H(IH)=¯β1ϕ1(IV)SHωm2β3ϕ3(IH)Λ(m2+mβ3ϕ3(IH))2(1+γμ)¯β1ϕ1(IV)ω+¯β2ϕ2(IH)SHω(1+γμ)¯β2ϕ2(IH)ω1.

    If R0>1, then

    H(0)=¯β2Kϕ2(0)ω+¯β1β3ΛKϕ1(0)ϕ3(0)m2ω1=R01>0,
    H(μKμ+γ)=(1+γμ)[¯β1ϕ1(IV(μKμ+γ))ω+¯β2ϕ2(μKμ+γ)ω]1<0,

    therefore, there exists a point ξ such that H(IH)>0 on [0,ξ) and H(IH)<0 on (ξ,μKμ+γ], i.e., H(IH) is monotonically increasing on [0,ξ) and monotonically decreasing on (ξ,μKμ+γ].

    Hence, there is a unique IH(0,μKμ+γ) such that H(IH)=0, which implies that when R0>1, system (1.4) has a unique endemic equilibrium E=(SH,IH,RH,SV,IV). This completes the proof.

    Proof. (i) Through calculation, Jacobian matrix of model (1.3) is obtained as follows

    J(SH,IH,RH,SVIV)=(μ¯β1ϕ1(IV)¯β2ϕ2(IH)d¯β2ϕ2(IH)SH00¯β1ϕ1(IV)SH¯β1ϕ1(IV)+¯β2ϕ2(IH)¯β2ϕ2(IH)SHω00¯β1ϕ2(IV)SH0γμ000β3Φ3(IH)SV0β3Φ3(IH)m00β3Φ3(IH)SV0β3Φ3(IH)m).

    Substituting E0 into the matrix J to get J0

    J0=(μd¯β2Kϕ2(0)00¯β1ϕ1(0)K0¯β2ϕ2(0)Kω00¯β1ϕ1(0)K0γμ000Λmβ3ϕ3(0)0m00Λmβ3ϕ3(0)00m),

    the corresponding characteristic polynomial is as follows

    ϕJ0(λ)=(λ+μ)(λ+μ)(λ+m)(λ+m)(λ+(ω¯β2Φ2(0))),

    obviously, if R0<1, then ω¯β2Φ2(0)>0, according to the Routh-Hurwitz criterion, there are only negative real part characteristic roots, so the disease-free equilibrium E0 is locally asymptotically stable.

    If R0>1, then ω¯β2Φ2(0)<0, this indicates that J0 has the eigenvalue of the positive real part, so the disease-free equilibrium E0 is unstable.

    Next we prove the global attractiveness of E0, define

    V=SHKKlogSHK+IH+SVΛmΛmlogmSVΛ+IV,

    Using Itˆo's formula for the above equation, we get

    LV=μK+dIHμSH¯β1Φ1(IV)SH¯β2Φ2(IH)SHμK2SHdIHKSH+¯β1Φ1(IV)K+¯β2Φ2(IH)K+μK+¯β1Φ1(IV)SH+¯β2Φ2(IH)SHωIH+Λβ3Φ3(IH)SVmSVΛ2mSV+β3Φ3(IH)Λm+Λ+β3Φ3(IH)SVmIV2μK+dIHμSHμK2SHdIHKSH+¯β1Φ1(0)IVK+¯β2Φ2(0)IHKωIHmSVΛ2mSV+β3Φ3(0)IHΛm+2ΛmIV,

    according to the method in [10], it is easy to see that SHK as t, SVΛm as t, and IH, IV0 as t. Then, LV0 and equal to 0 when it takes E0. Therefore, according to the LaSalle invariance principle, we conclude that when R0<1, E0 is globally asymptotically stable.

    (ii) According to the method in [32], the positive equilibrium point E of the model (1.3) satisfies the following system of equations

    {μK=μSH+(¯β1Φ1(IV)+¯β2Φ2(IH))SHdIH,(¯β1Φ1(IV)+¯β2Φ2(IH))SH=ωIH,γIH=μRH,Λ=β3Φ3(IH)SV+mSV,β3Φ3(IH)SV=mIV. (A.3)

    Define

    V1=SHSHSHlogSHSH+IHIHIHlogIHIH,

    after calculation

    d(SHlogSHSH)dt=SH(μKSH+μ+¯β1Φ1(IV)+¯β2Φ2(IH)dIHSH)=SH((μ+¯β1Φ1(IV)+¯β2Φ2(IH))SHdIHSH+μ+¯β1Φ1(IV)+¯β2Φ2(IH)dIHSH)=SH((μ+¯β1Φ1(IV)+¯β2Φ2(IH))SHSH+¯β1Φ1(IV)Φ1(IV)Φ1(IV)+¯β2Φ2(IH)Φ2(IH)Φ2(IH))+μSH(1SHSH)+SHSH(dIHdIH),
    d(IHlogIHIH)dt=¯β1Φ1(IV)SHSHΦ1(IV)IHSHΦ1(IV)IH¯β2Φ2(IH)SHSHΦ2(IH)IHSHΦ2(IH)IH+ωIH,
    d(SH+IH)dt=μKμSH+dIHωIH=(μ+¯β1Φ1(IV)+¯β2Φ2(IH))SHdIHμSH+dIH(¯β1Φ1(IV)+¯β2Φ2(IH))SHIHIH=μSH(1SHSH)+(¯β1Φ1(IV)+¯β2Φ2(IH))SH(1IHIH)d(IHIH),

    then we have

    dV1dt=μSH(1SHSH)+(¯β1Φ1(IV)+¯β2Φ2(IH))SH(1IHIH)d(IHIH)+ωIH+SH((μ+¯β1Φ1(IV)+¯β2Φ2(IH))SHSH+¯β1Φ1(IV)Φ1(IV)Φ1(IV)+¯β2Φ2(IH)Φ2(IH)Φ2(IH))+μSH(1SHSH)+SHSH(dIHdIH)¯β1Φ1(IV)SHSHΦ1(IV)IHSHΦ1(IV)IH¯β2Φ2(IH)SHSHΦ2(IH)IHSHΦ2(IH)IH=μSH(2SHSHSHSH)+¯β1Φ1(IV)SH[Φ1(IV)Φ1(IV)SHSHSHΦ1(IV)IHSHΦ1(IV)IHIHIH+3]+¯β2Φ2(IH)SH[Φ2(IH)Φ2(IH)SHSHSHΦ2(IH)IHSHΦ2(IH)IHIHIH+3]dIH(1SHSH)(1IHIH)μSH(2SHSHSHSH)+¯β1Φ1(IV)SH[Φ1(IV)Φ1(IV)logSHSHlogSHΦ1(IV)IHSHΦ1(IV)IHIHIH]+¯β2Φ2(IH)SH[Φ2(IH)Φ2(IH)logSHSHlogSHΦ2(IH)IHSHΦ2(IH)IHIHIH]=μ(SHSH)2SH+¯β1Φ1(IV)SH[Φ1(IV)Φ1(IV)logIHΦ1(IV)IHΦ1(IV)IHIH]+¯β2Φ2(IH)SH[Φ2(IH)Φ2(IH)logIHΦ2(IH)IHΦ2(IH)IHIH].

    Similarly, define

    V2=SVSVSVlogSVSV+IVIVIVlogIVIV,

    the same can be obtained

    dV2dt=mSV(1SVSV)+β3Φ3(IH)SV(1IVIV)+β3Φ3(IH)SV(SVSV+Φ3(IH)Φ3(IH))+mSV(1SVSV)β3SVΦ3(IH)IVSVΦ3(IH)IVSVΦ3(IH)+β3SVΦ3(IH)=mSV(1SVSVSVSV)+β3SVΦ3(IH)[Φ3(IH)Φ3(IH)SVSVIVSVΦ3(IH)IVSVΦ3(IH)IVIV+2]mSV(SVSV)2SV+β3SVΦ3(IH)[Φ3(IH)Φ3(IH)logIVΦ3(IH)IVΦ3(IH)IVIV].

    Next, define

    V3=V1+¯β1Φ1(IV)SH¯β3Φ3(IH)SVV2,

    we can get

    dV3dtμSH(SHSH)2SHm¯β1Φ1(IV)SH¯β3Φ3(IH)+¯β2Φ2(IH)SH[Φ2(IH)Φ2(IH)logIHΦ2(IH)IHΦ2(IH)IHIH]+¯β1Φ1(IV)SH[Φ1(IV)Φ1(IV)logIHΦ1(IV)IHΦ1(IV)IHIH+Φ3(IH)Φ3(IH)logIVΦ3(IH)IVΦ3(IH)IVIV]μSH(SHSH)2SHm¯β1Φ1(IV)SH¯β3Φ3(IH)+¯β2Φ2(IH)SH[Φ2(IH)Φ2(IH)+IHΦ2(IH)IHΦ2(IH)IHIH1]+¯β1Φ1(IV)SH[Φ1(IV)Φ1(IV)+IVΦ1(IV)IVΦ1(IV)IVIV1+IHΦ3(IH)IHΦ3(IH)+Φ3(IH)Φ3(IH)IHIH1],

    by the condition (A2) and (A3), we can know

    Φ1(IV)Φ1(IV)+IVΦ1(IV)IVΦ1(IV)IVIV1=IVΦ1(IV)Φ1(IV)[(Φ1(IV)Φ1(IV))(Φ1(IV)IVΦ1(IV)IV)]0,

    Φ2(IH) and Φ3(IH) similarly satisfy the structure of the above equation, combined with SH and SV are bounded, we finally get

    dV3dtμK(SHSH)2m2Λ¯β1Φ1(IV)SH¯β3Φ3(IH)SV(SVSV)2.

    Next, we define

    V4=(SHSH+IHIH)22, V5=(RHRH)22, V6=(SVSV+IVIV)22,

    similarly calculated

    dV4dt=μ(SHSH)2(d+γ)(IHIH)2ω(SHSH)(IHIH)μ(SHSH)2(d+γ)(IHIH)2+(d+γ)2(IHIH)2+ω22(d+γ)(SHSH)2ω22(d+γ)(SHSH)2(d+γ)2(IHIH)2,

    and

    dV5dt=γ(IHIH)(RHRH)μ(RHRH)2μ2(RHRH)2+γ22μ(IHIH)2,

    and

    dV6dt=m(SVSV)2m(IVIV)22m(SVSV)(IVIV)=m(SVSV)2m(IVIV)2+m2(IVIV)2+2m(SVSV)2m(SVSV)2m2(IVIV)2.

    Finally, define

    V=V3+A1V4+A2V5+A3V6,

    then

    dVdt(μKA1ω22(d+γ))(SHSH)2(d+γ2A1γ22μA2)(IHIH)2μ2A2(RHRH)2(m2Λ¯β1Φ1(IV)SH¯β3Φ3(IH)SVA3m)(SVSV)2m2A3(IVIV)2,

    take

    A1=μ(d+γ)Kω2, A2γ2μ=d+γ2A1, A3=m2Λ¯β1Φ1(IV)SH2m¯β3Φ3(IH)SV,

    we get

    dVdtμ2K(SHSH)2d+γ4A1(IHIH)2μ2A2(RHRH)2m2Λ¯β1Φ1(IV)SH2¯β3Φ3(IH)SV(SVSV)2m2A3(IVIV)2,

    the proof is done.

    Proof. It is obvious that the coefficients of the system is locally Lipschitz continuous, so there is a unique local solution (SH(t),IH(t),IV(t),β1(t),β2(t)) on t[0,τe), where τe represents explosion time. To show that the solution is global, according to the method in [33], we just need to verify that τe=+ a.s..

    Choose k0 be a sufficiently large integer for every component of (SH(0),IH(0),IV(0),β1(0),β2(0)) within the interval [1k0,k0]. For each integer kk0, define the stopping time as

    τk=inf{t[0,τe)|min(SH(t),IH(t),IV(t),β1(t),β2(t))1k or max(SH(t),IH(t),IV(t),β1(t),β2(t))k}.

    It can be seen that τk is monotonically increasing with respect to k. Then define inf{}=+ and τ=limt+τk. It is clearly visible that the solution is global due to the fact that τ<τe a.s., τ= leading to τe=. Next, we will prove τ= by the contradiction method. Assuming τ<+ a.s., then there are ε0(0,1) and T>0 such that P(τT)>ε, so there is a positive integer k1>k0 that makes

    P(τkT)ε, kk1.

    Define a non-negative C2-function V(SH,IH,IV,β1,β2) as follows

    V=SH1logSH+IH1logIH+IV1logIV+(KSHIH)1log(KSHIH)+(ΛmIV)1log(ΛmIV)+β11logβ1+β21logβ2.

    Applying Itˆo's formula to V, then we can get

    LV=μ(KSH)β1ϕ1(IV)SHβ2ϕ2(IH)SH+dIHμKSH+μ+β1ϕ1(IV)+β2ϕ2(IH)dIHSH+β1ϕ1(IV)SH+β2ϕ2(IH)SHωIHβ1ϕ1(IV)SHIHβ2ϕ2(IH)SHIH+ω+β3ϕ3(IH)(ΛmIV)mIVβ3ϕ3(IH)ΛmIV+β3ϕ3(IV)+m+1ΛmIV(β3ϕ3(IH)(ΛmIV)mIV)β3ϕ3(IH)(ΛmIV)+mIV+1K(SH+IH)(μ(KSH)+(dω)IH)μ(KSH)(dω)IH+β1(α1log¯β1α1logβ1+12σ21)α1(log¯β1logβ1)+β2(α2log¯β2α2logβ2+12σ22)α2(log¯β2logβ2)μK+dIH+μ+β1ϕ1(IV)+β2ϕ2(IH)+ω+Λmβ3ϕ3(IH)+m+2β3ϕ3(IH)+μ+γIH+β1(α1log¯β1α1logβ1+12σ21)α1(log¯β1logβ1)+β2(α2log¯β2α2logβ2+12σ22)α2(log¯β2logβ2).

    By the condition A3, we can know that

    ϕ1(IV)ϕ1(0)IV,ϕ2(IH)ϕ2(0)IH,ϕ3(IH)ϕ3(0)IH,

    then, we obtain

    LVμK+2μ+ω+m+β1ϕ1(0)IV+[β2ϕ2(0)+Λmβ3ϕ3(0)+2β3ϕ3(0)+γ]IH+β1(α1log¯β1α1logβ1+12σ21)α1(log¯β1logβ1)+β2(α2log¯β2α2logβ2+12σ22)α2(log¯β2logβ2)μK+2μ+ω+m+β1ϕ1(0)Λm+[β2ϕ2(0)+Λmβ3ϕ3(0)+2β3ϕ3(0)+γ]K+β1(α1log¯β1α1logβ1+12σ21)α1(log¯β1logβ1)+β2(α2log¯β2α2logβ2+12σ22)α2(log¯β2logβ2):=H(β1,β2).

    It's easy to see H(β1,β2) as β1+,β10,β2+,β20, so there's a positive constant H0 that makes LVH0. Integrating on both sides and taking the expectation, then

    0EW(SH(τkT),IH(τkT),IV(τkT),β1(τkT),β2(τkT))=EW(SH(0),IH(0),IV(0),β1(0),β2(0))+EτnT0LV(SH(τ),IH(τ),IV(τ),β1(τ),β2(τ))dτEV(SH(0),IH(0),IV(0),β1(0),β2(0))+H0T.

    One gets that for any ζGk, W(SH(τk,ζ),IH(τk,ζ),IV(τk,ζ),β1(τk,ζ),β2(τk,ζ)) will larger than (ek1k)(ek1+k), so

    EW(SH(0),IH(0),IV(0),β1(0),β2(0))+H0TEW(SH(τkT),IH(τkT),IV(τkT),β1(τkT),β2(τkT))E[IGk(ζ)W(SH(τkT),IH(τkT),IV(τkT),β1(τkT),β2(τkT))]P(Gk(ζ))W(SH(τk,ζ),IH(τk,ζ),IV(τk,ζ),β1(τk,ζ),β2(τk,ζ))ε0[(ek1k)(ek1+k)].

    Since k is an arbitrary constant, it can be contradictory by making k+

    +EV(SH(0),IH(0),IV(0),β1(0),β2(0))+H0T<+.

    Therefore τ=+ a.s., i.e., τe=+. Then system (1.6) has a unique global solution (SH(t),IH(t),IV(t),β1(t),β2(t)) on Γ.

    Proof. The theorem will be proved next in the following two steps.

    Step 1. Construct Lyapunov functions

    Using Itˆo's formula, we obtain

    L(logSH)=μKSH+μ+β1ϕ1(IV)+β2ϕ2(IH)dIHSH,
    L(logIH)=β1ϕ1(IV)SHIHβ2ϕ2(IH)SHIH+ω,
    L(logIV)=Λmβ3ϕ3(IH)1IV+β3ϕ3(IH)+m.

    Define a function V1=logIHc1logSHc2logSHc3logIV, where c1,c2,c3 are given in subsequent calculation. Using Itˆo's formula, then

    LV1=β1ϕ1(IV)SHIHβ2ϕ2(IH)SHIH+ωc1μKSH+c1μ+c1β1ϕ1(IV)+c1β2ϕ2(IH)c1dIHSHc2μKSH+c2μ+c2β1ϕ1(IV)+c2β2ϕ2(IH)c2dIHSHc3Λmβ3ϕ3(IH)1IV+c3β3ϕ3(IH)+c3mc4IVϕ1(IV)+c41ϕ1(0)+c4(IVϕ1(IV)1ϕ1(0))c5IHϕ3(IH)+c51ϕ3(0)+c5(IHϕ3(IH)1ϕ3(0))c6IHϕ2(IH)+c61ϕ2(0)+c6(IHϕ2(IH)1ϕ2(0)).

    By the condition A3, notice that

    (IVϕ1(IV))=IVϕ1(IV)1ϕ1(IV)IVm0,

    which can deduce

    IVϕ1(IV)1ϕ1(0)m0IV, (A.4)

    similarly, one gets

    IHϕ2(IH)1ϕ2(0)m0IH,IHϕ3(IH)1ϕ3(0)m0IH. (A.5)

    Combining (A.4) and (A.5), we have

    LV155c1c3c4c5β1β3μKΛm+c1μ+c3m+c41ϕ1(0)+c51ϕ3(0)33c2c6β2μK+c2μ+c61ϕ2(0)+ω+c1β1ϕ1(0)IV+c1β2ϕ2(0)IH+c3β3ϕ3(0)IH+c2β1ϕ1(0)IV+c2β2ϕ2(0)IH+c4m0IV+(c5+c6)m0IH=55c1c3c4c5~β1β3μKΛm+c1μ+c3m+c41ϕ1(0)+c51ϕ3(0)33c2c6~β2μK+c2μ+c61ϕ2(0)+ω+c1β1ϕ1(0)IV+c1β2ϕ2(0)IH+c3β3ϕ3(0)IH+c2β1ϕ1(0)IV+c2β2ϕ2(0)IH+c4m0IV+(c5+c6)m0IH+5(5c1c3c4c5~β1β3μKΛm5c1c3c4c5β1β3μKΛm)+3(3c2c6~β2μK3c2c6β2μK).

    Let c1,c2,c3,c4,c5 and c6 satisfy the following equalities

    c1μ=c3m=c41ϕ1(0)=c51ϕ3(0)=~β1β3ΛKϕ1(0)ϕ3(0)m2,
    c2μ=c61ϕ2(0)=~β2Kϕ2(0),

    then

    LV1~β1β3ΛKϕ1(0)ϕ3(0)m2~β2Kϕ2(0)+ω+(c1+c2)ϕ2(0)β2IH+c3ϕ3(0)β3IH+(c1+c2)ϕ1(0)β1IV+c4m0IV+(c5+c6)m0IH+5(5c1c3c4c5~β1β3μKΛm5c1c3c4c5β1β3μKΛm)+3(3c2c6~β2μK3c2c6β2μK)=ω(Rs01)+(c1+c2)ϕ2(0)β2IH+c3ϕ3(0)β3IH+(c1+c2)ϕ1(0)β1IV+c4m0IV+(c5+c6)m0IH+5(5c1c3c4c5~β1β3μKΛm5c1c3c4c5β1β3μKΛm)+3(3c2c6~β2μK3c2c6β2μK),

    where

    Rs0=~β1β3ΛKϕ1(0)ϕ3(0)m2ω+~β2Kϕ2(0)ω.

    By Holder inequality, for any positive constant δ, the following equations are true

    β1IV(δβ21+14δ)IVδβ21Λm+IV4δ=Λmδ¯β12eσ21α1+IV4δ+Λmδ(β21¯β12eσ21α1),β2IH(δβ22+14δ)IHδβ22K+IH4δ=Kδ¯β22eσ22α2+IH4δ+Kδ(β22¯β22eσ22α2).

    If take δ to be

    δ=ω2(Rs01)(¯β1β3μKϕ1(0)ϕ3(0)μm2+¯β2Kϕ2(0)μ)(ϕ2(0)K¯β22eσ22α2+ϕ1(0)Λm¯β12eσ21α1),

    then we can get

    LV1ω(Rs01)+(c1+c2)ϕ2(0)Kδ¯β22eσ22α2+(c1+c2)ϕ2(0)IH4δ+c3ϕ3(0)β3IH+(c1+c2)ϕ1(0)Λmδ¯β12eσ21α1+(c1+c2)ϕ1(0)IV4δ+c4m0IV+(c5+c6)m0IH+(c1+c2)ϕ2(0)Kδ(β22¯β22eσ22α2)+(c1+c2)ϕ1(0)Λmδ(β21¯β12eσ21α1)+5(5c1c3c4c5~β1β3μKΛm5c1c3c4c5β1β3μKΛm)+3(3c2c6~β2μK3c2c6β2μK):=ω2(Rs01)+[(c1+c2)ϕ2(0)14δ+c3ϕ3(0)β3+(c5+c6)m0]IH+[(c1+c2)ϕ1(0)14δ+c4m0]IV+F(β1,β2),

    where

    F(β1,β2)=(c1+c2)ϕ2(0)Kδ(β22¯β22eσ22α2)+(c1+c2)ϕ1(0)Λmδ(β21¯β12eσ21α1)+5(5c1c3c4c5~β1β3μKΛm5c1c3c4c5β1β3μKΛm)+3(3c2c6~β2μK3c2c6β2μK).

    Next we define

    V2=V1+(c1+c2)ϕ1(0)+4δc4m04mδIV,

    applying Itˆo's formula to V2, it leads to

    LV2ω2(Rs01)+[(c1+c2)ϕ2(0)14δ+c3ϕ3(0)β3+(c5+c6)m0]IH+(c1+c2)ϕ1(0)+4δc4m04δIV+F(β1,β2)+(c1+c2)ϕ1(0)+4δc4m04mδβ3ϕ3(IH)Λm(c1+c2)ϕ1(0)+4δc4m04mδβ3ϕ3(IH)IV(c1+c2)ϕ1(0)+4δc4m04δIVω2(Rs01)+AIH+F(β1,β2),

    where

    A=(c1+c2)ϕ2(0)14δ+c3ϕ3(0)β3+(c5+c6)m0+((c1+c2)ϕ1(0)+4δc4m0)β3ϕ3(0)Λ4m2δ.

    Next, define

    V3=logSHlogIVlog(K(SH+IH))log(ΛmIV)+(β11logβ1)+(β21logβ2),

    then, we have

    LV3=μKSH+μ+β1ϕ1(IV)+β2ϕ2(IH)dIHSHΛmβ3ϕ3(IH)1IV+β3ϕ3(IH)+mIHϕ3(IH)+1ϕ3(0)1K(SH+IH)[μ(K(SH+IH))+γIH]1ΛmIV[β3ϕ3(IH)(ΛmIV)+mIV]+IHϕ3(IH)1ϕ3(0)+β1(α1log¯β112α1logβ1+12σ21)α1(log¯β112logβ1)12β1α1logβ1+12α1logβ1+β2(α2log¯β212α2logβ2+12σ22)α2(log¯β212logβ2)12β2α2logβ2+12α2logβ2.μKSHΛβ3IHmIVγIHK(SH+IH)mIVΛmIV+W(β1,β2)12β1α1logβ1+12α1logβ112β2α2logβ2+12α2logβ2,

    where

    W(β1,β2)=2μ+m+1ϕ3(0)+β1ϕ1(0)Λm+β2ϕ2(0)K+2β3ϕ3(0)K+m0K+β1(α1log¯β112α1logβ1+12σ21)α1(log¯β112logβ1)+β2(α2log¯β212α2logβ2+12σ22)α2(log¯β212logβ2).

    Choose M is a large enough positive constant, let

    ˉV=MV2+V3,

    where M satisfying the following inequality

    Mω2(Rs01)+sup(β1,β2)R2+W(β1,β2)2.

    Notice that ˉV has a minimum value ˉVmin in the interior of Γ because ˉV+ as (SH,IH,IV,β1,β2) tends to the boundary of Γ. Ultimately, establish a non-negative C2-function V(SH,IH,IV,β1,β2):ΓR+ as follows

    V(SH,IH,IV,β1,β2)=ˉV(SH,IH,IV,β1,β2)ˉVmin,

    then we obtain

    LVMω2(Rs01)+MAIH+MF(β1,β2)μKSHΛβ3IHmIVγIK(SH+IH)mIVΛmIV12β1α1logβ1+12α1logβ112β2α2logβ2+12α2logβ2+W(β1,β2):=G(SH,IH,IV,β1,β2)+MF(β1,β2). (A.6)

    Step 2. Set up the closed set Uε

    Uε={(SH,IH,IV,β1,β2)Γ|IHε, SHε, IVε2, SH+IHKε2, IVΛmε3, εβ11ε, εβ21ε},

    where ε is a small enough constant and the complement of Uε can be divided into nine small sets as follows

    Uc1,ε={(SH,IH,IV,β1,β2)Γ|0<β1<ε},Uc2,ε={(SH,IH,IV,β1,β2)Γ|β1>1ε},
    Uc3,ε={(SH,IH,IV,β1,β2)Γ|0<β2<ε},Uc4,ε={(SH,IH,IV,β1,β2)Γ|β2>1ε},
    Uc5,ε={(SH,IH,IV,β1,β2)Γ|0<IH<ε},Uc6,ε={(SH,IH,IV,β1,β2)Γ|0<SH<ε},
    Uc7,ε={(SH,IH,IV,β1,β2)Γ|0<IV<ε2,IHε},
    Uc8,ε={(SH,IH,IV,β1,β2)Γ|SH+IH>Kε2,IHε},
    Uc9,ε={(SH,IH,IV,β1,β2)Γ|IV>Λmε3,IVε2},

    then the following results hold

    Case 1: (SH,IH,IV,β1,β2)Uc1,ε, then

    G(SH,IH,IV,β1,β2)=Mω2(Rs01)+MAIH+W(β1,β2)μKSHΛβ3IHmIVγIK(SH+IH)mIVΛmIV12β1α1logβ1+12α1logβ112β2α2logβ2+12α2logβ2Mω2(Rs01)+MAIH+W(β1,β2)+12α1logβ1+12α2logβ2Mω2(Rs01)+MAK+12α1logε+12α2logβ2+sup(β1,β2)R2+W(β1,β2)1.

    Case 2: (SH,IH,IV,β1,β2)Uc2,ε, then

    G(SH,IH,IV,β1,β2)Mω2(Rs01)+MAIH+W(β1,β2)+12α1logβ1+12α2logβ212β1α1logβ1Mω2(Rs01)+MAK+12α1logβ1+12α2logβ2α1log1ε2ε+sup(β1,β2)R2+W(β1,β2)1.

    Case 3: (SH,IH,IV,β1,β2)Uc3,ε, then

    G(SH,IH,IV,β1,β2)Mω2(Rs01)+MAIH+W(β1,β2)+12α1logβ1+12α2logβ2Mω2(Rs01)+MAK+12α2logε+12α1logβ1+sup(β1,β2)R2+W(β1,β2)1.

    Case 4: (SH,IH,IV,β1,β2)Uc4,ε, then

    G(SH,IH,IV,β1,β2)Mω2(Rs01)+MAIH+W(β1,β2)+12α1logβ1+12α2logβ212β2α2logβ2Mω2(Rs01)+MAKα2log1ε2ε+12α1logβ1+12α2logβ2+sup(β1,β2)R2+W(β1,β2)1.

    Case 5: (SH,IH,IV,β1,β2)Uc5,ε, then

    G(SH,IH,IV,β1,β2)Mω2(Rs01)+MAIH+W(β1,β2)+12α1logβ1+12α2logβ2Mω2(Rs01)+MAε+12α1logβ1+12α2logβ2+sup(β1,β2)R2+W(β1,β2)1.

    Case 6: (SH,IH,IV,β1,β2)Uc6,ε, then

    G(SH,IH,IV,β1,β2)Mω2(Rs01)+MAIH+W(β1,β2)μKSH+12α1logβ1+12α2logβ2Mω2(Rs01)+MAKμKε+12α1logβ1+12α2logβ2+sup(β1,β2)R2+W(β1,β2)1.

    Case 7: (SH,IH,IV,β1,β2)Uc7,ε, then

    G(SH,IH,IV,β1,β2)Mω2(Rs01)+MAIH+W(β1,β2)Λβ3IHmIV+12α1logβ1+12α2logβ2Mω2(Rs01)+MAKΛβ3mε+12α1logβ1+12α2logβ2+sup(β1,β2)R2+W(β1,β2)1.

    Case 8: (SH,IH,IV,β1,β2)Uc8,ε, then

    G(SH,IH,IV,β1,β2)Mω2(Rs01)+MAIH+W(β1,β2)γIHK(SH+IH)+12α1logβ1+12α2logβ2Mω2(Rs01)+MAKγε+12α1logβ1+12α2logβ2+sup(β1,β2)R2+W(β1,β2)1.

    Case 9: (SH,IH,IV,β1,β2)Uc9,ε, then

    G(SH,IH,IV,β1,β2)Mω2(Rs01)+MAIH+W(β1,β2)mIVΛmIV+12α1logβ1+12α2logβ2Mω2(Rs01)+MAKmε+12α1logβ1+12α2logβ2+sup(β1,β2)R2+W(β1,β2)1.

    According to the discussion of cases above, we can know that

    G(SH,IH,IV,β1,β2)1,(SH,IH,IV,β1,β2)ΓUε,

    in other words, let H is a positive constant that makes

    G(SH,IH,IV,β1,β2)H<+,(SH,IH,IV,β1,β2)Γ.

    For any initial value (SH(0),IH(0),IV(0),β1(0),β2(0))Γ, integrating the inequality (A.6) and taking the expectation, we get

    0E[V(SH(t),IH(t),IV(t),β1(t),β2(t))]t=E[V(SH(0),IH(0),IV(0),β1(0),β2(0))]t+1tt0E(LV(SH(τ),IH(τ),IV(τ),β1(τ),β2(τ)))dτEV(SH(0),IH(0),IV(0),β1(0),β2(0))t+1tt0E(G(SH(τ),IH(τ),IV(τ),β1(τ),β2(τ)))dτ+5M5c1c3c4c5β3μKΛm1tt0E(5~β15β1(τ))dτ+3M3c2c6μK1tt0E(3~β23β2(τ))dτ+(c1+c2)ϕ1(0)δΛm1tt0E(β21(τ)¯β12eσ21α1)dτ+(c1+c2)ϕ2(0)δK1tt0E(β22(τ)¯β22eσ22α2)dτ. (A.7)

    One gets that βi(i=1,2) is ergodic according to [34,35], then we can get that

    limt+1tt0βpi(τ)dτ=limt+1tt0eplogβi(τ)dτ=+epyiπ(yi)dyi=¯βipep2σ2i4αi,

    hence

    limt+1tt0β151(τ)dτ=¯β115eσ21100α1=~β115,limt+1tt0β132(τ)dτ=¯β213eσ2136α2=~β213,limt+1tt0β21(τ)dτ=¯β12eσ21α1,limt+1tt0β22(τ)dτ=¯β22eσ22α2.

    Then letting t+ and taking infimum to (A.7) it follows

    0lim inft+1tt0E(G(SH(τ),IH(τ),IV(τ),β1(τ),β2(τ)))dτ=lim inft+1tt0E(G(SH(τ),IH(τ),IV(τ),β1(τ),,β2(τ))I{(SH(τ),IH(τ),IV(τ),β1(τ),β2(τ)Uε})dτ+lim inft+1tt0E(G(SH(τ),IH(τ),IV(τ),β1(τ),β2(τ))I{(SH(τ),IH(τ),IV(τ),β1(τ),β2(τ)ΓUε})dτHlim inft+1tt0I{(SH(τ),IH(τ),IV(τ),β1(τ),β2(τ)Uε}dτlim inft+1tt0I{(SH(τ),IH(τ),IV(τ),β1(τ),β2(τ)ΓUε}dτ1+(H+1)lim inft+1tt0I{(SH(τ),IH(τ),IV(τ),β1(τ),β2(τ)Uε}dτ,

    which means

    lim inft+1tt0P{(SH(τ),IH(τ),IV(τ),β1(τ),β2(τ))Uε}dτ1H+1>0 a.s.,
    lim inft+1tt0P{τ,(SH(0),IH(0),IV(0),β1(0),β2(0)),Uε}dτ1H+1, (SH(0),IH(0),IV(0),β1(0),β2(0))Γ.

    According to the Lemma 2.4, we can conclude that when Rs0>1 system (1.6) has a stationary distribution on Γ.

    Proof. Define a C2-function G(IH,IV,β1,β2):ΓR by

    G(IH,IV,β1,β2)=v1IH+v2IV,

    where v1=R0,v2=¯β1ϕ1(0)Km. Applying Itˆo's formula to G(IH,IV,β1,β2), then we have

    L(logG)=1v1IH+v2IV[v1(β1ϕ1(IV)SH+β2ϕ2(IH)SHωIH)+v2(β3ϕ3(IH)(ΛmIV)mIV)]1v1IH+v2IV[v1β1ϕ1(0)KIV+v1β2ϕ2(0)KIHv1ωIH+v2β3ϕ3(0)(ΛmIV)IHv2mIV]=1v1IH+v2IV[v1¯β1ϕ1(0)KIV+v1¯β2ϕ2(0)KIHv1ωIH+v2β3ϕ3(0)(ΛmIV)IHv2mIV]+1v1IH+v2IV[v1(β1¯β1)ϕ1(0)KIV+v1(β2¯β2)ϕ2(0)KIH]1v1IH+v2IV[(v1¯β1ϕ1(0)Kv2m)IV+(v1¯β2ϕ2(0)K+v2β3ϕ3(0)Λmv1ω)IH]+mR0¯β1|β1¯β1|+ϕ2(0)K|β2¯β2|1v1IH+v2IV[¯β1ϕ1(0)K(R01)IV+¯β2ϕ2(0)K(R01)IH]+mR0¯β1|β1¯β1|+ϕ2(0)K|β2¯β2|min{m,¯β2ϕ2(0)KR0}(R01)+mR0¯β1|β1¯β1|+ϕ2(0)K|β2¯β2|.

    Integrating both sides of this equation from 0 to t and dividing by t, we get

    logG(t)logG(0)tmin{m,¯β2ϕ2(0)KR0}(R01)+mR0¯β11tt0|β1(τ)¯β1|dτ+ϕ2(0)K1tt0|β2(τ)¯β2|dτ. (A.8)

    According to the ergodicity of β1,β2, then

    limt1tt0|β1(τ)¯β1|dτlimt(1tt0(β1(τ)¯β1)2dτ)12=¯β1(eσ21α12eσ214α1+1)12,limt1tt0|β2(τ)¯β2|dτlimt(1tt0(β2(τ)¯β2)2dτ)12=¯β2(eσ22α22eσ224α2+1)12. (A.9)

    Letting t+, and submitting (A.9) into (A.8), then inequailty (A.8) becomes

    lim supt+logG(t)tmin{m,¯β2ϕ2(0)KR0}(R01)+mR0(eσ21α12eσ214α1+1)12+ϕ2(0)K¯β2(eσ22α22eσ224α2+1)12:=min{m,¯β2ϕ2(0)KR0}(RE01),

    where

    RE0=R0+mR0(eσ21α12eσ214α1+1)12+ϕ2(0)K¯β2(eσ22α22eσ224α2+1)12min{m,¯β2ϕ2(0)KR0}.

    If RE0<1,

    lim supt+logG(t)t<0

    will be true which indicates

    limt+IH(t)=0limt+IV(t)=0,

    this means the disease will die out exponentially.

    Proof. Step 1 Consider the following equation

    G21+AΣ1+Σ1AT=0, (A.10)

    where G1=diag(0,0,0,σ1,0).

    Let A1=J1AJ11, where

    J1=(0001010000010000010000001),

    then

    A1=(a440000a14a11a12a13a15a14a21a22a13a1500a32a3300000a55).

    Let A2=J2A1J12, where

    J2=(1000001000011000001000001),

    and

    A2=(a440000a14a12a11a12a13a150a12a11+a21+a22a12a22000a32a32a3300000a55).

    Due to a12a11+a21+a22=γ>0, let A3=J3A2J13, where

    J3=(10000010000010000a32γ1000001),

    and

    A3=(a440000a14a12a11a13a32γa12a13a150γa12a220000wa3300000a55).

    in which

    w=a32a32(a12+a22)γ+a32a33γ=mμ+β3ϕ3(I)>0.

    By using the methodology in [36,37], the standard transformation matrix of A3 has the following form

    M=(m1m2m3m4m50wγw(a12+a22+a33)a233000wa3300001000001),

    where m1=wγa14, m2=wγ(a33+a11+a22), m3=wa13a32wγa12+w(a12+a22+a33)(a12+a22)+wa233, m4=γwa13a333,m5=γwa15.

    Define A01=MA3M1, then we can get

    A01=(b1b2b3b4b51000001000001000000a55),

    in which

    b1=a11+a22+a33+a44,b2=a44(a11+a22+a33)+a11a22+a12a21+a11a33a13a32+a22a33,b3=a44(a11a22+a12a21+a11a33a13a32+a22a33)a11a13a32+a11a22a33+a12a21a33+a13a21a32,b4=a44(a11a22a33a11a13a32+a12a21a33+a13a21a32).

    Let J=J3J2J1, we can equivalently transform the Eq (A.10) into

    (MJ)G21(MJ)T+[(MJ)A(MJ)1][(MJ)Σ1(MJ)T]+[(MJ)Σ1(MJ)T][(MJ)A(MJ)1]T=0, (A.11)

    where (MJ)G21(MJ)T=diag((m1σ1)2,0,0,0,0), let ρ1=m1σ1, then (A.11) becomes

    G20+ρ21A01[(MJ)Σ1(MJ)T]+ρ21[(MJ)Σ1(MJ)T]AT01=0,

    then we obtian

    Σ01:=ρ21(MJ)Σ1(MJ)T=(b1b4b2b3b0b3b000b3b0b1b0b3b0b1b000b1b0b3b1b2b4b000000),

    where b=2[b4b21b1b2b3+b23]. We can obtain that the matrix Σ01 is a positive semi-definite matrix, the exact expression of Σ1 is as follows

    Σ1=ρ21(MJ)1Σ01[(MJ)1]T.

    Step 2 Consider the following equation

    G22+AΣ2+Σ2AT=0, (A.12)

    where G2=diag(0,0,0,0,σ2).

    Let B1=P1AP11, where

    P1=(0000110000010000010000010),

    then

    B1=(a550000a15a11a12a13a14a15a21a22a13a1400a32a3300000a44).

    Let B2=P2B1P12, where

    P2=(1000001000011000001000001),

    then

    B2=(a550000a15a12a11a12a13a140a12a11+a21+a22a12a22000a32a32a3300000a44).

    Similarly, due to a12a11+a21+a22=γ>0, let B3=P3B2P13, where

    P3=(10000010000010000a32γ1000001),

    and

    B3=(a550000a15a12a11a13a32γa12a13a140γa12a220000wa3300000a44),

    in which

    w=a32a32(a12+a22)γ+a32a33γ=mμ+β3ϕ3(I)>0.

    The standard transformation matrix of B3 has the following form

    N=(n1n2n3n400wγw(a33+a12+a22)a233000wa3300001000001),

    where n1=wγa15, n2=wγ(a33+a11+a22), n3=wa13a32wγa12+w(a12+a22+a33)(a12+a22)+wa233, n4=γwa13a333,n5=γwa14.

    Define B01=NB3N1, then we can get

    B01=(d1d2d3d4d51000001000001000000a44),

    in which

    d1=a11+a22+a33+a55,d2=a55(a11+a22+a33)+a11a22+a12a21+a11a33a13a32+a22a33,d3=a55(a11a22+a12a21+a11a33a13a32+a22a33)a11a13a32+a11a22a33+a12a21a33+a13a21a32,d4=a55(a11a22a33a11a13a32+a12a21a33+a13a21a32).

    Let P=P3P2P1, the equation (A.12) can be equivalently transformed into

    (NP)G22(NP)T+[(NP)A(NP)1][(NP)Σ2(NP)T]+[(NP)Σ2(NP)T][(NP)A(NP)1]T=0, (A.13)

    where (NP)G22(NP)T=diag((n1σ2)2,0,0,0,0), let ρ2=n1σ2, then (A.13) becomes

    G20+ρ22B01[(NP)Σ2(NP)T]+ρ23[(NP)Σ2(NP)T]BT01=0,

    then we obtian

    Σ02:=ρ22(NP)Σ2(NP)T=(d1d4d2d3d0d3d000d3d0d1d0d3d0d1d000d1d0d3d1d2d3d000000),

    where d=2(d4d21d1d2d3+d23), and we can obtain that the matrix Σ02 is a positive semi-definite matrix, the exact expression of Σ2 is as follows

    Σ2=ρ22(NP)1Σ02[(NP)1]T.

    Finally, Σ=Σ1+Σ2. Obviously, the matrix Σ is a positive definite matrix. The proof is complete.



    [1] H. Lee, S. Halverson, N. Ezinwa, Mosquito-borne diseases, Primary Care: Clin. Off. Pract., 45 (2018), 393–407. https://doi.org/10.1016/j.pop.2018.05.001
    [2] M. A. Tolle, Mosquito-borne diseases, Curr. Probl. Pediatr. Adolesc. Health Care, 39 (2009), 97–140. https://doi.org/10.1016/j.cppeds.2009.01.001 doi: 10.1016/j.cppeds.2009.01.001
    [3] J. Oliveira-Ferreira, M. V. G. Lacerda, P. Brasil, J. L. Ladislau, P. L. Tauil, C. T. Daniel-Ribeiro, Malaria in Brazil: an overview, Malar. J., 9 (2010), 1–15. https://doi.org/10.1186/1475-2875-9-115 doi: 10.1186/1475-2875-9-115
    [4] V. Wiwanitkit, Dengue fever: diagnosis and treatment, Expert Rev. Anti-Infect. Ther., 8 (2010), 841–845. https://doi.org/10.1586/eri.10.53 doi: 10.1586/eri.10.53
    [5] E. D. Barnett, Yellow fever: epidemiology and prevention, Clin. Infect. Dis., 44 (2007), 850–856. https://doi.org/10.1086/511869 doi: 10.1086/511869
    [6] E. B. Hayes, J. J. Sejvar, S. R. Zaki, R. S. Lanciotti, A. V. Bode, Virology, pathology, and clinical manifestations of West Nile virus disease, Emerging Infect. Dis., 11 (2005), 1174. https://doi.org/10.3201/eid1108.050289b doi: 10.3201/eid1108.050289b
    [7] L. Esteva, H. M. Yang, Mathematical model to assess the control of Aedes aegypti mosquitoes by the sterile insect technique, Math. Biosci., 198 (2005), 132–147. https://doi.org/10.1016/j.mbs.2005.06.004 doi: 10.1016/j.mbs.2005.06.004
    [8] E. A. Newton, P. Reiter, A model of the transmission of dengue fever with an evaluation of the impact of ultra-low volume (ULV) insecticide applications on dengue epidemics, Am. J. Trop. Med. Hyg., 47 (1992), 709–720. https://doi.org/10.4269/ajtmh.1992.47.709 doi: 10.4269/ajtmh.1992.47.709
    [9] H. R. Pandey, G. R. Phaijoo, D. B. Gurung, Analysis of dengue infection transmission dynamics in Nepal using fractional order mathematicalmodeling, Chaos, Solitons Fractals: X, 11 (2023), 100098. https://doi.org/10.1016/j.csfx.2023.100098 doi: 10.1016/j.csfx.2023.100098
    [10] R. Shi, H. Zhao, S. Tang, Global dynamic analysis of a vector-borne plant disease model, Adv. Differ. Equations, 59 (2014), 1–16. https://doi.org/10.1186/1687-1847-2014-59
    [11] N. S. Chong, J. M. Tchuenche, R. J. Smith, A mathematical model of avian influenza with half-saturated incidence, Theory Biosci., 133 (2014), 23–38. https://doi.org/10.1007/s12064-013-0183-6 doi: 10.1007/s12064-013-0183-6
    [12] Y. Li, F. Haq, K. Shah, G. Rahman, M. Shahzad, Numerical analysis of fractional order pine wilt disease model with bilinear incident rate, J. Math. Comput. Sci., 17 (2017), 420–428. https://doi.org/10.22436/jmcs.017.03.07 doi: 10.22436/jmcs.017.03.07
    [13] X. Yu, S. Yuan, T. Zhang, Persistence and ergodicity of a stochastic single species model with Allee effect under regime switching, Commun. Nonlinear Sci. Numer. Simul., 59 (2018), 359–374. https://doi.org/10.1016/j.cnsns.2017.11.028 doi: 10.1016/j.cnsns.2017.11.028
    [14] X. Lv, X. Meng, X. Wang, Extinction and stationary distribution of an impulsive stochastic chemostat model with nonlinear perturbation, Chaos, Solitons Fractals, 110 (2018), 273–279. https://doi.org/10.1016/j.chaos.2018.03.038 doi: 10.1016/j.chaos.2018.03.038
    [15] T. Su, Q. Yang, X. Zhang, D. Jiang, Stationary distribution, extinction and probability density function of a stochastic SEIV epidemic model with general incidence and Ornstein-Uhlenbeck process, Physica A, 615 (2023), 128605. https://doi.org/10.1016/j.physa.2023.128605 doi: 10.1016/j.physa.2023.128605
    [16] B. Zhou, D. Jiang, B. Han, T. Hayat, Threshold dynamics and density function of a stochastic epidemic model with media coverage and mean-reverting Ornstein-Uhlenbeck process, Math. Comput. Simul., 196 (2022), 15–44. https://doi.org/10.1016/j.matcom.2022.01.014 doi: 10.1016/j.matcom.2022.01.014
    [17] X. Mu, D. Jiang, T. Hayat, A. Alsaedi, Y. Liao, A stochastic turbidostat model with Ornstein-Uhlenbeck process: dynamics analysis and numerical simulations, Nonlinear Dyn., 107 (2022), 2805–2817. https://doi.org/10.1007/s11071-021-07093-9 doi: 10.1007/s11071-021-07093-9
    [18] E. Allen, Environmental variability and mean-reverting processes, Discrete Contin. Dyn. Syst. - Ser. B, 21 (2016), 2073–2089. https://doi.org/10.3934/dcdsb.2016037 doi: 10.3934/dcdsb.2016037
    [19] X. Mao, Stochastic Differential Equations and Applications, 2nd edition, Woodhead Publishing, 2011.
    [20] Z. Ma, Y. Zhou, C. Li, Qualitative and Stability Methods for Ordinary Differential Equations, Science Press, Beijing, 2001.
    [21] B. Zhou, D. Jiang, Y. Dai, T. Hayat, Stationary distribution and density function expression for a stochastic SIQRS epidemic model with temporary immunity, Nonlinear Dyn., 105 (2021), 931–955. https://doi.org/10.1007/s11071-020-06151-y doi: 10.1007/s11071-020-06151-y
    [22] S. P. Meyn, R. L. Tweedie, Stability of markovian processes III: foster-lyapunov criteria for continuous-time processes, Adv. Appl. Probab., 25 (1993), 518–548. https://doi.org/10.1017/s0001867800025532 doi: 10.1017/s0001867800025532
    [23] N. T. Dieu, Asymptotic properties of a stochastic SIR epidemic model with Beddington-DeAngelis incidence rate, J. Dyn. Partial Differ. Equations, 30 (2018), 93–106. https://doi.org/10.1007/s10884-016-9532-8 doi: 10.1007/s10884-016-9532-8
    [24] J. S. Muldowney, Compound matrices and ordinary differential equations, Rocky Mt. J. Math., 20 (1990), 857–872. https://doi.org/10.1216/rmjm/1181073047 doi: 10.1216/rmjm/1181073047
    [25] Z. Shi, D. Jiang, Dynamical behaviors of a stochastic HTLV-I infection model with general infection form and Ornstein-Uhlenbeck process, Chaos, Solitons Fractals, 165 (2022), 112789. https://doi.org/10.1016/j.chaos.2022.112789 doi: 10.1016/j.chaos.2022.112789
    [26] C. W. Gardiner, Handbook of Stochastic Methods for Physics, Chemistry and the Natural Sciences, Springer Berlin, 1985. https://doi.org/10.1007/978-3-662-02452-2
    [27] D. J. 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
    [28] M. Aguiar, V. Anam, K. B. Blyuss, C. D. S. Estadilla, B. V. Guerrero, D. Knopoff, et al., Mathematical models for dengue fever epidemiology: A 10-year systematic review, Phys. Life Rev., 40 (2022), 65–92. https://doi.org/10.1016/j.plrev.2022.02.001 doi: 10.1016/j.plrev.2022.02.001
    [29] J. K. K. Asamoah, E. Yankson, E. Okyere, G. Sun, Z. Jin, R. Jan, et al., Optimal control and cost-effectiveness analysis for dengue fever model with asymptomatic and partial immune individuals, Results Phys., 31 (2021), 104919. https://doi.org/10.1016/j.rinp.2021.104919 doi: 10.1016/j.rinp.2021.104919
    [30] J. Duan, An Introduction to Stochastic Dynamics, Cambridge University Press, 2015.
    [31] A. Yang, H. Wang, T. Zhang, S. Yuan, Stochastic switches of eutrophication and oligotrophication: Modeling extreme weather via non-Gaussian Lévy noise, Chaos, 32 (2022), 043116. https://doi.org/10.1063/5.0085560 doi: 10.1063/5.0085560
    [32] S. Khajanchi, S. Bera, T. K. Roy, Mathematical analysis of the global dynamics of a HTLV-I infection model, considering the role of cytotoxic T-lymphocytes, Math. Comput. Simul., 180 (2021), 354–378. https://doi.org/10.1016/j.matcom.2020.09.009 doi: 10.1016/j.matcom.2020.09.009
    [33] 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
    [34] X. Mu, D. Jiang, A. Alsaedi, Analysis of a stochastic phytoplankton–zooplankton model under non-degenerate and degenerate diffusions, J. Nonlinear Sci., 32 (2022), 35. https://doi.org/10.1007/s00332-022-09787-9 doi: 10.1007/s00332-022-09787-9
    [35] Y. Cai, J. Jiao, Z. Gui, Y. Liu, W. Wang, Environmental variability in a stochastic epidemic model, Appl. Math. Comput., 329 (2018), 210–226. https://doi.org/10.1016/j.amc.2018.02.009 doi: 10.1016/j.amc.2018.02.009
    [36] Q. Yang, X. Zhang, D. Jiang, Dynamical behaviors of a stochastic food chain system with Ornstein-Uhlenbeck process, J. Nonlinear Sci., 32 (2022), 34. https://doi.org/10.1007/s00332-022-09796-8 doi: 10.1007/s00332-022-09796-8
    [37] B. Zhou, D. Jiang, Y. Dai, T. Hayat, Threshold dynamics and probability density function of a stochastic avian influenza epidemic model with nonlinear incidence rate and psychological effect, J. Nonlinear Sci., 33 (2023), 1–52. https://doi.org/10.1007/s00332-022-09885-8 doi: 10.1007/s00332-022-09885-8
  • Reader Comments
  • © 2024 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(1150) PDF downloads(52) Cited by(0)

Figures and Tables

Figures(12)  /  Tables(1)

Other Articles By Authors

/

DownLoad:  Full-Size Img  PowerPoint
Return
Return

Catalog