The trifecta for curriculum sustainability in Australian universities

  • Commercialization and internationalization of tertiary education has opened a new way for universities to grow and make more profit. This in turn has supported sustainable growth of the higher education sector for the last three decades in developed countries. Curricula offered by accredited tertiary institutions must meet the quality standards set by both the governmental agencies and the professional accreditation bodies. These programs must also provide graduates with employment opportunities. Hence, quality and employment opportunities are the two key factors for sustainability of any degree program offered by tertiary institutions. However, changes in regulations and policies by the national government sometimes play a vital role in creating new programs, and maintaining or disestablishing some existing programs offered by institutions in the nation. These changes are not controlled by individual institutions, which has become the third unpredictable factor in curriculum creation and/or sustainability. Using the journey of a new master's program in information technology (IT) in an Australian university as a case study, we explore how this third factor impacted on the initialization, creation, and short life of this program primarily targeting international students in the mid-2000s. We then extend our discussion to the implications of the recent changes to tertiary tuitions from 2021 by the Australian Government on the sustainable future of the Australian tertiary education sector on a broad scale.

    Citation: William Guo, Wei Li, Roland Dodd, Ergun Gide. The trifecta for curriculum sustainability in Australian universities[J]. STEM Education, 2021, 1(1): 1-16. doi: 10.3934/steme.2021001

    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
  • Commercialization and internationalization of tertiary education has opened a new way for universities to grow and make more profit. This in turn has supported sustainable growth of the higher education sector for the last three decades in developed countries. Curricula offered by accredited tertiary institutions must meet the quality standards set by both the governmental agencies and the professional accreditation bodies. These programs must also provide graduates with employment opportunities. Hence, quality and employment opportunities are the two key factors for sustainability of any degree program offered by tertiary institutions. However, changes in regulations and policies by the national government sometimes play a vital role in creating new programs, and maintaining or disestablishing some existing programs offered by institutions in the nation. These changes are not controlled by individual institutions, which has become the third unpredictable factor in curriculum creation and/or sustainability. Using the journey of a new master's program in information technology (IT) in an Australian university as a case study, we explore how this third factor impacted on the initialization, creation, and short life of this program primarily targeting international students in the mid-2000s. We then extend our discussion to the implications of the recent changes to tertiary tuitions from 2021 by the Australian Government on the sustainable future of the Australian tertiary education sector on a broad scale.



    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 , notice that

    which can deduce

    (A.4)

    similarly, one gets

    (A.5)

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

    Let and satisfy the following equalities

    then

    where

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

    If take to be

    then we can get

    where

    Next we define

    applying It's formula to , it leads to

    where

    Next, define

    then, we have

    where

    Choose is a large enough positive constant, let

    where satisfying the following inequality

    Notice that has a minimum value in the interior of because as tends to the boundary of . Ultimately, establish a non-negative -function as follows

    then we obtain

    (A.6)

    Step 2. Set up the closed set

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

    then the following results hold

    Case 1: , then

    Case 2: , then

    Case 3: , then

    Case 4: , then

    Case 5: , then

    Case 6: , then

    Case 7: , then

    Case 8: , then

    Case 9: , then

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

    in other words, let is a positive constant that makes

    For any initial value , integrating the inequality (A.6) and taking the expectation, we get

    (A.7)

    One gets that is ergodic according to [34,35], then we can get that

    hence

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

    which means

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

    Proof. Define a -function by

    where . Applying It's formula to , then we have

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

    (A.8)

    According to the ergodicity of , then

    (A.9)

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

    where

    If ,

    will be true which indicates

    this means the disease will die out exponentially.

    Proof. Step 1 Consider the following equation

    (A.10)

    where .

    Let , where

    then

    Let , where

    and

    Due to , let , where

    and

    in which

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

    where

    Define , then we can get

    in which

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

    (A.11)

    where , let , then (A.11) becomes

    then we obtian

    where . We can obtain that the matrix is a positive semi-definite matrix, the exact expression of is as follows

    Step 2 Consider the following equation

    (A.12)

    where .

    Let , where

    then

    Let , where

    then

    Similarly, due to , let , where

    and

    in which

    The standard transformation matrix of has the following form

    where

    Define , then we can get

    in which

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

    (A.13)

    where , let , then (A.13) becomes

    then we obtian

    where , and we can obtain that the matrix is a positive semi-definite matrix, the exact expression of is as follows

    Finally, . Obviously, the matrix is a positive definite matrix. The proof is complete.



    [1] Developing a curriculum guarantee for overseas students. Higher Education Research & Development (1995) 14: 35-46.
    [2] Academic Standards, Internationalisation, and the Discursive Construction of "The International Student". Higher Education Research & Development (2003) 22: 155-166.
    [3] The writing on the wall: responses of Australian public universities to competition in global higher education. Higher Education Research & Development (2009) 28: 495-508.
    [4] Department of Education, Skills and Employment (DESE) (2020) Onshore higher education international students as a proportion of all students by university in 2019: Research Snapshot 2020: September.
    [5] Department of Education, Skills and Employment (DESE) (2020) International students in Australia up to 2019: Research Snapshot 2020: August.
    [6] Isserstedt, W., Schnitzer, K. (2005) Internationalization of higher education: foreign students in Germany and German students abroad, In: Federal Ministry of Education and Research, Berlin.
    [7] The internationalization of higher education: motivations and realities. Journal of Studies in International Education (2007) 11: 290-305.
    [8] Current trends in higher education: internationalization in Asia and Oceania. Comparative & International Higher Education (2009) 1: 3-4.
    [9] Internationalization of higher education: Challenges, strategies, policies and programs. US-China Education Review (2009) 6: 1-9.
    [10] The Development of international education in Australia: A framework for the future. Journal of Studies in International Education (2007) 11: 410-420.
    [11] Re-positioning Australia's international education in global knowledge economies: implications of shifts in skilled migration policies for universities. Journal of Higher Education Policy and Management (2012) 34: 341-354.
    [12] Employment outcomes of science graduates in Australia: Implications for choice and diversity in the curriculum. Tertiary Education and Management (2003) 9: 61-76.
    [13] Closing the gap: Preparing computing students for employment through embedding work-related learning in the taught curriculum. Industry and Higher Education (2012) 26: 315-320.
    [14] Surveying employment listings to inform curricula of environmental science degree programs. Journal of Environmental Studies and Sciences (2017) 7: 346-354.
    [15] Unpacking the learning-work nexus: 'priming' as lever for high-quality learning outcomes in work-integrated learning curricula. Studies in Higher Education (2015) 40: 22-42.
    [16] Assuring the quality of curriculum, pedagogy and assessment across satellite campuses. Journal of Further and Higher Education (2019) 43: 589-600.
    [17] Hammer, S., Ayriss, P., McCubbin, A. (2020) Style or substance: how Australian universities contextualise their graduate attributes for the curriculum quality space. Higher Education Research & Development, DOI: 10.1080/07294360.2020.1761304. doi: 10.1080/07294360.2020.1761304
    [18] Coulouris, G., Dollimore, J., Kindberg, T. (2001) Distributed systems: concepts and design (3rd ed.). Addison-Wesley.
    [19] Tanenbaum, A.S., van Steen, M. (2002) Distributed systems: principles and paradigms. Prentice Hall.
    [20] Incorporating statistical and neural network approaches for student course satisfaction analysis and prediction. Expert Systems with Applications (2010) 37: 3358-3365.
    [21] Analysis of student course evaluation data for an IT subject: Implications for improving STEM education. International Journal of Information and Education Technology (2017) 7: 635-640.
    [22] Department of Education, Skills and Employment (DESE) (2020) The 2019 export income by state and territory: Research Snapshot 2020: August.
    [23] Australian Government (2016) National Strategy for International Education 2025. Canberra
    [24] Chinese returnees from overseas study: An understanding of brain gain and brain circulation in the age of globalization. Frontiers of Education in China (2015) 10: 306-329.
    [25] Chinese Ministry of Education (2016) The Blue Book on the Return and Employment of Chinese Overseas Students. Beijing.
    [26] China CITIC Bank (2016) China overseas study 2016. Beijing
    [27] The Guardian, Australian universities to cut hundreds of courses as funding crisis deepens, 2020. Available from: https://www.theguardian.com/australia-news/2020/sep/30/australian-universities-to-cut-hundreds-of-courses-as-funding-crisis-deepens.
    [28] Financial Review, Second wave of university job cuts coming next year, 2020. Available from: https://www.afr.com/work-and-careers/education/second-wave-of-university-job-cuts-coming-next-year-20200917-p55wg2.
    [29] Business Insider Australia, This is how many jobs each Australian university has cut – or plans to – in 2020, 2020. Available from: https://www.businessinsider.com.au/australian-university-job-cuts-losses-tally-2020-9.
    [30] Coleborne, C. (2020) HASS experts are 'deeply disappointed' with the higher education reforms: opinion. Campus Review 2020: October.
    [31] Zaglas, W. (2020) Arts degree was 'probably the best thing that ever happened' to me: reforms may turn Indigenous Australians off university. Campus Review 2020: October.
  • Reader Comments
  • © 2021 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(3071) PDF downloads(279) Cited by(2)

Figures and Tables

Figures(6)  /  Tables(2)

Other Articles By Authors

/

DownLoad:  Full-Size Img  PowerPoint
Return
Return

Catalog