Loading [MathJax]/jax/output/SVG/jax.js
Research article Special Issues

Mathematical analysis on deterministic and stochastic lake ecosystem models

  • In this paper, we propose and study the deterministic and stochastic lake ecosystem models to investigate the impact of terrestrial organic matter upon persistence of the plankton populations. By constructing Lyapunov function and using the LaSalle's Invariance Principle, we establish global properties of the deterministic model. The dynamical behavior of solutions fits well with some experimental results. It is concluded that the terrestrial organic matter plays an important role in influencing interactions between phytoplankton and zooplankton. Based on the fluctuations of lake ecosystem, we further develop a stochastically perturbed model. Theoretic analysis implies that the stochastic model exists a stationary distribution which is ergodic. The key point of our analysis is to enhance our knowledge of the factors governing the dynamics of plankton population models.

    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

    Related Papers:

    [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
  • In this paper, we propose and study the deterministic and stochastic lake ecosystem models to investigate the impact of terrestrial organic matter upon persistence of the plankton populations. By constructing Lyapunov function and using the LaSalle's Invariance Principle, we establish global properties of the deterministic model. The dynamical behavior of solutions fits well with some experimental results. It is concluded that the terrestrial organic matter plays an important role in influencing interactions between phytoplankton and zooplankton. Based on the fluctuations of lake ecosystem, we further develop a stochastically perturbed model. Theoretic analysis implies that the stochastic model exists a stationary distribution which is ergodic. The key point of our analysis is to enhance our knowledge of the factors governing the dynamics of plankton population models.


    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)(r1a1P(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 cr1a1r2>0, there is a positive equilibrium E(P,Z), where P=r1a2+r2ba1a2+bc, Z=cr1a1r2a1a2+bc, which is globally asymptotically stable. If cr1a1r20, 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.

    Figure 1.  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 phytoplankton grows with growth rate r. γ is the natural mortality of zooplankton. K is the carrying capacity. m measures the strength of competition among zooplankton and the mortality of zooplankton caused by higher tropic. The biomass conversion between the phytoplankton and zooplankton is β2PZ. The biomass conversion between TOM and zooplankton is β3PZ.

    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)(1P(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)=C00,P(0)=P00,Z(0)=Z00. (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)et0β0Z(s)+δds0, P(t)=P(0)et0r(1P(s)K)β1Z(s)ds0 and Z(t)=Z(0)et0β2P(s)+β3C(s)γmZ(s)ds0. Hence, the solutions are non-negative. Clearly, dCdtδC+s yielding C(t)max{P0,sδ}:=Cmax. dPdtrP(1PK) yielding P(t)max{K,P0}:=Pmax. Therefore, Z(t)β2PmaxZ+β3CmaxZmZ2. 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 t0.

    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(1PK)β1PZ=0,β2PZ+β3CZγZmZ2=0. (2.4)

    If P=0, Z0, 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, Δ=b234a3c3>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+b234a3c32a3. Then C3=γ+mZ3β3.

    If P0,Z0, 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(1PK). 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 ss1, then E3(C3,0,Z3) is globally asymptotically stable.

    Proof. Consider the following Lyapunov function

    V1(t)=β32β0C3(CC3)2+ZZ3Z3lnZZ3+β2β1P. (3.2)

    Since the function h(z)=z1lnz 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(CC3)(β0CZδC+s)+β2rβ1Pβ2PZβ2rβ1KP2+(ZZ3)(β2P+β3CmZγ)=β3β0C3(CC3)(β0CZδC+β0C3Z3+δC3)m(Z3Z)2+β3(CC3)(ZZ3)+β2(rβ1Z3)Pβ2rβ1KP2=β3δβ0C3(CC3)2β3ZC3(CC3)2m(ZZ3)2β2rβ1KP2+β2(rβ1Z3)P.

    If ss1, then β2(rβ1Z3)0. It follows that dV1dt0. 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(CC4)2+ZZ4Z4lnZZ4+β2β1(PP4P4lnPP4). (3.3)

    Similarly, the derivative of V2(t) along the solution of (2.2) is given by

    dV2dt=β3β0C4(CC4)(β0CZδC+s)+β2β1(PP4)(rrPKβ1Z)+(ZZ4)(β2P+β3CmZγ)=β3β0C4(CC4)(β0CZδC+β0C4Z4+δC4)m(Z4Z)2+β3(CC4)(ZZ4)rβ2Kβ1(PP4)2=β3δβ0C4(CC4)2β3ZC4(CC4)2m(ZZ4)2rβ2Kβ1(PP4)2.

    It is easy to know dV2dt0. 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(CC2)2+Z+β2β1(PP2P2lnPP2). (3.4)

    Similarly, the derivative of V3(t) along the solution of (2.2) is given by

    dV3dt=β3β0C2(CC2)(β0CZδC+s)+β2β1(PP2)(rrPKβ1Z)+β2PZ+β3CZγZmZ2=β3β0C2(CC2)(β0CZδC+δC2)+β2β1(PP2)(rP2KrPKβ1Z)+β2PZ+β3CZγZmZ2=β3δβ0C2(CC2)2+β3β0C2(CC2)(β0CZ+β0C2Zβ0C2Z)rβ2Kβ1(PP2)2+β2P2Z+β3CZγZmZ2=β3δβ0C2(CC2)2β3ZC2(CC2)2rβ2Kβ1(PP4)2mZ2+(β3C2+β2P2γ)Z.

    If sδγβ3β2Kδβ3, then β3C2+β2P2γ=β3sδ+β2Kγ0. It is easy to know dV3dt0. 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.

    Table 1.  Stability of equilibria in (2.2).
    Kdβ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.

     | Show Table
    DownLoad: CSV

    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)(1P(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}t0,P) be a complete probability space with filtration {Ft}t0 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+={xR3:xi>0,1i3}.

    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{t0:XtD} 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 t0 with stochastic differential

    dx(t)=f(t)dt+g(t)dBt,

    where fL(R+;Rl) and gL2(R+;Rl×m). Let VC2,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=Vt,Vx=(Vx1,...,Vxl),Vxx=(2Vxixj)l×l.

    Suppose VC2,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 t0, 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 hh0, 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 t0. If not, then it exists constants T>0 and ϵ(0,1) such that P{τT}>ϵ. Hence, there is integer h1h0 such that

    P{τhT}>ϵ,hh1.

    Define a function V4: R3+R+ as

    V4(C,P,Z)=β3β0(C+1lnC)+β2β1(P+1lnP)+Z+1lnZ. (4.2)

    We can get the non-negativity of this function from f(u)=u+1lnu0,u>0. Applying Itˆo formula, it follows that

    dV4(C,P,Z)=(β3β0(11C)(β0CZδC+s)+β2β1(11P)(rP(1PK)β1PZ)+(11Z)(β2PZ+β3CZγZmZ2)+β3σ212β0+β2σ222β1+σ232)dt+β3σ1β0(C1)dB1(t)+β2σ2β1(P1)dB2(t)+σ3(Z1)dB3(t)=(β3β0(δC+s+β0Z+δsC)+β2β1(rP(1PK)r+rPK)+γ+Z(β2γmZ)β2Pβ3C+mZ+β3σ212β0+β2σ222β1+σ232)dt+β3σ1β0(C1)dB1(t)+β2σ2β1(P1)dB2(t)+σ3(Z1)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(C1)dB1(t)+β2σ2β1(P1)dB2(t)+σ3(Z1)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 ui2(ui+1lnui)[23]. Therefore, β3C+β2β1(r+rK+β2)P+q3Z2q2V4(C,P,Z), which implies

    dV4(C,P,Z)(q1+2q2V4(C,P,Z))dt+β3σ1β0(C1)dB1(t)+β2σ2β1(P1)dB2(t)+σ3(Z1)dB3(t),

    and hence,

    dV4(C,P,Z)q4(1+V4(C,P,Z))dt+β3σ1β0(C1)dB1(t)+β2σ2β1(P1)dB2(t)+σ3(Z1)dB3(t),

    where q4=max{q1,2q2}. Let us define ab=min{a,b}. If h1T, then

    τhh10dV4(C,P,Z)τhh10q4(1+V4(C,P,Z))dt+τhh10β3σ1β0(C1)dB1(t)+τhh10β2σ2β1(P1)dB2(t)+τhh10β3σ3β0(Z1)dB3(t).

    Taking expectation of both sides, we can get that

    EV4(C(τhh1),P(τhh1),Z(τhh1))V4(C(0),P(0),Z(0))+Eτhh10q4(1+V4(C,P,Z))dtV4(C(0),P(0),Z(0))+q4h1+Eτhh10V4(C,P,Z)dtV4(C(0),P(0),Z(0))+q4T+h10EV4(C(τhh1),P(τhh1),Z(τhh1))dt.

    By the Gronwall inequality,

    EV4(C(τhT),P(τhT),Z(τhT))(V4(C(0),P(0),Z(0))+q4T)eq4T.

    We set Ωh={τhT} for hh1, 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(τhT,ω),P(τhT,ω),Z(τhT,ω))min{β3β0(h+1lnh),β2β1(h+1lnh),h+1lnh,β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)eq4TE[1Ωh(ω)V4(C(τhT,ω),P(τhT,ω)Z(τhT,ω))]ϵmin{β3β0(h+1lnh),β2β1(h+1lnh),h+1lnh,β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(CC4)2+ZZ4Z4lnZZ4+β2β1(PP4P4lnPP4). (4.3)

    Then applying Itˆo's formula to (4.1), we obtain

    LV5=β3δβ0C4(CC4)2β3ZC4(CC4)2m(ZZ4)2rβ2Kβ1(PP4)2+β3σ212β0C4C2+β2σ222β1+σ232. (4.4)

    Since (a+b)22(a2+b2), we have that

    LV5β3β0C4(δσ21)(CC4)2β3ZC4(CC4)2m(ZZ4)2rβ2Kβ1(PP4)2+β3σ21β0C4+β2σ222β1+σ232β3β0C4(δσ21)(CC4)2m(ZZ4)2rβ2Kβ1(PP4)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)(CC4)2m(ZZ4)2rβ2Kβ1(PP4)2+η=0,

    lies entirely in R3+. Let U be a neighborhood of the ellipsoid with ˉUEl=R3+, so for xUEl, 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ξ23M|ξ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.

    Table 2.  Parameter definitions and the values used for model (2.2) and model (4.1).
    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.7d1 [9]
    β1 Zooplankton grazing rate on P 0.4μgcl1d1 -
    β2 Zooplankton grazing efficiency on P 0.25 -
    βmax Maximum grazing rate on P 6μgcl1d1 [9]
    β3 Zooplankton grazing efficiency on C 0.05 -
    γ Zooplankton mortality 0.3d1 [28]
    s Constant input rate of C 1 -
    m Zooplankton competitive mortality 0.1 -
    K Carrying capacity 100μgcl1 [28]

     | Show Table
    DownLoad: CSV

    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.

    Figure 2.  In the case where s=1<s1, the solutions of (2.2) that initially oscillate will dampen and tend to a stable steady state, which means the phytoplankton and zooplankton coexist in the end.

    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 (Ⅰ).

    Figure 3.  In the case where s=5>s1, the high input rate of TOM can result in the extinction of the phytoplankton.

    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.

    Figure 4.  Here, K=1,s=0.1, it can be seen that 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(PiPiK)β1PiZi)Δt+σ2Piξ2,iΔt+σ222Pi(ξ22,iΔtΔt),Zi+1=Zi+(β2PiZi+β3CiZiγZimZ2i)Δ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.

    Figure 5.  Here, σ1=σ2=σ3=0.1. The other parameter values are the same as those for Figure 2. The Figure 5a and Figure 5b are the trajectories of population densities of phytoplankton and zooplankton respectively. When the perturbation is small, the numerical simulations illustrate that mean value trajectory(the red trajectory) of (4.1) which is the average of ten trajectories(the blue trajectories) of (4.1) oscillates randomly surrounding the solution of (2.2) (the green trajectory). And the numerical simulations also show that the trajectory of (4.1) is totally different from that of (2.2).

    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.
  • This article has been cited by:

    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
  • Reader Comments
  • © 2019 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(4912) PDF downloads(519) Cited by(1)

Figures and Tables

Figures(5)  /  Tables(2)

Other Articles By Authors

/

DownLoad:  Full-Size Img  PowerPoint
Return
Return

Catalog