
Citation: Zhiwei Huang, Gang Huang. Mathematical analysis on deterministic and stochastic lake ecosystem models[J]. Mathematical Biosciences and Engineering, 2019, 16(5): 4723-4740. doi: 10.3934/mbe.2019237
[1] | Yansong Pei, Bing Liu, Haokun Qi . Extinction and stationary distribution of stochastic predator-prey model with group defense behavior. Mathematical Biosciences and Engineering, 2022, 19(12): 13062-13078. doi: 10.3934/mbe.2022610 |
[2] | Qi Zhou, Huaimin Yuan, Qimin Zhang . Dynamics and approximation of positive solution of the stochastic SIS model affected by air pollutants. Mathematical Biosciences and Engineering, 2022, 19(5): 4481-4505. doi: 10.3934/mbe.2022207 |
[3] | Yajuan Guo, Zijian Liu, Yuanshun Tan, Yawei Liu . Modeling and analysis of a stochastic giving-up-smoking model with quit smoking duration. Mathematical Biosciences and Engineering, 2023, 20(12): 20576-20598. doi: 10.3934/mbe.2023910 |
[4] | Dwi Lestari, Noorma Yulia Megawati, Nanang Susyanto, Fajar Adi-Kusumo . Qualitative behaviour of a stochastic hepatitis C epidemic model in cellular level. Mathematical Biosciences and Engineering, 2022, 19(2): 1515-1535. doi: 10.3934/mbe.2022070 |
[5] | Tingting Xue, Long Zhang, Xiaolin Fan . Dynamic modeling and analysis of Hepatitis B epidemic with general incidence. Mathematical Biosciences and Engineering, 2023, 20(6): 10883-10908. doi: 10.3934/mbe.2023483 |
[6] | Yiping Tan, Yongli Cai, Zhihang Peng, Kaifa Wang, Ruoxia Yao, Weiming Wang . Dynamics of a stochastic HBV infection model with drug therapy and immune response. Mathematical Biosciences and Engineering, 2022, 19(8): 7570-7585. doi: 10.3934/mbe.2022356 |
[7] | Shengnan Zhao, Sanling Yuan . A coral reef benthic system with grazing intensity and immigrated macroalgae in deterministic and stochastic environments. Mathematical Biosciences and Engineering, 2022, 19(4): 3449-3471. doi: 10.3934/mbe.2022159 |
[8] | 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 |
[9] | David F. Anderson, Tung D. Nguyen . Results on stochastic reaction networks with non-mass action kinetics. Mathematical Biosciences and Engineering, 2019, 16(4): 2118-2140. doi: 10.3934/mbe.2019103 |
[10] | Ying He, Yuting Wei, Junlong Tao, Bo Bi . Stationary distribution and probability density function analysis of a stochastic Microcystins degradation model with distributed delay. Mathematical Biosciences and Engineering, 2024, 21(1): 602-626. doi: 10.3934/mbe.2024026 |
Phytoplankton-derived carbon (internal primary production) is known to be essential for the somatic growth and reproduction of zooplankton and fish. For decades, aquatic food webs were taken as systems where carbon transfer was linear from phytoplankton to zooplankton to fish [1].
With the application of isotope labeling, substantial research shows that terrestrial organic matter (TOM) plays an important role in the lake food webs. TOM can not only affect the lake ecosystem in physical and chemical way but also can be exploited by consumer as a resource [2]. And particulate organic matter (POM) of terrestrial origin can be the key factor controlling whole-lake productivity in lakes where phytoplankton productivity is low [1,3]. Besides, in lakes, the inputs of TOM often equal to or exceed internal primary production [4].
Zooplankton composed of cladocerans, cyclopoids, and calanoids represents a critical food chain link between phytoplankton and larger fish in lakes. Lot of experimental results show that zooplankton can directly consume allochthonous POM that either entered the lake in particulate form or was formed through flocculation of allochthonous dissolved organic matter (DOM) [5]. Nevertheless, it is mentioned that Daphnia only supported by allochthonous carbon can survive and give birth to offspring [6]. For the reason that phytoplankton and zooplankton are important components of lakes, plankton mechanism considering TOM is widely studied. Lots of experiments have been done by isotope labeling [5,7,8].
Mathematical modelling is an important tool to investigate the plankton mechanism [9,10,11,12,13,14]. Nutrient-phytoplankton-zooplankton models were studied by many researchers [15,16,17]. In order to enhance the knowledge of plankton mechanism, we propose phytoplankton-zooplankton models considering TOM. We focus on the ecological function of TOM in the plankton mechanism. The results of the deterministic model can fit well with some experimental results and two hypotheses supported by experimental results in [18]. The related two hypotheses are given in the following part:
(Ⅰ) Catchment deposition hypothesis: Allochthony which is the portion of a zooplankter's body carbon content that is of terrestrial origin increases as more TOM is exported from the surrounding catchment.
(Ⅱ) Algal subtraction hypothesis: Allochthony increases with the availability of TOM, where algal production becomes limited by shading more than it benefits from the nutrients associated with TOM.
The rest of the paper is organized as follows: in the next section, by taking TOM into account, formulation of two deterministic models are discussed. In Section 3, global dynamical properties of the three-dimensional ODE model are completely established using Lyapunov function method. In Section 4, by including stochastic perturbation of the white noise type, we develop a stochastic model and show that there is a stationary distribution in the stochastic model. In Section 5, the numerical simulation is given to verify the theoretical predictions. In the biological sense, the numerical results are analyzed in detail. In the last section, the conclusion is given.
The research on the interaction between phytoplankton and zooplankton is important. The phytoplankton is not only the producer of the marine ecosystem but also the base of every food web. Phytoplankton is consumed by zooplankton, which is resource for consumers of higher tropical levels. By transferring organic matter and energy into higher tropical levels, zooplankton plays a critical role in lake ecosystem. Usually, the interaction between the phytoplankton and zooplankton can be described by the Prey-Predator system,
dP(t)dt=P(t)(r1−a1P(t)−bZ(t)),dZ(t)dt=Z(t)(−r2+cP(t)−a2Z(t)). | (2.1) |
Here, P(t) and Z(t) are the population densities of phytoplankton and zooplankton, respectively. The phytoplankton grows with rate r1. The interaction term taken here is Holling type-Ⅰ with b and c as the respective rates of grazing and biomass conversion by the zooplankton for its growth such that b>c. The natural mortality of zooplankton is r2. a1 reflects the density dependence of phytoplankton. The density dependence of predator species is wildly studied because of the environment factor [14]. In lakes where phytoplankton productivity is low, zooplankton will compete each other for the resources to survive, which means that it is more realistic to consider the density dependence of zooplankton population in system (2.1). Furthermore, the zooplankton is resource for higher tropical levels. The mortality caused by higher trophic levels is generally modeled by a nonlinear term [9]. Hence, a2 reflects the density dependence of zooplankton and the mortality of zooplankton caused by higher tropic. Here, it is easy to know that system (2.1) has two equilibria: O(0,0), A(r1a1,0). Besides, when cr1−a1r2>0, there is a positive equilibrium E∗(P∗,Z∗), where P∗=r1a2+r2ba1a2+bc, Z∗=cr1−a1r2a1a2+bc, which is globally asymptotically stable. If cr1−a1r2≤0, equilibrium A(r1a1,0) is globally asymptotically stable.
Model (2.1) ignores TOM in the lake ecosystem. Since TOM has critical influence in the lake food webs, it is more realistic to incorporate it into the phytoplankton-zooplankton model. Based on the important role of TOM in lake ecosystem, we propose the terrestrial organic matter-phytoplankton-zooplankton model. The meaning of the parameters and some assumptions are proposed in following part, firstly.
C(t) is the concentration of TOM. P(t), Z(t) are the population densities of phytoplankton and zooplankton, respectively. s is the constant input rate of TOM and δ is its sinking rate. The interaction term between the TOM and zooplankton is Holling type-Ⅰ. β0 is TOM uptake rate by zooplankton. TOM is converted for the growth of zooplankton with rate β3. The phytoplankton grows with growth rate r. K is the carrying capacity. The interaction term between the phytoplankton and zooplankton is also Holling type-Ⅰ with β1,β2 as the respective rates of grazing and biomass conversion by the zooplankton for their growth. γ is the natural mortality of zooplankton. m measures the strength of competition among zooplankton and the mortality of zooplankton caused by higher tropic. Above statements can be seen clearly in Figure 1.
Based on the above statements, the following ordinary differential equations can be derived,
dC(t)dt=−β0C(t)Z(t)−δC(t)+s,dP(t)dt=rP(t)(1−P(t)K)−β1P(t)Z(t),dZ(t)dt=β2P(t)Z(t)+β3C(t)Z(t)−γZ(t)−mZ2(t). | (2.2) |
In the biological context, the population densities of phytoplankton and zooplankton are nonnegative at any time t. The concentration of TOM is also nonnegative at any time t. We get the following initial conditions,
C(0)=C0≥0,P(0)=P0≥0,Z(0)=Z0≥0. | (2.3) |
Proposition 2.1. The solution (C(t),P(t),Z(t)) of (2.2) with the initial condition (2.3) is nonnegative and bounded in R3+.
Proof. From the equations of (2.2), it is easy to get that C(t)≥C(0)e−∫t0β0Z(s)+δds≥0, P(t)=P(0)e∫t0r(1−P(s)K)−β1Z(s)ds≥0 and Z(t)=Z(0)e∫t0β2P(s)+β3C(s)−γ−mZ(s)ds≥0. Hence, the solutions are non-negative. Clearly, dCdt≤−δC+s yielding C(t)≤max{P0,sδ}:=Cmax. dPdt≤rP(1−PK) yielding P(t)≤max{K,P0}:=Pmax. Therefore, Z(t)≤β2PmaxZ+β3CmaxZ−mZ2. Hence, Z(t)≤max{Z0,β2Pmax+β3Cmaxm}. It implies the solutions (C(t),P(t),Z(t)) of (2.2) are nonnegative and bounded for all t≥0.
The nonnegativity of the solution ensures that the model has practical significance. In the following part, we can get the equilibria of (2.2). By direct calculation, it can be known that (2.2) always has two equilibria,
E1=(sδ,0,0),E2=(C2,P2,0)=(sδ,K,0). |
Proposition 2.2. If s>δγβ3, there exists unique semi-trivial equilibrium E3(C3,0,Z3) of (2.2). If K≥γβ2 and s∈(0,s1), or K<γβ2 and s∈(s0,s1), where
s0=δγβ3−β2Kδβ3,s1=δγβ3+r2mβ0β3β21+r(γβ0+δm)β1β3, |
there exists a unique positive equilibrium E4(C4,P4,Z4).
Proof. The equilibria satisfy the following system of algebraic equations,
−β0CZ−δC+s=0,rP(1−PK)−β1PZ=0,β2PZ+β3CZ−γZ−mZ2=0. | (2.4) |
If P=0, Z≠0, the equilibrium is given by
C∗=γ+mZ∗β3,Z∗=−δC∗+sβ0C∗, | (2.5) |
where C∗ and Z∗ are constants. From the second equation of (2.5), we can get
a3Z2+b3Z+c3=0, | (2.6) |
where a3=−mβ0,b3=−(γβ0+δm),c3=β3s−γδ. If s>δγβ3, then c3>0. Thus, Δ=b23−4a3c3>0, Z31+Z32=−b3a3<0, Z31Z32=c3a3<0, which means that (2.6) exists only one positive root. It is easy to get that Z3=−b3+√b23−4a3c32a3. Then C3=γ+mZ3β3.
If P≠0,Z≠0, the equilibria of (2.4) are given by
H(P)=a4P2+b4P+c4, | (2.7) |
where
a4=−β0r(β2+rmβ1K),b4=(β0rK+δβ1K)(β2+rmβ1K+β0r(γ+rmβ1)),c4=K(β1β3s−(γ+rmβ1)(β0r+δβ1)). |
If s<s1, then c4<0. Hence, H(0)=c4<0. If K≥γβ2, or K<γβ2 and s>s0, H(K)=β3s+β2Kδ−γδ>0. And limP→+∞H(P)<0. Therefore, function H(P) has two positive roots, one of which is larger than K. However, Z=rβ1(1−PK). The root which is bigger than K should be discarded. Hence, there is a unique positive equilibrium.
Here, in biologic context, E1 means that both the phytoplankton and zooplankton go extinct eventually. E2 means that the zooplankton goes extinct, while the phytoplankton will exist at the density of K eventually. E3 means that the phytoplankton goes extinct, while the zooplankton will exist at the density of Z3 eventually. E4 means that the phytoplankton and zooplankton coexist. The existence of E3 and E4 is mainly determined by s, the constant input rate of TOM.
In the next section, we will analyse the stability of those equilibria. It should be pointed out that (2.2) is totally different from the nutrient-phytoplankton-zooplankton model, for the reason that s stands for TOM, which is ingested by zooplankton.
In this section, global dynamical properties of (2.2) are established.
Firstly, at the trivial equilibrium point E1(sδ,0,0) in (2.2), the corresponding characteristic equation is
(λ+δ)(λ−r)(λ+γ−β3sδ)=0. | (3.1) |
It is apparent that equation (3.1) has one positive eigenvalue, which implies that E1(sδ,0,0) is always unstable.
Theorem 3.1. If s≥s1, then E3(C3,0,Z3) is globally asymptotically stable.
Proof. Consider the following Lyapunov function
V1(t)=β32β0C3(C−C3)2+Z−Z3−Z3lnZZ3+β2β1P. | (3.2) |
Since the function h(z)=z−1−lnz is always nonnegative, it is easy to know that V1(t) is also nonnegative at any time t. Then the derivative of V1(t) along the solution of (2.2) is given by
dV1dt=β3β0C3(C−C3)(−β0CZ−δC+s)+β2rβ1P−β2PZ−β2rβ1KP2+(Z−Z3)(β2P+β3C−mZ−γ)=β3β0C3(C−C3)(−β0CZ−δC+β0C3Z3+δC3)−m(Z3−Z)2+β3(C−C3)(Z−Z3)+β2(rβ1−Z3)P−β2rβ1KP2=−β3δβ0C3(C−C3)2−β3ZC3(C−C3)2−m(Z−Z3)2−β2rβ1KP2+β2(rβ1−Z3)P. |
If s≥s1, then β2(rβ1−Z3)≤0. It follows that dV1dt≤0. Furthermore, dV1dt=0 if and only if C(t)=C3, P(t)=0 and Z(t)=Z3. Thus, the largest invariant set on which dV1dt is zero consists of just equilibrium E3. Therefore, by LaSalle's Invariance Principle, equilibrium E3 is globally asymptotically stable.
Theorem 3.2. If K≥γβ2 and s∈(0,s1), or K<γβ2 and s∈(s0,s1), then E4(C4,P4,Z4) is globally asymptotically stable.
Proof. Consider the following Lyapunov function
V2(t)=β32β0C4(C−C4)2+Z−Z4−Z4lnZZ4+β2β1(P−P4−P4lnPP4). | (3.3) |
Similarly, the derivative of V2(t) along the solution of (2.2) is given by
dV2dt=β3β0C4(C−C4)(−β0CZ−δC+s)+β2β1(P−P4)(r−rPK−β1Z)+(Z−Z4)(β2P+β3C−mZ−γ)=β3β0C4(C−C4)(−β0CZ−δC+β0C4Z4+δC4)−m(Z4−Z)2+β3(C−C4)(Z−Z4)−rβ2Kβ1(P−P4)2=−β3δβ0C4(C−C4)2−β3ZC4(C−C4)2−m(Z−Z4)2−rβ2Kβ1(P−P4)2. |
It is easy to know dV2dt≤0. dV2dt=0 if and only if C(t)=C4, P(t)=P4 and Z(t)=Z4, and hence the largest invariant set of (2.2) in the set {(C(t),P(t),Z(t))|dV2dt=0} is the singleton {E4}. Therefore, by the LaSalle's Invariance Principle, equilibrium E4 is globally asymptotically stable.
Theorem 3.3. If K<γβ2 and s∈(0,s0], then E2(sδ,K,0) is globally asymptotically stable.
Proof. Consider the following Lyapunov function
V3(t)=β32β0C2(C−C2)2+Z+β2β1(P−P2−P2lnPP2). | (3.4) |
Similarly, the derivative of V3(t) along the solution of (2.2) is given by
dV3dt=β3β0C2(C−C2)(−β0CZ−δC+s)+β2β1(P−P2)(r−rPK−β1Z)+β2PZ+β3CZ−γZ−mZ2=β3β0C2(C−C2)(−β0CZ−δC+δC2)+β2β1(P−P2)(rP2K−rPK−β1Z)+β2PZ+β3CZ−γZ−mZ2=−β3δβ0C2(C−C2)2+β3β0C2(C−C2)(−β0CZ+β0C2Z−β0C2Z)−rβ2Kβ1(P−P2)2+β2P2Z+β3CZ−γZ−mZ2=−β3δβ0C2(C−C2)2−β3ZC2(C−C2)2−rβ2Kβ1(P−P4)2−mZ2+(β3C2+β2P2−γ)Z. |
If s≤δγβ3−β2Kδβ3, then β3C2+β2P2−γ=β3sδ+β2K−γ≤0. It is easy to know dV3dt≤0. dV3dt=0 if and only if C(t)=sδ, P(t)=K and Z(t)=0. Similarly, by the LaSalle's Invariance Principle, every solution of (2.2) tends to M, where M={(C(t),P(t),Z(t))|dV3dt=0}={E2} is the largest invariant set of (2.2). It implies that E2 is globally asymptotically stable.
Global dynamical properties of (2.2) are shown clearly in Table 1. From biological viewpoint, Theorem 3.1 implies that if the constant input of TOM is larger than s1, it will result in the extinction of phytoplankton eventually which is impossible in the system (2.1). It also means that the increasing input of the TOM will inhibit the growth of phytoplankton, which fits well with the hypotheses(Ⅱ). In addition, in the extinct process of the phytoplankton, in order to survive, the zooplankton will uptake more TOM which is easy to get. That is the Allochthony of zooplankton will increase, which fits well with the hypotheses(Ⅰ). The coexistence of TOM and zooplankton is mentioned in [6]. It also means that the producers of catchment can support both the landscape and lake ecosystem, as told in [2]. Furthermore, from Theorem 3.2, it can be known that different from the system (2.1), with the enough constant input of the TOM(s>s0), zooplankton in the deterministic setting can produce and give offsprings at small carrying capacity(K). The result fits well with the experimental result in [19]. Besides, the threshold value s1 is positively related to the birth rate of phytoplankton, which means that phytoplankton with high birth rate is less possible to become extinct. The result of Theorem 3.3 implies that zooplankton goes extinct if both the carrying capacity and constant input rate of TOM are small.
K≥dβ2 | s∈[s1,+∞) | E3 is G.A.S. |
s∈(0,s1) | E4 is G.A.S. | |
K<dβ2 | s∈[s1,+∞) | E3 is G.A.S. |
s∈(s0,s1) | E4 is G.A.S. | |
s∈(0,s0] | E2 is G.A.S. |
In lake ecosystem, the weather, temperature and other physical factors are hard to predict. TOM, phytoplankton and zooplankton are easily affected by these factors. The deterministic model studied in the previous section does not take environmental randomness into consideration. Taking the randomness into account is more realistic. Recently, many researchers studied the stochastic biological models. Different types of stochastic perturbations are introduced into the models. Prey-predator models with fluctuations around the positive equilibrium were studied [20,21]. The stochastic perturbations which are proportional to the variables were discussed [22,23]. Fluctuations manifesting in the transmission coefficient rate were discussed [24,25].
In this section, by considering role of unpredictable environmental factors in plankton dynamics, we get following stochastic equations,
{dC(t)=(−β0C(t)Z(t)−δC(t)+s)dt+σ1C(t)dB1(t),dP(t)=(rP(t)(1−P(t)K)−β1P(t)Z(t))dt+σ2P(t)dB2(t),dZ(t)=(β2P(t)Z(t)−γZ(t)+β3C(t)Z(t)−mZ2(t))dt+σ3Z(t)dB3(t), | (4.1) |
where Bi(t),i=1,2,3, are independent Brownian motions and σi, i=1,2,3 are the corresponding intensities of stochastic perturbations. Obviously, all the equilibria of (2.2) are no longer the equilibria of the stochastic model (4.1).
In this paper, let (Ω,F,{Ft}t≥0,P) be a complete probability space with filtration {Ft}t≥0 satisfying the usual conditions(i.e. it is right continuous and F0 contains all P-null sets). Let Bi(t), i=1,2,3, be a standard one-dimensional Brownian motion defined on this complete probability space. Let R3+={x∈R3:xi>0,1≤i≤3}.
Definition 4.1 ([26]). A non-negative random variable τ(ω), which is allowed to the value ∞, is called stopping time(with respect to filtration Ft), if for each t, the event {ω:τ≤t}∈Ft.
Lemma 4.2 ([26]). Let {Xt}t be right-continuous on Rn and adapted to Ft. Defineinf∅=∞. If D is an open interval on Rn, τD=inf{t≥0:Xt∉D} is a Ft stopping time.
Lemma 4.3 ([26]). Let x(t)=(x1(t),...,xn(t)) be a regular adapted process n-dimensional vector process. We can have any number n of process driven by a d-dimensional Brownian motion
dx(t)=b(t)dt+σ(t)dBt, |
where σ(t) is n×d matrix valued function, Bt is d-dimensional Brownian motion, x(t),b(t) are n-dim vector-valued functions, the integrals with respect to Brownian motion are Itˆo integrals. Then x(t) is called an Itˆo process. The only restriction is that this dependence results in:
(i) for any i=1,2...n, bi(t) is adapted and ∫T0|bi(t)|<∞ a.s..
(ii) for any i=1,2,...n, σij(t) is adapted and ∫T0σ2ij(t)<∞ a.s.
Lemma 4.4 ([27]). Let x(t) be a l-dimensional Itˆo process on t≥0 with stochastic differential
dx(t)=f(t)dt+g(t)dBt, |
where f∈L(R+;Rl) and g∈L2(R+;Rl×m). Let V∈C2,1(Rl×R+;R). Then V(x(t),t) is till an Itˆo process with the stochastic differential given by
dV(x,t)=[Vt+Vx(t,x)f(t)+12Tr[gT(t)Vxxg(t)]]dt+Vx(x,t)g(t)dB(t)a.s., |
where
Vt=∂V∂t,Vx=(∂V∂x1,...,∂V∂xl),Vxx=(∂2V∂xi∂xj)l×l. |
Suppose V∈C2,1(Rl×R+;R), define a operator LV from Rl×R+ to R by
LV(t,x)=Vt+Vx(t,x)f(t)+12Tr[gT(t)Vxxg(t)]. |
Theorem 4.5. For any given initial value (C(0),P(0),Z(0))∈R3+, (4.1) has almost surely (a.s.) a unique positive solution (C(t),P(t),Z(t)) for t⩾0, and the solution remains in R3+ with probability 1.
Proof. Since the coefficients of the equation are locally Lipschitz continuous, for any given initial value (C(0),P(0),Z(0))∈R3+, it can be obtained that there is a unique local solution (C(t),P(t),Z(t)) on t∈[0,τe) where τe is the explosion time. Explosion refers to the situation when the process reaches infinite values in finite time. Solution can be considered until the time of explosion[26]. If we can check that τe=∞ a.s., then the solution is global. Let h0 be sufficient large, such that every coordinate of (C(0), P(0), Z(0)) lies within the interval [1h0,h0]. For each integer h≥h0, stopping time is defined by
τh=inf{t∈[0,τe):P(t)∉(1h,h)orC(t)∉(1h,h)orZ(t)∉(1h,h)}, |
where we set inf∅=∞. Clearly, τh is increasing as h→∞. Let τ∞=limh→∞τh, where τ∞≤τe a.s. If τ∞=∞ a.s., then τe=∞ and (C(t),P(t),Z(t))∈R3+ a.s. for all t≥0. If not, then it exists constants T>0 and ϵ∈(0,1) such that P{τ∞≤T}>ϵ. Hence, there is integer h1≥h0 such that
P{τh≤T}>ϵ,h≥h1. |
Define a function V4: R3+→R+ as
V4(C,P,Z)=β3β0(C+1−lnC)+β2β1(P+1−lnP)+Z+1−lnZ. | (4.2) |
We can get the non-negativity of this function from f(u)=u+1−lnu≥0,∀u>0. Applying Itˆo formula, it follows that
dV4(C,P,Z)=(β3β0(1−1C)(−β0CZ−δC+s)+β2β1(1−1P)(rP(1−PK)−β1PZ)+(1−1Z)(β2PZ+β3CZ−γZ−mZ2)+β3σ212β0+β2σ222β1+σ232)dt+β3σ1β0(C−1)dB1(t)+β2σ2β1(P−1)dB2(t)+σ3(Z−1)dB3(t)=(β3β0(−δC+s+β0Z+δ−sC)+β2β1(rP(1−PK)−r+rPK)+γ+Z(β2−γ−mZ)−β2P−β3C+mZ+β3σ212β0+β2σ222β1+σ232)dt+β3σ1β0(C−1)dB1(t)+β2σ2β1(P−1)dB2(t)+σ3(Z−1)dB3(t). |
Hence,
dV4(C,P,Z)≤(β3β0(s+δ)+γ+β3C+β2β1(r+rK+β2)P+(β2+β3+m)Z+β3σ212β0+β2σ222β1+σ232)dt+β3σ1β0(C−1)dB1(t)+β2σ2β1(P−1)dB2(t)+σ3(Z−1)dB3(t). |
Let
q1=β3β0(s+δ)+γ+β3σ212β0+β2σ222β1+σ232,q2=max{β0,r+rK+β2,β2+β3+m},q3=β2+β3+m. |
Since the function f is positive, it yields ui≤2(ui+1−lnui)[23]. Therefore, β3C+β2β1(r+rK+β2)P+q3Z≤2q2V4(C,P,Z), which implies
dV4(C,P,Z)≤(q1+2q2V4(C,P,Z))dt+β3σ1β0(C−1)dB1(t)+β2σ2β1(P−1)dB2(t)+σ3(Z−1)dB3(t), |
and hence,
dV4(C,P,Z)≤q4(1+V4(C,P,Z))dt+β3σ1β0(C−1)dB1(t)+β2σ2β1(P−1)dB2(t)+σ3(Z−1)dB3(t), |
where q4=max{q1,2q2}. Let us define a∧b=min{a,b}. If h1≤T, then
∫τh∧h10dV4(C,P,Z)≤∫τh∧h10q4(1+V4(C,P,Z))dt+∫τh∧h10β3σ1β0(C−1)dB1(t)+∫τh∧h10β2σ2β1(P−1)dB2(t)+∫τh∧h10β3σ3β0(Z−1)dB3(t). |
Taking expectation of both sides, we can get that
EV4(C(τh∧h1),P(τh∧h1),Z(τh∧h1))≤V4(C(0),P(0),Z(0))+E∫τh∧h10q4(1+V4(C,P,Z))dt≤V4(C(0),P(0),Z(0))+q4h1+E∫τh∧h10V4(C,P,Z)dt≤V4(C(0),P(0),Z(0))+q4T+∫h10EV4(C(τh∧h1),P(τh∧h1),Z(τh∧h1))dt. |
By the Gronwall inequality,
EV4(C(τh∧T),P(τh∧T),Z(τh∧T))≤(V4(C(0),P(0),Z(0))+q4T)eq4T. |
We set Ωh={τh≤T} for h≥h1, which implies P(Ωh)≥ϵ. Note that for every ω∈Ωh, we can have that C(τh,ω) equals either h or 1h or P(τh,ω) equals either h or 1h or Z(τh,ω) equals either h or 1h, and hence,
V4(C(τh∧T,ω),P(τh∧T,ω),Z(τh∧T,ω))≥min{β3β0(h+1−lnh),β2β1(h+1−lnh),h+1−lnh,β3β0(1h+1+lnh),β2β1(1h+1+lnh),1h+1+lnh}. |
Then it implies that
(V4(C(0),P(0),Z(0))+q4T)eq4T≥E[1Ωh(ω)V4(C(τh∧T,ω),P(τh∧T,ω)Z(τh∧T,ω))]≥ϵmin{β3β0(h+1−lnh),β2β1(h+1−lnh),h+1−lnh,β3β0(1h+1+lnh),β2β1(1h+1+lnh),1h+1+lnh}, |
where 1Ωh is the indicator function of Ωh. Then, we let h→∞, which leads to the contradiction ∞>(V4(C(0),P(0),Z(0))+q4T)eq4T=∞. Therefore, τh=∞ a.s.
Now, we give the main result of this section in the following part.
Theorem 4.6. Assume σ21<δ, σ2, σ3>0 such that
η<min{β3β0C4(δ−σ21)C24,mZ24,rβ2Kβ1P24}. |
Then there is a stationary distribution μ(⋅) for (4.1) with initial value (P0,C0,Z0)∈R3+, which has ergodic property.
Proof. Define a Lyapunov functional
V5(C,P,Z)=β32β0C4(C−C4)2+Z−Z4−Z4lnZZ4+β2β1(P−P4−P4lnPP4). | (4.3) |
Then applying Itˆo's formula to (4.1), we obtain
LV5=−β3δβ0C4(C−C4)2−β3ZC4(C−C4)2−m(Z−Z4)2−rβ2Kβ1(P−P4)2+β3σ212β0C4C2+β2σ222β1+σ232. | (4.4) |
Since (a+b)2≤2(a2+b2), we have that
LV5≤−β3β0C4(δ−σ21)(C−C4)2−β3ZC4(C−C4)2−m(Z−Z4)2−rβ2Kβ1(P−P4)2+β3σ21β0C4+β2σ222β1+σ232≤−β3β0C4(δ−σ21)(C−C4)2−m(Z−Z4)2−rβ2Kβ1(P−P4)2+η, |
where η=β3σ21β0C4+β2σ222β1+σ232. When η<min{β3β0(δ−σ21)C4,mZ24,rβ2Kβ1P24}, we can get that the ellipsoid
−β3β0C4(δ−σ21)(C−C4)2−m(Z−Z4)2−rβ2Kβ1(P−P4)2+η=0, |
lies entirely in R3+. Let U be a neighborhood of the ellipsoid with ˉU⊂El=R3+, so for x∈U∖El, LV5<−ζ(ζ is a positive constant). Therefore, we have that condition (B.2) in Lemma 2.1 of [23] is satisfied(El denotes euclidean l-space). Besides, there is M>0 such that
σ21C2ξ21+σ22P2ξ22+σ23Z2ξ23≥M|ξ2|all(C,P,Z)∈ˉU,ξ∈R3+. |
Thus, we have that the condition (B.1) of [23] is also satisfied. Therefore, there exists a stable stationary distribution μ(⋅) which is ergodic in (4.1).
In this section, the numerical simulations of (2.2) and (4.1) are given to study the proposed models in details. Different simulated results can fit with different hypotheses mentioned in the first section. Unless otherwise stated, parameter values given in Table 2 are used for the simulations. In order to get biologically plausible results, many values are taken from parameter ranges found in the literature. And the numerical simulations of (2.2) are prepared by tool kit ode45 of Matlab.
Parameters | Definition | Values | Reference |
β0 | Zooplankton grazing rate on C | 0.15 | - |
δ | Sinking rate of C | 0.2 | - |
r | Intrinsic growth rate of P | 0.7d−1 | [9] |
β1 | Zooplankton grazing rate on P | 0.4μgcl−1d−1 | - |
β2 | Zooplankton grazing efficiency on P | 0.25 | - |
βmax | Maximum grazing rate on P | 6μgcl−1d−1 | [9] |
β3 | Zooplankton grazing efficiency on C | 0.05 | - |
γ | Zooplankton mortality | 0.3d−1 | [28] |
s | Constant input rate of C | 1 | - |
m | Zooplankton competitive mortality | 0.1 | - |
K | Carrying capacity | 100μgcl−1 | [28] |
It is easy to get s1=4.4, which is the threshold value of model (2.2). It should be noticed that since it is impossible that the absorbed energy is fully used for reproduction, the values of our parameters should satisfy β0>β3 and β1>β2. In the experimental results of [13], it points out that for zooplankton, the growth efficiency of assimilating TOM is lower than that of ingesting phytoplankton. In following analysis, hence, the values of our parameters should satisfy β1>β0.
In the case where s=1<s1, Figure 2(a) shows that solutions dampen and tend to a stable steady state. The carrying capacity varies with the geographical location of the lake. Especially, the carrying capacity of poor nutrient lake is very small. If we decrease K to 1, from Theorem 3.1, we can know dynamical behavior of (2.2) is similar. However, the solution of (2.1) with the same parameter values will tend to A(r1a1,0), which means the extinction of zooplankton. Above statements can fit well with the experimental results of [28] that Daphnia magna(zooplankton) can use terrestrial-derived dissolved organic matter (t-DOM) to support growth and reproduction when alternative food sources are limiting.
The catchment areas around different lakes are different. The constant input of TOM is mainly determined by the catchment areas. It is reasonable to consider the variation of parameter value s. The values of other parameters are kept invariant.
In the case where s=5>s1, we see from Figure 3(a) that the phytoplankton goes extinct eventually. In particular, because of the extinction of phytoplankton, the growth of zooplankton is determined by TOM. That is, the zooplankton shows a preference for TOM, which fits well with the hypothesis (Ⅰ).
In addition, it also means the increasing input rate of TOM will hinder the growth of phytoplankton, which fits well with the hypothesis (Ⅱ).
In this part, we decrease K=100 to 1 and s=1 to 0.1. The other parameter values are the same as before. The result of numerical simulation is shown in Figure 4(a). We can notice that the density of zooplankton decreases to zero, which means if both environmental capacity and constant input of TOM are small, the zooplankton goes extinct in the end.
Finally, we consider the stochastic model. the values of the parameters are same as those for Figure 2. And σ1=σ2=σ3=0.1. The following discretization equations in Milsteins type[29] are used to iteratively calculate the approximate solutions of stochastic system (4.1) in Matlab programs:
Ci+1=Ci+(−δCi−β0CiZi+s)Δt+σ1Ciξ1,i√Δt+σ212Ci(ξ21,iΔt−Δt),Pi+1=Pi+(r(Pi−PiK)−β1PiZi)Δt+σ2Piξ2,i√Δt+σ222Pi(ξ22,iΔt−Δt),Zi+1=Zi+(β2PiZi+β3CiZi−γZi−mZ2i)Δt+σ3Ziξ3,i√Δt+σ232Zi(ξ23,iΔt−Δt), | (5.1) |
where ξ1,i,ξ2,i, and ξ3,i are N(0,1)-distributed independent Gaussian random variable, σ1,σ2 and σ3 are intensities of white noise and time increment Δt>0. The simulation results are shown in Figure 5.
In this paper, firstly, phytoplankton-zooplankton model including TOM is proposed to investigate effects of TOM upon planktonic dynamics. By constructing Lyapunov functions and using LaSalle's Invariance Principle, global stability of the equilibria which is determined by threshold values s0 and s1 is established. The analytical results of (2.2) indicate that TOM has significant effects on the lake ecosystem stability. Different constant input rates of TOM in the model can result in different results, which represents different biological meanings.
Plankton populations which reside in lake ecosystem are persistently influenced by the environmental fluctuations. In order to make more precise biological findings, we develop a stochastic model. We get that when the perturbation is small, the plankton populations reach a stationary distribution, which is ergodic. Then we simulate model (2.2) with many parameter values taken from the literature by varying the constant input rate of TOM along with the carrying capacity value. The simulations indicate that the magnitude of constant input rate of TOM and carrying capacity affect the persistence of phytoplankton populations greatly. We also simulate the stochastic model with the same parameter values. We can notice that mean value trajectory of (4.1) which is the average of ten trajectories oscillates randomly surrounding the solution of system (2.2), rather than reach the stable steady state (Figure 5). And the simulations also show that the trajectory of (4.1) is totally different from that of (2.2).
In summary, we set up exploratory models to consider the potential role of TOM that mediates interactions between trophic levels in a simple plankton food-chain. We have shown that, in principle, TOM plays an important role in influencing interactions between phytoplankton and zooplankton. The theoretical results of our deterministic model fit well with some experimental results. Further, when the perturbation is small, the stochastic model exists a stationary distribution which is ergodic. We emphasize that the qualitative behaviour of our models that we are interested in and this study are only meant as an initial exploratory attempt to consider the influence of TOM on the lake ecosystem.
This research is supported by the National Natural Sciences Foundation of China (Grant No. 11571326).
The authors declare there is no conflict of interest.
[1] | M. Rautio, H. Mariash and L. Forsstrom, Seasonal shifts between autochthonous and allochthon-ous carbon contributions to zooplankton diets in a subarctic lake, Limnol. Oceanogr., 56 (2011), 1513–1524. |
[2] | S. E. Jones, C. T. Solomon and B. C. Weidel, Subsidy or subtraction: how do terrestrial inputs influence consumer production in lakes, Freshwater Rev., 5 (2012), 37–49. |
[3] | K. A. Emery, G. M. Wilkinson, M. L. Pace, et al., Use of allochthonous resources by zooplankton in reservoirs, Hydrobiologia, 758 (2015), 257–269. |
[4] | J. Cole, S. R. Carpenter, J. F. Kitchell, et al., Differential support of lake food webs by three types of terrestrial organic carbon, Ecol. Lett., 9 (2006), 558–568. |
[5] | M. Berggren, S. E. Ziegler, N. F. Gelais, et al., Contrasting patterns of allochthony among three major groups of crustacean zooplankton in boreal and temperate lakes, Ecology, 95 (2014), 1947–1959. |
[6] | M. T. Brett, M. Kainz, S. J. Taipale, et al., Phytoplankton, not allochthonous carbon, sustains herbivorous zooplankton production, Proc. Natl. Acad. Sci. USA, 21 (2009), 197–201. |
[7] | M. Berggren, A. Bergstrom and J. Karlsson, Intraspecific autochthonous and allochthonous resource use by zooplankton in a humic lake during the transitions between winter, summer and fall, Plos One, 10 (2014), e0120575. |
[8] | A. D. Persaud and P. J. Dillon, Differences in zooplankton feeding rates and isotopic signatures from three temperate lakes, Aquat. SCI., 73 (2011), 261–273. |
[9] | N. D. Lewis, M. N. Breckels, M. Steinke, et al., Role of infochemical mediated zooplankton grazing in a phytoplankton competition model, Ecol. Complex., 16 (2013), 41–50. |
[10] | Y. Lv and Y. Pei, Harvesting of a phytoplankton-zooplankton model, Nonlinear Anal. Real. World. Appl., 11 (2010), 3608–3619. |
[11] | A. Sharma, A. K. Sharma and K. Agnihotri, Analysis of a toxin producing phytoplankton zooplankton interaction with Holling IV type scheme and time delay, Nonlinear Dynam., 81 (2015), 13–25. |
[12] | J. E. Truscott and J. Brindley, Ocean plankton population as excitable media, B. Math. Biol., 56 (1994), 981–998. |
[13] | Y. Pei, Y. Lv and C. Li, Evolutionary consequences of harvesting for a two-zooplankton one-phytoplankton system, Appl. Math. Model, 36 (2012), 1752–1765. |
[14] | L. Qian, Q. Lu, Q. Meng et al., Dynamical behaviors of a prey-predator system with impulsive control, J. Math. Anal. Appl., 1 (2010), 345–356. |
[15] | S. R. Jang and E. J. Allen, Deterministic and stochastic nutrient phytoplankton zooplankton models with periodic toxin producing phytoplankton, Appl. Math. Comput., 271 (2015), 52–67. |
[16] | S. Ruan, Persistence and coexistence in zooplankton-phytoplankton nutrient models with instantaneous nutrient recycling, J. Math. Biol., 31 (1993), 633–654. |
[17] | J. S. Wroblewski, S. L. Sarmiento and G. R. Flierl, An ocean basin scale model of plankton dynamics in the North Atlantic, solutions for the climatological oceanographic condition in May, Global Biogeochem. CY., 2 (1988), 199–218. |
[18] | A. J. Tanentzap, B. W. Kielstra, G. M. Wilkinson, et al., Terrestrial support of lake food webs: synthesis reveals controls over cross-ecosystem resource use, Sci. Adv., 3 (2017), e1601765. |
[19] | B. C. Mcmeans, A. Koussoroplis, M. T. Arts, et al., Terrestrial dissolved organic matter supports growth and reproduction of Daphnia magna when algae are limiting, J. Plankton. Res., 0 (2015), 1–9. |
[20] | U. Forys, M. Qiao and A. Liu, Asymptotic dynamics of a deterministic and stochastic predator- prey model with disease in the prey species, Math. Method. Appl. Sci., 37 (2014), 306–320. |
[21] | G. Ge, H. Wang and J. Xu, A stochastic analysis for a phytoplankton zooplankton model, JPCS, 96 (2008), 012168. |
[22] | P. Ghosh, P. Das and D. Mukherjee, Persistence and stability of a seasonally perturbed three species stochastic model of salmonoid aquaculture, Differ. Equ. Dyn. Syst., 2016. Available from: https://doi.org/10.1007/s12591-016-0283-0. |
[23] | C. Ji, D. Jiang and N. Shi, A note on a predator-prey model with modified Leslie-Gower and Holling-type II schemes with stochastic perturbation, J. Math. Anal. Appl., 377 (2011), 435–440. |
[24] | J. Djordjevic, C. J. Silva and D. F. Torres, A stochastic SICA epidemic model for HIV transmission, Appl. Math. Lett., 84 (2018), 168–175. |
[25] | R. Rajaji and M. Pitchaimani, Analysis of stochastic viral infection model with immune impairment, Int. J. Ap. Mat. Com-Pol., 3 (2017), 3561–3574. |
[26] | F. C. Klebaner, Introduction to Stochastic Calculus with Applications, 2nd edition, World Scientific Publishing, London, 2005. |
[27] | X. Mao, Stochastic Diferential Equations and Applications, 1nd edition, Horwood Publishing, Chichester, UK, 2008. |
[28] | N. D. Lewis, M. N. Breckels, S. D. Archer, et al., Grazing-induced production of DMS can stabilize food-web dynamics and promote the formation of phytoplankton blooms in a multitrophic plankton model, Biogeochem., 110 (2012), 303–313. |
[29] | D. J. Higham, An algorithmic introduction to numerical simulation of stochastic differential equations, SIAM Rev., 43 (2001), 525–546. |
[30] | N. Dalal, D. Greenhalgh and X. Mao, A stochastic model for internal HIV dynamics, J. Math. Anal. Appl., 341 (2008), 1084–1101. |
[31] | L. Shaikheta and A. Korobeinikov, Stability of a stochastic model for HIV-1 dynamics within a host, Appl. Anal., 6 (2016), 1228–1238. |
1. | Tiancai Liao, Hengguo Yu, Chuanjun Dai, Min Zhao, Impact of Cell Size Effect on Nutrient-Phytoplankton Dynamics, 2019, 2019, 1076-2787, 1, 10.1155/2019/8205696 |
K≥dβ2 | s∈[s1,+∞) | E3 is G.A.S. |
s∈(0,s1) | E4 is G.A.S. | |
K<dβ2 | s∈[s1,+∞) | E3 is G.A.S. |
s∈(s0,s1) | E4 is G.A.S. | |
s∈(0,s0] | E2 is G.A.S. |
Parameters | Definition | Values | Reference |
β0 | Zooplankton grazing rate on C | 0.15 | - |
δ | Sinking rate of C | 0.2 | - |
r | Intrinsic growth rate of P | 0.7d−1 | [9] |
β1 | Zooplankton grazing rate on P | 0.4μgcl−1d−1 | - |
β2 | Zooplankton grazing efficiency on P | 0.25 | - |
βmax | Maximum grazing rate on P | 6μgcl−1d−1 | [9] |
β3 | Zooplankton grazing efficiency on C | 0.05 | - |
γ | Zooplankton mortality | 0.3d−1 | [28] |
s | Constant input rate of C | 1 | - |
m | Zooplankton competitive mortality | 0.1 | - |
K | Carrying capacity | 100μgcl−1 | [28] |
K≥dβ2 | s∈[s1,+∞) | E3 is G.A.S. |
s∈(0,s1) | E4 is G.A.S. | |
K<dβ2 | s∈[s1,+∞) | E3 is G.A.S. |
s∈(s0,s1) | E4 is G.A.S. | |
s∈(0,s0] | E2 is G.A.S. |
Parameters | Definition | Values | Reference |
β0 | Zooplankton grazing rate on C | 0.15 | - |
δ | Sinking rate of C | 0.2 | - |
r | Intrinsic growth rate of P | 0.7d−1 | [9] |
β1 | Zooplankton grazing rate on P | 0.4μgcl−1d−1 | - |
β2 | Zooplankton grazing efficiency on P | 0.25 | - |
βmax | Maximum grazing rate on P | 6μgcl−1d−1 | [9] |
β3 | Zooplankton grazing efficiency on C | 0.05 | - |
γ | Zooplankton mortality | 0.3d−1 | [28] |
s | Constant input rate of C | 1 | - |
m | Zooplankton competitive mortality | 0.1 | - |
K | Carrying capacity | 100μgcl−1 | [28] |