Abstract
A toxin producing phytoplankton-zooplankton model with inhibitory exponential substrate and time delay has been formulated and analyzed. Since the liberation of toxic substances by phytoplankton species is not an instantaneous process but is mediated by some time lag required for maturity of the species and the zooplankton mortality due to the toxic phytoplankton bloom occurs after some time laps of the bloom of toxic phytoplankton, we induced a discrete time delay to both of the consume response function and distribution of toxic substance term. Furthermore, based on the fact that the predation rate decreases at large toxic-phytoplankton density, the system is modelled via a Tissiet type functional response. We study the dynamical behaviour and investigate the conditions to guarantee the coexistence of two species. Analytical methods and numerical simulations are used to obtain information about the qualitative behaviour of the models.
1. Introduction
Phytoplankton are one of the most important components of the marine ecosystem. They not only form a basis for all aquatic food chains but also perform a very useful service by producing a huge amount of oxygen for human and other living animals after absorbing carbon-dioxide from surrounding environments [1]. Zooplankton are microscopic animals that eat other plankton and serve as a most favorable food source for fish and other aquatic animals. During the recent years, many authors have studied the system between zooplankton and phytoplankton. Authors in [2] have dealt with a nutrient-plankton model in an aquatic environment in the context of phytoplankton bloom. In [3] the effect of seasonality and periodicity on plankton dynamics is investigated. In [4], two plankton ecosystem models with explicit representation of viruses and virally infected phytoplankton are presented.
The most common features of the phytoplankton population is rapid increase of biomass due to rapid cell proliferation and almost equally rapid decrease in populations, separated by some fixed time period. This type of rapid change in phytoplankton population density is known as “bloom” [5]. Due to the accumulation of high biomass or to the presence of toxicity, some of these blooms, more adequately called “harmful algal blooms” [6], are noxious to marine ecosystems or to human health and can produce great socioeconomic damage. There has been a global increase in harmful plankton blooms in last two decades [7–9].
Because of the difficulty in measuring plankton biomass, mathematical modeling of plankton population is an important alternative method of improving our knowledge of the physical and biological processes relating to plankton ecology [10–15]. In [16], nutrient-plankton-zooplankton interaction models with a toxic substance which inhibits either the growth of phytoplankton, zooplankton, or both trophic levels are proposed and studied. In [17], the authors have constructed a mathematical model for describing the interaction between a nontoxic and toxic phytoplankton with a single nutrient.
Based on the fact that the release of toxin from phytoplankton species is not an instantaneous process but is mediated by some time lag required for maturity of the species and the zooplankton may die after some time lapse of the bloom of toxic phytoplankton (see http://www.mote.org/, http://www.mdsg.umd.edu/), models incorporating time delay in diverse biological models are extensively reviewed by Beretta and Kuang [18], Gopalsamy [19], Cooke and Grossman [20], and Cushing [21]. The discrete time delay has potential to change the qualitative behavior of the dynamical systems [22–29]. Chattopadhyay et al. proposed a delay model incorporating time lag in toxin liberation by phytoplankton to avoid predation by zooplankton [28]. In [28], the authors introduced distribution delay and discrete to toxin liberation term. Due to discrete time delay in toxin liberation, the local existence of periodic solution through Hopf bifurcation has been obtained in [28].
In [5, 10, 28–31], the following plausible toxic-phytoplankton-zooplankton system has been studied: where is the density of the toxin-producing phytoplankton (TPP) population and is the density of zooplankton population at any instant of time . In model (1), is the intrinsic growth rate and is the environmental carrying capacity of TPP population. The term describes the functional response for the grazing of phytoplankton by zooplankton, and is the maximum uptake rate for zooplankton. denotes the ratio of biomass conversion and is the natural death rate of zooplankton. The function represents the distribution of toxic substance which ultimately contributes to the death of zooplankton populations, and the parameter denotes the rate of toxic substances produced by per unit biomass of phytoplankton.
Model (1) has been studied for the following cases: when is linear but is Holling type II [10, 28, 29] or Holling type III [10]; when and are both Holling type II [5, 10] or Holling type III [10]; when is Holling type II while is Holling type III [10].
In above cases, and are both increasing functions of over the entire interval . However, in some cases, very high substrate concentrations in the lakes actually inhabit the growth of phytoplankton cells. Moreover, with the substrate concentrations increasing unlimitedly, some kind of microorganism will die eventually [32]. To describe the above phenomenon accurately, we consider from a different point of view. Adopting the idea used in [32], we assume that there exists a constant such that is increasing over the interval and is decreasing on the interval . More precisely, we use the so-called Tissiet functional response of the form of (see, e.g., [32]). This type of functional response takes care of the fact that the predation rate decreases at large toxic-phytoplankton density.
The remainder of this paper is organized as follows. The model is described in Section 2. In Section 3, we state and prove the positivity and boundedness of the solutions. Then, in Section 4, equilibria, their existence, and local asymptotic stability are considered. Under the aids of numerical simulation method, we further analyse the model and determine if there is a parameter range for the delay parameter , where more complicated dynamics occur. In Section 5, the permanence of the system is discussed by some analytic techniques on limit sets of differential dynamical systems. Finally, a brief discussion is presented in Section 6.
2. State of the Model
In a real ecological context, the interaction between phytoplankton and zooplankton will not be essentially instantaneous. Instead, the response of zooplankton to contacts with phytoplankton is likely to be delayed due to a gestation period. Another fact, during the interaction between phytoplankton and zooplankton, is that the liberation of toxic substances by phytoplankton must be mediated by some time lag which is required for the maturity of toxic-phytoplankton. Let stand for the time delay in conversion of food to viable biomass for the species, and is the discrete time period required for the maturity of phytoplankton cells to liberate toxic substances. Based on model 3 in [5], we intend to study a model system with the assumption that and are described by same type of function, namely, Tissiet functional response. Then we arrive at the following model: where the parameter is the intrinsic growth rate and is the environmental carrying capacity of TPP population. The constant is the maximum per capita grazing rate, denotes the ratio of biomass conversion, denotes the rate of toxic substances produced by per unit biomass, is the natural death rate of the zooplankton. All the parameters in system (2) are positive constants with their usual ecological meanings.
Here we observe that if there is no delay (i.e., ) and , then . Hence, throughout our analysis, we assume that .
For the sake of simplicity of mathematical analysis, in this paper, we consider model (2) for the special case: the time delay in conversion of food to viable biomass for the species equal to the discrete time period required for the maturity of phytoplankton cells to liberate toxic substances; that is, .
By performing the following scaling for model (2): we obtain a dimensionless system in the state variables , , which can be written, by removing the stars, in the form with initial conditions: where .
3. Positive and Boundedness
In this section, we consider the positive and boundedness of the solutions of model (4) with initial condition (5). One has the following theorem.
Theorem 1. Let , on and , . Then(a)all the solutions of (4) with initial condition (5) exist on for some constant and are unique and positive for ,(b), , where ,(c)if , then , where . Further, the subset is positively invariant with respect to (4).
Proof. By Theorems 2.1 and 2.3 in Hale and Lunel [33], solutions of (4) with initial data exist on for some and are unique. Suppose is a solution of (4) for . Without loss of generality, assume that is the maximum internal of the solution and if the solution exists for any . Integrating the first equation of (4) gives
To prove the for any , use the method of contradiction. Suppose there exists a such that
From the second equation of system (4), we have
This contradicts . Hence, for all . This completes the proof of conclusion (a) in Theorem 1.
It follows from the first equation of (4) that , which implies that . Define , . Then from (4) we obtain
Applying the theorem of differential inequality we obtain that
Therefore, . Thus, there is a constant , such that for all . This completes the proof of conclusion (b) in Theorem 1.
From the first equation of system (4) we get
which implies that .
For any , let be the solution of (4) with the initial function . If there is a such that , then for some and . Hence, it follows from the first equation of (4) that
which is a contradiction to . So, for all .
It is easy to prove that if there is a such that , then for all one has . Which implies that for all . This completes the proof of conclusion (c) in Theorem 1.
Remark. The conclusions (b) and (c) in Theorem 1 indicate that if the ratio of biomass conversion is less than certain value, then phytoplankton population will be persistent.
4. Equilibrium, Stability, and Hopf Bifurcation
4.1. Existence of Equilibria
It is easy to see that model (4) has two boundary equilibria and . To discuss the existence of the positive equilibria, we work on the equation Denote the left-hand side of the second equation in (14): Then , , Let , we have , where From this we know that the function is monotone increasing in the interval and monotone decreasing in the interval . Hence, reaches its maximum on the interval . Therefore, we have the following results.
Theorem 2. The following statements hold.(1)If , then equation has no roots in the interval (Figure 1(a)). In this case, for system (4), there is no positive equilibrium.(2)If and , then equation has only one root, , in the interval , where (Figure 1(c)). In this case, for system (4), there exists a single positive equilibrium .(3)If and , then equation has two distinct roots, and , in the interval , where (Figure 1(d)). In this case, for system (4), there are two distinct positive equilibria and .(4)If , then equation has a unique root, , in the interval , (Figure 1(b)). In this case, for system (4), there is a unique positive equilibrium .(5)If and , then equation has two distinct roots, and , in the interval , where (Figure 1(e)). In this case, for system (4), there is also unique positive equilibrium ,where , , are determined by the first equation in (14).
(a)
(b)
(c)
(d)
(e)
4.2. Local Stability and Hopf Bifurcation
Let denote any one of the equilibrium points of system (4). Linearizing (4) about we obtain In what follows, existence of the interior equilibria, dynamical properties of the zooplankton-free equilibrium, and the interior equilibrium are investigated.
Theorem 3. For any time delay , is locally asymptotically stable if and is unstable if ; the trivial solution of the linearized system of (4) about is stable for .
Proof. From (18) we know that, at any equilibrium point , the characteristic equation is
where
in which .
We first consider local stability of . For , (19) becomes
One of the roots of (21) is , and other roots are given by solution of the following equality:
It follows from [34] that (i)if , then all roots of (22) have negative parts for any time delay . Hence, is locally asymptotically stable for any time delay ;(ii)if , then (22) has roots which have positive real parts for any time delay . Hence, is unstable for any time delay ;(iii)if , this is a critical case and (22) is equivalent to
From [34], it has that except , any root of (23) has negative real part for any time delay . Hence, the trivial solution of the linearized system of (4) about is stable for any time delay .
Theorem 4. Suppose that the positive equilibria , , of system (4) exist. For system (4), the following statements hold. (i)If and , then are locally asymptotically stable for , are unstable for , and there is a periodic solution around for . If , then is unstable for any , where is determined by (30), .(ii) is unstable for any .(iii)If , then the trivial solution of the linearized system of (4) about is stable for any . If , then the trivial solution of the linearized system of (4) about is unstable for any ,where , .
Proof. When , the transcendental equation (19) becomes
where
if , and if , .
By the argument in Section 4.1, as long as exists, it must be , ; thus, ; as long as exists, it must be ; thus, ; as long as exists, it must be ; thus, ; hence, from Routh-Hurwitz theorem, we have that if , , then and are locally asymptotically stable. If , , then and are unstable; if , , then and are nonhyperbolic equilibria; is unstable; is a critical case.
Suppose that , , is a root of (19) for some . Substituting in the characteristic equation (19), and separating real and imaginary parts, then we get
Eliminating from (26), we get a biquadratic equation in as
Its roots are
By (20) we obtain
We consider the stability of , .
Since , if , , then there is only one root, , , such that
That is, (19) has one imaginary solution, , . From (26), we obtain the following set of for which there are imaginary root:
where and
If , then , when , , are locally asymptotically stable. Hence, if and , then are locally asymptotically stable for (, ).
Differentiating (19) with respect to , we obtain
It follows that
From (35), we have
Thus, we have
Hence, there is a Hopf bifurcation at , . Therefore, if and , then is locally asymptotically stable for , is unstable for , and there is a periodic solution around for . If and or , then is unstable for any , .
We consider the stability of .
Since , if , then (19) has only one imaginary solution. From proof of , . Hence, is unstable for any . If and , then (19) has only one imaginary solution. If and , then (19) does not have imaginary solution. Wnen , is unstable. Therefore, is unstable for any .
We consider the stability of .
Consider the characteristic equation (19). is a root of (19) for all . Suppose is a root of (19), then substituting in the characteristic equation, we have
and separating real and imaginary parts, then we get
From (39), we obtain
(i) Assume . Suppose (19) has a root , for all . Then, from (40), we have
The left-hand size of (41) is equal to
Obviously,
and by assumption,
Therefore,
According to (41), this is impossible, since we are assuming . Hence, all roots of (19) have nonpositive real parts, and this implies that the trivial solutions of the linearized system of (4) about are stable.
(ii) Assume . Consider the following real function:
It is clear that and . There exists a such that if , , we also have
Since yields , which gives , hence, there exists a such that when , . Therefore, there exists at least one , , such that , that is, (19) has at least a positive root. This indicates that trivial solution of the linearized system of (18) about is unstable.
4.3. Numerical Simulation
For numerical simulation, let us first consider case (i) in Theorem 4. We choose parameter values as , , , and . For this choice, , . The stability of the interior equilibrium point depends upon the magnitude of delay . From (30) and (31) we obtain , . Consider , . When , Figures 2(a) and 2(b) show that the system is approaching the equilibrium point . When increases to , becomes unstable (Figure 2(c)). Stable oscillations appear when and the Hopf bifurcation periodic solutions are depicted in Figure 3.
(a)
(b)
(c)
For case (ii) in Theorem 4, consider the following set of parametric values: , , , and . Then we have , . By Theorem 4, the system (4) has two equilibria and . From direct calculation, we get that , , , , , , and . For , the numerical simulation shows that equilibrium point is stable (Figure 4). When increases to , the equilibrium becomes unstable (Figure 6). There is a periodic solution around the equilibrium when (Figure 5).
For case (iii) in Theorem 4, let us take the parametric values , , , and . For these set of values, by direct calculation, we get , , and . By Theorem 1, system (4) has only one equilibrium , and from Theorem 4 we know that the equilibrium is stable (Figure 7(a)). Take parametric values as , , , and . Then we have , , , and . By Theorem 2, system (4) has only one equilibrium . The numerical simulation shows is unstable (Figure 7(b)).
(a)
(b)
5. Permanence
In this section, we will use the same methods as [35] to prove the permanence of system (4).
Theorem 5. Suppose that , , and . Then, for any time delay , system (4) is permanent.
Proof. From Theorem 1, we see that it is enough to consider the solution with initial condition . Theorem 1 implies that the omega limit set of is nonempty, compact, and invariant and . It follows from definition of permanence and Theorem 1, that we only need to show
Here is some positive constant which dose not depend on the initial function .
Let us first show
In fact, if (49) is not true, then, from , we see that there exists a positive time sequence such that
Note that the solution is bounded on by Theorem 1. It follows from (4) that is uniformly continuous on . Hence, from Ascoli’s theorem there is a subsequence of , still denoted by , such that
holds uniformly on in the wide sense. From Theorem 1, we have that for any and that, for any , the function of is the solution of (4) with the initial function . Here we note that and for any .
We claim that for any . From Theorem 1, we know that if , then the solution of (4) exists and and . Thus, from , we have that for any . Then, it follows from (4) that for any , and for any . Hence,
Since is bounded for , we must have , which implies that for . It follows from (4) and the invariant of that for any . This shows that above claim holds. In particular, we have that
For sufficiently small and sufficiently large , , we have
Hence, for large , from (50) and above inequality, we obtain
which is a contradiction to . This completes the proof of .
Next, let us show
For any initial functions sequence , let be the solution of (4) with the initial function . Let be the omega limit set of . We have that there exists some compact and invariant set such that as . Here, means Hausdor distance [35, 36].
If (48) does not hold, for some initial function sequence such that , we have that there is some such that for some . Now, let be the solution of (4) with the initial function . Then, by the invariant of , we have that for all and . Note and the positivity of all solutions, we easily have that for all . Hence, it follows from (4) that , , and for all . This implies that , , for all . If , we see that the negative semiorbit is unbounded. This is contradiction.
If , we have that , for all . This shows that . Let us show is factually isolated [35, 36]. that is, there exists some neighborhood of in such that is the largest invariant set in . In fact, let us choose
for some sufficiently small positive constant and . We will show that is the largest invariant set in for some .
If not, for any sufficiently small , there exists some invariant set such that is not empty. Let and be the solution of (4) with the initial function . Then, for all .
If , by the invariance of and Theorem 1, we also have the contradiction that or that the negative semiorbit of (4) through is unbounded.
If , from the Theorem 1, we see that for all . Now, let us consider the continuous function:
for some constant . Because of , we have . The time derivative of along the solution satisfies
Since and , we have that . We can choose , such that . From (49), there exists a constant such that for some constant and . Hence, it follows from (59)
for all . Thus, as . This contradicts Theorem 1 and shows that is isolated.
We easily see that the semigroup defined by the solution of (4) satisfies the conditions of Lemma 4.3 in [36] with . Thus, by Lemma 4.3 in [36], we have that there is some , such that . Here, denotes the stable set of .
If , again by the invariant of and Theorem 1, we also have the contradiction that or that the negative semiorbit of (4) through is unbounded.
If , from Theorem 1, we see that , for all . It follows from that , , which contradicts (49). This shows that (48) holds. Thus, (4) is permanent. This completes the proof of Theorem 5.
6. Discussion
In a real ecological context, the interaction between phytoplankton and zooplankton will not be essentially instantaneous. Instead, the response of zooplankton to contact with phytoplankton is likely to be delayed due to a gestation period. Another fact, during the interaction between phytoplankton and zooplankton, is that the liberation of toxic substances by phytoplankton must be mediated by some time lag which is required for the maturity of toxic-phytoplankton. And in some cases, very high substrate concentrations in the lakes actually inhabit the growth of phytoplankton cells. Moreover, with the substrate concentrations increasing unlimitedly, some kind of microorganism will die eventually [32]. Based on the above fact, in this paper, we introduce time delay to a phytoplankton-zooplankton interaction model with exponential substrate uptake and exponential distribution of toxic substance term. The model (4) accounts for some natural phenomenon. By using comparison principle for functional differential equations and traditional analysis technique for transcendental equations [34], we give a detailed analysis on boundedness of solutions of system (4) and local asymptotic stability of the equilibria of system (4). Our results show that time delay is factually harmless for the local asymptotic stability of the zooplankton free equilibrium of (4), but it is not always harmless for the positive equilibrium; that is to say, because of the time delay the positive equilibrium becomes unstable (Theorem 4). Based on some known techniques on limit sets of differential dynamical systems, we show that, for any time delay, the phytoplankton-zooplankton interaction model is permanent if and only if one positive equilibrium exists.
Under the aids of the numerical simulation we further investigate the delayed model system (4). Figures 1–6 illustrate that the delay plays crucial role in determining the asymptotic behavior of solutions of model (4). The stability or oscillatory coexistence depends upon not only the parametric restriction but also upon the gestation delay (liberation delay). Our system exhibits stable behavior when , where is the threshold value of the parameter . This threshold value of is determined by our numerical simulations. This value of is 3.5 for case (i) of Theorem 4 and is 5.1 for case (ii). These stable solutions are shown in Figures 2(a), 2(b), and 4. When crosses , there is a delay induced instability demonstrated in Figures 2(c) and 5 for . Stable oscillations appear when and the Hopf bifurcation periodic solution occurs (see Figures 3 and 5). Thus, there is a range of gestation delay (liberation delay) which initially imparts stability, then induces instability, and ultimately leads to periodic behavior.
Conflict of Interests
All authors declare that there is no conflict of interests regarding the publication of this paper.
Acknowledgments
The authors are grateful to the anonymous reviewers for their careful reading and valuable suggestions that greatly improved the quality of this paper. This work is supported by the National Natural Science Foundation of China (Grants no. 11261058, 11261056, 11271312, and 11361059).