Processing math: 100%

A note for the global stability of a delay differential equation of hepatitis B virus infection

  • The global stability for a delayed HIV-1 infection model is investigated. It is shown that the global dynamics of the system can be completely determined by the reproduction number, and the chronic infected equilibrium of the system is globally asymptotically stable whenever it exists. This improves the related results presented in [S. A. Gourley,Y. Kuang and J.D.Nagy, Dynamics of a delay differential equation model of hepatitis B virus infection, Journal of Biological Dynamics, 2(2008), 140-153].

    Citation: Bao-Zhu Guo, Li-Ming Cai. A note for the global stability of a delay differential equation of hepatitis B virus infection[J]. Mathematical Biosciences and Engineering, 2011, 8(3): 689-694. doi: 10.3934/mbe.2011.8.689

    Related Papers:

    [1] 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
    [2] 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
    [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] 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
    [5] 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
    [6] 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
    [7] 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
    [8] 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
    [9] Jia Chen, Yuming Chen, Qimin Zhang . Ergodic stationary distribution and extinction of stochastic pertussis model with immune and Markov switching. Electronic Research Archive, 2025, 33(5): 2736-2761. doi: 10.3934/era.2025121
    [10] 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
  • The global stability for a delayed HIV-1 infection model is investigated. It is shown that the global dynamics of the system can be completely determined by the reproduction number, and the chronic infected equilibrium of the system is globally asymptotically stable whenever it exists. This improves the related results presented in [S. A. Gourley,Y. Kuang and J.D.Nagy, Dynamics of a delay differential equation model of hepatitis B virus infection, Journal of Biological Dynamics, 2(2008), 140-153].


    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, $ S_V $ and $ I_V $, and $ N = S_V+I_V $, while human host population is divided into three categories, $ S_H, I_H $ and $ R_H $. According to [10], it is reasonable to assume that the total number of human $ K = S_H+I_H+R_H $ is a positive constant.

    Table 1.  Variables and parameters in model (1.1).
    Variables and Parameters Description
    $ S_H $ number of the susceptible human host
    $ I_H $ number of the infected human host
    $ R_H $ number of the recovered human host
    K sum of the total human host
    $ S_V $ density of the susceptible mosquito vectors
    $ I_V $ density of the infected mosquito vectors
    N sum of the total mosquito vectors density
    $ \bar{\beta_1} $ biting rate of an infected vector on the susceptible human
    $ \bar{\beta_2} $ infection incidence between infected and susceptible hosts
    $ \beta_3 $ infection ratio between infected hosts and susceptible vectors
    $ \eta_1 $ determines the level at which the force of infection saturates
    $ \eta_2 $ determines the level at which the force of infection saturates
    $ \eta_3 $ determines the level at which the force of infection saturates
    $ \gamma $ the conversion rate of infected hosts to recovered hosts
    $ \mu $ natural death rate of human
    $ \Lambda $ 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 $ \frac{S_bI_b}{H_b+I_b} $, $ \frac{S_HI_a}{H_a+I_a} $ and $ \frac{S_hI_m}{H_m+I_m} $. In addition, Li et al. [12] also carried out numerical analysis of friteral order pine wilt disease model with bilinear incidence rate $ S_hI_v $ and $ I_hS_v $. 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 $ \phi_1(I_V) $, $ \phi_2(I_H) $ and $ \phi_3(I_H) $, 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

    $ (\mathcal{A}_1)\ \phi_1(0) = \phi_2(0) = \phi_3(0) = 0 $,

    $ (\mathcal{A}_2)\ \phi_1^{'}(I_V) \geq 0, \phi_2^{'}(I_H) \geq 0, \phi_3^{'}(I_H) \geq 0, \forall \ I_V, I_H \geq 0 $,

    $ (\mathcal{A}_3)\ 0\leq(\frac{I_V}{\phi_1(I_V)})^{'} \leq m_0, 0\leq(\frac{I_H}{\phi_2(I_H)})^{'} \leq m_0, 0\leq(\frac{I_H}{\phi_3(I_H)})^{'} \leq m_0 $, where $ m_0 $ is a positive constant.

    By looking at the model (1.3), the following equation is valid, $ \frac{dN}{dt} = \Lambda-mN $, that indicates, $ N(t) \rightarrow \frac{\Lambda}{m} \ as \ t \rightarrow \infty $. Note that $ S_H+I_H+R_H = K $, $ S_V+I_V = \frac{\Lambda}{m} $, this means that $ R_H = K-S_H-I_H $, $ S_V = \frac{\Lambda}{m}-I_V $, let $ \omega = d+\mu+\gamma $, 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 $ \beta_i(i = 1, 2) $ take the following form

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

    where $ \alpha_i $ represent the speed of reversal and $ \sigma_i $ represent the intensity of fluctuation. Solve the Eq (1.5), we can get

    $ \beta_i(t) = e^{-\alpha_it}\beta_i(0)+\bar{\beta_i}(1-e^{-\alpha_it})+\sigma_i\int_{0}^{t}e^{-\alpha_i(t-s)}dB_i(s), $

    where $ \beta_i(0) $ is the initial value of $ \beta_i(t) $. For arbitrary initial value $ \beta_i(0) $, $ \beta_i(t) $ follows a Gaussian distribution $ \beta_i(t) \sim \mathcal{N}(\bar{\beta_i}, \frac{\sigma_i^{2}}{2\alpha_i}) \; (t\rightarrow \infty) $. Furthermore, setting $ \beta_i(0) = \bar{\beta_i} $, then the average value of $ \beta_i(t) $ satisfies

    $ \bar{\beta_i}(t) = \frac{1}{t}\int_{0}^{t}\beta_i(s)ds = \bar{\beta_i}+\frac{1}{t}\int_{0}^{t}\frac{\sigma_i}{\alpha_i}(1-e^{\alpha_i(s-t)})dB_i(s), $

    it is known that the mathematical expectation and variance of $ \beta_i(t) $ are $ \bar{\beta_i} $ and $ \frac{\sigma_i^{2}t}{3}+\mathcal{O}(t^{2}) $, respectively, where $ \mathcal{O}(t^{2}) $ is the higher order infinitesimal of $ t^2 $. Obviously, the variance becomes zero instead of infinity as $ t\rightarrow 0 $. 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 $ \beta_1 $ and $ \beta_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 $ \phi_1(I_V) $, $ \phi_2(I_H) $ and $ \phi_3(I_H) $ 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 $ \beta_1 $ and $ \beta_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 $ \beta_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 $ R_+^n = \left\{(y_1, y_2, ...y_n) \in R_n | y_j > 0, 1 \leq j \leq n\right\} $. $ I_n $ represents the $ n $-dimensional unit matrix. $ I_A $ denotes the indicator function of set $ A $, and it means that when $ x \in A, I_A = 1 $, otherwise, $ I_A = 0 $. If $ A $ is a matrix or vector, then $ A^T $ stands for its inverse matrix, and $ A^{-1} $ stands for its inverse matrix.

    Lemma 2.1. (It$ \hat{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) = (B_1(t), B_2(t), ..., B_n(t)) $ and it represents n-dimensional Brownian motion defined on a complete probability space, let $ \mathcal{L} $ act on a function $ V \in C^{2, 1}(R_n \times R^+; R) $, then we have

    $ dV(x(t),t) = \mathcal{L}V(x(t),t)dt+V_x(x(t),t)g(t)dB(t),a.s., $

    where

    $ \mathcal{L}V(x(t),t) = V_t(x(t),t)+V_x(x(t),t)f(t)+\frac{1}{2}trace(g^T(t)V_{xx}(x(t),t)g(t)), $

    it represents the differential operator, and

    $ V_t = \frac{\partial V}{\partial t}, V_x = (\frac{\partial V}{\partial x_1},...,\frac{\partial V}{\partial x_n}),V_{xx} = (\frac{\partial^2 V}{\partial x_ix_j})_{n \times n}. $

    Lemma 2.2. (Ma et al. [20]) Letting $ \phi(\lambda) = \lambda^{n}+a_{1} \lambda^{n-1}+a_{2} \lambda^{n-2}+\cdots+a_{n-1} \lambda+a_{n} $ 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

    $ H_{k} = \left|a1a3a5a2k11a2a4a2k20a1a3a2k301a2a2k4000ak
    \right| > 0, $

    $ k = 1, 2, \cdots, n $, among them $ j > n $, replenishing definition $ a_{j} = 0 $.

    Lemma 2.3. ([21]) For five-dimension algebraic equation $ G_0^2+L\Theta+\Theta L^T = 0 $, where $ \Theta $ is a symmetric matrix, $ G_0 = diag(1, 0, 0, 0, 0) $ and

    $ L = \left(l1l2l3l4l51000001000001000000l6
    \right). $

    If $ l_1 > 0, l_2 > 0, l_3 > 0, l_4 > 0 \ and\ l_1l_2l_3-l_3^2-l_1^2l_4 > 0 $, then the symmetric matrix $ \Theta $ is a positive semi-definite matrix. Thus, we have

    $ \Theta = \left(l1l4l2l3l0l3l000l3l0l1l0l3l0l1l000l1l0l3l1l2ll4000000
    \right), $

    where $ l = 2(l_4l_1^2-l_1l_2l_3+l_3^2) $.

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

    $ \liminf\limits_{t \rightarrow +\infty} \frac{1}{t} \int_{0}^{t} P(\tau, X(0), U)d\tau > 0,a.s., $

    in which $ P(\tau, X(0), U) $ represents the transition probability of $ X(t) $, then $ X(t) $ has an invariant probability measure on $ R^n $, 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 = \left(¯β2ϕ2(0)K¯β1ϕ1(0)K00
    \right), V = \left(ω0β3Λϕ3(0)mm
    \right), $

    therefore, the next generation matrix is

    $ FV^{-1} = \left(¯β2Kϕ2(0)ω+¯β1β3ΛKϕ1(0)ϕ3(0)m2ω¯β1ϕ1(0)Km00
    \right), $

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

    $ R_0 = \rho(FV^{-1}) = \frac{\bar{\beta_2}K\phi_2^{'}(0)}{\omega}+\frac{\bar{\beta_1} \beta_3\Lambda K\phi_1^{'}(0)\phi_3^{'}(0)}{m^2\omega}. $

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

    Theorem 3.1. The disease-free equilibrium $ E_0 $ of model (1.3) is $ E_0 = (S_{H0}, \ 0, \ 0, \ S_{V0}, \ 0) = (K, \ 0, \ 0, \ \frac{\Lambda}{m}, \ 0) $ which always exists. If $ R_0>1 $, there is a unique local equilibrium $ E^{*} = (S_H^{*}, \ I_H^{*}, \ R_H^{*}, \ S_V^{*}, \ I_V^{*}) $.

    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 $ R_0 < 1 $, the disease-free equilibrium point $ E_0 $ is globally asymptotically stable. If $ R_0 > 1 $, $ E_0 $ is unstable. (ii) If $ R_0 > 1 $, the endemic equilibrium point $ E^* $ is globally asymptotically stable.

    Remark 3.2 The global stability of $ E_0 $ 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

    $ \Gamma = \left\{(S_H, I_H, I_V, \beta_1, \beta_2) \in R_{+}^{5}: S_H+I_H < K , I_V < \frac{\Lambda}{m} \right\} $

    is the positive invariant set, and the subsequent research will be discussed in $ \Gamma $.

    Theorem 4.1. For any initial value $ \left(S_H(0), I_H(0), I_V(0), \beta_1(0), \beta_2(0) \right) \in \Gamma $, there exists a unique solution $ (S_H(t), I_H(t), I_V(t), \beta_1(t), \beta_2(t)) $ of system (1.6) and the solution will remain in $ \Gamma $ 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

    $ R_{0}^{s} = \frac{\tilde{\beta_2}K\phi_2^{'}(0)}{\omega}+\frac{\tilde{\beta_1} \beta_3\Lambda K\phi_1^{'}(0)\phi_3^{'}(0)}{m^2\omega}, $

    where

    $ \tilde{\beta_1} = \bar{\beta_1}e^{\frac{\sigma_1^2}{20\alpha_1}}, \\ \tilde{\beta_2} = \bar{\beta_2}e^{\frac{\sigma_2^2}{12\alpha_2}}. $

    Theorem 4.2. If $ R_{0}^{s} > 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 $ R_0^s > 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

    $ R_{0}^{E} = R_{0}+\frac{mR_0(e^\frac{\sigma_1^2}{\alpha_1}-2e^\frac{\sigma_1^2}{4\alpha_1}+1)^\frac{1}{2}+\phi_2^{'}(0)K\bar{\beta_2}(e^\frac{\sigma_2^2}{\alpha_2}-2e^\frac{\sigma_2^2}{4\alpha_2}+1)^\frac{1}{2}}{\min\{m,\frac{\bar{\beta_2}\phi_2^{'}(0)K}{R_0}\}}. $

    Theorem 4.3. If $ R_{0}^{E} < 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 $ I_H, \ I_V $. From the expressions of $ R_0^s $ and $ R_0^E $, the relationships $ R_0^s\geq R_0 $ and $ R_0^E\geq R_0 $ are deduced, and the equal sign holds if and only if $ \sigma_1 = \sigma_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^* = (S_H^*, I_H^*, I_V^*, \log\beta_1^*, \log\beta_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 $ (z_{1}, z_{2}, z_{3}, x_{1}, x_2)^{T} = (S_H-S_H^{*}, I_H-I_H^{*}, I_V-I_V^{*}, \log\beta_1-\log\beta_1^*, \log\beta_2-\log\beta_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 $ a_{11} = \mu+\bar{\beta_1}\phi_1(I_V^*)+\bar{\beta_2}\phi_2(I_H^*), \ a_{12} = \bar{\beta_2}\phi_2^{'}(I_H^*)S_H^*-d, \ a_{13} = \bar{\beta_1}\phi_1^{'}(I_V^*)S_H^*, \ a_{14} = \bar{\beta_1}\phi_1(I_V^*)S_H^*, $ $ a_{15} = \bar{\beta_2}\phi_2(I_H^*)S_H^*, \ a_{21} = \bar{\beta_1}\phi_1(I_V^*)+\bar{\beta_2}\phi_2(I_H^*), \ a_{22} = \omega-\bar{\beta_2}\phi_2^{'}(I_H^*)S_H^*, $ $ a_{32} = \beta_3\phi_3^{'}(I_H^*)(\frac{\Lambda}{m}-I_V^*), \ a_{33} = \beta_3\phi_3(I_H^*)+m, \ a_{44} = \alpha_1, \ a_{55} = \alpha_2. $ Apparently, $ a_{ij} > 0 $.

    By denoting

    $ A = \left(a11a12a13a14a15a21a22a13a14a150a32a3300000a4400000a55
    \right),\ G = \left(000σ1σ2
    \right), $
    $ \boldsymbol{B(t)} = (0,0,0,B_1(t),B_2(t))^{T},\ \boldsymbol{Z}(t) = (z_{1},z_{2},z_{3},x_{1},x_{2})^{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) = e^{At}X(0)+\int_{0}^{t}e^{A(t-s)}GdB(t), $

    since $ \int_{0}^{t}e^{A(t-s)}GdB(t) $ obeys a normal distribution $ \mathcal{N}(0, \hat{\Sigma}(t)) $ at time t, where $ \hat{\Sigma}(t) = \int_{0}^{t}e^{A(t-s)}GG^Te^{A^{T}(t-s)}ds, $ then, we can get $ X(t) \sim \mathcal{N}(e^{At}X(0), \hat{\Sigma}(t)) $.

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

    $ \varphi_{A}(\lambda) = \left|λ+a11a12a13a14a15a21λ+a22a13a14a150a32λ+a3300000λ+a4400000λ+a55
    \right| = \left(\lambda+a_{44}\right)\left(\lambda+a_{55}\right)\left|λ+a11a12a13a21λ+a22a230a32λ+a33
    \right|. $

    Obviously there are $ \lambda_1 = -a_{44}, \lambda_2 = -a_{55} $, 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

    $ \lim\limits_{t \rightarrow +\infty}e^{At}X(0) = 0, $
    $ \Sigma \triangleq \lim\limits_{t\rightarrow +\infty}\hat{\Sigma} = \lim\limits_{t \rightarrow +\infty}\int_{0}^{t}e^{A(t-s)}GG^Te^{A^{T}(t-s)}ds = \int_{0}^{+\infty}e^{A^Tt}G^2e^{At}dt. $

    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 $ (S_H(0), I_H(0), I_V(0), \beta_1(0), \beta_2(0)) \in \Gamma $, if $ R_0 > 1 $ and $ m-\mu+\beta_3\phi_3(I_H^*) > 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 $ \Sigma $ 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 $ \phi_1(I_V) = I_V $, $ \phi_2(I_H) = I_H $, $ \phi_3(I_H) = I_H $, letting $ x_i(t) = \log \beta_i(t)-\log \bar{\beta}_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 $ \xi_{i, j} $ be random variables that follow a Gaussian distribution $ \mathcal{N}(0, 1) $ for $ i = 1, 2 $ and $ j = 1, 2, ..., n $. The time interval is denoted by $ \Delta t > 0 $. The values $ (S_H^{j}, I_H^{j}, I_V^{j}, x_{1}^{j}, x_{2}^{j}) $ 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

    $ \mu = 0.05, K = 100, \bar{\beta_1} = 0.15, \bar{\beta_2} = 0.1, \beta_3 = 0.1, \omega = 0.8, d = 0.5, $
    $ \gamma = 0.25, \Lambda = 5, m = 0.5, \alpha_1 = 0.8, \alpha_2 = 0.8, \sigma_1 = 0.1, \sigma_2 = 0.1, $

    after calculation, we can get $ m-\mu+\beta_3I_H^* = 2.0391 > 0 $ and the indexes of deterministic system and stochastic system can be obtained respectively, as shown below

    $ R_0 = \frac{\bar{\beta_2}K}{\omega}+\frac{\bar{\beta_1}\beta_3\Lambda k}{m^2\omega} = 50 > 1,\ R_0^s = \frac{\tilde{\beta_2}K}{\omega}+\frac{\tilde{\beta_1}\beta_3\Lambda k}{m^2\omega} = 50.0365 > 1, $

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

    $ (S_H^*,\ I_H^*,\ I_V^*,\ \log\beta_1^*,\ \log\beta_2^*) = (4.6565,\ 15.8906,\ 7.6066,\ \log0.15,\ \log0.1), $
    $ \Sigma = \left(0.05290.03770.00380.00930.01300.03770.03520.00310.00720.01000.00380.00310.00040.00060.00080.00930.00720.00060.006200.01300.01000.000800.0062
    \right). $

    Therefore, we can get that the solution $ (S_H(t), I_H(t), I_V(t), \beta_1(t), \beta_2(t)) $ obeys the normal density function

    $ \Phi(S_H,I_H,I_V,\beta_1,\beta_2) \sim \mathcal{N}((4.6565,\ 15.8906,\ 7.6066,0.15,0.1)^{T},\Sigma). $

    The marginal density functions are as follows

    $ \Phi_{S_H} = 1.7345e^{-9.4518(S_H-4.6565)^{2}},\ \Phi_{I_H} = 2.1264e^{-14.2045(I_H-15.8906)^{2}},\ \Phi_{I_V} = 19.9471e^{-1250(I_V-7.6066)^{2}}. $

    Using the above parameters, we can get trajectories of $ S_H(t), \ I_H(t) $ and $ I_V(t) $ respectively, which are present in Figure 1. It is used to represent the variation of the solution $ (S_H(t), I_H(t), I_V(t)) $ in the deterministic model (1.4) and the stochastic model (1.6). Frequency histograms and marginal density function curves for $ S_H(t), I_H(t) $ and $ I_V(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 $ (S_H(t), I_H(t), I_V(t)) $ of the stochastic and deterministic systems under perturbations $ \sigma_1 = 0.1, \sigma_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 $ S_H(t), I_H(t) $ and $ I_V(t) $ are given in Figure 2, respectively, which are highly consistent. Therefore, we deduce that the solutions $ (S_H(t), I_H(t), I_V(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 $ S_H(t), I_H(t) $ and $ I_V(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, $ \bar{\beta_1} = 0.015, \bar{\beta_2} = 0.001, \beta_3 = 0.015. $ These are denoted as Number 2. The crucial value $ R_0^E $ takes the form

    $ R_{0}^{E} = R_{0}+(e^{\frac{\sigma_1^2}{\alpha_1}}-2e^{\frac{\sigma_1^2}{4\alpha_1}}+1)^\frac{1}{2}+\frac{K\bar{\beta_2}}{m}(e^{\frac{\sigma_2^2}{\alpha_2}}-2e^{\frac{\sigma_2^2}{4\alpha_2}}+1)^\frac{1}{2} = 0.7986 < 1 $

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

    Figure 3.  The trajectory of solution $ (S_H(t), I_H(t), I_V(t)) $ under the condition $ R_0^{E} < 1 $.

    Further, by choosing the parameters in Numbers 1 and 2, respectively, the left panel in Figure 4 satisfies the condition $ R_0^s > 1 $ and the right panel satisfies $ R_0^E < 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 $ (S_H(t), I_H(t), I_V(t)) $ in the stochastic model (1.6) for $ R_0^s > 1 $, and the right panel shows the trajectories of the solutions of the stochastic model (1.6) for $ R_0^E < 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 $ \alpha $ and volatility intensity $ \sigma $ to plot the graphs, respectively. Taking $ \alpha_1 = \alpha_2 = 0.8 $ and different volatility intensity as shown in the Figure 5, the icons are red line $ \sigma_i = 0.05 $, blue line $ \sigma_i = 0.1 $ and green line $ \sigma_i = 0.15, \ i = 1, 2 $, the trends of the solution $ (S_H(t), I_H(t), I_V(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 $ \sigma_1 = \sigma_2 = 0.1 $, the reversion speed $ \alpha_i = 0.1 $ as shown by the red line in the figure, the blue line shows $ \alpha_i = 1.0 $, and similarly, the green line indicates $ \alpha_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 $ (S_H(t), I_H(t), I_V(t)) $ of the stochastic model (1.6) at the reversion speed $ \alpha_i = 0.8, \ i = 1, 2 $ and different volatility intensity is shown in the icon with the red line $ \sigma_1 = \sigma_2 = 0.05 $, the blue line $ \sigma_1 = \sigma_2 = 0.1 $ and the green line $ \sigma_1 = \sigma_2 = 0.15 $.
    Figure 6.  Trajectory plots of the solution $ (S_H(t), I_H(t), I_V(t)) $ of the stochastic model (1.6) at the volatility intensity $ \sigma_i = 0.1, \ i = 1, 2 $ and different reversion is shown in the icon with the red line $ \alpha_1 = \alpha_2 = 0.1 $, the blue line $ \alpha_1 = \alpha_2 = 1 $ and the green line $ \alpha_1 = \alpha_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 $ R_0, R_0^s, $ and $ R_0^E $ under different volatility intensity and different reversion speed, respectively, and the range of the two variables we choose is $ \left[0, 1\right] $. Combining the information in the two figures, it can be concluded that higher reversion speed and lower volatility intensity can make $ R_0^E $ and $ R_0^s $ smaller.

    Figure 7.  Trend plots of $ R_0, R_0^s, $ and $ R_0^E $ at fixed reversion speed $ \alpha_1 = \alpha_2 = 0.8 $ and volatility intensity $ \sigma_1 \in \left[0.01, 1\right] $, $ \sigma_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 $ R_0, R_0^s, $ and $ R_0^E $ at fixed volatility intensity $ \sigma_1 = \sigma_2 = 0.1 $ and reversion speed $ \alpha_1 \in \left[0.01, 1\right] $, $ \alpha_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 $ \alpha $ and volatility intensity $ \sigma $ on $ R_0^s $ and $ R_0^E $ and their magnitude relationships under different conditions.

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

    Figure 9.  Color plot of the trend of $ R_0^E $ and $ R_0^s $ with variables $ (\alpha_1, \sigma_1) \in [0.5, 0.6] \times [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 $ R_0^E $ increases with the increase of $ \sigma_1 $ and decreases with the increase of $ \alpha_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 $ R_0^s $ and $ R_0^E $ are less than 1, and $ R_0^E $ is greater than $ R_0^s $.

    Figure 10.  Color plot of the trend of $ R_0^E $ and $ R_0^s $ with variables $ (\alpha_1, \sigma_1) \in [0.5, 0.6] \times [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 $ (S_H(0), I_H(0), I_V(0)) $, we aim to examine the time it takes for the system to evolve from this initial state to either a stationary state ($ MFPT_1 $) or to an extinction state ($ MFPT_2 $).

    Then we define $ \tau_1 $ as the $ FPT $ from the initial state to the persistent state, and $ \tau_2 $ as the $ FPT $ from the initial state to the extinct state

    $ \tau_1 = \inf\left\{t:S_H < S_H^*, I_H > I_H^*, I_V > I_V^*\right\}, $
    $ \tau_2 = \inf\left\{t:I_{H0} < 0.0001,I_{V0} < 0.0001\right\}. $

    Then we have

    $ MFPT_1 = E(\tau_1),\ MFPT_2 = E(\tau_2). $

    Using Monte Carlo numerical simulation method, if $ S_H(n\Delta t) < S_H^* $, $ I_H(n\Delta t) > I_H^* $, $ I_V(n\Delta t) > I_V^* $, then $ \tau_1 = n\Delta t $, assuming that the number of simulations is $ N $, then

    $ MFPT_1 = \frac{\sum_ {i = 1}^ {N} n_i\Delta t}{N}. $

    Similarly, if $ I_{H0} < 0.0001 $ and $ I_{V0} < 0.0001 $, then $ \tau_2 = m\Delta t $ and

    $ MFPT_2 = \frac{\sum_ {i = 1}^ {N} n_i\Delta t}{N}. $

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

    Figure 11.  The mean first passage time for transitioning from the initial state values $ (S_H(0), I_H(0), I_V(0)) $ = (3, 1, 1) to the state of stationary with $ \sigma_i \in \left[0.01, 0.1\right] $, $ \alpha_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 $ (S_H(0), I_H(0), I_V(0)) $ = (5,200, 10) to the state of extinction with $ \sigma_i \in \left[0.01, 0.1\right] $, $ \alpha_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 $ R_0 $ are obtained, and their global asymptotic stability are deduced. Specificaly, the endemic equilibrium is globally asymptotically stable if $ R_0 > 1 $, and the disease free equilibrium is globally asymptotically stable if $ R_0 < 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.

    $ \bullet $ When $ R_0^s > 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.

    $ \bullet $ Also, from Figure 3, we can further find that when $ R_0^E < 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.

    $ \bullet $ 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.

    $ \bullet $ Then, correlation plots of $ R_0^s $, $ R_0^E $ with reversion speed and volatility intensity are obtained in Figures 710. We conculde that higher reversion speed and lower volatility intensity can make $ R_0^E $ and $ R_0^s $ more smaller.

    $ \bullet $ Finally, Figures 11 and 12 visually demonstrate the relationship between $ MFPT $ and the $ \alpha_i $ and $ \sigma_i, \ i = 1, 2 $ in a stochastic system (1.6) with bilinear incidence. We can see that if voltility intensity $ \sigma_i $ is much bigger (or the reversion speed $ \alpha_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

    $ R_H = K-S_H-I_H,\ S_V = \frac{\Lambda}{m}-I_V, $

    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 $ E_0 = (S_{H0}, \ 0, \ 0, \ S_{V0}, \ 0) = (K, \ 0, \ 0, \ \frac{\Lambda}{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, $ I_H \in [0, \frac{\mu K}{\mu+\gamma}] $, let

    $ H(I_H) = \frac{\bar{\beta_1}\phi_1(I_V)S_H}{\omega}+\frac{\bar{\beta_2}\phi_2(I_H)S_H}{\omega}-I_H, $

    then by calculation, $ H(0) = 0, H(\frac{\mu K}{\mu+\gamma}) = -\frac{\mu K}{\mu+\gamma} < 0 $, and

    $ H^{'}(I_H) = \frac{\bar{\beta_1}\phi_1^{'}(I_V)S_H}{\omega}\frac{m^2\beta_3\phi_3^{'}(I_H)\Lambda}{(m^2+m\beta_3\phi_3(I_H))^2}-(1+\frac{\gamma}{\mu})\frac{\bar{\beta_1}\phi_1(I_V)}{\omega}+\frac{\bar{\beta_2}\phi_2^{'}(I_H)S_H}{\omega}-(1+\frac{\gamma}{\mu})\frac{\bar{\beta_2}\phi_2(I_H)}{\omega}-1. $

    If $ R_0 > 1 $, then

    $ H^{'}(0) = \frac{\bar{\beta_2}K\phi_2^{'}(0)}{\omega}+\frac{\bar{\beta_1} \beta_3\Lambda K\phi_1^{'}(0)\phi_3^{'}(0)}{m^2\omega}-1 = R_0-1 > 0, $
    $ H^{'}(\frac{\mu K}{\mu+\gamma}) = -(1+\frac{\gamma}{\mu})[\frac{\bar{\beta_1}\phi_1(I_V(\frac{\mu K}{\mu+\gamma}))}{\omega}+\frac{\bar{\beta_2}\phi_2(\frac{\mu K}{\mu+\gamma})}{\omega}]-1 < 0, $

    therefore, there exists a point $ \xi $ such that $ H^{'}(I_H) > 0 $ on $ \left[0, \xi\right) $ and $ H^{'}(I_H) < 0 $ on $ \left(\xi, \frac{\mu K}{\mu+\gamma}\right] $, i.e., $ H(I_H) $ is monotonically increasing on $ \left[0, \xi\right) $ and monotonically decreasing on $ \left(\xi, \frac{\mu K}{\mu+\gamma}\right] $.

    Hence, there is a unique $ I_H^* \in (0, \frac{\mu K}{\mu+\gamma}) $ such that $ H(I_H^*) = 0 $, which implies that when $ R_0 > 1 $, system (1.4) has a unique endemic equilibrium $ E^* = (S_H^*, I_H^*, R_H^*, S_V^*, I_V^*) $. 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 $ E_0 $ into the matrix $ J $ to get $ J_0 $

    $ J_0 = \left(μ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
    \right), $

    the corresponding characteristic polynomial is as follows

    $ \phi_{J_0}(\lambda) = (\lambda+\mu)(\lambda+\mu)(\lambda+m)(\lambda+m)(\lambda+(\omega-\bar{\beta_2}\Phi_2^{'}(0))), $

    obviously, if $ R_0 < 1 $, then $ \omega-\bar{\beta_2}\Phi_2^{'}(0) > 0 $, according to the Routh-Hurwitz criterion, there are only negative real part characteristic roots, so the disease-free equilibrium $ E_0 $ is locally asymptotically stable.

    If $ R_0 > 1 $, then $ \omega-\bar{\beta_2}\Phi_2^{'}(0) < 0 $, this indicates that $ J_0 $ has the eigenvalue of the positive real part, so the disease-free equilibrium $ E_0 $ is unstable.

    Next we prove the global attractiveness of $ E_0 $, define

    $ V = S_H-K-K\log\frac{S_H}{K}+I_H+S_V-\frac{\Lambda}{m}-\frac{\Lambda}{m}\log\frac{mS_V}{\Lambda}+I_V, $

    Using It$ \hat{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 $ S_H \rightarrow K \ as \ t \rightarrow \infty, \ S_V \rightarrow \frac{\Lambda}{m} \ as\ t \rightarrow \infty $, and $ I_H, \ I_V \rightarrow 0 \ as\ t \rightarrow \infty $. Then, $ \mathcal{L} V \leq 0 $ and equal to 0 when it takes $ E_0 $. Therefore, according to the LaSalle invariance principle, we conclude that when $ R_0 < 1 $, $ E_0 $ 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

    $ V_1 = S_H-S_H^*-S_H^*log\frac{S_H}{S_H^*}+I_H-I_H^*-I_H^*log\frac{I_H}{I_H^*}, $

    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

    $ V_2 = S_V-S_V^*-S_V^*log\frac{S_V}{S_V^*}+I_V-I_V^*-I_V^*log\frac{I_V}{I_V^*}, $

    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

    $ V_3 = V_1+\frac{\bar{\beta_1}\Phi_1(I_V^*)S_H^*}{\bar{\beta_3}\Phi_3(I_H^*)S_V^*}V_2, $

    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 ($ \mathcal{A}_2 $) and ($ \mathcal{A}_3 $), we can know

    $ \frac{\Phi_1({I_V})}{\Phi_1({I_V^*})}+\frac{I_V\Phi_1({I_V^*})}{I_V^*\Phi_1({I_V})}-\frac{I_V}{I_V^*}-1 = \frac{I_V}{\Phi_1(I_V)\Phi_1(I_V^*)}\left[(\Phi_1(I_V)-\Phi_1(I_V^*))(\frac{\Phi_1(I_V)}{I_V}-\frac{\Phi_1(I_V^*)}{I_V^*})\right]\leq 0, $

    $ \Phi_2(I_H) $ and $ \Phi_3(I_H) $ similarly satisfy the structure of the above equation, combined with $ S_H $ and $ S_V $ are bounded, we finally get

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

    Next, we define

    $ V_4 = \frac{(S_H-S_H^*+I_H-I_H^*)^2}{2},\ V_5 = \frac{(R_H-R_H^*)^2}{2},\ V_6 = \frac{(S_V-S_V^*+I_V-I_V^*)^2}{2}, $

    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 = V_3+A_1V_4+A_2V_5+A_3V_6, $

    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

    $ A_1 = \frac{\mu(d+\gamma)}{K\omega^2},\ A_2\frac{\gamma^2}{\mu} = \frac{d+\gamma}{2}A_1,\ A_3 = \frac{m^2}{\Lambda}\frac{\bar{\beta_1}\Phi_1(I_V^*)S_H^*}{2m\bar{\beta_3}\Phi_3(I_H^*)S_V^*}, $

    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 $ (S_H(t), I_H(t), I_V(t), \beta_1(t), \beta_2(t)) $ on $ t \in \left[0, \tau_{e}\right) $, where $ \tau_e $ represents explosion time. To show that the solution is global, according to the method in [33], we just need to verify that $ \tau_e = +\infty \ a.s. $.

    Choose $ k_0 $ be a sufficiently large integer for every component of $ (S_H(0), I_H(0), I_V(0), \beta_1(0), \beta_2(0)) $ within the interval $ [\frac{1}{k_0}, k_0] $. For each integer $ k \geq k_0 $, define the stopping time as

    $ \tau_{k} = \inf \left\{t \in\left[0, \tau_{e}\right)| min(S_H(t), I_H(t), I_V(t), \beta_1(t), \beta_2(t))\leq \frac{1}{k}\ or \ max(S_H(t), I_H(t), I_V(t), \beta_1(t), \beta_2(t)) \geq k \right\}. $

    It can be seen that $ \tau_k $ is monotonically increasing with respect to $ k $. Then define $ \inf \{\varnothing\} = +\infty $ and $ \tau_{\infty} = \lim _{t \rightarrow +\infty} \tau_{k} $. It is clearly visible that the solution is global due to the fact that $ \tau_\infty < \tau_e \ a.s. $, $ \tau_\infty = \infty $ leading to $ \tau_e = \infty $. Next, we will prove $ \tau_\infty = \infty $ by the contradiction method. Assuming $ \tau_{\infty} < +\infty $ a.s., then there are $ \varepsilon_{0} \in \left(0, 1\right) $ and $ T > 0 $ such that $ P(\tau_\infty\leq T) > \varepsilon $, so there is a positive integer $ k_1 > k_0 $ that makes

    $ P(\tau_k \leq T) \geq \varepsilon,\ \forall k \geq k_1. $

    Define a non-negative $ C^2 $-function $ V(S_H, I_H, I_V, \beta_1, \beta_2) $ as follows

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

    Applying It$ \hat{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 $ \mathcal{A}_3 $, we can know that

    $ \phi_1(I_V) \leq \phi_1^{'}(0)I_V,\phi_2(I_H) \leq \phi_2^{'}(0)I_H, \phi_3(I_H) \leq \phi_3^{'}(0)I_H, $

    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(\beta_1, \beta_2) \rightarrow -\infty \ as \ \beta_1 \rightarrow +\infty, \beta_1 \rightarrow 0, \beta_2 \rightarrow +\infty, \beta_2 \rightarrow 0 $, so there's a positive constant $ H_0 $ that makes $ \mathcal{L}V\leq H_0 $. 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 $ \zeta \in G_k $, $ W\left(S_H\left(\tau_{k}, \zeta\right), I_H\left(\tau_{k}, \zeta\right), I_V\left(\tau_{k}, \zeta\right), \beta_1\left(\tau_{k}, \zeta\right), \beta_2\left(\tau_{k}, \zeta\right)\right) $ will larger than $ \left(e^k-1-k\right) \wedge\left(e^{-k}-1+k\right) $, 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 \rightarrow +\infty $

    $ +\infty\leq E V\left(S_H(0), I_H(0), I_V(0), \beta_1(0), \beta_2(0)\right)+H_0 T < +\infty. $

    Therefore $ \tau_{\infty} = +\infty $ a.s., i.e., $ \tau_{e} = +\infty $. Then system (1.6) has a unique global solution $ \left(S_H(t), I_H(t), I_V(t), \beta_1(t), \beta_2(t)\right) $ on $ \Gamma $.

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

    Step 1. Construct Lyapunov functions

    Using It$ \hat{o} $'s formula, we obtain

    $ \mathcal{L}(-\log S_H) = -\frac{\mu K}{S_H}+\mu+\beta_1\phi_1(I_V)+\beta_2\phi_2(I_H)-\frac{dI_H}{S_H}, $
    $ \mathcal{L}(-\log I_H) = -\frac{\beta_1\phi_1(I_V)S_H}{I_H}-\frac{\beta_2\phi_2(I_H)S_H}{I_H}+\omega, $
    $ \mathcal{L}(-\log I_V) = -\frac{\Lambda}{m}\beta_3\phi_3(I_H)\frac{1}{I_V}+\beta_3\phi_3(I_H)+m. $

    Define a function $ V_1 = -\log I_H-c_1\log S_H-c_2\log S_H-c_3\log I_V $, where $ c_1, c_2, c_3 $ are given in subsequent calculation. Using It$ \hat{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 $ \mathcal{A}_3 $, notice that

    $ (\frac{I_V}{\phi_1(I_V)})^{'} = \frac{\frac{I_V}{\phi_1(I_V)}-\frac{1}{\phi_1(I_V)}}{I_V}\leq m_0, $

    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 $ c_{1}, c_{2}, c_3, c_4, c_5 $ and $ c_6 $ satisfy the following equalities

    $ c_1\mu = c_3m = c_4\frac{1}{\phi_1^{'}(0)} = c_5\frac{1}{\phi_3^{'}(0)} = \frac{\tilde{\beta_1}\beta_3\Lambda K\phi_1^{'}(0)\phi_3^{'}(0)}{m^2}, $
    $ c_2\mu = c_6\frac{1}{\phi_2^{'}(0)} = \tilde{\beta_2}K\phi_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

    $ R_0^s = \frac{\tilde{\beta_1}\beta_3\Lambda K\phi_1^{'}(0)\phi_3^{'}(0)}{m^2\omega}+\frac{\tilde{\beta_2}K\phi_2^{'}(0)}{\omega}. $

    By Holder inequality, for any positive constant $ \delta $, 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 $ \delta $ to be

    $ \delta = \frac{\frac{\omega}{2}(R_0^s-1)}{(\frac{\bar{\beta_1}\beta_3\mu K\phi_1^{'}(0)\phi_3^{'}(0)}{\mu m^2}+\frac{\bar{\beta_2}K\phi_2^{'}(0)}{\mu})(\phi_2^{'}(0)K\bar{\beta_2}^2e^{\frac{\sigma_2^2}{\alpha_2}}+\phi_1^{'}(0)\frac{\Lambda}{m}\bar{\beta_1}^2e^{\frac{\sigma_1^2}{\alpha_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

    $ V_2 = V_1+\frac{(c_1+c_2)\phi_1^{'}(0)+4\delta c_4m_0}{4m\delta}I_V, $

    applying It$ \hat{o} $'s formula to $ V_{2} $, 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 = (c_1+c_2)\phi_2^{'}(0)\frac{1}{4\delta}+c_3\phi_3^{'}(0)\beta_3+(c_5+c_6)m_0+\frac{((c_1+c_2)\phi_1^{'}(0)+4\delta c_4m_0)\beta_3\phi_3^{'}(0)\Lambda}{4m^2\delta}. $

    Next, define

    $ V_3 = -\log S_H-\log I_V-\log(K-(S_H+I_H))-\log(\frac{\Lambda}{m}-I_V)+(\beta_1-1-\log\beta_1)+(\beta_2-1-\log\beta_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

    $ \bar{V} = MV_2+V_3, $

    where $ M $ satisfying the following inequality

    $ -\frac{M\omega}{2}(R_0^s-1)+\sup\limits_{(\beta_1, \beta_2)\in R_{+}^{2}}W(\beta_1, \beta_2)\leq-2. $

    Notice that $ \bar{V} $ has a minimum value $ \bar{V}_{min} $ in the interior of $ \Gamma $ because $ \bar{V} \rightarrow +\infty $ as $ (S_H, I_H, I_V, \beta_1, \beta_2) $ tends to the boundary of $ \Gamma $. Ultimately, establish a non-negative $ C^2 $-function $ V(S_H, I_H, I_V, \beta_1, \beta_2): \Gamma \rightarrow R_+ $ as follows

    $ V(S_H, I_H, I_V, \beta_1, \beta_2) = \bar{V}(S_H, I_H, I_V, \beta_1, \beta_2) - \bar{V}_{min}, $

    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_{\varepsilon} $

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

    where $ \varepsilon $ is a small enough constant and the complement of $ U_{\varepsilon} $ can be divided into nine small sets as follows

    $ U_{1,\varepsilon}^{c} = \{(S_H,I_H,I_V,\beta_1,\beta_2)\in \Gamma|0 < \beta_1 < \varepsilon\}, U_{2,\varepsilon}^{c} = \{(S_H,I_H,I_V,\beta_1,\beta_2)\in \Gamma|\beta_1 > \frac{1}{\varepsilon}\}, $
    $ U_{3,\varepsilon}^{c} = \{(S_H,I_H,I_V,\beta_1,\beta_2)\in \Gamma|0 < \beta_2 < \varepsilon\}, U_{4,\varepsilon}^{c} = \{(S_H,I_H,I_V,\beta_1,\beta_2)\in \Gamma|\beta_2 > \frac{1}{\varepsilon}\}, $
    $ U_{5,\varepsilon}^{c} = \{(S_H,I_H,I_V,\beta_1,\beta_2)\in \Gamma|0 < I_H < \varepsilon\}, U_{6,\varepsilon}^{c} = \{(S_H,I_H,I_V,\beta_1,\beta_2)\in \Gamma|0 < S_H < \varepsilon\}, $
    $ U_{7,\varepsilon}^{c} = \{(S_H,I_H,I_V,\beta_1,\beta_2)\in \Gamma|0 < I_V < \varepsilon^2,I_H\geq\varepsilon \}, $
    $ U_{8,\varepsilon}^{c} = \{(S_H,I_H,I_V,\beta_1,\beta_2)\in \Gamma|S_H+I_H > K-\varepsilon^2,I_H\geq\varepsilon\}, $
    $ U_{9,\varepsilon}^{c} = \{(S_H,I_H,I_V,\beta_1,\beta_2)\in \Gamma|I_V > \frac{\Lambda}{m}-\varepsilon^3,I_V\geq\varepsilon^2\}, $

    then the following results hold

    Case 1: $ (S_H, I_H, I_V, \beta_1, \beta_2)\in U_{1, \varepsilon}^{c} $, 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: $ (S_H, I_H, I_V, \beta_1, \beta_2)\in U_{2, \varepsilon}^{c} $, 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: $ (S_H, I_H, I_V, \beta_1, \beta_2)\in U_{3, \varepsilon}^{c} $, 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: $ (S_H, I_H, I_V, \beta_1, \beta_2)\in U_{4, \varepsilon}^{c} $, 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: $ (S_H, I_H, I_V, \beta_1, \beta_2)\in U_{5, \varepsilon}^{c} $, 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: $ (S_H, I_H, I_V, \beta_1, \beta_2)\in U_{6, \varepsilon}^{c} $, 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: $ (S_H, I_H, I_V, \beta_1, \beta_2)\in U_{7, \varepsilon}^{c} $, 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: $ (S_H, I_H, I_V, \beta_1, \beta_2)\in U_{8, \varepsilon}^{c} $, 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: $ (S_H, I_H, I_V, \beta_1, \beta_2)\in U_{9, \varepsilon}^{c} $, 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\left(S_H, I_H, I_V, \beta_1, \beta_2\right) \leq-1, \quad \forall\left(S_H, I_H, I_V, \beta_1, \beta_2 \right) \in \Gamma \backslash U_{\varepsilon}, $

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

    $ G\left(S_H, I_H, I_V, \beta_1, \beta_2 \right) \leq H < +\infty, \quad \forall\left(S_H, I_H, I_V, \beta_1, \beta_2 \right) \in \Gamma. $

    For any initial value $ (S_H(0), I_H(0), I_V(0), \beta_1(0), \beta_2(0))\in\Gamma $, 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 $ \beta_i(i = 1, 2) $ is ergodic according to [34,35], then we can get that

    $ \lim _{t \rightarrow +\infty} \frac{1}{t} \int_{0}^{t} \beta_i^p(\tau) d\tau = \lim _{t \rightarrow +\infty} \frac{1}{t} \int_{0}^{t}e^{plog\beta_i(\tau)}d\tau = \int_{-\infty}^{+\infty}e^{py_i}\pi(y_i)dy_i = \bar{\beta_i}^pe^{\frac{p^2\sigma_i^2}{4\alpha_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 \rightarrow +\infty $ 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

    $ \liminf\limits_{t \rightarrow +\infty} \frac{1}{t} \int_{0}^{t} P\left\{\left(S_H(\tau), I_H(\tau), I_V(\tau), \beta_1(\tau), \beta_2(\tau)\right) \in U_{\varepsilon}\right\} d \tau \geq \frac{1}{H+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 $ R_0^s > 1 $ system (1.6) has a stationary distribution on $ \Gamma $.

    Proof. Define a $ C^2 $-function $ G(I_H, I_V, \beta_1, \beta_2): \Gamma \rightarrow R $ by

    $ G(I_H, I_V, \beta_1, \beta_2) = v_1I_H+v_2I_V, $

    where $ v_1 = R_0, v_2 = \frac{\bar{\beta_1}\phi_1^{'}(0)K}{m} $. Applying It$ \hat{o} $'s formula to $ G(I_H, I_V, \beta_1, \beta_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 $ \beta_1, \beta_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 \rightarrow +\infty $, 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

    $ R_{0}^{E} = R_{0}+\frac{mR_0(e^\frac{\sigma_1^2}{\alpha_1}-2e^\frac{\sigma_1^2}{4\alpha_1}+1)^\frac{1}{2}+\phi_2^{'}(0)K\bar{\beta_2}(e^\frac{\sigma_2^2}{\alpha_2}-2e^\frac{\sigma_2^2}{4\alpha_2}+1)^\frac{1}{2}}{\min\{m,\frac{\bar{\beta_2}\phi_2^{'}(0)K}{R_0}\}}. $

    If $ R_{0}^{E} < 1 $,

    $ \limsup _{t \rightarrow +\infty} \frac{\log G(t)}{t} < 0 $

    will be true which indicates

    $ \lim\limits_{t \rightarrow +\infty}I_H(t) = 0 \lim\limits_{t \rightarrow +\infty}I_V(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 $ G_1 = diag(0, 0, 0, \sigma_1, 0) $.

    Let $ A_{1} = J_{1}AJ_{1}^{-1} $, where

    $ J_{1} = \left(0001010000010000010000001
    \right), $

    then

    $ A_{1} = \left(a440000a14a11a12a13a15a14a21a22a13a1500a32a3300000a55
    \right). $

    Let $ A_{2} = J_{2}A_{1}J_{2}^{-1} $, where

    $ J_{2} = \left(1000001000011000001000001
    \right), $

    and

    $ A_{2} = \left(a440000a14a12a11a12a13a150a12a11+a21+a22a12a22000a32a32a3300000a55
    \right). $

    Due to $ a_{12}-a_{11}+a_{21}+a_{22} = \gamma > 0 $, let $ A_{3} = J_{3}A_{2}J_{3}^{-1} $, where

    $ J_{3} = \left(10000010000010000a32γ1000001
    \right), $

    and

    $ A_{3} = \left(a440000a14a12a11a13a32γa12a13a150γa12a220000wa3300000a55
    \right). $

    in which

    $ w = a_{32}-\frac{a_{32}(a_{12}+a_{22})}{\gamma}+\frac{a_{32}a_{33}}{\gamma} = m-\mu+\beta_3\phi_3(I^*) > 0. $

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

    $ M = \left(m1m2m3m4m50wγw(a12+a22+a33)a233000wa3300001000001
    \right), $

    where $ m_{1} = -w\gamma a_{14}, \ m_{2} = -w\gamma(a_{33}+a_{11}+a_{22}), \ m_{3} = wa_{13}a_{32}-w\gamma a_{12}+w(a_{12}+a_{22}+a_{33})(a_{12}+a_{22})+wa_{33}^2, \ m_{4} = -\gamma wa_{13}-a_{33}^3, m_{5} = -\gamma wa_{15}. $

    Define $ A_{01} = MA_3M^{-1} $, then we can get

    $ A_{01} = \left(b1b2b3b4b51000001000001000000a55
    \right), $

    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 = J_3J_2J_1 $, 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)G_1^2(MJ)^T = diag((m_{1}\sigma_1)^2, 0, 0, 0, 0) $, let $ \rho_1 = m_{1}\sigma_1 $, then (A.11) becomes

    $ G_0^2+\rho_1^{-2}A_{01}[(MJ)\Sigma_1(MJ)^T]+\rho_1^{-2}[(MJ)\Sigma_1(MJ)^T]A_{01}^T = 0, $

    then we obtian

    $ \Sigma_{01}: = \rho_1^{-2}(MJ)\Sigma_1(MJ)^T = \left(b1b4b2b3b0b3b000b3b0b1b0b3b0b1b000b1b0b3b1b2b4b000000
    \right) , $

    where $ b = 2[b_4b_1^2-b_1b_2b_3+b_3^2] $. We can obtain that the matrix $ \Sigma_{01} $ is a positive semi-definite matrix, the exact expression of $ \Sigma_1 $ is as follows

    $ \Sigma_1 = \rho_1^{2}(MJ)^{-1}\Sigma_{01}[(MJ)^{-1}]^{T}. $

    Step 2 Consider the following equation

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

    where $ G_2 = diag(0, 0, 0, 0, \sigma_2) $.

    Let $ B_{1} = P_{1}AP_{1}^{-1} $, where

    $ P_{1} = \left(0000110000010000010000010
    \right), $

    then

    $ B_{1} = \left(a550000a15a11a12a13a14a15a21a22a13a1400a32a3300000a44
    \right). $

    Let $ B_{2} = P_{2}B_1P_{2}^{-1} $, where

    $ P_{2} = \left(1000001000011000001000001
    \right), $

    then

    $ B_{2} = \left(a550000a15a12a11a12a13a140a12a11+a21+a22a12a22000a32a32a3300000a44
    \right). $

    Similarly, due to $ a_{12}-a_{11}+a_{21}+a_{22} = \gamma > 0 $, let $ B_{3} = P_{3}B_{2}P_{3}^{-1} $, where

    $ P_{3} = \left(10000010000010000a32γ1000001
    \right), $

    and

    $ B_{3} = \left(a550000a15a12a11a13a32γa12a13a140γa12a220000wa3300000a44
    \right), $

    in which

    $ w = a_{32}-\frac{a_{32}(a_{12}+a_{22})}{\gamma}+\frac{a_{32}a_{33}}{\gamma} = m-\mu+\beta_3\phi_3(I^*) > 0. $

    The standard transformation matrix of $ B_3 $ has the following form

    $ N = \left(n1n2n3n400wγw(a33+a12+a22)a233000wa3300001000001
    \right), $

    where $ n_{1} = -w\gamma a_{15}, \ n_{2} = -w\gamma(a_{33}+a_{11}+a_{22}), \ n_{3} = wa_{13}a_{32}-w\gamma a_{12}+w(a_{12}+a_{22}+a_{33})(a_{12}+a_{22})+wa_{33}^2, \ n_{4} = -\gamma wa_{13}-a_{33}^3, n_{5} = -\gamma wa_{14}. $

    Define $ B_{01} = NB_3N^{-1} $, then we can get

    $ B_{01} = \left(d1d2d3d4d51000001000001000000a44
    \right), $

    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 = P_3P_2P_1 $, 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)G_2^2(NP)^T = diag((n_{1}\sigma_2)^2, 0, 0, 0, 0) $, let $ \rho_2 = n_{1}\sigma_2 $, then (A.13) becomes

    $ G_0^2+\rho_2^{-2}B_{01}[(NP)\Sigma_2(NP)^T]+\rho_3^{-2}[(NP)\Sigma_2(NP)^T]B_{01}^T = 0, $

    then we obtian

    $ \Sigma_{02}: = \rho_2^{-2}(NP)\Sigma_2(NP)^T = \left(d1d4d2d3d0d3d000d3d0d1d0d3d0d1d000d1d0d3d1d2d3d000000
    \right) , $

    where $ d = 2(d_4d_1^2-d_1d_2d_3+d_3^2) $, and we can obtain that the matrix $ \Sigma_{02} $ is a positive semi-definite matrix, the exact expression of $ \Sigma_2 $ is as follows

    $ \Sigma_2 = \rho_2^{2}(NP)^{-1}\Sigma_{02}[(NP)^{-1}]^{T}. $

    Finally, $ \Sigma = \Sigma_1+\Sigma_2 $. Obviously, the matrix $ \Sigma $ is a positive definite matrix. The proof is complete.

  • Reader Comments
  • © 2011 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(3156) PDF downloads(648) Cited by(17)

Other Articles By Authors

/

DownLoad:  Full-Size Img  PowerPoint
Return
Return

Catalog