Research article

Stationary distribution and probability density function analysis of a stochastic Microcystins degradation model with distributed delay

  • Received: 08 September 2023 Revised: 09 November 2023 Accepted: 10 December 2023 Published: 18 December 2023
  • A stochastic Microcystins degradation model with distributed delay is studied in this paper. We first demonstrate the existence and uniqueness of a global positive solution to the stochastic system. Second, we derive a stochastic critical value Rs0 related to the basic reproduction number R0. By constructing suitable Lyapunov function types, we obtain the existence of an ergodic stationary distribution of the stochastic system if Rs0>1. Next, by means of the method developed to solve the general four-dimensional Fokker-Planck equation, the exact expression of the probability density function of the stochastic model around the quasi-endemic equilibrium is derived, which is the key aim of the present paper. In the analysis of statistical significance, the explicit density function can reflect all dynamical properties of a chemostat model. To validate our theoretical conclusions, we present examples and numerical simulations.

    Citation: Ying He, Yuting Wei, Junlong Tao, Bo Bi. Stationary distribution and probability density function analysis of a stochastic Microcystins degradation model with distributed delay[J]. Mathematical Biosciences and Engineering, 2024, 21(1): 602-626. doi: 10.3934/mbe.2024026

    Related Papers:

    [1] Buyu Wen, Bing Liu, Qianqian Cui . Analysis of a stochastic SIB cholera model with saturation recovery rate and Ornstein-Uhlenbeck process. Mathematical Biosciences and Engineering, 2023, 20(7): 11644-11655. doi: 10.3934/mbe.2023517
    [2] Sanling Yuan, Xuehui Ji, Huaiping Zhu . Asymptotic behavior of a delayed stochastic logistic model with impulsive perturbations. Mathematical Biosciences and Engineering, 2017, 14(5&6): 1477-1498. doi: 10.3934/mbe.2017077
    [3] Huili Wei, Wenhe Li . Dynamical behaviors of a Lotka-Volterra competition system with the Ornstein-Uhlenbeck process. Mathematical Biosciences and Engineering, 2023, 20(5): 7882-7904. doi: 10.3934/mbe.2023341
    [4] Helong Liu, Xinyu Song . Stationary distribution and extinction of a stochastic HIV/AIDS model with nonlinear incidence rate. Mathematical Biosciences and Engineering, 2024, 21(1): 1650-1671. doi: 10.3934/mbe.2024072
    [5] Ying He, Junlong Tao, Bo Bi . Stationary distribution for a three-dimensional stochastic viral infection model with general distributed delay. Mathematical Biosciences and Engineering, 2023, 20(10): 18018-18029. doi: 10.3934/mbe.2023800
    [6] Ke Qi, Zhijun Liu, Lianwen Wang, Qinglong Wang . Survival and stationary distribution of a stochastic facultative mutualism model with distributed delays and strong kernels. Mathematical Biosciences and Engineering, 2021, 18(4): 3160-3179. doi: 10.3934/mbe.2021157
    [7] Fathalla A. Rihan, Hebatallah J. Alsakaji . Analysis of a stochastic HBV infection model with delayed immune response. Mathematical Biosciences and Engineering, 2021, 18(5): 5194-5220. doi: 10.3934/mbe.2021264
    [8] H. J. Alsakaji, F. A. Rihan, K. Udhayakumar, F. El Ktaibi . Stochastic tumor-immune interaction model with external treatments and time delays: An optimal control problem. Mathematical Biosciences and Engineering, 2023, 20(11): 19270-19299. doi: 10.3934/mbe.2023852
    [9] Jing Hu, Zhijun Liu, Lianwen Wang, Ronghua Tan . Extinction and stationary distribution of a competition system with distributed delays and higher order coupled noises. Mathematical Biosciences and Engineering, 2020, 17(4): 3240-3251. doi: 10.3934/mbe.2020184
    [10] Ming-Zhen Xin, Bin-Guo Wang, Yashi Wang . Stationary distribution and extinction of a stochastic influenza virus model with disease resistance. Mathematical Biosciences and Engineering, 2022, 19(9): 9125-9146. doi: 10.3934/mbe.2022424
  • A stochastic Microcystins degradation model with distributed delay is studied in this paper. We first demonstrate the existence and uniqueness of a global positive solution to the stochastic system. Second, we derive a stochastic critical value Rs0 related to the basic reproduction number R0. By constructing suitable Lyapunov function types, we obtain the existence of an ergodic stationary distribution of the stochastic system if Rs0>1. Next, by means of the method developed to solve the general four-dimensional Fokker-Planck equation, the exact expression of the probability density function of the stochastic model around the quasi-endemic equilibrium is derived, which is the key aim of the present paper. In the analysis of statistical significance, the explicit density function can reflect all dynamical properties of a chemostat model. To validate our theoretical conclusions, we present examples and numerical simulations.



    Microcystins (MCs) represent a highly prevalent form of cyanotoxins derived from harmful cyanobacterial blooms. These toxins can introduce toxicity into vital organs including liver, kidney, heart etc., posing major threats to human health [1]. Biodegradation has emerged as an economically and environmentally favorable approach for MCs degradation, garnering considerable research attention [2,3]. Recently, Song et al. [4] proposed a degradation model as follows:

    {x1(t)=Da10a12x1(t)x2(t)a13x1(t)x3(t)(D+d1)x1(t),x2(t)=a21x2(t)x1(t)a20x2(t)(D+d2)x2(t),x3(t)=a30x2(t)a31x1(t)x3(t)(D+d3)x3(t), (1.1)

    where x1(t),x2(t) and x3(t) denote the concentrations of MCs, MC-degrading bacteria and degrading enzymes produced by MC-degrading bacteria at time t, respectively. a10>0 is the input concentration of MCs; D is the washout rate; a120 is the consumption rate of MCs; a210 is the maximal growth rate of MC-degrading bacteria; a200 is the consumption rate of MC-degrading bacteria; a300 is the maximal growth rate of degrading enzymes; a130 and a310 are the degradation rate of MCs and the consumption rate of degrading enzymes, respectively; di>0(i=1,2,3) represents the death rates of MCs, MC-degrading bacteria and degrading enzymes, as distinguished by, i = 1, 2, 3, respectively. Re-scaling model (1.1) by

    x1=a10x,x2=y,x3=z,t=TD,D1=D+d1D,D2=a20+D+d2D,D3=D+d3D,a1=a12D,a2=a13D,b1=a21a10D,c1=a30D,c2=a31a10D,

    the authors of [5] obtained the following dimensionless model:

    {x(t)=1a1x(t)y(t)a2x(t)z(t)D1x(t),y(t)=b1x(t)y(t)D2y(t),z(t)=c1y(t)c2x(t)z(t)D3z3(t). (1.2)

    Time delays cannot be ignored in models for MCs degradation, especially the delay between MC-degrading bacteria and the production of new degrading enzymes particles. There has been extensive discussion of the presence of such delays and their impact on microbial population dynamics (see [6,7,8]). Song et al. [9] considered the following model with time delay:

    {x(t)=1a1x(t)y(t)a2x(t)z(t)D1x(t),y(t)=b1eδτx(tτ)y(tτ)D2y(t),z(t)=c1y(t)c2x(t)z(t)D3z3(t).

    Here the constant τ0 represents the time consumed by MC-degrading bacteria to convert MCs, once consumed, into viable biomass and eδτ represents the survival probability of degrading bacteria during the conversion process. However, a more realistic delayed MC-degrading model should be an equation with distributed delay, which describes the cumulative effect of the past history of a system. In this paper, we will mainly consider the following degradation model with distributed delay

    dz(t)=[c1tK(ts)y(s)dsc2x(t)z(t)D3z3(t)]dt.

    According to MacDonald [10], the distributed delay function f(t) can be described by the Gamma distribution

    K(t)=αn+1tneαtn!,

    as a kernel, where α>0 and n is a nonnegative integer. In this article, we mainly consider the weak kernel, that is K(t)=αeαt(n=0) with α>0.

    On the other hand, in the real world, population systems are inevitably influenced by the environmental noise (see [11,12,13,14,15]). According to May [16], the birth rate, carrying capacity, competition coefficient, and other parameters should exhibit random fluctuations as a result of environmental noise. In order to investigate the impact of environmental noise on population dynamics, various scholars have added white noise into deterministic systems (see [17,18,19,20,21]).

    Therefore, based on the model (1.2), we further consider the following stochastic MCs degradation model with distribution delay

    {dx(t)=[1a1x(t)y(t)a2x(t)z(t)D1x(t)]dt+σ1x(t)dB1(t),dy(t)=[b1x(t)y(t)D2y(t)]dt+σ2y(t)dB2(t),dz(t)=[c1tK(ts)y(s)dsc2x(t)z(t)D3z3(t)]dt+σ3z(t)dB3(t), (1.3)

    where Bi(t) denotes mutually independent standard Brownian motions and σ2i>0 denotes the intensity of white noises i=1,2,3. Set

    w(t)=tK(ts)y(s)ds=tαeα(ts)y(s)ds.

    Then, by using the linear chain technique, system (1.3) is transformed into the following equivalent system:

    {dx(t)=[1a1x(t)y(t)a2x(t)z(t)D1x(t)]dt+σ1x(t)dB1(t),dy(t)=[b1x(t)y(t)D2y(t)]dt+σ2y(t)dB2(t),dz(t)=[c1w(t)c2x(t)z(t)D3z(t)]dt+σ3z(t)dB3(t),dw(t)=[αy(t)αw(t)]dt. (1.4)

    In this paper, we propose a Microcystins degradation model with distributed delay and stochastic perturbation. As is known to us, there is no work concerning the existence of stationary distribution in system (1.3). The main difficulty is to find appropriate Lyapunov functions. It should be emphasized that there are relatively few studies focusing on the explicit expression of probability density functions due to the difficulty of solving the corresponding Fokker-Planck equation.

    The paper is organized as follows. In Section 2, we prove that the stochastic Microcystins degradation system (1.4) with distributed delay has a unique global positive solution. In Section 3, a unique global solution of system (1.4) is proved to be stationary by constructing several suitable Lyapunov functions. Numerical simulations are presented to illustrate our analytical results in Section 4.

    To study the dynamical behavior of Microcystins degradation model, we first need to consider whether the solution is global and positive. In this section, we will prove the existence and uniqueness of a global solution to system (1.4) with any positive initial value.

    Theorem 2.1: For any initial value (x(0),y(0),z(0),w(0))R4+, there is a unique solution of the system (1.4) on t0 and the solution will remain in R4+ with probability one.

    Proof. The coefficients of system (1.3) satisfy the local Lipschitz condition, so there is a unique local solution (x(t),y(t),z(t),w(t)) on t[0,τe), where τe is the explosion time. Now we shall show that this solution is global, i.e., prove τe= a.s. Let k01 be sufficiently large so that x(0),y(0),z(0)andw(0) all lie within the interval [1k0,k0]. For each integer kk0, define the stopping time:

    τk=inf{t[0,τe):min{x(t),y(t),z(t),w(t)}1kormax{x(t),y(t),z(t),w(t)}k},

    where throughout this paper we set inf=. Evidently, τk is increasing as k. Denote τ=limkτk, where ττe a.s. If we can prove that τ= a.s., then τe= and (x(t),y(t),z(t),w(t))R4+ a.s. for all t0. In other words, to complete the proof all we need to prove is that τ= a.s. If this statement is false, then there is a pair of constants T>0 and ε(0,1) such that

    P{τT}>ε.

    Hence there is an integer k1k0 such that

    P{τkT}>εfor all kk1. (2.1)

    Define the nonnegative C2 Lyapunov function V0 as follows

    V0(x,y,z,w)=x1lnx+aylny(1+lna)+czdlnz(d+lncd)+bwlnw(1+lnb),

    where a,b,c,d are positive constants which will be determined later.

    Based on the basic inequality u1lnu0, we have

    czdlnz(d+lncd)=d(czd1lnczd)0.

    Applying Itô's formula, we have

    dV0=LV0dt+σ1(x1)dB1(t)+σ2(ay1)dB2(t)+σ3(czd)dB3(t), (2.2)
    LV0=1a1xya2xzD1x+a(b1xyD2y)+c(c1wc2xzD3z)+b(αyαw)1x+a1y+a2z+D1+12σ21b1x+D2+12σ22+d(c1wz+c2x+D3+12σ23)αyw+α1+D1+D2+dD3+α+12(σ21+σ22+dσ23)+(ab1a1)xy+(a2cD3)z+(D1b1+dc2)x+(aD2+a1+bα)y+(cc1bα)w.

    Choose b=a2c1αD3,c=a2D3,d=b+D1c and a=min{a1b1,a1D3+c1a2D2D3} such that ab1a10and a1aD2+αb0, then, we have

    LV01+D1+D2+dD3+α+12(σ21+σ22+dσ23):=K1,

    where K1 is a positive constant. Integrating equation (2.2) from 0to τkT and then taking the expectation on both sides, we have

    EV0(x(τkT),y(τkT),z(τkT),w(τkT)))V0(x(0),y(0),z(0),w(0))+K1E(τkT)V0(x(0),y(0),z(0),w(0))+K1T (2.3)

    Set Ωk={τkT}for kk1 and in view of Eq (2.1), we have P(Ωk)ϵ. Note that for every ωΩk, there is x(τk,ω)or y(τk,ω)or z(τk,ω)or w(τk,ω) that is equal to either kor1k. Thus V0(x(τk,ω),y(τk,ω),z(τk,ω),w(τk,ω)) is no less than either

    (k1lnk)(aklnk(1+lna))d(cdk1lnckd)(bklnk(1+lnb))or(1k1ln1k)(a1kln1k(1+lna))d(cd1k1lncdk)(b1kln1k(1+lnb))=(k1lnk)(aklnk(1+lna))d(cdk1ln(ck)+lnd)(bklnk(1+lnb))or(1k1+lnk)(a1k+lnk(1+lna))d(cd1k1lnc+ln(dk))(b1k+lnk(1+lnb))

    Therefore, we obtain

    V0(x(τk,ω),y(τk,ω),z(τk,ω),w(τk,ω)))(k1lnk)(aklnk(1+lna))d(cdk1ln(ck)+lnd)(bklnk(1+lnb))(1k1+lnk)(a1k+lnk(1+lna))d(cd1k1lnc+ln(dk))(b1k+lnk(1+lnb))

    It follows from Eq (2.3) that

    V0(x(0),y(0),z(0),w(0))+K1TE[IΩk(ω)V0(x(τk,ω),y(τk,ω),z(τk,ω),w(τk,ω)))]ϵ(k1lnk)(aklnk(1+lna))d(cdk1ln(ck)+lnd)(bklnk(1+lnb))(1k1+lnk)(a1k+lnk(1+lna))d(cd1k1lnc+ln(dk))(b1k+lnk(1+lnb))

    where IΩk(ω) is the indicator function of Ωk. Let k, we obtain

    >V0(x(0),y(0),z(0),w(0))+K1T=

    which is a contradiction and so we must have that τ=. This completes the proof.

    For the stochastic model, we mainly focus on the existence of a stationary distribution. Define a critical value

    Rs0=b1(D1+12σ21)(D2+12σ22).

    Consider the integral equation

    X(t)=X(t0)+tt0b(s,X(s))ds+kr=1tt0σr(s,X(s))dB(s). (3.1)

    Lemma 3.1. ([22]) Suppose that the coefficient of (3.1) is independent of t and satisfies the following condition for some constant B:

    |b(s,x)b(s,y)|+kr=1|σr(s,x)σr(s,y)|B|xy|,|b(s,x)|+kr=1|σr(s,x)|B(1+|x|), (3.2)

    in URRd+ for every R>0, also suppose that there exists a nonnegative C2-function V(x)Rd+ such that

    LV(x)1outside some compact set.

    Then system (3.1) has a solution, which is a stationary distribution.

    Remark 3.1. According to the proof of Lemma 2.1 in [23], we know that condition (3.2) in Lemma 3.1 can be replaced by the global existence of the solutions. Hence, Theorem 2.1 indicates that (3.2) in Lemma 3.1 is satisfied.

    Theorem 3.1. For any initial value (x(0),y(0),z(0),w(0))R4+, if Rs0>1, then system (1.4) has at least one stationary solution (x(t),y(t),z(t),w(t))R4+.

    Proof. By Lemma 3.1, we only need to construct a nonnegative C2-function ˜V(x,y,z,w) and a closed set UR4+ such that

    L˜V(x,y,z,w)1,for any (x,y,t,w)R4+U.

    Define

    V1(x,y,z,w)=lnxe1lny+a2D3z+c1a2αD3w,

    where e1 is a positive constant to be chosen later. Then applying Itô's formula to V1, we have

    LV11xe1b1x+e1(D2+12σ22)+D1+12σ21+a1y+a2z+a2D3(c1wc2xzD3z)+c1a2D3α(αyαw)2b1e1+e1(D2+12σ22)+D1+12σ21+(a1+c1a2D3)y (3.3)

    Taking e1=b1(D2+12σ22)2, then in view of (3.3) one can see that

    LV1b1D2+12σ22+D1+12σ21+(a1+c1a2D3)y=(D1+12σ21)(Rs01)+(a1+c1a2D3)y=:λ+(a1+c1a2D3)y, (3.4)

    where

    Rs0=b1(D1+12σ21)(D2+12σ22),λ=(D1+12σ21)(Rs01).

    Define

    V2=lnxlnwlnz.

    Then applying Itô's formula to V2, one gets

    LV21xαywc1wz+a1y+a2z+D1+12σ21+c2x+D3+12σ23+α. (3.5)

    Select

    V3(x,y,z,w)=1θ+2(x+a1b1y+a1D24b1c1z+a1D22αb1w)θ+2,

    where 0<θ<min{D1σ212D1+σ212,D2σ22D2+σ22,D3σ232D3+σ232}. Applying the same method to V3, we have

    LV3(x+a1b1y+a1D24b1c1z+a1D22αb1w)θ+1{1D1xa1D22b1ya1D2D34b1c1za1D24b1w}+θ+12(x+a1b1y+a1D24b1c1z+a1D22αb1w)θ[σ21x2+σ22(a1b1)2y2+σ23(a1D24c1b1)2z2]}(x+a1b1y+a1D24b1c1z+a1D22αb1w)θ+1D1xθ+2D22(a1b1)θ+2yθ+2D3(a1D24b1c1)θ+2zθ+2α2(a1D22αb1)θ+2wθ+2+θ+12(x+a1b1y+a1D24b1c1z+a1D22αb1w)θ[σ21x2+σ22(a1b1)2y2+σ23(a1D24c1b1)2z2]D1θxθ+2D2θ2(a1b1)θ+2yθ+2D3θ(a1D24b1c1)θ+2zθ+2α4(a1D22αb1)θ+2wθ+2+B0, (3.6)

    where

    B0=sup(x,y,z,w)R4+{(x+a1b1y+a1D24b1c1z+a1D22αb1w)θ+1D1(1θ)xθ+2D2(1θ)2(a1b1)θ+2yθ+2D3(1θ)(a1D24b1c1)θ+2zθ+2α4(a1D22αb1)θ+2wθ+2+θ+12(x+a1b1y+a1D24b1c1z+a1D22αb1w)θ[σ21x2+σ22(a1b1)2y2+σ23(a1D24c1b1)2z2]}<.

    Consequently, we have

    V=M0V1+V2+V3, (3.7)

    where M0>0 satisfies that M0Rs0+c02. Also,

    c0=sup(x,y,z,w)R4+{D1θxθ+2D2θ2(a1b1)θ+2yθ+2D3θ(a1D24b1c1)θ+2zθ+2α4(a1D22αb1)θ+2wθ+2+a2z+c2x+D1+D3+12σ21+12σ32+α+B0}.

    It is obvious that

    lim infn,(x,y,z,w)R4+KnV(x,y,z,w)=+,

    and Kn=(1n,n)×(1n,n)×(1n,n)×(1n,n). It is evident that V(x,y,z,w) is a continuous functions, so there exists a minimum point (x0,y0,z0,w0)R4+. Thereafter, we formulate the nonnegative C2 function ˜V as follows:

    ˜V(x,y,z,w)=V(x,y,z,w)V(x0,y0,z0,w0).

    Combining with (3.4), (3.5), (3.6) and (3.7), we obtain that

    L˜VM0λ+[M0(a1+a2c1D3)+a1]y+a2z+c2x+D1+12σ21+D3+12σ23+αD1θxθ+2D2θ2(a1b1)θ+2yθ+2D3θ(a1D24b1c1)θ+2zθ+2α4(a1D22αb1)θ+2wθ+2+B01xαywc1wz.

    Next, we define the following bounded closed set

    Uε={εx1ε,εy1ε,ε2w1ε3,ε3z1ε3}.

    For convenience, we can divide R4+Uε into eight regions:

    U1={(x,y,z,w)R4+,0<y<ε},U2={(x,y,z,w)R4+,0<x<ε},U3={(x,y,z,w)R4+,0<w<ε2,y>ε},U4={(x,y,z,w)R4+,0<z<ε3,wε2},U5={(x,y,z,w)R4+,y>1ε},U6={(x,y,z,w)R4+,x>1ε},U7={(x,y,z,w)R4+,z>1ε3},U8={(x,y,z,w)R4+,w>1ε2}.

    Case 1. If (x,y,z,w)U1, in view of (3.8) we get

    L˜VM0λ+[M0(a1+a2c1D3)+a1]y+a2z+c2xD1θxθ+2D2θ2(a1b1)θ+2yθ+2D3θ(a1D24b1c1)θ+2zθ+2α4(a1D22αb1)θ+2wθ+2+B0+D1+12σ21+D3+12σ23+αM0λ+{M0(a1+a2c1D3)+a1]ε+c01;

    here, ε is a positive constant small enough to satisfy the following conditions:

    M0λ+[M0(a1+a2c1D3)+a1]ε+c01, (3.8)
    1ε+E1, (3.9)
    αε+E1, (3.10)
    c1ε+E1, (3.11)
    D2θ4(a1b1ε)θ+2+E1, (3.12)
    D1θ2(1ε)θ+2+E1, (3.13)
    D3θ2(a1D24b1c1ε3)θ+2+E1, (3.14)
    α8(a1D22αb1ε2)θ+2+E1, (3.15)

    where

    E=sup(x,y,z,w)R4+{D1θ2xθ+2D2θ4(a1b1)θ+2yθ+2D3θ2(a1D24b1c1)θ+2zθ+2α8(a1D22αb1)θ+2wθ+2+a2z+c2x+[M(a1+a2c1D3)+a1]y+D1+12σ21+D3+12σ23+α+B0}.

    Case 2 : If (x,y,z,w)U2 according to (3.9), it implies that

    L˜V1x+[M0(a1+a2c1D3)+a1]y+a2z+c2xD1θxθ+2D2θ2(a1b1)θ+2yθ+2D3θ(a1D24b1c1)θ+2zθ+2α4(a1D22αb1)θ+2wθ+2+B0+D1+12σ21+D3+12σ23+α1ε+E1.

    Case 3 : If (x,y,z,w)U3, by (3.10) we have

    L˜Vαyw+Eαε+E1.

    Case 4 : If (x,y,z,w)U4, according to (3.11), we deduce that

    L˜Vc1wz+Ec1ε+E1.

    Case 5 : If (x,y,z,w)U5, by condition (3.12), we conclude that

    L˜VD2θ4(a1b1)θ+2yθ+2+ED2θ4(a1b1ε)θ+2+E1.

    Case 6 : If (x,y,z,w)U6, and the condition (3.13) is satisfied, it follows that

    L˜VD1θ2xθ+2+ED1θ2(1ε)θ+2+E1.

    Case 7 : If (x,y,z,w)U7, in view of (3.14) we get

    L˜VD3θ2(a1D24b1c1)θ+2zθ+2+ED3θ2(a1D24b1c1ε3)θ+2+E1.

    Case 8 : If (x,y,z,w)U8, from (3.15), it is deduced that

    L˜Vα8(a1D22αb1)θ+2wθ+2+Eα8(a1D22αb1ε2)θ+2+E1.

    Clearly, there exists a small enough constant ε to make the conclusion

    L˜V1,for all (x,y,z,w)UCε.

    Remark 3.2. In the proof of the above theorem, the construction of V3 means to obtain the constant E and high-order terms of x,yandz.

    Let u1(t)=lnx,u2(t)=lny,u3(t)=lnzandu4(t)=lnw. Using the Itô's formula, an equivalent equation of system (1.4) is given as follows

    {du1=[eu1a1eu2a2eu3(D1+12σ21)]dt+σ1dB1(t),du2=[b1eu1(D2+12σ22)]dt+σ2dB2(t),du3=[c1eu4u3c2eu1(D3+12σ23)]dt+σ3dB3(t),du4=(αeu2u4α)dt. (4.1)

    Assume that Rs0>1, then, there is the quasi-endemic equilibrium E=(x,y,z,w)=(eu1,eu2,eu3,eu4) where

    x=D2+12σ22b1,y=[b1(D3+12σ23)+c2(D2+12σ22)][b1(D1+12σ21)(D2+12σ22)]a2b1c1(D2+12σ22)+a1b1(D2+12σ22)(D3+12σ23)+a1c2(D2+12σ22)2,z=b1c1[b1(D1+12σ21)(D2+12σ22)]a2b1c1(D2+12σ22)+a1b1(D2+12σ22)(D3+12σ23)+a1c2(D2+12σ22)2,w=y=[b1(D3+12σ23)+c2(D2+12σ22)][b1(D1+12σ21)(D2+12σ22)]a2b1c1(D2+12σ22)+a1b1(D2+12σ22)(D3+12σ23)+a1c2(D2+12σ22)2.

    Before proving Theorem 4.1, we need to firstly introduce two important standard matrices in the following Lemmas 4.1 and 4.2.

    Lemma 4.1 ([24]). For the algebraic equation H20+B0Σ0+Σ0BT0=0, where H0=diag(1,0,0,0) and Σ0 is a real symmetric matrix, we have the standard matrix

    B0=(γ1γ2γ3γ4100001000010).

    If γ1>0,γ3>0,γ4>0 and γ1γ2γ3γ23γ21γ4>0, then Σ0 is a positive definite matrix, where

    Σ0=(γ2γ3γ1γ42(γ1γ2γ3γ23γ21γ4)0γ32(γ1γ2γ3γ23γ21γ4)00γ32(γ1γ2γ3γ23γ21γ4)0γ12(γ1γ2γ3γ23γ21γ4)γ32(γ1γ2γ3γ23γ21γ4)0γ12(γ1γ2γ3γ23γ21γ4)00γ12(γ1γ2γ3γ23γ21γ4)0γ1γ2γ32γ4(γ1γ2γ3γ23γ21γ4)).

    Here B0 in this form is called the standard R1 matrix.

    Lemma 4.2 ([24]). For the algebraic equation H20+E0Ω0+Ω0ET0=0, where H0=diag(1,0,0,0) and Ω0 is a real symmetric matrix, we have the standard matrix

    E0=(τ1τ2τ3τ410000100000τ5).

    If τ1>0,τ3>0andτ1τ2τ3>0, then the matrix Ω0 is a semi-positive definite matrix, which means that

    Ω0=(τ22(τ1τ2τ3)012(τ1τ2τ3)0012(τ1τ2τ3)0012(τ1τ2τ3)0τ12τ3(τ1τ2τ3)00000).

    Denote P=(p1,p2,p3,p4)T=(u1u1,u2u2,u3u3,u4u4)T(k=1,2,3,4). Therefore, the linearized equation system of Eq (4.1) is obtained as follows

    {dp1=(b11p1b12p2b13p3)dt+σ1dB1(t),dp2=b21p1dt+σ2dB2(t),dp3=(b31p1b33p3+b33p4)dt+σ3dB3(t),dp4=(αp2αp4)dt, (4.2)

    where

    b11=1x>b12+b13,b12=a1y,b13=a2z,b21=b1x,b31=c2x,b33=c1wz>b31.

    Theorem 4.1. Assume that Rs0>1; for any initial value (x(0),y(0),z(0),w(0))R4+, then the solution (x,y,z,w) of system (1.4) will have a log-normal probability density Φ(x,y,z,w) around E, which is given by

    Φ(x,y,z,w)=(2π)2Σ12e12(lnxx,lnyy,lnzz,lnww)Σ1(lnxx,lnyy,lnzz,lnww)T

    with Σ=Σ1+Σ2+Σ3 as a positive definite matrix; Σ1,Σ2andΣ3 are defined as follows

    Σ1={(b31b33σ1)2(T3T2T1)1¯Ω2[(T3T2T1)1]T,if q1=0,(b31b33q1σ1)2(T4T2T1)1Ω1[(T4T2T1)1]T,if q10,
    Σ2={α2b213(b12b31+αb33)2σ22b212(T7T6T5)1Ω1[(T7T6T5)1]T,if q2=0,(b12b31+b33α)2σ22(T9T8T6T5)1˜Ω2[(T9T8T6T5)1]T,if q20,q3=0,(b12b31+b33α)2q23σ22(T10T8T6T5)1Ω1[(T10T8T6T5)1]T,if q20,q30,
    Σ3=(b13b21ασ3)2(T12T11)1Ω1[(T12T11)1]T,

    where q1=b21α[b13b33(b21+b33)α]b231b33b33,q2=α(b11+α)b12,q3=b13αb12b12b33q2b12b31+b33α+b12q2(b12b33q2+b12b31α+b33α2)(b12b31+b33α)2, and the matrices T1,T12,Ω1,¯Ω2,˜Ω3 are defined in the following proof.

    Proof. For convenience and simplicity, let B(t)=(B1(t),B2(t),B3(t),0)τ,

    G=diag(σ1,σ2,σ3,0) and

    M=(b11b12b130b21000b310b33b330α0α),

    the linearized system (4.2) can be equivalently rewritten as

    dP=MP(t)dt+GdB(t). (4.3)

    Since G is a constant matrix, according to Gardiner [27], we have

    G2+MΣ+ΣMT=0.

    In view of the independence of Brownian motions B1(t),B2(t) and B3(t), by the principle of finite independent superposition, it is easy to know that Eq (4.4) can be equivalently developed into the sum of the solution to the following three algebraic sub-equations:

    G2k+MΣk+ΣkMT=0,k=1,2,3,

    where G21=diag(σ21,0,0,0),G22=diag(0,σ22,0,0),G23=diag(0,0,σ23,0) and Σ=Σ1+Σ2+Σ3.

    Next it will be proved that M has all negative real-part eigenvalues. The characteristic polynomial of M is defined as

    ΨM(λ)=λ4+r1λ3+r2λ2+r3λ+r4,

    where

    r1=b11+b33+α>0,r2=b11b33b13b31+b12b21+α(b11+b33),r3=b12b21b33+α(b11b33b13b31+b12b21),r4=b21b31α(b12+b13)>0.

    By calculation, we can obtain

    b11b33b13b31>b13b31b13b31=0,

    which yields that r2 and r3>0. Furthermore, we can verify that r1r2r3>0andr1r2r3r23r21r4>0

    So, the matrix M is a Hurwitz matrix. Next, we will prove the definiteness of Σ by following three steps.

    Step 1. Consider the algebraic equation

    G21+MΣ1+Σ1MT=0. (4.4)

    Let M1=T1MT11, where

    T1=(100001000b31b21100001),M1=(b11b12+b13b31b21b130b210000b31b33b21b33b330α0α).

    Take M2=T2M1T12, where

    T2=(10000100001000b21αb31b331),M2=(b11b12+b13b31b21b130b210000b31b33b21b33+b21αb31b3300q1(b21+b31)αb31),

    where q1=b21α[b31b33(b21+b31)α]a231b33.

    Case 1.1: If q1=0, the standard transformation matrix T3 is given by

    T3=(b31b33b33(b31b33b21+α)(b33b21αb31)2b33(b33+α)0b31b33b21b33+b21αb31b3300100001).

    Then, we calculate

    Q1=T3M2T13=(τ11τ12τ13τ1410000100001(b21+b31)αb31),

    where

    τ11=b11+b33b21αb31,τ12=b12b21b13b31+b11b33b11b21αb31,τ13=b21(b12b33b13α+b12b21αb31),τ14=b33{b12b21b13b31+(b21+b31)α[b11b31+(b21+b31)α]b231}.

    Therefore the equivalent equation of (4.4) can be written as (T3T2T1)G21(T3T2T1)T+Q1[(T3T2T1)Σ1(T3T2T1)T]+[(T3T2T1)Σ1(T3T2T1)T]QT1=0. From Lemma 4.2, we have

    (T3T2T1)Σ1(T3T2T1)T=(b31b33σ21)¯Ω2,

    where

    ¯Ω2=(τ122(τ11τ12τ13)012(τ11τ12τ13)0012(τ11τ12τ13)0012(τ11τ12τ13)0τ112τ13(τ11τ12τ13)00000)

    is a symmetric positive semi-definite matrix. Hence,

    Σ1=(b31b33σ1)2(T3T2T1)1¯Ω2[(T3T2T1)1]T.

    Case 1.2: If q10, we can find standardized transformation T4 such that Q2=T4M2T14, where

    T4=(b31b33q1j1j2j30b31b33q1b21q1(b33+α)b33q1+(b21+b31)2α2b23100q1(b21+b31)αb310001),

    where j1=b31b33q1(b33+α)b21, j2=q1[b231b33(b33+q1)+b31(b21+b31)b33α+(b221+b21b31+b231)α2]b231, j3=b33q1(b33+α)(b21+b31)α[b231b33q1+(b21+b31)2α2]b331.

    Then, we have

    Q2=(r1r2r3r4100001000010).

    Likewise, (4.4) can be transformed into the following form: (T4T2T1)G21(T4T2T1)T+Q2[(T4T2T1)Σ1(T4T2T1)T]+[(T4T2T1)Σ1(T4T2T1)T]QT2=0. From Lemma 4.1, we have

    (T4T2T1)Σ1(T4T2T1)T=(b31b33q1σ1)2Ω1,

    where

    Ω1=(r2r3r1r42(r1r2r3r23r21r4)0r32(r1r2r3r23r21r4)00r32(r1r2r3r23r21r4)0r12(r1r2r3r23r21r4)r32(r1r2r3r23r21r4)0r12(r1r2r3r23r21r4)00r12(r1r2r3r23r21r4)0r1r2r32r4(r1r2r3r23r21r4)).

    is a positive definite symmetric matrix; hence,

    Σ1=(b31b33q1σ1)2(T4T2T1)1Ω1[(T4T2T1)1]T,

    is also a positive definite matrix.

    Step 2. Consider

    G22+MΣ2+Σ2MT=0. (4.5)

    Let M3=T5MT15, where

    T5=(0100100000100001),M3=(0b2100b12b11b1300b31b33b33α00α).

    Take M4=T6M3T16, where

    T6=(1000010000100αb1201),M4=(0b2100b12b11b1300b12b31+b33αb12b33b330q2b13αb12α),

    with q2=α(b11+α)b12.

    Case 2.1: If q2=0, we let Q3=T7M4T17, where the standard transformation matrix T7 is given by

    T7=(b13α(b12b31+b33α)b12j4j5j60b13α(b12b31+b33α)b212b13α(b33+α)b12α(b13b33b12+α)00b13αb12α0001),
    Q3=(r1r2r3r4100001000010),

    where j4=b13α(b11+b33+α)(b12b31+b33α)b212,j5=b13α(b13b31+b233+b33α+α2)b12,j6=α3+b13b33α(b33+2α)b12. Therefore, we have

    (T7T6T5)G22(T7T6T5)T+M[(T7T6T5)Σ2(T7T6T5)T]+[(T7T6T5)Σ2(T7T6T5)T]MT=0,

    where

    (T7T6T5)Σ2(T7T6T5)T=[b13α(b12b31+b33α)b12]2Ω1.

    Hence

    Σ2=α2b213(b12b31+b33α)2σ22b212(T7T6T5)1Ω1[(T7T6T5)1)]T,

    is a positive define matrix.

    Case 2.2: Consider that q20. Let M5=T8M4T18, where

    T8=(10000100001000b12q2b12b31+b33α1),M5=(0b2100b12b11b1300b12b31+b33αb12b33(1b12q2b12b31+b33α)b3300q3α+b12b33q2b12b31+b33α),

    where q3=b13αb12b12b33q2b12b31+b33α+b12q2(b12b33q2+b12b31α+b33α2)(b12b31+b33α)2.

    Case 2.2.1. If q20andq3=0, let Q4=T9M5T19, where

    T9=(b12b31+b33αj7j8j90b12b31+b33αb12b33(1b12q2b12b31+b33α)b3300100001),
    Q4=(τ21τ22τ23τ2410000100001α+b12b33q2b12b31+b33α),

    where j7=b11b31+b33(b31+q2)+b33α(b11+b33)b12, j8=b13(b31+b33αb12)+(b33+b12b33q2b12b31+b33α)2, j9=b33(b33+α),τ21=b11+b33(1+b12q2b12b31+b33α), τ22=(b12b21+b13b31b11b33+b13b33αb12b11b12b33q2b12b31+b33α), τ23=b12b21b33[b12(b31+q2)+b33α]b12b31+b33α, τ24=b33{α(b11+α)+b12[b21+b33q2(b11b12b31+b12b33q22b12b31α+b11b33α2b33α2)(b12b31+b33α)2]}. Therefore, we conclude that

    (T9T8T6T5)G22(T9T8T6T5)T+Q4[(T9T8T6T5)Σ2(T9T8T6T5)T]+[(T9T8T6T5)Σ2(T9T8T6T5)T]QT4=0.

    where

    (T9T8T6T5)Σ2(T9T8T6T5)T=[(b12b31+b33α)σ2]2~Ω2.

    Beside, it follows from Lemmas 4.2 that the specific form of the positive semi-definite matrix ~Ω2 is

    ~Ω2=(τ222(τ21τ22τ23)012(τ21τ23τ23)0012(τ21τ22τ23)0012(τ21τ22τ23)0τ212τ23(τ21τ22τ23)00000).

    Therefore, we conclude that

    Σ2=(b12b31+b33α)2σ22(T9T8T6T5)1~Ω2[(T9T8T6T5)1]T,

    is semi-positive definite.

    Case 2.2.2. If q20,q30, we let Q5=T10M5T110, where

    T10=(q3(b12b31+b33α)j10j11j120q3(b12b31+b33α)b12q3(b33+α)b33q3+(αb12b33q2b12b31+b33α)200q3α+b12b33q2b12b31+b33α0001),
    Q5=(r1r2r3r4100001000010),

    with

    j10=q3(b11+b33+α)(b12b31+b33α)b12,j11=q3{α2+b33(b33+q3+α)+b13(b31+b33αb12)+b12b33q2[b12b33(b31+q2)b12b31α+b33(b33α)α](b12b31+b33α)2},j12=b33q3(b33+α)+(α+b12b33q2b12b31+b33α)[b33q3+(αb12b33q2b12b31+b33α)2].

    Therefore, we have

    (T10T8T6T5)G22(T10T8T6T5)T+Q5[(T10T8T6T5)Σ2(T10T8T6T5)T]+[(T10T8T6T5)Σ2(T10T8T6T5)T]QT5=0,

    where

    (T10T8T6T5)Σ2(T10T8T6T5)T=q23(b12b31+b33α)2σ22Ω1.

    Hence

    Σ2=q23(b12b31+b33α)2σ22(T10T8T6T5)1Ω1[(T10T8T6T5)1]T,

    is also a positive definite matrix.

    Step 3. Consider

    G23+MΣ3+Σ3MT=0. (4.6)

    We let M6=T11MT111, where

    T11=(0010100001000001),M6=(b33b310b33b13b11b1200b210000αα).

    Let Q6=T12M6T112, where the standard transformation matrix T12 is given by

    T12=(b13b21αb21α(b11+α)b12b21α+α3α30b21αα2α200αα0001).

    By direct calculation, we obtain that

    Q6=(r1r2r3r4100001000010).

    Thus (4.6) can be transformed into the following from:

    (T12T11)G23(T12T11)T+Q6[(T12T11)Σ3(T12T11)T]+[(T12T11)Σ3(T12T11)T]QT6=0,

    where

    (T12T11)Σ3(T12T11)T=(b13b21ασ3)2Ω1.

    Thus,

    Σ3=(b13b21ασ3)2(T12T11)1Ω1[(T12T11)1]T,

    is a positive definite matrix.

    Summing up the above steps, we have Σ=Σ1+Σ2+Σ3 is a positive definite matrix. The solution (x(t),y(t),z(t),w(t)) of system (1.4) has a log-normal probability density function

    Φ(x,y,z,w)=(2π)2Σ12e12(lnxx,lnyy,lnzz,lnww)Σ1(lnxx,lnyy,lnzz,lnww)T.

    Theorem 5.1. Let (x(t),y(t),z(t),w(t)) be the solution of the system (1.4) with any initial value (x(0),y(0),z(0),w(0))R4+. If RE0=b1D1(D2+12σ22)<1, then

    lim supt+lny(t)t(D2+12σ22)(RE01)<0,a.s.,

    which means the extinction of y(t).

    Proof. Consider the following SDE:

    dˆx(t)=[1D1ˆx(t)]dt+σ1ˆx(t)dB1(t),

    with the same initial value ˆx(0)=x(0)>0, making use of the stochastic comparison theorem, we have

    x(t)ˆx(t),a.s.

    Moreover,

    limt+1tt0ˆx(s)ds=+0sh(s)ds=1D1,

    where h(s)=[(2σ21)1+2D1σ21Γ1(1+2D1σ21)]s2(1+2D1σ21)e2σ21s,s>0. Applying Itô's formula to lny(t), we have

    dlny=[b1x(D2+12σ22)]dt+σ2dB2(t)][b1ˆx(D2+12σ22)]dt+σ2dB2(t)], (5.1)

    Integrating (5.1) from 0 to t on both sides, one can see that

    lny(t)tlny(0)tb11tt0ˆx(s)ds(D2+12σ22)+t0σ2dB2(s)t. (5.2)

    Applying the strong law of large numbers yields

    limt+t0σ2dB2(s)t=0.

    Next, we take the superior limit on both sides of (5.2)

    lim supt+lny(t)tb1lim supt+1tt0ˆx(s)ds(D2+12σ22)=b1D1(D2+12σ22)=(D2+12σ22)(RE01)<0.

    which implies that

    limt+y(t)=0,a.s.

    In this section, employing the Milstein higher-order method, we will present numerical simulations to verify the theoretical results. The discretization equations of model (1.4) are given by

    {xk+1=xk+(1a1xkyka2xkzkD1xk)t+σ1xktη1,k+σ21xk2(η21,k1)t,yk+1=yk+(b1xkykD2yk)t+σ2yktη2,k+σ22yk2(η22,k1)t,zk+1=zk+(c1wkc2xkzkD3zk)t+σ3zktη3,k+σ23zk2(η23,k1)t,wk+1=wk+(αykαwk)t, (6.1)

    where ηi,k(i=1,2,3;k=1,2,n) denotes independent Gaussian random variables which follow the distribution N(0,1).

    Example 6.1. Let us choose a1=1,a2=4,D1=1.01,D2=8,D3=1.02,b1=20,c1=5,c2=5,α=0.09,σ1=σ2=0.01,σ3=0.05 and the initial value (x(0),y(0),z(0),w(0))=(0.9,0.1,0.1,0.1). We conclude that R0=b1D1D2=2.4752>1 and Rs0=b1(D1+12σ21)(D2+12σ22)=2.4751>1. Theorem 4.1 shows that there exists an ergodic stationary distribution of stochastic model (1.4). Noting that (x,y,z,w)=(0.4,0.1955,0.3236,0.1955), we have following the covariance matrix

    Σ=(1.3643e044.7814e062.3417e041.0541e044.7814e068.7590e036.2162e045.2528e052.3417e046.2162e047.5113e042.3141e041.0541e045.2528e052.3141e042.4051e04).

    Thus, we have the corresponding probability density function

    Φ(x,y,z,w)=1.0882×105×e12(lnx2.5,lny5.1142,lnz3.09,lnw5.1142)Σ1(lnx2.5,lny5.1142,lnz3.09,lnw5.1142)T.

    As a result, Φ(x,y,z,w) has the following four marginal density functions:

    Φx=1x2πw11e(lnxlnx)22w11=134.1553xe(lnx+0.9163)22.7286e04,Φy=1y2πw22e(lnylny)22w22=10.2346ye(lny+1.632)20.0175,Φz=1z2πw33e(lnzlnz)22w33=10.0687ze(lnz+1.1282)20.0015,Φw=1w2πw44e(lnwlnw)22w44=10.0389we(lnw+1.632)24.8101e04;

    we can conclude that system (1.2) admits a global positive stationary solution on R4+; see Figure 1 and Figure 2.

    Figure 1.  Left-hand column presents the numbers of x,y,z and w for system (6.1) with (σ1,σ2,σ3)=(0.01,0.01,0.05), and its deterministic system, respectively. Right-hand columns shows the frequency histograms and corresponding marginal density function curves of x,y,z and w, respectively.
    Figure 2.  The marginal density functions and frequency fitting curves for x,y and z in system (1.4).

    Example 6.2. Let us choose b1=20,D1=2.3,D2=8andσ2=1.8. Then the condition RE0=0.693<1 is satisfied. Theorem 5.1 shows that MC-degrading bacteria of system (1.4) will go to extinction with probability one, which is numerically confirmed by Figure 3.

    Figure 3.  Corresponding numbers for solution (x(t),y(t),z(t)) to system (1.4) with random perturbations (σ1,σ2,σ3)=(0.01,1.8,0.05) and main parameters (b1,D1,D2)=(20,2.3,8).

    In the current paper, we consider a stochastic Microcystins degradation model with distributed delay. We have established sufficient conditions for the existence of an ergodic stationary distribution of the positive solutions to system (1.4) by constructing a suitable stochastic Lyapunov function. The result shows that a small amount of white noise can guarantee the existence of an ergodic stationary distribution of the positive solutions to system (1.4). In addition, we have obtained the exact probability density function around a quai-equilibrium point.

    Some interesting topics deserve further consideration. On the one hand, the paper focuses on the dynamical evolution and stability of the system (1.4). It should be noted that model (1.4) can potentially be solved analytically by using the Lie algebra method [25,26]. On the other hand, because our model is autonomous and only disturbed by white noise, it would be interesting to introduce telegraph noise, such as a continuous-time Markov chain, into system (1.3). Moreover, it would also be interesting to study more complicated MCs degradation models, such as multi-group MCs degradation models. We will leave these problems as our future work.

    The authors declare that they have not used artificial intelligence tools in the creation of this article.

    Thanks to the reviewers for their helpful comments and suggestions. This work was supported by the Hainan Provincial Natural Science Foundation of China (No. 122RC679,121RC554), Talent Program of Hainan Medical University (No. XRC2020030) and Northeast Petroleum University Special Research Team Project (No. 2022TSTD-05).

    The authors declare that there is no conflict of interest.



    [1] J. Li, R. Li, J. Li, Current research scenario for microcystins biodegradation—a review on fundamental knowledge application prospects and challenges, Sci. Total Environ., 595 (2017), 615–632. https://doi.org/10.1016/j.scitotenv.2017.03.285 doi: 10.1016/j.scitotenv.2017.03.285
    [2] J. Li, K. Shimizu, H. Maseda, Z. Lu, M. Utsumi, Z. Zhang, et al., Investigations into the biodegradation of microcystin-LR mediated by the biofilm in wintertime from a biological treatment facility in a drinking-water treatment plant, Bioresour. Technol., 106 (2012), 27–35. https://doi.org/10.1016/j.biortech.2011.11.099 doi: 10.1016/j.biortech.2011.11.099
    [3] K. Shimizu, H. Maseda, K. Okano, T. Kurashima, Y. Kawauchi, Q. Xue, et al., Enzymatic pathway for biodegrading microcystin LR in sphingopyxis sp. C-1, J. Biosci. Bioeng., 114 (2012), 630–634. https://doi.org/10.1016/j.jbiosc.2012.07.004 doi: 10.1016/j.jbiosc.2012.07.004
    [4] K. Song, W. Ma, K. Guo, Global behavior of a dynamic model with biodegradation of microcystins, J. Appl. Anal. Comput., 9 (2019), 1261–1276. https://doi.org/10.11948/2156-907X.20180215 doi: 10.11948/2156-907X.20180215
    [5] K. Ge, K. Song, W. Ma, Existence of positive periodic solutions of a delayed periodic Microcystins degradation model with nonlinear functional responses, Appl. Math. Lett., 131 (2022), 108056. https://doi.org/10.1016/j.aml.2022.108056 doi: 10.1016/j.aml.2022.108056
    [6] J. Caperon, Time lag in population growth response of isochrysis galbana to a variable nitrate environment, Ecology, 50 (1969), 188–192. https://doi.org/10.2307/1934845 doi: 10.2307/1934845
    [7] M. Droop, Vitamin B12 and marine ecology. IV. The kinetics of uptake, growth and inhibition in monochrysis lutheri, J. Mar. Biol. Assoc. UK, 48 (1968), 689–733. https://doi.org/10.1017/S0025315400019238 doi: 10.1017/S0025315400019238
    [8] B. Li, G. Wolkowicz, Y. Kuang, Global asymptotic behavior of a chemostat model with two perfectly complementary resources and distributed delay, SIAM J. Appl. Math., 60 (2000), 2058–2086. https://doi.org/10.1137/S0036139999359756 doi: 10.1137/S0036139999359756
    [9] K. Song, W. Ma, Z. Jiang, Bifurcation analysis of modeling biodegradation of microcystins, Int. J. Biomath., 12 (2019), 1950028. https://doi.org/10.1142/S1793524519500281 doi: 10.1142/S1793524519500281
    [10] N. Macdonald, Time lags in biological models, in Lecture Notes in Biomathematics, Springer-Verlag, Heidelberg, 1978. https://doi.org/10.1007/978-3-642-93107-9
    [11] T. C. Gard, Persistence in stochastic food web models, Bull. Math. Biol., 46 (1984), 357–370. https://doi.org/10.1007/bf02462011 doi: 10.1007/bf02462011
    [12] M. Bandyopadhyay, J. Chattopadhyay, Ratio-dependent predator–prey model: Effect of environmental fluctuation and stability, Nonlinearity, 18 (2005), 913–936. https://doi.org/10.1088/0951-7715/18/2/022 doi: 10.1088/0951-7715/18/2/022
    [13] I. Addition, Extinction and persistence in mean of a novel delay impulsive stochastic infected predator–prey system with jumps, Complexity, 2017 (2017), 1950970. https://doi.org/10.1155/2017/1950970 doi: 10.1155/2017/1950970
    [14] X. Meng, S. Zhao, T. Feng, T. Zhang, Dynamics of a novel nonlinear stochastic SIS epidemic model with double epidemic hypothesis, J. Math. Anal. Appl., 433 (2016), 227–242. https://doi.org/10.1016/j.jmaa.2015.07.056 doi: 10.1016/j.jmaa.2015.07.056
    [15] X. Leng, T. Feng, X. Meng, Stochastic inequalities and applications to dynamics analysis of a novel SIVS epidemic model with jumps, J. Inequal. Appl., 2017 (2017), 138. https://doi.org/10.1186/s13660-017-1418-8 doi: 10.1186/s13660-017-1418-8
    [16] R. M. May, Stability and Complexity in Model Ecosystems, Princeton University Press, 2001.
    [17] B. Han, D. Jiang, T. Hayat, A. Alsaedi, B. Ahmad, Stationary distribution and extinction of a stochastic staged progression AIDS model with staged treatment and second-order perturbation, Chaos Solitons Fractals, 140 (2020), 110238. https://doi.org/10.1016/j.chaos.2020.110238 doi: 10.1016/j.chaos.2020.110238
    [18] B. Zhou, B. Han, D. Jiang, Ergodic property, extinction and density function of a stochastic SIR epidemic model with nonlinear incidence and general stochastic perturbations, Chaos Solitons Fractals, 152 (2021), 111338. https://doi.org/10.1016/j.chaos.2021.111338 doi: 10.1016/j.chaos.2021.111338
    [19] 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
    [20] S. Zhang, X. Meng, X. Wang, Application of stochastic inequalities to global analysis of a nonlinear stochastic SIRS epidemic model with saturated treatment function, Adv. Differ. Equations, 2018 (2018), 50. https://doi.org/10.1186/s13662-018-1508-z doi: 10.1186/s13662-018-1508-z
    [21] 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
    [22] R. Z. Khas'minskii, Stochastic Stability of Differential Equations, Sijthoff Noordhoff, Alphen ana den Rijn, Netherlands, 1980. https://doi.org/10.1007/978-94-009-9121-7
    [23] D. Xu, Y. Huang, Z. Yang, Existence theorems for periodic Markov process and stochastic functional differential equations, Discrete Contin. Dyn. Syst., 24 (2009), 1005–1023. https://doi.org/10.3934/dcds.2009.24.1005 doi: 10.3934/dcds.2009.24.1005
    [24] J. Ge, W. Zuo, D. Jiang, Stationary distribution and density function analysis of a stochastic epidemic HBV model, Math. Comput. Simul., 191 (2022), 232–255. https://doi.org/10.1016/j.matcom.2021.08.003 doi: 10.1016/j.matcom.2021.08.003
    [25] Y. Shang, Lie algebraic discussion for affinity based information diffusion in social networks, Open Phys., 15 (2017), 705–711. https://doi.org/10.1515/phys-2017-0083 doi: 10.1515/phys-2017-0083
    [26] Y. Shang, Analytical solution for an in-host viral infection model with time-inhomogeneous rates, Acta Phys. Pol. B, 46 (2015), 1567–1577. https://doi.org/10.5506/aphyspolb.46.1567 doi: 10.5506/aphyspolb.46.1567
    [27] C. W. Gardiner, Handbook of Stochastic Methods for Physics, Springer, 1983.
  • Reader Comments
  • © 2024 the Author(s), licensee AIMS Press. This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0)
通讯作者: 陈斌, bchen63@163.com
  • 1. 

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

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

Metrics

Article views(1442) PDF downloads(53) Cited by(0)

Figures and Tables

Figures(3)

Other Articles By Authors

/

DownLoad:  Full-Size Img  PowerPoint
Return
Return

Catalog