
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
[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:
{x′1(t)=Da10−a12x1(t)x2(t)−a13x1(t)x3(t)−(D+d1)x1(t),x′2(t)=a21x2(t)x1(t)−a20x2(t)−(D+d2)x2(t),x′3(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; a12≥0 is the consumption rate of MCs; a21≥0 is the maximal growth rate of MC-degrading bacteria; a20≥0 is the consumption rate of MC-degrading bacteria; a30≥0 is the maximal growth rate of degrading enzymes; a13≥0 and a31≥0 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)=1−a1x(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)=1−a1x(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)=[c1∫t−∞K(t−s)y(s)ds−c2x(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)=[1−a1x(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)=[c1∫t−∞K(t−s)y(s)ds−c2x(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)=∫t−∞K(t−s)y(s)ds=∫t−∞αe−α(t−s)y(s)ds. |
Then, by using the linear chain technique, system (1.3) is transformed into the following equivalent system:
{dx(t)=[1−a1x(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 t≥0 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 k0≥1 be sufficiently large so that x(0),y(0),z(0)andw(0) all lie within the interval [1k0,k0]. For each integer k≥k0, 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 t≥0. 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 k1≥k0 such that
P{τk≤T}>εfor all k≥k1. | (2.1) |
Define the nonnegative C2− Lyapunov function V0 as follows
V0(x,y,z,w)=x−1−lnx+ay−lny−(1+lna)+cz−dlnz−(d+lncd)+bw−lnw−(1+lnb), |
where a,b,c,d are positive constants which will be determined later.
Based on the basic inequality u−1−lnu≥0, we have
cz−dlnz−(d+lncd)=d(czd−1−lnczd)≥0. |
Applying Itô's formula, we have
dV0=LV0dt+σ1(x−1)dB1(t)+σ2(ay−1)dB2(t)+σ3(cz−d)dB3(t), | (2.2) |
LV0=1−a1xy−a2xz−D1x+a(b1xy−D2y)+c(c1w−c2xz−D3z)+b(αy−αw)−1x+a1y+a2z+D1+12σ21−b1x+D2+12σ22+d(−c1wz+c2x+D3+12σ23)−αyw+α≤1+D1+D2+dD3+α+12(σ21+σ22+dσ23)+(ab1−a1)xy+(a2−cD3)z+(−D1−b1+dc2)x+(−aD2+a1+bα)y+(cc1−bα)w. |
Choose b=a2c1αD3,c=a2D3,d=b+D1c and a=min{a1b1,a1D3+c1a2D2D3} such that ab1−a1≤0and a1−aD2+αb≤0, then, we have
LV0≤1+D1+D2+dD3+α+12(σ21+σ22+dσ23):=K1, |
where K1 is a positive constant. Integrating equation (2.2) from 0to τk∧T and then taking the expectation on both sides, we have
EV0(x(τk∧T),y(τk∧T),z(τk∧T),w(τk∧T)))≤V0(x(0),y(0),z(0),w(0))+K1E(τk∧T)≤V0(x(0),y(0),z(0),w(0))+K1T | (2.3) |
Set Ωk={τk≤T}for k≥k1 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
(k−1−lnk)∧(ak−lnk−(1+lna))∧d(cdk−1−lnckd)∧(bk−lnk−(1+lnb))or(1k−1−ln1k)∧(a1k−ln1k−(1+lna))∧d(cd1k−1−lncdk)∧(b1k−ln1k−(1+lnb))=(k−1−lnk)∧(ak−lnk−(1+lna))∧d(cdk−1−ln(ck)+lnd)∧(bk−lnk−(1+lnb))or(1k−1+lnk)∧(a1k+lnk−(1+lna))∧d(cd1k−1−lnc+ln(dk))∧(b1k+lnk−(1+lnb)) |
Therefore, we obtain
V0(x(τk,ω),y(τk,ω),z(τk,ω),w(τk,ω)))≥(k−1−lnk)∧(ak−lnk−(1+lna))∧d(cdk−1−ln(ck)+lnd)∧(bk−lnk−(1+lnb))∧(1k−1+lnk)∧(a1k+lnk−(1+lna))∧d(cd1k−1−lnc+ln(dk))∧(b1k+lnk−(1+lnb)) |
It follows from Eq (2.3) that
V0(x(0),y(0),z(0),w(0))+K1T≥E[IΩk(ω)V0(x(τk,ω),y(τk,ω),z(τk,ω),w(τk,ω)))]≥ϵ(k−1−lnk)∧(ak−lnk−(1+lna))∧d(cdk−1−ln(ck)+lnd)∧(bk−lnk−(1+lnb))∧(1k−1+lnk)∧(a1k+lnk−(1+lna))∧d(cd1k−1−lnc+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+k∑r=1∫tt0σ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)|+k∑r=1|σr(s,x)−σr(s,y)|≤B|x−y|,|b(s,x)|+k∑r=1|σr(s,x)|≤B(1+|x|), | (3.2) |
in UR⊂Rd+ 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 U⊂R4+ such that
L˜V(x,y,z,w)≤−1,for any (x,y,t,w)∈R4+∖U. |
Define
V1(x,y,z,w)=−lnx−e1lny+a2D3z+c1a2αD3w, |
where e1 is a positive constant to be chosen later. Then applying Itô's formula to V1, we have
LV1≤−1x−e1b1x+e1(D2+12σ22)+D1+12σ21+a1y+a2z+a2D3(c1w−c2xz−D3z)+c1a2D3α(αy−αw)≤−2√b1e1+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
LV1≤−b1D2+12σ22+D1+12σ21+(a1+c1a2D3)y=−(D1+12σ21)(Rs0−1)+(a1+c1a2D3)y=:−λ+(a1+c1a2D3)y, | (3.4) |
where
Rs0=b1(D1+12σ21)(D2+12σ22),λ=(D1+12σ21)(Rs0−1). |
Define
V2=−lnx−lnw−lnz. |
Then applying Itô's formula to V2, one gets
LV2≤−1x−αyw−c1wz+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{1−D1x−a1D22b1y−a1D2D34b1c1z−a1D24b1w}+θ+12(x+a1b1y+a1D24b1c1z+a1D22αb1w)θ[σ21x2+σ22(a1b1)2y2+σ23(a1D24c1b1)2z2]}≤(x+a1b1y+a1D24b1c1z+a1D22αb1w)θ+1−D1xθ+2−D22(a1b1)θ+2yθ+2−D3(a1D24b1c1)θ+2zθ+2−α2(a1D22αb1)θ+2wθ+2+θ+12(x+a1b1y+a1D24b1c1z+a1D22αb1w)θ[σ21x2+σ22(a1b1)2y2+σ23(a1D24c1b1)2z2]≤−D1θxθ+2−D2θ2(a1b1)θ+2yθ+2−D3θ(a1D24b1c1)θ+2zθ+2−α4(a1D22αb1)θ+2wθ+2+B0, | (3.6) |
where
B0=sup(x,y,z,w)∈R4+{(x+a1b1y+a1D24b1c1z+a1D22αb1w)θ+1−D1(1−θ)xθ+2−D2(1−θ)2(a1b1)θ+2yθ+2−D3(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+c0≤−2. Also,
c0=sup(x,y,z,w)∈R4+{−D1θxθ+2−D2θ2(a1b1)θ+2yθ+2−D3θ(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˜V≤−M0λ+[M0(a1+a2c1D3)+a1]y+a2z+c2x+D1+12σ21+D3+12σ23+α−D1θxθ+2−D2θ2(a1b1)θ+2yθ+2−D3θ(a1D24b1c1)θ+2zθ+2−α4(a1D22αb1)θ+2wθ+2+B0−1x−αyw−c1wz. |
Next, we define the following bounded closed set
Uε={ε≤x≤1ε,ε≤y≤1ε,ε2≤w≤1ε3,ε3≤z≤1ε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˜V≤−M0λ+[M0(a1+a2c1D3)+a1]y+a2z+c2x−D1θxθ+2−D2θ2(a1b1)θ+2yθ+2−D3θ(a1D24b1c1)θ+2zθ+2−α4(a1D22αb1)θ+2wθ+2+B0+D1+12σ21+D3+12σ23+α≤−M0λ+{M0(a1+a2c1D3)+a1]ε+c0≤−1; |
here, ε is a positive constant small enough to satisfy the following conditions:
−M0λ+[M0(a1+a2c1D3)+a1]ε+c0≤−1, | (3.8) |
−1ε+E≤−1, | (3.9) |
−αε+E≤−1, | (3.10) |
−c1ε+E≤−1, | (3.11) |
−D2θ4(a1b1ε)θ+2+E≤−1, | (3.12) |
−D1θ2(1ε)θ+2+E≤−1, | (3.13) |
−D3θ2(a1D24b1c1ε3)θ+2+E≤−1, | (3.14) |
−α8(a1D22αb1ε2)θ+2+E≤−1, | (3.15) |
where
E=sup(x,y,z,w)∈R4+{−D1θ2xθ+2−D2θ4(a1b1)θ+2yθ+2−D3θ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˜V≤−1x+[M0(a1+a2c1D3)+a1]y+a2z+c2x−D1θxθ+2−D2θ2(a1b1)θ+2yθ+2−D3θ(a1D24b1c1)θ+2zθ+2−α4(a1D22αb1)θ+2wθ+2+B0+D1+12σ21+D3+12σ23+α≤−1ε+E≤−1. |
Case 3 : If (x,y,z,w)∈U3, by (3.10) we have
L˜V≤−αyw+E≤−αε+E≤−1. |
Case 4 : If (x,y,z,w)∈U4, according to (3.11), we deduce that
L˜V≤−c1wz+E≤−c1ε+E≤−1. |
Case 5 : If (x,y,z,w)∈U5, by condition (3.12), we conclude that
L˜V≤−D2θ4(a1b1)θ+2yθ+2+E≤−D2θ4(a1b1ε)θ+2+E≤−1. |
Case 6 : If (x,y,z,w)∈U6, and the condition (3.13) is satisfied, it follows that
L˜V≤−D1θ2xθ+2+E≤−D1θ2(1ε)θ+2+E≤−1. |
Case 7 : If (x,y,z,w)∈U7, in view of (3.14) we get
L˜V≤−D3θ2(a1D24b1c1)θ+2zθ+2+E≤−D3θ2(a1D24b1c1ε3)θ+2+E≤−1. |
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+E≤−1. |
Clearly, there exists a small enough constant ε to make the conclusion
L˜V≤−1,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=[e−u1−a1eu2−a2eu3−(D1+12σ21)]dt+σ1dB1(t),du2=[b1eu1−(D2+12σ22)]dt+σ2dB2(t),du3=[c1eu4−u3−c2eu1−(D3+12σ23)]dt+σ3dB3(t),du4=(αeu2−u4−α)dt. | (4.1) |
Assume that Rs0>1, then, there is the quasi-endemic equilibrium E∗=(x∗,y∗,z∗,w∗)=(eu∗1,eu∗2,eu∗3,eu∗4) 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)0−12(τ1τ2−τ3)0012(τ1τ2−τ3)00−12(τ1τ2−τ3)0τ12τ3(τ1τ2−τ3)00000). |
Denote P=(p1,p2,p3,p4)T=(u1−u∗1,u2−u∗2,u3−u∗3,u4−u∗4)T(k=1,2,3,4). Therefore, the linearized equation system of Eq (4.1) is obtained as follows
{dp1=(−b11p1−b12p2−b13p3)dt+σ1dB1(t),dp2=b21p1dt+σ2dB2(t),dp3=(−b31p1−b33p3+b33p4)dt+σ3dB3(t),dp4=(αp2−αp4)dt, | (4.2) |
where
b11=1x∗>b12+b13,b12=a1y∗,b13=a2z∗,b21=b1x∗,b31=c2x∗,b33=c1w∗z∗>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⏐Σ⏐−12e−12(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 q1≠0, |
Σ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 q2≠0,q3=0,(b12b31+b33α)2q23σ22(T10T8T6T5)−1Ω1[(T10T8T6T5)−1]T,if q2≠0,q3≠0, |
Σ3=(−b13b21ασ3)2(T12T11)−1Ω1[(T12T11)−1]T, |
where q1=b21α[b13b33−(b21+b33)α]b231b33b33,q2=α(−b11+α)b12,q3=−b13αb12−b12b33q2b12b31+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=(−b11−b12−b130b21000−b310−b33b330α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=b11b33−b13b31+b12b21+α(b11+b33),r3=b12b21b33+α(b11b33−b13b31+b12b21),r4=b21b31α(b12+b13)>0. |
By calculation, we can obtain
b11b33−b13b31>b13b31−b13b31=0, |
which yields that r2 and r3>0. Furthermore, we can verify that r1r2−r3>0andr1r2r3−r23−r21r4>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=T1MT−11, where
T1=(100001000b31b21100001),M1=(−b11−b12+b13b31b21b130b210000b31b33b21−b33b330α0−α). |
Take M2=T2M1T−12, where
T2=(10000100001000−b21αb31b331),M2=(−b11−b12+b13b31b21−b130b210000b31b33b21−b33+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+α)(b33−b21αb31)2−b33(b33+α)0b31b33b21−b33+b21αb31b3300100001). |
Then, we calculate
Q1=T3M2T−13=(−τ11−τ12−τ13−τ1410000100001−(b21+b31)αb31), |
where
τ11=b11+b33−b21αb31,τ12=b12b21−b13b31+b11b33−b11b21αb31,τ13=−b21(−b12b33−b13α+b12b21αb31),τ14=−b33{b12b21−b13b31+(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)0−12(τ11τ12−τ13)0012(τ11τ12−τ13)00−12(τ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 q1≠0, we can find standardized transformation T4 such that Q2=T4M2T−14, where
T4=(b31b33q1j1j2j30b31b33q1b21−q1(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=(−r1−r2−r3−r4100001000010). |
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=(r2r3−r1r42(r1r2r3−r23−r21r4)0−r32(r1r2r3−r23−r21r4)00r32(r1r2r3−r23−r21r4)0−r12(r1r2r3−r23−r21r4)−r32(r1r2r3−r23−r21r4)0r12(r1r2r3−r23−r21r4)00−r12(r1r2r3−r23−r21r4)0r1r2−r32r4(r1r2r3−r23−r21r4)). |
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=T5MT−15, where
T5=(0100100000100001),M3=(0b2100−b12−b11−b1300−b31−b33b33α00−α). |
Take M4=T6M3T−16, where
T6=(1000010000100αb1201),M4=(0b2100−b12−b11−b1300−b12b31+b33αb12−b33b330q2−b13αb12−α), |
with q2=α(−b11+α)b12.
Case 2.1: If q2=0, we let Q3=T7M4T−17, where the standard transformation matrix T7 is given by
T7=(−b13α(b12b31+b33α)b12j4j5j60b13α(b12b31+b33α)b212b13α(b33+α)b12α(−b13b33b12+α)00−b13αb12−α0001), |
Q3=(−r1−r2−r3−r4100001000010), |
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 q2≠0. Let M5=T8M4T−18, where
T8=(10000100001000b12q2b12b31+b33α1),M5=(0b2100−b12−b11−b1300−b12b31+b33αb12b33(−1−b12q2b12b31+b33α)b3300q3−α+b12b33q2b12b31+b33α), |
where q3=−b13αb12−b12b33q2b12b31+b33α+b12q2(−b12b33q2+b12b31α+b33α2)(b12b31+b33α)2.
Case 2.2.1. If q2≠0andq3=0, let Q4=T9M5T−19, where
T9=(b12b31+b33αj7j8j90−b12b31+b33αb12b33(−1−b12q2b12b31+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+b13b31−b11b33+b13b33αb12−b11b12b33q2b12b31+b33α), τ23=b12b21b33[b12(b31+q2)+b33α]b12b31+b33α, τ24=−b33{α(−b11+α)+b12[b21+b33q2(b11b12b31+b12b33q2−2b12b31α+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)0−12(τ21τ23−τ23)0012(τ21τ22−τ23)00−12(τ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 q2≠0,q3≠0, we let Q5=T10M5T−110, where
T10=(q3(b12b31+b33α)j10j11j120−q3(b12b31+b33α)b12−q3(b33+α)b33q3+(α−b12b33q2b12b31+b33α)200q3−α+b12b33q2b12b31+b33α0001), |
Q5=(−r1−r2−r3−r4100001000010), |
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=T11MT−111, where
T11=(0010100001000001),M6=(−b33−b310b33−b13−b11−b1200b210000α−α). |
Let Q6=T12M6T−112, 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=(−r1−r2−r3−r4100001000010). |
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⏐Σ⏐−12e−12(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)(RE0−1)<0,a.s., |
which means the extinction of y(t).
Proof. Consider the following SDE:
dˆx(t)=[1−D1ˆ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→+∞1t∫t0ˆx(s)ds=∫+∞0s⋅h(s)ds=1D1, |
where h(s)=[(2σ21)1+2D1σ21Γ−1(1+2D1σ21)]s−2(1+2D1σ21)e−2σ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)t≤lny(0)t−b11t∫t0ˆ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)t≤b1lim supt→+∞1t∫t0ˆx(s)ds−(D2+12σ22)=b1D1−(D2+12σ22)=(D2+12σ22)(RE0−1)<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+(1−a1xkyk−a2xkzk−D1xk)△t+σ1xk√△tη1,k+σ21xk2(η21,k−1)△t,yk+1=yk+(b1xkyk−D2yk)△t+σ2yk√△tη2,k+σ22yk2(η22,k−1)△t,zk+1=zk+(c1wk−c2xkzk−D3zk)△t+σ3zk√△tη3,k+σ23zk2(η23,k−1)△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.3643e−04−4.7814e−06−2.3417e−04−1.0541e−04−4.7814e−068.7590e−03−6.2162e−045.2528e−05−2.3417e−04−6.2162e−047.5113e−042.3141e−04−1.0541e−045.2528e−052.3141e−042.4051e−04). |
Thus, we have the corresponding probability density function
Φ(x,y,z,w)=1.0882×105×e−12(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=1x√2πw11e−(lnx−lnx∗)22w11=134.1553xe−(lnx+0.9163)22.7286e−04,∂Φ∂y=1y√2πw22e−(lny−lny∗)22w22=10.2346ye−(lny+1.632)20.0175,∂Φ∂z=1z√2πw33e−(lnz−lnz∗)22w33=10.0687ze−(lnz+1.1282)20.0015,∂Φ∂w=1w√2πw44e−(lnw−lnw∗)22w44=10.0389we−(lnw+1.632)24.8101e−04; |
we can conclude that system (1.2) admits a global positive stationary solution on R4+; see Figure 1 and Figure 2.
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.
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. |