Loading [MathJax]/jax/element/mml/optable/Latin1Supplement.js
Research article Special Issues

Threshold dynamics of a stochastic general SIRS epidemic model with migration


  • In this study, a stochastic SIRS epidemic model that features constant immigration and general incidence rate is investigated. Our findings show that the dynamical behaviors of the stochastic system can be predicted using the stochastic threshold RS0. If RS0<1, the disease will become extinct with certainty, given additional conditions. Conversely, if RS0>1, the disease has the potential to persist. Moreover, the necessary conditions for the existence of the stationary distribution of positive solution in the event of disease persistence is determined. Our theoretical findings are validated through numerical simulations.

    Citation: Zhongwei Cao, Jian Zhang, Huishuang Su, Li Zu. Threshold dynamics of a stochastic general SIRS epidemic model with migration[J]. Mathematical Biosciences and Engineering, 2023, 20(6): 11212-11237. doi: 10.3934/mbe.2023497

    Related Papers:

    [1] Yanmei Wang, Guirong Liu . Dynamics analysis of a stochastic SIRS epidemic model with nonlinear incidence rate and transfer from infectious to susceptible. Mathematical Biosciences and Engineering, 2019, 16(5): 6047-6070. doi: 10.3934/mbe.2019303
    [2] Yoichi Enatsu, Yukihiko Nakata, Yoshiaki Muroya . Global stability for a class of discrete SIR epidemic models. Mathematical Biosciences and Engineering, 2010, 7(2): 347-361. doi: 10.3934/mbe.2010.7.347
    [3] Abdennasser Chekroun, Mohammed Nor Frioui, Toshikazu Kuniya, Tarik Mohammed Touaoula . Global stability of an age-structured epidemic model with general Lyapunov functional. Mathematical Biosciences and Engineering, 2019, 16(3): 1525-1553. doi: 10.3934/mbe.2019073
    [4] Baoxiang Zhang, Yongli Cai, Bingxian Wang, Weiming Wang . Dynamics and asymptotic profiles of steady states of an SIRS epidemic model in spatially heterogenous environment. Mathematical Biosciences and Engineering, 2020, 17(1): 893-909. doi: 10.3934/mbe.2020047
    [5] Tingting Xue, Xiaolin Fan, Zhiguo Chang . Dynamics of a stochastic SIRS epidemic model with standard incidence and vaccination. Mathematical Biosciences and Engineering, 2022, 19(10): 10618-10636. doi: 10.3934/mbe.2022496
    [6] Yoichi Enatsu, Yukihiko Nakata . Stability and bifurcation analysis of epidemic models with saturated incidence rates: An application to a nonmonotone incidence rate. Mathematical Biosciences and Engineering, 2014, 11(4): 785-805. doi: 10.3934/mbe.2014.11.785
    [7] Fang Zhang, Wenzhe Cui, Yanfei Dai, Yulin Zhao . Bifurcations of an SIRS epidemic model with a general saturated incidence rate. Mathematical Biosciences and Engineering, 2022, 19(11): 10710-10730. doi: 10.3934/mbe.2022501
    [8] Toshikazu Kuniya, Yoshiaki Muroya, Yoichi Enatsu . Threshold dynamics of an SIR epidemic model with hybrid of multigroup and patch structures. Mathematical Biosciences and Engineering, 2014, 11(6): 1375-1393. doi: 10.3934/mbe.2014.11.1375
    [9] Han Ma, Qimin Zhang . Threshold dynamics and optimal control on an age-structured SIRS epidemic model with vaccination. Mathematical Biosciences and Engineering, 2021, 18(6): 9474-9495. doi: 10.3934/mbe.2021465
    [10] Yilei Tang, Dongmei Xiao, Weinian Zhang, Di Zhu . Dynamics of epidemic models with asymptomatic infection and seasonal succession. Mathematical Biosciences and Engineering, 2017, 14(5&6): 1407-1424. doi: 10.3934/mbe.2017073
  • In this study, a stochastic SIRS epidemic model that features constant immigration and general incidence rate is investigated. Our findings show that the dynamical behaviors of the stochastic system can be predicted using the stochastic threshold RS0. If RS0<1, the disease will become extinct with certainty, given additional conditions. Conversely, if RS0>1, the disease has the potential to persist. Moreover, the necessary conditions for the existence of the stationary distribution of positive solution in the event of disease persistence is determined. Our theoretical findings are validated through numerical simulations.



    Many well-known epidemic models [1,2,3,4,5,6,7] have been proposed and discussed over the years. For instance, De la Sen et al. [8] in their study analyzed an epidemic model that incorporates delayed, distributed disease transmission and a general vaccination policy. Weera et al. conducted a numerical investigation of a nonlinear computer virus epidemic model with time delay effects [9]. Li et al. [3] examined an SIRS epidemic model with a general incidence rate and constant immigration, which took the following form

    {˙S=aAβf(N)SIμS+δR,˙I=bA+βf(N)SI(μ+γ+α)I,˙R=cA+γI(μ+δ)R, (1.1)

    where N==S+I+R and the biological implications are shown in Table 1, and the infectious force βf(N) is a continuous and twice differentiable function of total population and β>0 is adequate contact rate. Furthermore, f satisfies the following hypotheses

    Table 1.  Variables in model (1.1).
    Variables Biological implications
    S Numbers of susceptible individuals
    I Numbers of infectious individuals
    R Numbers of removed individuals
    N Total population
    A Rate of input to the total population
    a Fraction of input to susceptible class
    b Fraction of input to infectious class
    c Fraction of input to removed class
    μ Natural death rate
    γ Recovery rate
    α Mortality due to virulence
    δ Rate of losing immunity

     | Show Table
    DownLoad: CSV

    1)fC2((0,);(0,)).

    2)f(N)0 for any N>0.

    3)[f(N)N]0 for any N>0.

    Their research [3] found that

    R0=βf(Aμ)A(δ+(1c)μ)(μ+γ+α)(μ+δ)μ

    is the basic reproduction number. Furthermore, one gets

    ● If b=0 and R0<1, then system (1.1) has a disease-free equilibrium E0=(S0,I0,R0)=(AμcAμ+δ,0,cAμ+δ), which is globally asymptotically stable (GAS).

    ● If R0>1 and b=0, there exists a unique endemic equilibrium E=(S,I,R) which is GAS.

    ● Otherwise if b>0, there is no disease-free equilibrium in system (1.1) and there exists a unique endemic equilibrium P=(S1,I1,R1) which is locally asymptotically stable. In addition, when αμ+2δ, the endemic equilibrium P is GAS.

    However, in reality, variations in environmental factors affect the transmission coefficients of infectious diseases. As a result, stochastic modelling is an appropriate way to model epidemics in a variety of situations. For example, stochastic models can account for the randomness of infectious contacts that may occur during potential and infectious periods [10]. In comparison to deterministic models, stochastic epidemic models can provide more realism. A growing number of authors have recently focused on stochastic epidemic models [4,5,6,7,11,12,13,14,15,16,17,18,19,20,21,22]. Cai et al. [7] discovered that the global dynamics of a general SIRS epidemic model can determine the existence of either a unique stationary distribution free of disease or a unique stationary distribution with endemic disease. Liu et al. [18] found that in a stochastic SIRS epidemic model with standard incidence, in which two threshold parameters RS0 and ^RS0 exist.

    Inspired by Mao et al. [23], this paper posits that fluctuations in the environment primarily manifest as fluctuations in the transmission coefficient,

    ββ+σ˙B(t),

    where B(t) is a standard Brownian motion and σ2>0 indicates its intensity. Then we have

    {dS(t)=[aAβf(N)S(t)I(t)μS(t)+δR(t)]dtσf(N)S(t)I(t)dB(t),dI(t)=[bA+βf(N)S(t)I(t)(μ+γ+α)I(t)]dt+σf(N)S(t)I(t)dB(t),dR(t)=[cA+γI(t)(μ+δ)R(t)]dt. (1.2)

    Our study is based on the deterministic SIRS epidemic model, which has proven to be an effective tool for investigating the spread of infectious diseases. Our approach incorporates two crucial elements: constant immigration and a general incidence rate, which are essential for understanding the impact of environmental fluctuations on disease dynamics.

    One of the main strengths of our study lies in the fact that we have integrated these essential components into a stochastic framework. This has enabled us to analyze the effects of random fluctuations in disease transmission and immigration rates, which are significant factors that can profoundly influence the dynamics of infectious diseases. By examining these effects, we can obtain a more comprehensive understanding of the factors that contribute to the spread and persistence of diseases. Furthermore, our research has established the necessary conditions for the existence of a stationary distribution of positive solutions in the case of disease persistence. This novel contribution to the field has significant implications for the development of effective strategies for managing and controlling infectious diseases. Ultimately, our study provides valuable insights that can inform public health policies and initiatives aimed at reducing the impact of infectious diseases on global health.

    The purpose of this paper is to explore the impact of environmental fluctuations on disease dynamics by analyzing the global dynamics of the stochastic SIRS epidemic model (1.2). The paper is structured as follows: In Section 2, we provide some preliminaries. Section 3 outlines the necessary conditions for disease extinction and persistence. We determine sufficient conditions for the existence of stationary distributions for persistent solutions of the model in Section 4. The paper concludes with numerical simulations and conclusions.

    In this paper, unless specified otherwise, let (Ω,F,{Ft}t0,P) denote a complete probability space with a filtration {Ft}t0 that satisfies the regular conditions. Let B(t) be defined on this complete probability space. Denote ab=max{a,b}for anya,bR, and X={(x1,x2,x3)R3:x1>0,x2>0,x3>0}.

    Lemma 1. [24] (Strong Law of Large Numbers) Let M={M}t0 be a real-valued continuous local martingale vanishing at t=0. Then

    limtM,Mt=,a.s.limtMtM,Mt=0a.s.,

    and

    lim suptM,Mtt<a.s.limtMtt=0a.s.

    Theorem 1. For any (S(0),I(0),R(0))X, there is a unique solution (S(t),I(t),R(t)) of system (1.2) that remain in X with probability one.

    The proof is standard and hence is omitted here.

    Remark 1. From Theorem 2.1, we have

    [A(α+μ)N]dtdN[AμN]dt,dR[cA(μ+δ)R]dt,t[0,),a.s.

    This implies that

    Γ={(S,I,R)X:Aα+μ<N<Aμ,R>cAμ+δ}

    is a positively invariant set of system (1.2). Hence throughout this paper we always assume that the initial value (S(0),I(0),R(0))Γ.

    In contrast to the deterministic system (1.1), the purpose of this section is to study the dynamics of the system (1.2) when b=0 holds. Denote

    Rs0:=βf(Aμ)A(δ+(1c)μ)(μ+γ+α)(μ+δ)μσ2f2(Aμ)A22(μ+γ+α)μ2=R0σ2f2(Aμ)A22(μ+γ+α)μ2.

    Theorem 2. Let b=0 and (S(t),I(t),R(t)) be a solution of system (1.2). If

    σ2>max{β22(μ+γ+α),βμf(Aμ)A} (3.1)

    or

    Rs0<1andσ2<βμf(Aμ)A, (3.2)

    then

    lim suptlnI(t)t<0,limtS(t)=AμcAμ+δ,limtR(t)=cAμ+δa.s.

    Proof. Making the use of Itô's formula [24] to lnI, we have

    dlnI=(βf(N)S(μ+γ+α)σ22f2(N)S2)dt+σf(N)SdB(t).

    Integrating the above equality from 0 to t and then dividing by t on both sides, one obtains

    lnI(t)lnI(0)t=t0ϕ(τ)dτt+G(t)t, (3.3)

    where

    ϕ(τ)=βf(N(τ))S(τ)(μ+γ+α)σ22f2(N(τ))S2(τ),G(t)=t0σf(N(τ))S(τ)dB(τ).

    Noting that G(t) is a local martingale (since it is a right continuous adapted process defined on (Ω,F,{Ft}t0,P)) whose quadratic variation is

    G,Gt=t0σ2f2(N(τ))S2(τ)dτσ2f2(Aμ+α)A2μ2t.

    Making the use of Lemma 2.1 leads to limtG(t)t=0 a.s. Combining (3.1), we have

    ϕ(τ)=σ22(f(N(τ))S(τ)βσ2)2+β22σ2(μ+γ+α)β22σ2(μ+γ+α).

    Substituting the above inequality into (3.3) and taking the limit on both sides, we obtain

    limtlnI(t)tβ22σ2(μ+γ+α)<0a.s. (3.4)

    Consider the case σ2<βμf(Aμ)A, we get

    ϕ(τ)=βf(N(τ))N(τ)S(τ)N(τ)(μ+γ+α)σ22f2(N(τ))S2(τ)βf(Aμ)Aμ(1R(τ)N(τ))σ22f2(Aμ)A2μ2(μ+γ+α). (3.5)

    Noting that Aα+μ<N<Aμ, R>cAμ+δ and substituting them into (3.5), we have

    ϕ(τ)βf(Aμ)A(δ+(1c)μ)μ(μ+δ)σ22f2(Aμ)A2μ2(μ+γ+α)=(Rs01)(μ+γ+α).

    From (3.2) and (3.3), we get

    limtlnI(t)t(μ+γ+α)(Rs01)<0a.s. (3.6)

    Then we have

    limtI(t)=0,a.s., (3.7)

    which means that for arbitrary small ε>0 there are t0 and Ωε such that P(Ωε)1ε and αIε for tt0 and ωΩε. In view of system (1.2), we have

    AεμlimtN(t)Aμa.s.

    Due to the arbitrariness of ε, one has

    limtN(t)=Aμa.s. (3.8)

    Similarly as getting equality (3.8), we have

    limtR(t)=cAμ+δa.s. (3.9)

    In view of (3.7)–(3.9), we have

    limtS(t)=AμcAμ+δa.s.

    Remark 2. According to Theorem 3.1, if Rs0<1 and σ is not large, the disease will inevitably die out. It is worth noting that the expressions Rs0 and R0 reveal that Rs0<R0. Furthermore, if σ=0, Rs0=R0. In simpler terms, the conditions for the disease to die out in system (1.2) are considerably easier than those in the corresponding deterministic system (1.1).

    In this section, we will prove that if b=0 and Rs0>1 or b>0, the densities of the distributions of the solutions to system (1.2) can converge in L1 to an invariant density.

    Theorem 3. The distribution of (S(t),I(t),R(t)) has a density U(t,x,y,z) for t>0. If b=0 and Rs0>1 or b>0, then there is a unique density U(x,y,z) such that

    limtΓ|U(t,x,y,z)U(x,y,z)|dxdydz=0.

    The following steps constitute the proof of Theorem 4.1 above:

    ● First, the kernel function of (S(t),I(t),R(t)) is absolutely continuous.

    ● We demonstrate that the kernel function is positive on X.

    ● The Markov semigroup is either sweeping with respect to compact sets or asymptotically stable.

    ● Due to the presence of Khasminskiˇi function, we exclude sweeping.

    For definitions related to Markov semigroups and their asymptotic properties, the reader is referred to the papers [25,26,27,28,29,30,31]. We will show this by Lemmas 4.1–4.5.

    Lemma 2. For t>0 and any initial value (x0,y0,z0)X, the transition probability function P(t,x0,y0,z0,B) has a continuous density k(t,x,y,z;x0,y0,z0).

    Proof. Similar to the proof method in [31], the Lie bracket is given by

    [a,b]j(u)=3i=1(aibjui(u)biajui(u)),j=1,2,3.

    Let a0(S,I,R)=(aAβf(N)SIμS+δRbA+βf(N)SI(μ+γ+α)IcA+γI(μ+δ)R) and a1(S,I,R)=(σf(N)SIσf(N)SI0). Direct calculation leads to

    a2=[a0,a1]=(a21a22σγf(N)SI),

    with

    a21=(AμNαI)σSIf(N)σf(N)(βf(N)S2I2+(aA+δR)I+(bA(μ+γ+α)I)S),a22=σf(N)SI(AμNαI)+σf(N)(I(aA(2μ+γ+α)S+δR)+bAS),

    and

    a3=[a1,a2]=(a31a32σ2γSIf2(N)(IS)),

    where

    a31=σ2βS2I2f2(N)f(N)(I2S)+σ2f2(N)(μS2I(aA+δR)I2+bAS2)σ2f(N)f(N)SI(βf(N)S2I2+μSI),a32=σ2f2(N)(μSI+βf(N)S2I3+(aA+δR)I2+μS2IbAS2)+σ2f(N)f(N)SI(βf(N)S2I2+μSI).

    Therefore, we have

    |a1a2a3|=σ{σ2γSIf2(N)(IS)(a21+a22)+σγf(N)SI(a31+a32)}<0.

    According to H¨ormander Theorem [23], one obtains that P(t,x0,y0,z0,B) has a continuous density k(t,x,y,z;x0,y0,z0).

    Next, fixing a point (x0,y0,z0)X and a function ψL2([0,T];R), we have

    {xψ(t)=x0+t0(f1(xψ(τ),yψ(τ),zψ(τ))σψxψ(τ)yψ(τ)f(Nψ(τ)))dτ,yψ(t)=y0+t0(f2(xψ(τ),yψ(τ),zψ(τ))σψxψ(τ)yψ(τ)f(Nψ(τ)))dτ,zψ(t)=z0+t0f3(xψ(τ),yψ(τ),zψ(τ))dτ, (4.1)

    where

    Nψ=xψ+yψ+zψ,f1=aAβf(x+y+z)xyμx+δz,f2=bA+βf(x+y+z)xy(μ+γ+α)y,f3=cA+γy(μ+δ)z.

    Let DX0;ψ be the Freˊchet derivative. If for some \psi\in L^2([0, T]; \mathbb{R}) , the rank of D_{X_0;\psi} is 3, then k(T, x, y, x; x_0, y_0, z_0) > 0 for X = X_{\psi}(T) . Let

    \Psi (t) = f'(X_{\psi }(t))+\psi g'(X_{\psi }(t)),

    where f' and g' are the Jacobians of

    f = \begin{pmatrix} f_1 \\ f_2 \\ f_3 \end{pmatrix}, \quad g = \begin{pmatrix} -\sigma xyf(x+y+z)\\ \sigma xyf(x+y+z)\\ 0 \end{pmatrix}.

    For T\geq t\geq t_0\geq 0 , let Q(t, t_0) be a matrix function such that Q(t_0, t_0) = Id , \frac{\partial Q(t, t_0)}{\partial t} = \Psi (t)Q(t, t_0). Then

    D_{X_0;\psi }h = \int_{0}^{T}Q(T, \tau )g(\tau )h(\tau )d\tau .

    Lemma 3. For each (x_0, y_0, z_0), (x, y, z)\in \Gamma , there is T > 0 satistying k(T, x, y, z; x_0, y_0, z_0) > 0 .

    Proof. Since we only need to find a continuous control function \psi , system (4.1) can be rewritten as follows

    \begin{align} \left\{\begin{matrix} x'_{\psi }(t) = f_1(x_{\psi }(t), y_{\psi }(t), z_{\psi }(t))-\sigma \psi x_{\psi }(t)y_{\psi }(t)f(N_{\psi }(t)) , \\ y'_{\psi }(t) = f_2(x_{\psi }(t), y_{\psi }(t), z_{\psi }(t))-\sigma \psi x_{\psi }(t)y_{\psi }(t)f(N_{\psi }(t)) , \\ z'_{\psi }(t) = f_3(x_{\psi }(t), y_{\psi }(t), z_{\psi }(t)) , \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \end{matrix}\right. \end{align} (4.2)

    First, we verify that the rank of D_{X_0;\psi} is 3. Let \varepsilon \in (0, T) and

    h(t) = \frac{\chi _{[T-\varepsilon , T]}}{x_{\psi}(t)y_{\psi}(t)f(N_{\psi}(t))} , \; t\in [0, T] ,

    where \chi denotes the indicator function of the interval [T-\varepsilon, T] . Since

    Q(T, \tau ) = Id+\Psi (T)(\tau -T)+\frac{1}{2}\Psi ^2(T)(\tau -T)^2+\circ ((\tau -T)^2),

    we have

    D_{X_0;\psi }h = \varepsilon v-\frac{\varepsilon ^2}{2}\Psi (T)v+\frac{\varepsilon ^3}{6}\Psi ^2(T)v+\circ (\varepsilon ^3),

    where v = \begin{pmatrix} -\sigma \\ \sigma \\ 0 \end{pmatrix} . Direct calculation leads to

    \Psi (T)v = \sigma \begin{pmatrix} (\beta +\psi \sigma )f(N)(I-S)+\mu \\ (\beta +\psi \sigma )f(N)(I-S)-(\mu +\gamma +\alpha )\\ \gamma \end{pmatrix},
    \Psi ^2(T)v = \sigma \begin{pmatrix} c_{11} \\ c_{21}\\ \sigma \gamma (\beta +\psi \sigma )(S-I)f(N)-\sigma \gamma (2\mu +\gamma +\alpha +\delta ) \end{pmatrix},

    where

    \begin{aligned} c_{11} = &-\sigma (S-I)^2(\beta +\psi \sigma )^2f^2(N)+\sigma (\beta +\psi \sigma )(2\mu (S-I)+(\alpha +\gamma )S)f(N)\\ &+\sigma (\gamma \delta -\mu ^2)+\sigma (\beta +\psi \sigma )(\alpha +\gamma )SIf'(N), \\ c_{21} = &\sigma (S-I)^2(\beta +\psi \sigma )^2f^2(N)+\sigma (\beta +\psi \sigma )(-2(\alpha S-\mu I)+(\alpha +\gamma )I-2S(\mu +\gamma ))f(N) \\ &-\sigma (\beta +\psi \sigma )(\alpha +\gamma )SIf'(N)+\sigma (\mu +\gamma +\alpha )^2. \end{aligned}

    Thus the rank of D_{X_0;\psi } is 3.

    Then let w_{\psi} = x_{\psi}+y_{\psi}+z_{\psi} , (4.2) becomes

    \begin{align} \left\{\begin{aligned} x'_{\psi }(\tau) = &g_1(x_{\psi }(\tau), w_{\psi }(\tau), z_{\psi }(\tau))-\sigma \psi x_{\psi }(\tau)(w_{\psi }(\tau)-x_{\psi }(\tau)-z_{\psi }(\tau))f(w_{\psi }(\tau)) , \\ w'_{\psi }(\tau) = &g_2(x_{\psi }(\tau), w_{\psi }(\tau), z_{\psi }(\tau)), \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \\ z'_{\psi }(\tau) = &g_3(x_{\psi }(\tau), w_{\psi }(\tau), z_{\psi }(\tau)), \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \end{aligned}\right. \end{align} (4.3)

    where

    \begin{align} \begin{matrix} g_1 = aA-\beta f(w)x(w-x-z)-\mu x+\delta z, \\ g_2 = bA-(\mu +\alpha )w+\alpha (x+z), \\ g_3 = cA+\gamma (w-x)-(\gamma +\mu +\delta )z. \end{matrix} \end{align} (4.4)

    Let

    \begin{align} \Gamma _0 = \left \{ (x, w, z)\in \mathbb{X}:0 < x < w, \frac{cA}{\mu +\delta } < z < w\; {\rm and}\; \frac{A}{\alpha +\mu } < w < \frac{A}{\mu } \right \}. \end{align} (4.5)

    First, we find a positive constant T and a differentiable function

    w_\psi:[0, T]\to \left(\frac{A}{\alpha +\mu }, \frac{A}{\mu }\right)

    such that w_\psi(0) = w_0 , w_\psi(T) = w_1 , w_\psi'(0) = g_2(x_0, w_0, z_0) = w_0^d , w_\psi'(T) = g_2(x_1, w_1, z_1) = w_T^d and

    A-(\alpha+\mu) w_\psi(t) < w_\psi'(t) < A-\mu w_\psi(t), \quad t\in[0, T].

    We split the construction of the function w_\psi on three intervals [0, \tau] , [\tau, T-\tau] and [T-\tau, T] , where 0 < \tau < T /2 . Let

    \xi = \frac{1}{2}\min\left\{w_0-\frac{A}{\alpha +\mu }, w_1-\frac{A}{\alpha +\mu }, \frac{A}{\mu }-w_0, \frac{A}{\mu }-w_1 \right\}.

    If w_\psi \in \left(\frac{A}{\alpha +\mu }+\theta, \frac{A}{\mu }-\theta\right) , we have

    A-(\alpha +\mu)w_\psi(t) < -(\alpha +\mu)\theta < 0, \; 0 < \mu \theta < A-\mu w_\psi(t), \quad t\in[0, T].

    Then we construct a C^2 -function w_\psi : [0, \tau] \to \left(\frac{A}{\alpha +\mu }+\theta, \frac{A}{\mu }-\theta\right) such that

    w_\psi(0) = w_0, \; w_\psi'(0) = w_0^d, \; w_\psi'(\tau) = 0,

    and for t\in[0, \tau] , w_\psi satisfies

    A-(\alpha+\mu) w_\psi(t) < w_\psi'(t) < A-\mu w_\psi(t) .

    Analogously, we can construct a C^2 -function w_\psi : [T-\tau, T] \to \left(\frac{A}{\alpha +\mu }+\theta, \frac{A}{\mu }-\theta\right) such that

    w_\psi(T) = w_1, \; w_\psi'(T) = w_T^d, \; w_\psi'(T-\tau) = 0,

    and for t\in[T-\tau, T] , w_\psi satisfies

    A-(\alpha+\mu) w_\psi(t) < w_\psi'(t) < A-\mu w_\psi(t) .

    Taking T sufficiently large, then we can extend the function w_\psi : [0, \tau]\bigcup [T-\tau, T] \to \left(\frac{A}{\alpha +\mu }+\theta, \frac{A}{\mu }-\theta\right) to a C^2 -function w_\psi defined on the whole interval [0, T] such that

    A-(\alpha+\mu) w_\psi(t) < w_\psi'(t) < A-\mu w_\psi(t).

    Thus, we can find C^1 -functions x_\psi and z_\psi that satisfy (4.3). Finally we can determine a continuous function \psi . and T > 0 such that x_{\psi}(0) = x_0 , w_{\psi}(0) = w_0 , z_{\psi}(0) = z_0 , x_{\psi}(T) = x , w_{\psi}(T) = w , z_{\psi}(T) = z . This completes the proof.

    Lemma 4. If b = 0 and R_0^{s} > 1 or b > 0 . For \left \{ P(t) \right \}_{t\geq 0} and every density g , one has

    \lim\limits_{t\rightarrow \infty }{\iiint}_{\Gamma }P(t)g(x, y, z)dxdydz = 1.

    Proof. System (1.2) can be rewriten as

    \begin{align} \left\{\begin{matrix} dS = g_1(S, N, R)dt-\sigma S(N-S-R)f(N)dB(t), \\ dN = g_2(S, N, R)dt, \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \\ dR = g_3(S, N, R)dt.\; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \end{matrix}\right. \end{align} (4.6)

    From Remark 2.1, we get

    \begin{align} A-(\alpha +\mu )N < \frac{dN}{dt} < A-\mu N\; {\rm and}\; \frac{dR}{dt} > cA-(\mu +\delta )R\; {\rm for}\; t\in (0, \infty )\; {\rm a.s.} \end{align} (4.7)

    For almost every w\in \Omega , there is t_0\in t_0(w) such that

    \frac{A}{\alpha +\mu } < N(t, w) < \frac{A}{\mu }\; {\rm and}\; R(t, w) > \frac{cA}{\mu +\delta }\; {\rm for}\; t > t_0.

    As a matter of fact, there exist three possible cases:

    1) N(0, w)\in \left (\frac{A}{\alpha +\mu }, \frac{A}{\mu }\right) . In this case, our statement is obvious from (4.7).

    2) N(0, w)\in \left (0, \frac{A}{\alpha +\mu }\right) . Assume that our claim is not satisfied. Then there is \Omega '\subset \Omega with \mathbb{P}(\Omega') > 0 such that N(t, w)\in (0, \frac{A}{\alpha +\mu }), w\in \Omega' . By (4.7), we obtain that for any w\in \Omega' , N(t, w) is strictly increasing on [0, \infty) and

    \lim\limits_{t\rightarrow \infty }N(t, w) = \frac{A}{\alpha +\mu }, \; w\in \Omega '.

    According to system (4.6), we get that \lim_{t\rightarrow \infty }S(t, w) = \lim_{t\rightarrow \infty }R(t, w) = 0 , w\in \Omega ' and thus, \lim_{t\rightarrow \infty }I(t, w) = \frac{A}{\alpha +\mu} , w\in \Omega ' .

    Consider the case b = 0 , making the use of It \hat{o} 's formula, we have

    d\ln I = \left ( \beta f(N)S-(\mu +\gamma +\alpha )-\frac{\sigma ^2}{2}S^2f^2(N) \right )dt+\sigma Sf(N)dB(t).

    Thus

    \begin{align*} \frac{\ln I(t)-\ln I(0)}{t}& = \frac{1}{t}\int_{0}^{t}\left ( \beta f(N(\tau ))S(\tau )-(\mu +\gamma +\alpha )-\frac{\sigma ^2}{2}S^2(\tau )f^2(N(\tau )) \right )d\tau +\frac{1}{t}\int_{0}^{t}\sigma S(\tau )f(N(\tau ))dB(\tau )\\ & = \frac{1}{t}\int_{0}^{t}\left ( \beta f(N(\tau ))S(\tau )-(\mu +\gamma +\alpha )-\frac{\sigma ^2}{2}S^2(\tau )f^2(N(\tau )) \right )d\tau +\frac{G(t)}{t}, \end{align*}

    where G(t): = \frac{1}{t}\int_{0}^{t}\sigma S(\tau)f(N(\tau))dB(\tau) . Applying Lemma 2.1, we have

    \lim\limits_{t\rightarrow \infty }\frac{G(t)}{t} = 0\; \; {\rm a.s.}

    Thus, due to S(t) , I(t) , f(N(t)) are continuous,

    \lim\limits_{t\rightarrow \infty }\frac{1}{t}\int_{0}^{t}\left ( \beta f(N(\tau ))S(\tau )-(\mu +\gamma +\alpha )-\frac{\sigma ^2}{2}S^2(\tau )f^2(N(\tau )) \right )d\tau +\lim\limits_{t\rightarrow \infty }\frac{G(t)}{t} = -(\mu +\gamma +\alpha ).

    This contradicts the assumption

    \lim\limits_{t\rightarrow \infty }\frac{\ln I(t)-\ln I(0)}{t} = 0\; a.s.

    Then let us consider the case b > 0 . Since \lim_{t\rightarrow \infty }N(t, w) = \frac{A}{\alpha+\mu} and \lim_{t\rightarrow \infty }S(t, w) = \lim_{t\rightarrow \infty }R(t, w) = 0 for w\in \Omega' , which contradicts that R(t, w) > 0 for w\in \Omega' , t\in (0, \infty) and the claim follows.

    3) N(0, w)\in (\frac{A}{\mu}, \infty) . We suppose, by contradiction, and analogous arguments to 2), that there is \Omega'\subset \Omega with \mathbb{P}(\Omega') > 0 such that

    \lim\limits_{t\rightarrow \infty }N(t, w) = \frac{A}{\mu}, w\in \Omega'.

    Firstly, consider the case b = 0 , by the second and third equations of (4.6), for any w\in \Omega' , one gets

    N(t, w) = e^{-(\mu +\alpha )t}\left ( N(0, w)+\int_{0}^{t}e^{(\mu +\alpha )\tau }(A+\alpha (S(\tau , w)+R(\tau , w)))d\tau \right ),
    R(t, w) = e^{-(\mu +\delta )t}\left ( R(0, w)+\int_{0}^{t}e^{(\mu +\delta )\tau }(cA+\gamma I(\tau , w))d\tau \right ).

    For all w\in \Omega' , one has

    \lim\limits_{t\rightarrow \infty }S(t, w) = \frac{A}{\mu }-\frac{cA}{\mu +\delta }, \; \; \lim\limits_{t\rightarrow \infty }I(t, w) = 0, \; \; \lim\limits_{t\rightarrow \infty }R(t, w) = \frac{cA}{\mu +\delta }\; \; {\rm a.s.}

    Therefore

    \begin{align*} \lim\limits_{t\rightarrow \infty }\frac{\ln I(t)-\ln I(0)}{t}& = \lim\limits_{t\rightarrow \infty }\frac{1}{t}\int_{0}^{t}\left ( \beta f(N(\tau ))S(\tau )-(\mu +\gamma +\alpha )-\frac{\sigma ^2}{2}S^2(\tau )f^2(N(\tau )) \right )d\tau +\lim\limits_{t\rightarrow \infty }\frac{G(t)}{t}\\ & = \lim\limits_{t\rightarrow \infty }\frac{1}{t}\int_{0}^{t}\left ( \beta f(N(\tau ))S(\tau )-(\mu +\gamma +\alpha )-\frac{\sigma ^2}{2}S^2(\tau )f^2(N(\tau )) \right )d\tau \\ & = \beta f(\frac{A}{\mu })\frac{A(\delta +(1-c)\mu )}{\mu (\mu +\delta )}-(\mu +\gamma +\alpha )-\frac{\sigma ^2}{2}\left ( \frac{A}{\mu}-\frac{cA}{\mu +\delta } \right )^2f^2(\frac{A}{\mu })\\ & > \beta f(\frac{A}{\mu })\frac{A(\delta +(1-c)\mu )}{\mu (\mu +\delta )}-(\mu +\gamma +\alpha )-\frac{\sigma ^2}{2}\frac{A^2}{\mu ^2}f^2(\frac{A}{\mu })\\ & = (\mu +\gamma +\alpha )(R_0^s-1)\\ & > 0\; \; {\rm a.s.\; on\; } \Omega'. \end{align*}

    This contradicts the assumption \lim_{t\rightarrow \infty }I(t) = 0 a.s. In other words, for almost all w\in \Omega , there is t_0 = t_0(w) such that

    \frac{A}{\alpha +\mu } < N(t, w) < \frac{A}{\mu }\; {\rm for\; }t > t_0.

    When b > 0 , we get that I(t, w) > 0 for t\in (0, \infty) and w\in \Omega' . This contradicts the assumption \lim_{t\rightarrow \infty }N(t, w) = \frac{A}{\mu} , w\in \Omega ' and the claim holds.

    Similar to the proof of 2) and 3), one obtains that for almost all w\in \Omega , there is t_1 = t_1(w) such that

    R(t, w) > \frac{cA}{\mu +\delta }\; {\rm for\; }t > t_1.

    Lemma 5. \left \{ P(t) \right \}_{t\geq 0} is asymptotically stable or is sweeping with respect to compact sets.

    Proof. In view of Lemma 4.1, \left \{ P(t) \right \}_{t\geq 0} is an integral Markov semigroup with kernel k(t, x, y, z; x_0, y_0, z_0) . According to Lemma 4.3, it suffices to consider the restriction of \left \{ P(t) \right \}_{t\geq 0} to the space L^1(\Gamma) . By Lemma 4.2, one gets

    \int_{0}^{\infty }P(t)gdt > 0\; a.s.

    on \Gamma , for every g\in \mathbb{D} . Then \left \{ P(t) \right \}_{t\geq 0} is asymptotically stable or is sweeping with respect to compact sets.

    Lemma 6. Assume that b = 0 and R_0^s > 1 or b > 0 , then \{P(t)\}_{t\geq0} is asymptotically stable.

    Proof. From Lemma 4.4, \{P(t)\}_{t\geq0} satisfies the Foguel alternative. In order to exclude sweeping it is sufficient to construct a nonnegative C^2 -Khasminskiǐ function V and a closed set D_\epsilon\in\Sigma such that

    \sup\limits_{(S, I, R)\in\mathbb{X}\setminus D_\epsilon}\mathscr{A}^{*}V(S, I, R) < 0.

    First of all, we consider the case b = 0 and R_0^s > 1 . Define

    \begin{aligned} H = &M(-\ln I-\ell_1N+\ell_2R)-\ln S-\ln (\frac{A}{\mu}-N )-\ln (N-\frac{A}{\mu+\alpha} )-\ln (R-\frac{cA}{\mu+\delta} )\\ : = &MV_1+V_2+V_3+V_4+V_5, \end{aligned}

    where V_1 = -\ln I-\ell_1N+\ell_2R , V_2 = -\ln S , V_3 = -\ln(\frac{A}{\mu}-N) , V_4 = -\ln(N-\frac{A}{\mu+\alpha}) , V_5 = -\ln(R-\frac{cA}{\mu+\delta}) , \ell_1 = \frac{\beta f(\frac{A}{\mu})}{\mu} , \ell_2 = \frac{\beta f(\frac{A}{\mu})}{\mu+\delta} and M is a positive constant satisfying

    \begin{align} \begin{aligned} \frac{\beta A}{\mu}f (\frac{A}{\alpha+\mu} )+4\mu+\alpha+\delta+\frac{\sigma^2}{2}\frac{A^2}{\mu^2}f^2 (\frac{A}{\alpha+\mu} )\leq M(\mu+\gamma+\alpha)(R_0^s-1)-2. \end{aligned} \end{align} (4.8)

    It is easy to find that H reaches a minimum at (S_*, I_*, R_*) . Then we define

    \begin{aligned} V = MV_1+V_2+V_3+V_4+V_5-H(S_*, I_*, R_*). \end{aligned}

    Thus we have

    \begin{aligned} \mathscr{A}^{*}V_1 = &-\beta Sf(N)+(\mu+\gamma+\alpha)+\frac{\sigma^2}{2}S^2f^2(N)-\ell_1A+\ell_1\mu N+\ell_1\alpha I+\ell_2cA+\ell_2\gamma I-\ell_2(\mu+\delta)R\\ \leq&-\beta f (\frac{A}{\mu} )S+\ell_1\mu S+\ell_1\mu R-\ell_2(\mu+\delta)R+(\mu+\gamma+\alpha)+\frac{\sigma^2}{2}\frac{A^2}{\mu^2}f^2 (\frac{A}{\mu} )\\ &-\ell_1A+\ell_2cA+(\ell_1\mu+\ell_1\alpha+\ell_2\gamma)I\\ = &(\mu+\gamma+\alpha)+\frac{\sigma^2}{2}\frac{A^2}{\mu^2}f^2 (\frac{A}{\mu} )-\ell_1A+\ell_2cA+(\ell_1\mu+\ell_1\alpha+\ell_2\gamma)I\\ = &-\beta f (\frac{A}{\mu} )\frac{A(\delta+(1-c)\mu)}{\mu(\mu+\delta)}+(\mu+\gamma+\alpha)+\frac{\sigma^2}{2}\frac{A^2}{\mu^2}f^2 (\frac{A}{\mu} )+(\ell_1\mu+\ell_1\alpha+\ell_2 \gamma)I\\ = &-(\mu+\gamma+\alpha)(R_0^s-1)+(\ell_1\mu+\ell_1\alpha+\ell_2\gamma)I. \end{aligned}

    Similarly, we obtain

    \begin{aligned} \mathscr{A}^{*}V_2 = &- (\frac{aA}{S}-\beta If(N)-\mu+\frac{\delta R}{S}-\frac{\sigma^2}{2}I^2f^2(N) )\\ \leq&-\frac{aA}{S}+\beta\frac{A}{\mu}f (\frac{A}{\alpha+\mu} )+\mu+\frac{\sigma^2}{2}\frac{A^2}{\mu^2}f^2 (\frac{A}{\alpha+\mu} ), \end{aligned}
    \begin{aligned} \mathscr{A}^{*}V_3 = \mu-\frac{\alpha I}{\frac{A}{\mu}-N}, \end{aligned}
    \begin{aligned} \mathscr{A}^{*}V_4 = &-\frac{A-\mu N-\alpha I}{N-\frac{A}{\mu+\alpha}}\leq\mu+\alpha-\frac{\alpha I}{N-\frac{A}{\mu+\alpha}} \end{aligned}

    and

    \begin{aligned} \mathscr{A}^{*}V_5 = &-\frac{cA+\gamma I-(\mu+\delta)R}{R-\frac{cA}{\mu+\delta}} = \mu+\delta-\frac{\gamma I}{R-\frac{cA}{\mu+\delta}}. \end{aligned}

    Therefore

    \begin{aligned} \mathscr{A}^{*}V\leq&-M(\mu+\gamma+\alpha)(R_0^s-1)+M(\ell_1\mu+\ell_1\alpha+\ell_2\gamma)I-\frac{aA}{S}-\frac{\alpha I}{\frac{A}{\mu}-N}-\frac{\alpha I}{N-\frac{A}{\mu+\alpha}}-\frac{\gamma I}{R-\frac{cA}{\mu+\delta}} \\ &+\frac{\beta A}{\mu}f (\frac{A}{\alpha+\mu} )+4\mu+\alpha+\delta+\frac{\sigma^2}{2}\frac{A^2}{\mu^2}f^2 (\frac{A}{\alpha+\mu} ). \end{aligned}

    Define

    D_\epsilon = \{(S, I, R)\in\Gamma:\epsilon\leq S, \epsilon\leq I, \frac{cA}{\mu+\delta}+\epsilon^2\leq R, \frac{A}{\mu+\alpha}+\epsilon^2\leq N\leq\frac{A}{\mu}-\epsilon^2 \},

    where \epsilon\in(0, 1) is sufficiently small satisfying

    \begin{align} -\frac{aA}{\epsilon}+M(\ell_1\mu+\ell_1\alpha+\ell_2\gamma)\frac{A}{\mu}+\beta\frac{A}{\mu}f (\frac{A}{\alpha+\mu} )+4\mu+\alpha+\delta+\frac{\sigma^2}{2}\frac{A^2}{\mu^2}f^ 2 (\frac{A}{\alpha+\mu} ) < -1, \end{align} (4.9)
    \begin{align} \epsilon < \frac{1}{M(\ell_1\mu+\ell_1\alpha+\ell_2\gamma)}, \end{align} (4.10)
    \begin{align} -\frac{\gamma}{\epsilon}+M(\ell_1\mu+\ell_1\alpha+\ell_2\gamma)\frac{A}{\mu}+\beta\frac{A}{\mu}f (\frac{A}{\alpha+\mu} )+4\mu+\alpha+\delta+\frac{\sigma^2}{2}\frac{A^2}{\mu^ 2}f^2 (\frac{A}{\alpha+\mu} ) < -1. \end{align} (4.11)
    \begin{align} -\frac{\alpha}{\epsilon}+M(\ell_1\mu+\ell_1\alpha+\ell_2\gamma)\frac{A}{\mu}+\beta\frac{A}{\mu}f (\frac{A}{\alpha+\mu} )+4\mu+\alpha+\delta+\frac{\sigma^2}{2}\frac{A^2}{\mu^ 2}f^2 (\frac{A}{\alpha+\mu} ) < -1. \end{align} (4.12)

    Denote

    D_1 = \{(S, I, R)\in\Gamma:S < \epsilon\}, \; D_2 = \{(S, I, R)\in\Gamma:I < \epsilon\}, \; D_3 = \{(S, I, R)\in\Gamma:I\geq\epsilon, R < \frac{cA}{\mu+\delta}+\epsilon^2 \},
    D_4 = \{(S, I, R)\in\Gamma:I\geq\epsilon, \frac{A}{\mu}-\epsilon^2 < N \}, \; D_5 = \{(S, I, R)\in\Gamma:I\geq\epsilon, N < \frac{A}{\mu+\alpha}+\epsilon^2 \}.

    Then we prove that \mathscr{A}^{*}V(S, I, R) < -1 for any (S, I, R)\in\Gamma\setminus D_\epsilon = D_1\bigcup D_2\bigcup D_3\bigcup D_4\bigcup D_5 .

    Case 1. For any (S, I, R)\in D_{1} , from (4.9),

    \begin{aligned} \mathscr{A}^{*}V\leq& -\frac{aA}{S}+M(\ell_1\mu+\ell_1\alpha+\ell_2\gamma)\frac{A}{\mu}+\beta\frac{A}{\mu}f (\frac{A}{\alpha+\mu} )+4\mu+\alpha+\delta+\frac{\sigma^2}{2}\frac{A^2}{\mu^2}f^2 (\frac{A}{\alpha+\mu} )\\ < &-\frac{aA}{\epsilon}+M(\ell_1\mu+\ell_1\alpha+\ell_2\gamma)\frac{A}{\mu}+\beta\frac{A}{\mu}f (\frac{A}{\alpha+\mu} )+4\mu+\alpha+\delta+ \frac{\sigma^2}{2}\frac{A^2}{\mu^2}f^2 (\frac{A}{\alpha+\mu} )\\ < &-1. \end{aligned}

    Thus

    \mathscr{A}^{*}V < -1\; \text{for any}\; (S, I, R)\in D_{1}.

    Case 2. For any (S, I, R)\in D_{2} , from (4.8) and (4.10),

    \begin{aligned} \mathscr{A}^{*}V\leq &-M(\mu+\gamma+\alpha)(R_0^s-1)+M(\ell_1\mu+\ell_1\alpha+\ell_2\gamma)\epsilon+\beta\frac{A}{\mu}f (\frac{A}{\alpha+\mu} )+4\mu+\alpha+\delta+\frac{\sigma^2}{2}\frac{A^2}{\mu^2}f^2 (\frac{A}{\alpha+\mu} )\\ < &-2+1\\ = &-1. \end{aligned}

    Therefore

    \mathscr{A}^{*}V < -1\; \text{for any}\; (S, I, R)\in D_{2}.

    Case 3. For any (S, I, R)\in D_{3} , from (4.11),

    \begin{aligned} \mathscr{A}^{*}V\leq&-\frac{\gamma I}{R-\frac{cA}{\mu+\delta}}+M(\ell_1\mu+\ell_1\alpha+\ell_2\gamma)I+\beta\frac{A}{\mu}f (\frac{A}{\alpha+\mu} )+4\mu+\alpha+\delta+\frac{\sigma^2}{2}\frac{A^2}{\mu^2}f^2 ( \frac{A}{\alpha+\mu} )\\ < &-\frac{\gamma\epsilon}{\epsilon^{2}}+M(\ell_1\mu+\ell_1\alpha+\ell_2\gamma)\frac{A}{\mu}+\beta\frac{A}{\mu}f (\frac{A}{\alpha+\mu} )+ 4\mu+\alpha+\delta+\frac{\sigma^2}{2}\frac{A^2}{\mu^2}f^2 (\frac{A}{\alpha+\mu} )\\ = &-\frac{\gamma}{\epsilon}+M(\ell_1\mu+\ell_1\alpha+\ell_2\gamma)\frac{A}{\mu}+\beta\frac{A} {\mu}f (\frac{A}{\alpha+\mu} )+4\mu+\alpha+\delta+\frac{\sigma^2}{2}\frac{A^2}{\mu^2}f^2 (\frac{A}{\alpha+\mu} )\\ < &-1. \end{aligned}

    Hence

    \mathscr{A}^{*}V < -1\; \text{for any}\; (S, I, R)\in D_{3}.

    Case 4. For any (S, I, R)\in D_{4} , from (4.12),

    \begin{aligned} \mathscr{A}^{*}V\leq&-\frac{\alpha I}{\frac{A}{\mu}-N}+M(\ell_1\mu+\ell_1\alpha+\ell_2\gamma)I+\beta\frac{A}{\mu}f (\frac{A}{\alpha+\mu} )+4\mu+\alpha+\delta+\frac{\sigma^2}{2}\frac{A^2}{\mu^2}f^2 (\frac{A} {\alpha+\mu} )\\ \leq &-\frac{\alpha}{\epsilon}+M(\ell_1\mu+\ell_1\alpha+\ell_2\gamma)\frac{A}{\mu}+\beta\frac{A}{\mu} f (\frac{A}{\alpha+\mu} )+4\mu+\alpha+\delta+\frac{\sigma^2}{2}\frac{A^2}{\mu^2}f^2 (\frac{A}{\alpha+\mu} )\\ < &-1, \end{aligned}

    Then

    \mathscr{A}^{*}V < -1\; \text{for any}\; (S, I, R)\in D_{4}.

    Case 5. For any (S, I, R)\in D_{5} , from (4.12),

    \begin{aligned} \mathscr{A}^{*}V\leq&-\frac{\alpha I}{N-\frac{A}{\mu+\alpha}}+M(\ell_1\mu+\ell_1\alpha+\ell_2\gamma)I+\beta\frac{A}{\mu}f (\frac{A}{\alpha+\mu} )+4\mu+\alpha+\delta+\frac{\sigma^2}{2}\frac{A^2}{\mu^2}f^2 ( \frac{A}{\alpha+\mu} \\ \leq &-\frac{\alpha}{\epsilon}+M(\ell_1\mu+\ell_1\alpha+\ell_2\gamma)\frac{A}{\mu}+\beta\frac{A}{ \mu}f (\frac{A}{\alpha+\mu} )+4\mu+\alpha+\delta+\frac{\sigma^2}{2}\frac{A^2}{\mu^2}f^2 (\frac{A}{\alpha+\mu} )\\ < &-1. \end{aligned}

    Thus

    \mathscr{A}^{*}V < -1\; \text{for any}\; (S, I, R)\in D_{5}.

    In summary,

    \sup\limits_{(S, I, R)\in\Gamma\setminus D_\epsilon}\mathscr{A}^{*}V(S, I, R) < -1.

    Using similar arguments to those in [27], we can obtain that \{P(t)\}_{t\geq0} is asymptotically stable.

    Next, we consider the case b > 0 , define

    E = -\ln S-\ln I-\ln (R-\frac{cA}{\mu+\delta} )-\ln (N-\frac{A}{\mu+\alpha} )-\ln (\frac{A}{\mu}-N ).

    Obviously, E has a minimum point (S_{1*}, I_{1*}, R_{1*}) in the interior of \Gamma. Then we define

    W = -\ln S-\ln I-\ln (R-\frac{cA}{\mu+\delta} )-\ln (N-\frac{A}{\mu+\alpha} )-\ln (\frac{A}{\mu}-N ) -E(S_{1*}, I_{1*}, R_{1*}).

    Then we have

    \begin{aligned} \mathscr{A}^{*}W\leq &-\frac{aA}{S}-\frac{bA}{I}-\frac{\gamma I}{R-\frac{cA}{\mu+\delta}}-\frac{\alpha I}{\frac{A}{\mu}-N}-\frac{\alpha I}{N-\frac{A}{\mu+\alpha}}+5\mu+2\alpha+\gamma+\delta+\frac{\beta A}{\mu}f (\frac{A}{\mu+\alpha} )+\frac{\sigma^{2}A^{2}}{\mu^{2}}f^2 (\frac{A}{\mu+\alpha} ).\\ \end{aligned}

    Similarly, define

    U_{\epsilon_1} = \{(S, I, R)\in\Gamma:\epsilon_1\leq S, \epsilon_1\leq I, \frac{cA}{\mu+\delta}+\epsilon_1^2\leq R, \frac{A}{\mu+\alpha}+\epsilon_1^2\leq N\leq\frac{A}{\mu}-\epsilon_1^2 \},

    where \epsilon_1\in (0, 1) is sufficiently small satisfying

    -\frac{aA}{\epsilon_1}+5\mu+2\alpha+\gamma+\delta+\frac{\beta A}{\mu}f (\frac{A}{\mu+\alpha} )+\frac{\sigma^{2}A^{2}}{\mu^{2}}f^2 (\frac{A}{\mu+\alpha} ) < -1,
    -\frac{bA}{\epsilon_1}+5\mu+2\alpha+\gamma+\delta+\frac{\beta A}{\mu}f (\frac{A}{\mu+\alpha} )+\frac{\sigma^{2}A^{2}}{\mu^{2}}f^2 (\frac{A}{\mu+\alpha} ) < -1,
    -\frac{\gamma}{\epsilon_1}+5\mu+2\alpha+\gamma+\delta+\frac{\beta A}{\mu}f (\frac{A}{\mu+\alpha} )+\frac{\sigma^{2}A^{2}}{\mu^{2}}f^2 (\frac{A}{\mu+\alpha} ) < -1,
    -\frac{\alpha}{\epsilon_1}+5\mu+2\alpha+\gamma+\delta+\frac{\beta A}{\mu}f (\frac{A}{\mu+\alpha} )+\frac{\sigma^{2}A^{2}}{\mu^{2}}f^2 (\frac{A}{\mu+\alpha} ) < -1.

    For convenience, we divide \Gamma\setminus U_{\epsilon_1} as

    U_1 = \{(S, I, R)\in\Gamma:S < \epsilon_{1}\}, \; U_2 = \{(S, I, R)\in\Gamma:I < \epsilon_{1}\}, \; U_3 = \{(S, I, R)\in\Gamma:I\geq\epsilon_{1}, R < \frac{cA}{\mu+\delta}+\epsilon_1^2 \},
    U_4 = \{(S, I, R)\in\Gamma:I\geq\epsilon_{1}, \frac{A}{\mu}-\epsilon_1^2 < N \}, \; U_5 = \{(S, I, R)\in\Gamma:I\geq\epsilon_{1}, N < \frac{A}{\mu+\alpha}+\epsilon_1^2 \}.

    The rest of the proof is omitted here due to it is similar to the case of b = 0 . This completes the proof.

    Remark 4.1. The stationary distribution of the correct solution refers to the long-term behavior of a stochastic system when the probability of the disease persisting is not zero. In other words, if the random threshold R_0^s is greater than 1, the disease may not be eradicated and will persist in the population. In this case, the stable distribution of the correct solution refers to the probability distribution of infected individuals in the population over time once the system has reached a steady state. This distribution is said to be stationary because it does not change over time, while the correct solution refers to the non-zero probability of individuals being infected.

    Remark 4.2. According to Theorems 3.1 and 4.1, if R^s_0 < 1 , the disease will become extinct under mild additional conditions, whereas if R^s_0 > 1 , the disease will be stochastically persistent. The value of R^s_0 can determine the extinction of the disease or not, and thus it can be considered as a threshold for the stochastic system (1.2).

    In this section, we give several numerical examples to support our results. Employing Milstein's high-order method [32], the discretized system is

    \begin{equation} \left\{ \begin{aligned} S^{k+1} = &S^k+[aA-\beta f(N^k)S^k I^k-\mu S^k+\delta R^k] \Delta t-\sigma f(N^k) S^k I^k \sqrt{\Delta t}\varrho_k\\ &+\frac{1}{2}\sigma^2f(N^k)S^k(I^k)^2(f'(N^k)S^k+f(N^k))(\varrho_k^2-1)\Delta t, \\ I^{k+1} = &I^k+[bA+\beta f(N^k)S^k I^k-(\mu+\gamma+\alpha)I^k] \Delta t+\sigma f(N^k) S^k I^k \sqrt{\Delta t}\varrho_k\\ &+\frac{1}{2}\sigma^2f(N^k)I^k(S^k)^2(f'(N^k)I^k+f(N^k))(\varrho_k^2-1)\Delta t, \\ R^{k+1} = &R^k+[cA+\gamma I-(\mu+\delta)R^k] \Delta t, \end{aligned} \right. \end{equation} (5.1)

    where the time increment \Delta t > 0 , \varrho_k for k = 1, 2, ..., n are Gaussian random variables following the standard normal distribution.

    In this part, we focus on the dynamical behavior of system (1.2) with standard incidence. Let

    f(N) = \frac{\lambda}{N}.

    Assume

    \begin{equation} \begin{aligned} &A = 6, \ \ a = 0.9, \ \ \beta = 0.1, \ \ \alpha = 0.2, \ \ \mu = 0.02, \ \ \delta = 0.1, \\ & \lambda = 10, \ \ {\gamma = 0.5}, \ \ S(0) = 500, \ \ I(0) = 1, \ \ R(0) = 1, \end{aligned} \end{equation} (5.2)

    Parameters b , c and \sigma will take different values in different examples.

    Example 1. (Stationary distribution) Let b = 0 and c = 0.1 , then we obtain R_0 = 1.3657 > 1 . From [3], the disease of the deterministic system (1.1) will persist in a long term (Figure 1).

    Figure 1.  The pictures on the left present the numbers of S , I and R of system (1.2) with b = 0 and R_0^s = 1.3588 , and its deterministic system (1.1) with R_0 = 1.3657 . The pictures on the right show the corresponding frequency histogram of S , I and R with 50,000 iteration points, respectively. The run time of our code is about 1.6488 seconds on a standard computer with a 2.0 GHz processor and 8 GB of RAM.

    For system (1.2), let \sigma = 0.01 and one obtains

    R_0^s = R_0-\frac{\sigma^2f^2(\frac{A}{\mu})A^2}{2(\mu+\gamma+\alpha)\mu^2} = 1.3588 > 1.

    From Theorem 4.1, system (1.2) admits an ergodic stationary distribution (Figure 1).

    For the case with b = 0.1 and c = 0 , we choose \sigma = 0.075 such that R_0^s = R_0-\frac{\sigma^2f^2(\frac{A}{\mu})A^2}{2(\mu+\gamma+\alpha)\mu^2} = 0.9983 < 1 . From Theorem 4.1, system (1.2) admits an ergodic stationary distribution (Figure 2).

    Figure 2.  The pictures on the left present the numbers of S , I and R of system (1.2) with b = 0.1 and R_0^s = 0.9983 , and its deterministic system (1.1) with R_0 = 1.3657 . The pictures on the right show the corresponding frequency histogram of S , I and R with 50,000 iteration points, respectively. The run time of our code is about 1.7667 seconds.

    Example 2. (Extinction) Let b = 0 , c = 0.1 , \sigma = 0.1 , and the other parameters are shown in (5.2) such that

    \sigma^2-\max \{\frac{\beta^2}{2(\mu+\gamma+\alpha)}, \frac{\beta \mu}{f(\frac{A}{\mu})A} \} = 1.7347\times 10^{-18} > 0,

    then from Theorem 3.1, the disease of system (1.2) will become extinct, see Figure 3.

    Figure 3.  The pictures present the numbers of S , I and R of system (1.2) with b = 0 and \sigma = 0.1 , and its deterministic system (1.1) with R_0 = 1.3763 > 1 .

    Let b = 0 , c = 0.1 and \sigma = 0.08 and the other parameters are shown in (5.2) such that

    R_0^s = R_0-\frac{\sigma^2f^2(\frac{A}{ \mu})A^2}{2(\mu+\gamma+\alpha)\mu^2} = 0.9318 < 1,

    and \sigma^2-\frac{\beta\mu}{f(\frac{A}{\mu})A} = -0.0036 < 0 . According to Theorem 3.1, the disease of system (1.2) will be extinct (Figure 4).

    Figure 4.  The phase diagram presents the numbers of S , I and R of system (1.2) with b = 0 , \sigma = 0.08 and R_0^s = 0.9318 , and its deterministic system (1.1) with R_0 = 1.3763 > 1 .

    In this part, we investigate the threshold dynamics of deterministic system (1.1) and stochastic system (1.2) with mass action incidence. Let

    f(N) = \lambda.

    where \lambda is a positive constant. Assume

    \begin{equation} \begin{aligned} A = 10, \ \ a = 0.9, \ \ \alpha = 0.2, \ \ \mu = 0.02, \ \ \delta = 0.2, \ \ \lambda = 1, \ \ S(0) = 500, \ \ I(0) = 1, \ \ R(0) = 1. \end{aligned} \end{equation} (5.3)

    Parameters \beta , b , c and \sigma will take different values in different examples.

    Example 3. (Stationary distribution) First, consider the persistence of the disease of system (1.2) with \beta = 0.002 , b = 0 and c = 0.1 . Then we obtain R_0 = 1.3763 > 1 . From [3], the disease of the deterministic system (1.1) will persist in a long term, see Figure 1.

    For the stochastic system (1.2), let \sigma = 0.0002 and we obtain

    R_0^s = R_0-\frac{\sigma^2f^2(\frac{A}{\mu})A^2}{2(\mu+\gamma+\alpha)\mu^2} = 1.37626 > 1.

    From Theorem 4.1, the stochastic system (1.2) admits an ergodic stationary distribution. See Figure 5.

    Figure 5.  The pictures on the left present the numbers of S , I and R of system (1.2) with b = 0 and R_0^s = 1.37626 , and its deterministic system (1.1) with R_0 = 1.3763 . The pictures on the right show the corresponding frequency histogram of S , I and R with 50,000 iteration points, respectively. The run time of our code is about 1.6803 seconds.

    For the case with b = 0.1 and c = 0 , we choose \beta = 0.001 and \sigma = 0.0002 such that R_0^s = R_0-\frac{\sigma^2f^2(\frac{A}{\mu})A^2}{2(\mu+\gamma+\alpha)\mu^2} = 0.6944 < 1 . From Theorem 4.1, system (1.2) admits an ergodic stationary distribution. See Figure 6.

    Figure 6.  The pictures on the left present the numbers of S , I and R of the stochastic system (1.2) with b = 0.1 and R_0^s = 0.6944 , and its deterministic system (1.1) with R_0 = 1.3763 . The pictures on the right show the corresponding frequency histogram of S , I and R with 50,000 iteration points, respectively. The run time of our code is about 1.7259 seconds.

    Example 4. (Extinction) Let b = 0 , c = 0.1 , \beta = 0.0015 and \sigma = 0.002 , and the other parameters are shown in (5.3) such that

    \sigma^2-\max \{\frac{\beta^2}{2(\mu+\gamma+\alpha)}, \frac{\beta \mu}{f(\frac{A}{\mu})A} \} = 10^{-6} > 0,

    thus from Theorem 3.1, the disease of system (1.2) will be extinct exponentially in a long term (Figure 7).

    Figure 7.  The phase diagram presents the numbers of S , I and R of system (1.2) with b = 0 , c = 0.1 , \beta = 0.0015 and \sigma = 0.002 , and its deterministic system (1.1) with R_0 = 1.0322 .

    Let b = 0 , c = 0.1 , \beta = 0.0014 and \sigma = 0.0015 and the other parameters are shown in (5.3) such that

    R_0^s = R_0-\frac{\sigma^2f^2(\frac{A}{ \mu})A^2}{2(\mu+\gamma+\alpha)\mu^2} = 0.9634 < 1,

    and \sigma^2-\frac{\beta\mu}{f(\frac{A}{\mu})A} = -5.5\times 10^{-7} < 0 . According to Theorem 3.1, the disease of system (1.2) will be extinct exponentially in a long term (Figure 8).

    Figure 8.  The phase diagram presents the numbers of S , I and R of system (1.2) with b = 0 , c = 0.1 , \beta = 0.0014 , \sigma = 0.0015 and R_0^s = 0.9634 , and its deterministic system (1.1) with R_0 = 0.9634 .

    In this study, we present a stochastic SIRS epidemic model with constant immigration and general incidence rate. Our results show that the threshold parameter

    R_0^s = R_0-\frac{\sigma^2f^2(\frac{A}{\mu})A^2}{2(\mu+\gamma+\alpha)\mu^2}

    for this model is lower than its deterministic counterpart ( R_0^s < 1 < R_0 ). In this scenario, the deterministic system may have an endemic state, while the stochastic system leads to disease extinction with probability one (Theorem 3.1). On the other hand, if R_0^s > 1 , the distribution of solution converge in L^{1} to an invariant density (Theorem 4.1), indicating that environmental fluctuations can positively impact the control of infectious diseases. Moreover, if there is a constant influx of infected population, i.e. b > 0 , the stationary distribution will always exist and the disease will persist. We contend that conducting a comprehensive analysis of the influence of migration on the dynamics of our model will yield valuable insights into the intricate interplay between migration and disease transmission.

    This work was supported by Department of Science and Technology of Jilin Province (No. 20210509040RQ), National Natural Science Foundation of China (No. 12271201), Innovation and Entrepreneurship Talent Funds of Jilin Province (No. 2022ZY22) and the Research Funds of Jilin University of Finance and Economics (No. 2022YB025).

    The authors declare there is no conflict of interest.



    [1] W. O. Kermack, A. G. McKendrick, Contributions to the mathematical theory of epidemics-I, Bltn. Mathcal. Biology, 53 (1991), 33–55. https://doi.org/10.1007/bf02464423 doi: 10.1007/bf02464423
    [2] Z. Ma, Y. Zhou, J. Wu, Modeling and dynamics of infectious diseases, World Scientific Publishing, New Jersey, 2009. https://doi.org/10.1142/7223
    [3] J. Li, J. Zhang, Z. Ma, Global analysis of some epidemic models with general contact rate and constant immigration, Appl. Math. Mech., 25 (2004), 396–404. https://doi.org/10.1007/bf02437523 doi: 10.1007/bf02437523
    [4] Y. Cai, Y. Kang, M. Banerjee, W. Wang, A stochastic SIRS epidemic model with infectious force under intervention strategies, J. Differ. Equations, 259 (2015), 7463–7502. https://doi.org/10.1016/j.jde.2015.08.024 doi: 10.1016/j.jde.2015.08.024
    [5] Y. Cai, Y. Kang, M. Banerjee, W. Wang, A stochastic epidemic model incorporating media coverage, Commun. Math. Sci., 14 (2016), 893–910. https://doi.org/10.4310/CMS.2016.v14.n4.a1 doi: 10.4310/CMS.2016.v14.n4.a1
    [6] B. Cao, M. Shan, Q. Zhang, W. Wang, A stochastic SIS epidemic model with vaccination, Phys. A, 485 (2017), 127–143. https://doi.org/10.1016/j.physa.2017.05.083 doi: 10.1016/j.physa.2017.05.083
    [7] Y. Cai, Y. Kang, W. Wang, A stochastic SIRS epidemic model with nonlinear incidence rate, Appl. Math. Comput., 305 (2017), 221–240. https://doi.org/10.1016/j.amc.2017.02.003 doi: 10.1016/j.amc.2017.02.003
    [8] M. De la Sen, S. Alonso-Quesada, A. Ibeas, On the stability of an SEIR epidemic model with distributed time-delay and a general class of feedback vaccination rules, Appl. Math. Comput., 270 (2015), 953–976. https://doi.org/10.1016/j.amc.2015.08.099 doi: 10.1016/j.amc.2015.08.099
    [9] W. Weera, T. Botmart, T. La-inchua, Z. Sabir, R. Núñez, M. Abukhaled, et al., A stochastic computational scheme for the computer epidemic virus with delay effects, AIMS Math., 8 (2023), 148–163. https://doi.org/10.3934/math.2023007 doi: 10.3934/math.2023007
    [10] T. Britton, D. Lindenstrand, Epidemic modelling: aspects where stochasticity matters, Math. Biosci., 222 (2009), 109–116. https://doi.org/10.1016/j.mbs.2009.10.001 doi: 10.1016/j.mbs.2009.10.001
    [11] F. Wang, X. Wang, S. Zhang, C. Ding, On pulse vaccine strategy in a periodic stochastic SIR epidemic model, Chaos, Solitons Fractals, 66 (2014), 127–135. https://doi.org/10.1016/j.chaos.2014.06.003 doi: 10.1016/j.chaos.2014.06.003
    [12] Z. Shi, D. Jiang, N. Shi, T. Hayat, A. Alsaedi, The impact of nonlinear perturbation to the dynamics of HIV model, Math. Methods Appl. Sci., 45 (2022), 2542–2562. https://doi.org/10.1002/mma.7939 doi: 10.1002/mma.7939
    [13] Z. Shi, D. Jiang, N. Shi, A. Alsaedi, Virus infection model under nonlinear perturbation: Ergodic stationary distribution and extinction, J. Franklin Inst., 359 (2022), 11039–11067. https://doi.org/10.1016/j.jfranklin.2022.03.035 doi: 10.1016/j.jfranklin.2022.03.035
    [14] L. Wang, H. Huang, A. Xu, W. Wang, Stochastic extinction in an SIRS epidemic model incorporating media coverage, Abstr. Appl. Anal., 2013 (2013), 891765. https://doi.org/10.1155/2013/891765 doi: 10.1155/2013/891765
    [15] Y. Lin, D. Jiang, Long-time behaviour of a perturbed SIR model by white noise, Discrete Contin. Dyn. Syst. Ser. B, 18 (2013), 1873–1887. https://doi.org/10.3934/dcdsb.2013.18.1873 doi: 10.3934/dcdsb.2013.18.1873
    [16] Q. Liu, D. Jiang, N. Shi, T. Hayat, A. Alsaedi, Stationarity and periodicity of positive solutions to stochastic SEIR epidemic models with distributed delay, Discrete Contin. Dyn. Syst. Ser. B, 22 (2017), 2479–2500. https://doi.org/10.3934/dcdsb.2017127 doi: 10.3934/dcdsb.2017127
    [17] Y. Zhao, D. Jiang, The threshold of a stochastic SIS epidemic model with vaccination, Appl. Math. Comput., 243 (2014), 718–727. https://doi.org/10.1016/j.amc.2014.05.124 doi: 10.1016/j.amc.2014.05.124
    [18] Q. Liu, D. Jiang, N. Shi, T. Hayat, A. Alsaedi, Stationary distribution and extinction of a stochastic SIRS epidemic model with standard incidence, Phys. A, 469 (2017), 510–517. https://doi.org/10.1016/j.physa.2016.11.077 doi: 10.1016/j.physa.2016.11.077
    [19] Z. Shi, X. Zhang, D. Jiang, Dynamics of an avian influenza model with half-saturated incidence, Appl. Math. Comput., 355 (2019), 399–416. https://doi.org/10.1016/j.amc.2019.02.070 doi: 10.1016/j.amc.2019.02.070
    [20] X. Zhang, D. Jiang, T. Hayat, B. Ahmad, Dynamical behavior of a stochastic SVIR epidemic model with vaccination, Phys. A, 483 (2017), 94–108. https://doi.org/10.1016/j.physa.2017.04.173 doi: 10.1016/j.physa.2017.04.173
    [21] Z. Shi, D. Jiang, Dynamical behaviors of a stochastic HTLV-I infection model with general infection form and Ornstein-Uhlenbeck process, Chaos, Solitons Fractals, 165 (2022), 112789. https://doi.org/10.1016/j.chaos.2022.112789 doi: 10.1016/j.chaos.2022.112789
    [22] Y. Lin, D. Jiang, Threshold behavior in a stochastic SIS epidemic model with standard incidence, J. Dyn. Differ. Equations, 26 (2014), 1079–1094. https://doi.org/10.1007/s10884-014-9408-8 doi: 10.1007/s10884-014-9408-8
    [23] X. Mao, G. Marion, E. Renshaw, Environmental Brownian noise suppresses explosions in population dynamics, Stoch. Processes Appl., 97 (2002), 95–110. https://doi.org/10.1016/S0304-4149(01)00126-0 doi: 10.1016/S0304-4149(01)00126-0
    [24] A. Friedman, Stochastic differential equations and applications, in Stochastic Differential Equations, Springer, Berlin, 2010. https://doi.org/10.1007/978-3-642-11079-5_2
    [25] D. W. Stroock, S. R. S. Varadhan, On the support of diffusion processes with applications to the strong maximum principle, in Proceedings of the Sixth Berkeley Symposium on Mathematical Statistics and Probability, 3 (1972), 333–359.
    [26] G. B. Arous, R. Léandre, Décroissance exponentielle du noyau de la chaleur sur la diagonale (II), Probab. Th. Rel. Fields, 90 (1991), 377–402. https://doi.org/10.1007/BF01193751 doi: 10.1007/BF01193751
    [27] K. Pichór, R. Rudnicki, Stability of Markov semigroups and applications to parabolic systems, J. Math. Anal. Appl., 215 (1997), 56–74. https://doi.org/10.1006/jmaa.1997.5609 doi: 10.1006/jmaa.1997.5609
    [28] R. Rudnicki, Long-time behaviour of a stochastic prey-predator model, Stoch. Processes Appl., 108 (2003), 93–107. https://doi.org/10.1016/S0304-4149(03)00090-5 doi: 10.1016/S0304-4149(03)00090-5
    [29] R. Rudnicki, K. Pichór, Influence of stochastic perturbation on prey-predator systems, Math. Biosci., 206 (2007), 108–119. https://doi.org/10.1016/j.mbs.2006.03.006 doi: 10.1016/j.mbs.2006.03.006
    [30] R. Rudnicki, K. Pichór, M. Tyran-Kamińska, Markov semigroups and their applications, in Dynamics of Dissipation, Springer, Berlin, 2022. https://doi.org/10.1007/3-540-46122-1_9
    [31] X. Mu, D. Jiang, A. Alsaedi, Analysis of a stochastic phytoplankton Czooplankton model under non-degenerate and degenerate diffusions, J. Nonlinear Sci., 32 (2022), 35. https://doi.org/10.1007/s00332-022-09787-9 doi: 10.1007/s00332-022-09787-9
    [32] D. J. Higham, An algorithmic introduction to numerical simulation of stochastic differential equations, SIAM Rev., 43 (2001), 525–546. https://doi.org/10.1137/s0036144500378302 doi: 10.1137/s0036144500378302
  • This article has been cited by:

    1. Nabeela Anwar, Iftikhar Ahmad, Adiqa Kausar Kiani, Muhammad Shoaib, Muhammad Asif Zahoor Raja, Novel neuro-stochastic adaptive supervised learning for numerical treatment of nonlinear epidemic delay differential system with impact of double diseases, 2024, 0228-6203, 1, 10.1080/02286203.2024.2303577
  • Reader Comments
  • © 2023 the Author(s), licensee AIMS Press. This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0)
通讯作者: 陈斌, bchen63@163.com
  • 1. 

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

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

Metrics

Article views(1722) PDF downloads(79) Cited by(1)

Other Articles By Authors

/

DownLoad:  Full-Size Img  PowerPoint
Return
Return

Catalog