Mathematical analysis of a nutrient-plankton system with delay.

A mathematical model describing the interaction of nutrient-plankton is investigated in this paper. In order to account for the time needed by the phytoplankton to mature after which they can release toxins, a discrete time delay is incorporated into the system. Moreover, it is also taken into account discrete time delays which indicates the partially recycled nutrient decomposed by bacteria after the death of biomass. In the first part of our analysis the sufficient conditions ensuring local and global asymptotic stability of the model are obtained. Next, the existence of the Hopf bifurcation as time delay crosses a threshold value is established and, meanwhile, the phenomenon of stability switches is found under certain conditions. Numerical simulations are presented to illustrate the analytical results.

Background

Plankton refers to organisms that are living in water bodies (such as oceans, lakes, rivers and ponds) freely drifting and weakly mobile (Abdllaoui et al. 2002; Odum 1971). Plant forms of plankton community are known as phytoplankton, they serve as the basic food source and occupy the first trophic level of all aquatic food chains. Animals in the plankton community are known as zooplankton. They consume phytoplankton which are their most favourable food source. Phytoplankton are not only the basis for all aquatic food chains, but also they do huge services for our earth by supplying the essential oxygen and absorbing the harmful carbon dioxide which contributes to global warming (Odum 1971). In addition to these benefits phytoplankton act as the biological indicators of water quality. Excess blooming of the phytoplankton will deteriorate the water quality. For example, increase of phytoplankton population in lakes (reservoirs), especially the extension of the growing season and over growing of the cyanobacteria are important causes for eutrophication in lakes. Eutrophication refers to the enrichment of an ecosystem with chemical nutrients such as nitrogen, phosphate and so on, leading to the over growth of biomass and their rapid reproduction in water bodies. This leads to the decrease of dissolved oxygen in water, which in turn results in the death of aquatic organisms. These dead aquatic organisms get settled at the bottom of the lake and then decomposed by microorganisms which once again consume a large amount of dissolved oxygen. Consequently, the dissolved oxygen content of the water body is further reduced and the water quality deteriorates further. This process affects the survival of aquatic organism and greatly accelerates the process of eutrophication in water bodies. The occurrence of the eutrophication, because of a large amount of reproduction of plankton, often makes the water bodies appear in different colors such as blue, red, brown, white, and so on. This phenomenon occurring in water bodies is called “algae bloom” and “red tide” in sea. These algae are foul smelling, poisonous and can’t be eaten by fish. And also they prevent sunlight from reaching the submerged plants and leading to their death by hindering their photosynthesis. These dead submerged plants releasing nitrogen, phosphorus and other nutrients after decaying and then the algae use these nutrients. Because of the high biomass accumulation or the presence of toxicity, some of these blooms, more adequately called “harmful algal blooms” (Smayda 1997), are noxious to marine ecosystems or to human health and can produce great socioeconomic damage. Therefore, the study of marine plankton ecology is an important consideration for the survival of our earth.

Due to the difficulty of 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 (Edwards and Brindley 1999). The problems of zooplankton–phytoplankton systems have been discussed by many authors (Rose 2012; Saha and Bandyopadhyay 2009; Chakraborty and Dasb 2015; Yunfei et al. 2014; Rehim and Imran 2012; Ruan 1995) in resent years. These systems can exhibit rich dynamic behavior, such as stability of equilibria, Hopf bifurcation, global stability, global Hopf bifurcation and so on. However, the importance of nutrients to the growth of plankton leads to explicit incorporation of nutrients concentrations in the phytoplankton–zooplankton models. Therefore, a better understanding of mechanisms that determine the plankton is to consider plankton–nutrient interaction models. Recently, a nutrient–plankton model system for a water ecosystem is proposed by Fan et al. (2013) and its global dynamics behavior under different levels of nutrient has been studied. He and Ruan (1998), Zhang and Wang (2012), Pardo (2000) studied nutrient–phytoplankton interaction model and observed the global behavior of the system. Huppert et al. (2004) studied a simple nutrient–phytoplankton model to explore the dynamics of phytoplankton bloom. Huppert et al. (2004) provided a full mathematical investigation of the effects of three different features in an excitable system framework.

The understanding of the dynamic of plankton–nutrient system becomes complex when additional effects of toxicity (caused due to the release of toxin substances by some phytoplankton species known as harmful phytoplankton) are considered. The role of toxin and nutrient on the plankton system have been discussed by many researchers (Chakarborty et al. 2008; Pal et al. 2007; Khare et al. 2010; Jang et al. 2006; Chowdhury et al. 2008; Upadhyay and Chattopadhyay 2005; Chatterjee et al. 2011).

Time delays of one type or another have been incorporated into biological models by many researchers (Aiello and Freedman 1990; Chen et al. 2007; Cooke and Grossman 1982; Hassard et al. 1981; Song et al. 2004). In general, delay differential equations exhibit much more complicated dynamics than ordinary differential equations since a time delay could cause a stable equilibrium to become unstable and induce oscillations and periodic solutions. Therefore, more realistic models of population interactions should take into account the effect of time delays. The interaction of plankton–nutrient model with delay due to gestation and nutrient recycling has been studied by Ruan (1995) and Das and Ray (2008). Chattopadhyay et al. (2002) proposed and analysed a mathematical model of toxic phytoplankton–zooplankton interaction and assumed that the liberation of toxic substances by the phytoplankton species is not an instantaneous process but is mediated by some time lag required for maturity of species. Extending the work of Chattopadhyay et al. (2002), Saha and Bandyopadhyay (2009) and Rehim and Imran (2012) have studied the global stability of the toxin producing phytoplankton–zooplankton system.

The effect of nutrient recycling on food chain dynamics has been extensively studied. Nisbet et al. (1983), Ruan (1993), Angelis (1992) and Ghosh and Sarkar (1998) studied the effect of nutrient recycling for ecosystem. In their model the nutrient recycling is considered as an instantaneous process and the time required to regenerate nutrient from dead organic is neglected. Beretta et al. (1990), Bischi (1992) and Ruan (2001) studied nutrient recycling model with time delay. They have performed a stability and bifurcation analysis of the system and estimated an interval of recycling delay that preserves the stability switch for the model.

In the present paper, motivated by the above work, a model for the nutrient–plankton consists of dissolved nutrient (N), phytoplankton (p) and herbivorous zooplankton (z) is considered. We assume that the functional form of biomass conversion by the zooplankton is Holling type-II and the predator is obligated that is they dose not take nutrient directly. The toxic substance term which induces extra mortality in zooplankton is also expressed by Holling type II functional form.

In order to account for the time needed by the phytoplankton to mature after which they can release toxins, a discrete time delay is incorporated into the system. Moreover, the discrete delays also indicate the partially recycled nutrient decomposed by bacteria after the death of biomass. The models in Fan et al. (2013)and Das and Ray (2008), time required to regenerate nutrient from dead organisms is neglected. Also the term of toxin liberation has not take into their model. But these are one of the most important features in the real ecosystem (Sarkara et al. 2005; Chattopadhayay et al. 2002; Mukhopadhyay and Bhattacharyya 2010). In comparison with literature (Fan et al. 2013; Das and Ray 2008), the model proposed in this paper is more general and realistic.

The organization of the paper is as follows. In next section, a nutrient–plankton delay differential equations with delay will be proposed and its boundedness criteria will be established. In “Equilibria and its stability” section, we analyze the dynamical properties such as existence of the equilibria and its stability, possible bifurcations with variation of the parameters. In “Numerical simulation” section, numerical studies of the models are performed to support our analytical results. Discussion are drawn in the final section.

The model

Let N(t), p(t) and z(t) are the concentration of nutrient, phytoplankton and zooplankton population at time t, respectively. Let be the constant input of nutrient concentration and D be the washout rates for nutrient, phytoplankton and zooplankton, respectively. The constant delay parameter , and are considered in the decomposition of phytoplankton population, zooplankton population and the discrete time period required for the maturity of toxic-phytoplankton, respectively. With these assumptions, we write the following system of delay differential equations describing nutrient–plankton interaction

We assume that all parameters are non-negative and are interpreted as follows:Here we observe that, if there is no delay (i.e., ) and , then . If and , then we also get . Hence, throughout our analysis, we assume that

—nutrient uptake rate for the phytoplankton

—the maximum zooplankton ingestion rate

—the natural death rate of phytoplankton

—the natural death rate of zooplankton

—the nutrient recycle rate after the death of phytoplankton population

—the nutrient recycle rate after the death of zooplankton population

—the conversion factor from nutrient to phytoplankton

—the conversion factor from phytoplankton to zooplankton

—the rate of toxic substances produced by per unit biomass of phytoplankton

a—the half saturation constant

—the intra-specific competition coefficient or the density dependent mortality rate of phytoplankton population.

As pointed out by Holling (1965), Ma (1996) and Das and Ray (2008), the functional response of Holling type I is applied to lower organisms, for example, alga and unicellular organism. Therefore, in this paper we let the functional response of phytoplankton to nutrient be Holling type I.

As phytoplankton is the most favorable food source for zooplankton within aquatic environments and the Holling type-II functional form is a reasonable assumption to describe the law of predation (Chattopadhyay et al. 2002; Ludwig et al. 1978; Das and Ray 2008). It is quite reasonable to assume that the law of grazing must be same whether it contributes toward the growth of zooplankton species or it suppresses the rate of grazing due to presence of toxic substances.

In fact, the liberation of toxic substance by phytoplankton is not an instantaneous phenomenon, since it must be mediated by some time lag which is required for the maturity of toxic-phytoplankton. However, the liberation of toxic substance at the time t depends on the population size of toxic phytoplankton species at time . So, the zooplankton mortality due to the toxic phytoplankton is described by the term . In model (1), the term describe the distribution of toxic substance which ultimately contributes to the death of zooplankton populations.

Our model consider that the phytoplankton have competition among themselves for their survival (Barton and Dutkiewicz 2010; Jana et al. 2012; Ruan et al. 2007; Wang et al. 2014). is the reduction term for the phytoplankton population.

From the standpoint of biology, we are only interested in the dynamics of model (1) in the closed first octant .

In accordance with the biological meaning, the initial conditions of the system (1) are taken as followswhere .

Regarding the positivity and boundedness of the solution for the system (1) we state the following theorem.

Theorem 1

All solutions of system (1) with initial conditions (2) are positive and bounded.

Proof

The proof of positivity of the solutions of system (1) is easy, so we omit it here. As for boundedness of the solutions of (1), we define function

Derivative of X(t) with respect to system (1), we obtain

Therefore, for all large t, where is an arbitrarily small positive constant. Thus, N(t), p(t) and z(t) are ultimately bounded by some positive constant.

Equilibria and its stability

Equilibria

System (1) possesses three possible nonnegative equilibria, namely the extinction equilibrium , the zooplankton-eradication equilibrium and the coexistence equilibrium . For the zooplankton-eradication equilibrium , the and satisfy the following equation:

From first equation of system (3) we have

On substituting (4) into second equation of (3) we derive that

If , then Eq. (5) has exactly one positive real rootwhere .

For the coexistence equilibrium , the , and satisfy the following equation:

From third equation of system (6) we have

Again from first and second equations we have

Let . If , and , then and . Thus

From above analysis we obtain the following theorem.

Theorem 2

The extinction equilibriumalways exists. Furthermore, suppose that. Then the zooplankton-eradication equilibriumexists, and the unique coexistence equilibriumexists only ifand.

In what follows, we will analysis the stability of the system (1) around different equilibria.

Model (1) without delay

In this subsection, we give the basic dynamical behavior of system (1) without delay.

Theorem 3

If, then the extinction equilibriumis locally asymptotically stable andunstable if.

Suppose that. If, then thezooplankton-eradication equilibriumis locally asymptotically stable andunstable if.

Suppose that the coexistence equilibriumexists. Then it is locally asymptotically stable if the following inequality hold

Proof

The characteristic equation about is given byIt is clear that Eq. (8) has negative root and . So, if , then

From this we have that the extinction equilibrium is locally asymptotically stable. If , then is unstable.

The characteristic equation about is . If , then Further, and . Which implies that and . Therefore, if , then the zooplankton-eradication equilibrium is locally asymptotically stable and unstable if .

The characteristic equation about is given bywhere

If , and , then

Therefore, all roots of (9) have negative real parts. By the Routh–Hurwitz criterion we obtain that the coexistence equilibrium is locally asymptotically stable.

Remark 1

From above analysis we see that the input concentration of the nutrient, density dependent mortality rate of phytoplankton population and the death rate of the plankton play an important role in controlling the dynamics of the system.

After studying the local stability behavior we perform a global analysis around the equilibrium point.

Theorem 4

If, then the extinction equilibriumis globally asymptotically stable.

Define a positive definite functionCalculating the derivative of along the positive solution of system (1) we haveSince N(t), p(t) and z(t) are positive, if , then . if and only if . Thus is globally asymptotically stable by Lyapunov–LaSalle invariance principle.

Remark 2

Theorem 4 shows that too low of the conversion rate of the plankton will cause species extinction. This is consistent with the real ecosystem.

For the globally asymptotically stability of the equilibrium , we first consider the transformations , , . With these transformations, the model (1) reduces to

Then, (0, 0, 0) is an equilibrium point of (10). Define a positive functionwhere , are to be chosen. Now, calculating the derivative of along the positive solution of system (10) we haveUsing the inequalitywe haveSetwith , . By choosing , properly it is possible to set and such that , that is , we can choose and such thatSo, if (12) holds, . if and only if . Thus by Lyapunov–LaSalle invariance principle we obtain the following theorem

Theorem 5

Suppose that the equilibrium pointof system (1) exists. Then it is globally asymptotically stable if (12) holds, where, , , , are given by (11).

Let us consider the transformations , , . With these transformations, the model system (1) reduces to

Then, (0, 0, 0) is an equilibrium point of (13). Define a positive functionwhere , are to be chosen. Now, calculating the derivative of along the positive solution of system (13) we haveUsing the inequalitywe haveSetwith , . By choosing , properly it is possible to set and such that , that is, we can choose and such that

So, if (15) holds, . if and only if . Thus by Lyapunov–LaSalle invariance principle we obtain the following theorem.

Theorem 6

Suppose that the equilibrium pointof system (1) exists. Then it is globally asymptotically stable if condition (15) hold, where, , , , are given by (14).

Model (1) with delay

In this section, we discuss the asymptotic stability of coexistence equilibrium and the existence of Hopf bifurcations of the delayed model (1). To simplify the analysis, it is assumed that all the delays are of equal magnitude, i.e. and , namely reconversion of dead zooplankton biomass into nutrient is ignored.

We need the following result which was proved in Ruan and Wei (2003) by using Rouches theorem and it is a generalization of the lemma in Dieudonne (1960).

Lemma 1

Consider the exponential polynomialwhereandare constants. Asvary, the sum of the orders of the zeros ofon the open right half plane can change only if a zero appears on or crosses the imaginary axis.

From “Model (1) without delay” section 3.2 we know that the coexistence equilibrium is locally asymptotically stable for if (7) holds. For , the linearization of system (1) at is

Then the associated characteristic equation of (16) iswhere , , , , .

In the following, we study the Hopf bifurcation of the coexistence equilibrium. Now for , if is a root of , then we have

Separating the real and imaginary parts, we haveAdding up the squares of both equations, we obtainDenot , then (19) becomeswhere , and . Let

By the ideal of Li and Wei (2005), Ruan and Wei (2001), Song and Wei (2004), in what follows, we study the distribution of the zeros of (20). From we can easily get the following lemma.

Lemma 2

Equation (20) has at least one negative real root.

Lemma 3

If, then Eq. (20) has no positive roots.

From (21) we have . SetThen the roots of Eq. (22) can be expressed as

If , then (22) has no real roots or exists one root. So the function g(r) is monotone increasing with r. Therefore, Eq. (20) has no positive real roots due to .

Obviously, if , then is the local minimum of g(r). Thus, we get the following result.

Lemma 4

Equation (20) has positive roots if and only ifand.

The sufficiency is obvious. We only need to prove the necessity. Otherwise, we assume that either or and . Since g(r) is increasing for and , g(r) has no positive real zeros for . If and , since is the local maximum value, it gives that . Hence, by we obtain that g(r) has no positive real zeros. This completes the proof.

From above discussion, we get the following lammas.

Lemma 5

Suppose thatis defined by (23).

Eq. (20) has at least one negative real root.

If, then Eq. (20) has no positive roots.

Eq. (20) has positive roots if and only ifand.

Assume that Eq. (20) has positive roots. Without loss of generality, we suppose that it has two positive roots, denoted by , , respectively. Then (19) has two positive roots, sayBy (18) we haveLetwhere

Then is a pair of purely imaginary roots of (17), . define

Therefore, applying Lemmas 1 and 5 to (17), we obtain the following lemma.

Lemma 6

Suppose that the inequality (7) holds. Then we have

If, then all roots of equation (17) have negative real parts for all.

If, and, then Eq. (17) has a pair of imaginary roots. Furthermore, if, then all roots of equation (17) have negative real parts.

Let be the root of (17) near satisfying , . Let . Then we have the following transversality condition.

Lemma 7

Suppose. If the conditions of Lemma 6(b) are satisfied, then, andhave the same sign.

Differentiating (17) with respect to , we obtainIt follows thatThenFrom (19), we have

Thus, we haveSince , we conclude that , and have the same sign. This completes the proof.

Theorem 7

Suppose that the inequality (7) holds.

If, then the coexistence equilibriumof system (1) is asymptotically stable for all.

If, andhold, then system (1) at the equilibriumis asymptotically stable for, and unstable when. System (1) undergoes a Hopf bifurcation atwhen.

Numerical simulation

To substantiate analytical findings a set of hypothetical parameter values have been considered for numerical simulation (see Table 1). Most of the parameters in Table 1 used by authors in Chattopadhayay et al. (2002) and Fan et al. (2013).

Table 1

Parameter values used in numerical simulation

Parameters	Symbols	Values
Dilution rate	D	0.4(1day^-1)
Constant input of nutrient concentration	N_0	40 (mg dm^-1 )
Nutrient uptake rate for the phytoplankton	α	0.7(1day^-1)
Maximum zooplankton ingestion rate	β	0.6(1day^-1)
Conversion factor from death phytoplankton	m_1	0.8
Conversion factor from death zooplankton	m_2	0.5
Natural death rate of phytoplankton	d_1	0.025(day^-1)
Natural death rate of zooplankton	d_2	0.02(day^-1)
Conversion factor from nutrient to phytoplankton	k_1	0.9677
Conversion factor from phytoplankton to zooplankton	k_2	0.9661
Toxin-production rate	k_3	0.0186(day^-1)
Half-saturation coefficient	a	2(mgdm^-1)
Intra-specific competition coefficient	υ	0.1

First, we consider the special case of system (1), that is, and . In order to verify the results of Theorem 7, we consider as bifurcation parameter and for case (a) taken parameters in Table 1. It is easy to compute that Our numerical simulations show that for all , interior equilibrium is stable. Figure 1 shows the simulation result for the system (1) with . For case (b) of Theorem 7, we take parameters as , , and other parameters the same as that in Table 1. A direct computation gives , and holds. After calculations we find the minimum value of the delay parameter ‘’ for system (1) for which the stability behaviour changes and the first critical values are given by , such that is locally stable for and is unstable for . From our analytical findings we have seen that is locally asymptotically stable for . Figure 2 shows the simulation result for system (1) with . Interior equilibrium point looses its stability as passes through its critical value and a Hopf bifurcation occurs. A periodic solution is depicted in Fig. 2d, e.

Fig. 1

a, b The asymptotical stability of the coexistence equilibrium with . Here, initial value is (5, 3, 15 )

Fig. 2

a, b The asymptotical stability of the coexistence equilibrium with . c–e Coexistence equilibrium loses its stability when . Stable periodic solution arising from Hopf bifurcation at . Here

Next, We present some numerical results on the case of system (1) that and . Take , , , , and other parameters the same as that in Table 1. With the help of this parameter set we obtain the interior equilibrium as . Let us fix and gradually increase the value of . After some calculations one can find the minimum value of the delay parameter “” for the model system (1) for which the stability behaviour changes and the first critical values are given by , such that is stable for and unstable for . Figure 3 shows the simulation result for the model system (1) with . Interior equilibrium point looses its stability as passes through its critical value and a Hopf bifurcation occurs, a stable Hopf-bifurcating periodic solution is depicted in Fig. 3d, e. The equilibrium point remains locally asymptotically stable whenever the delay parameter lies in the range (6.759, 10.6156). again switches from stability to instability as passes through and an unstable solution for the model system (1) is shown in Fig. 4. The numerical simulations we have done here illustrate the stable periodic solution arising from Hopf bifurcation at and , respectively, and the switching of stability that occurs as the magnitude of the delay parameter increases gradually. For the above set of parameter values, when fixing and varying the value of or fixing and varying , the dynamical behavior of the system (1) explored by numerical simulation are the same as above, so we omit it here.

Fig. 3

a, b The asymptotical stability of the coexistence equilibrium for with . c Coexistence equilibrium loses its stability at with . Stable periodic solution arising from Hopf bifurcation at with . d, e Stable limit cycle is observed at with . Here and

Fig. 4

a, b The coexistence equilibrium becomes stable for with . c Coexistence equilibrium loses its stability at with . Stable periodic solution arising from Hopf bifurcation at . d, e Stable limit cycle is observed at with . Here

Conclusions and discussion

In the present analysis, we have proposed and analyzed a three component model consisting of nutrient, phytoplankton and zooplankton. It is assumed that the grazing on phytoplankton , zooplankton growth rate and the zooplankton mortality due to the toxin phytoplankton are Holling type II forms. According to the facts that reconversion of dead biomass into nutrient is not an instantaneous process, but is mediated by some time lag required, and the toxin liberation by the phytoplankton species also need time period, our model in present paper incorporate delayed nutrient recycling and delayed toxic liberation. In comparison with literatures (Fan et al. 2013; Das and Ray 2008), the model (1) in this paper is more general and realistic.

In the absence of the time delay, the dynamical behavior of system (1) was studied extensively around all feasible equilibria. Conditions were also derived both for the local and global stability of the system at all possible equilibria. Theorem 4 indicates that if the conversion rate from nutrient to phytoplankton and phytoplankton to zooplankton lower than certain values, then plankton will extinct. This result is consistent with real ecosystem. Theorem 5 shows that a high concentration of the input nutrient (together with a high mortality rate of the zooplankton population) will cause eradication of the zooplankton. Theorem 6 reveals that low values of mortality rate of the both phytoplankton and zooplankton population ensures coexistence of the plankton. Thus, the concentration of the input nutrient, the mortality rate of the plankton plays a major role in controlling the local and global dynamics of the basic model around the various stationary states.

Next we have studied the model with discrete delay in the term modeling plankton recycling and the term of toxin liberation. Numerically it is shown that the behavior of the system around the interior equilibrium depends on the time delay. When we fix time delay and gradually increase the value of , the numerical simulations which we have performed show that there are threshold limit such that as the time delay crosses the threshold value , the delayed nutrient–plankton system enters into a Hopf bifurcation and we have a periodic orbit around the coexisting equilibrium point . The interior equilibrium point is stable whenever and unstable for . This phenomenon is known as switching of stability which arises for our model system. The most interesting as well as mathematically important results we have presented in this paper is the stability criteria for the Hopf-bifurcating periodic solution by considering the discrete time lag as bifurcation parameter. These findings demonstrate the delayed effect of plankton and the cyclic nature of blooms in this nutrient–plankton system. This is one of the most important findings of our analysis. In their analysis, Fan et al. (2013) and Das and Ray (2008) were unable to exhibit the periodic nature of blooms by considering non-delayed nutrient–plankton system. From Figs. 3c and 4c we find the solutions oscillate around . Figures 3c and 4c show that the plankton system can occurs the peak phenomenon, which corresponds to blooms, and also occurs the valley effect, which corresponds the low values of phytoplankton. Our mathematical and numerical results provide certain threshold values for the delay parameters for which we can maintain a stable situation for all the species and can control bloom dynamics.

Now let us make a comparison with result of previous studies and present study. Fan et al. (2013) investigated a nutrient–plankton system with nutrient recycling from dead plankton, but the time required to regenerate nutrient from dead organic is neglected. Besides, the effects of the toxin which produced by phytoplankton did not take account into their model. Sharma et al. (2014) studied a nutrient–toxin phytoplankton–zooplankton model with nutrient recycling, but there is only nutrient recycling from dead phytoplankton and the recycling assumed to be instantaneous. The model studied by Fan et al. (2013) has an unique interior equilibrium under the condition , which also ensures local asymptotic stability of the interior equilibrium . Our results shows that conditions and ensure the existence of the interior equilibrium . But for locally asymptotically stability of the , we need condition (7). Therefore, here the coexistence equilibrium in this setting possesses more restrictive existence and stability condition, since they involve the intra-specific competition parameter and toxin liberation parameter , see Theorems 2 and 3. The local and global stability of the interior equilibrium have not studied by Sharma et al. (2014). While the globally asymptotically stability of the interior equilibrium have been proved by Fan et al. (2013) only for the special case of model (1) (with . In this paper, we obtained sufficient conditions which ensures for the interior equilibrium of the model (1.2) to be globally asymptotically stable. This can be seen as one of the novelty of this paper. Moreover, the results obtained by Fan et al. (2013) indicate that the concentration of the input nutrient and the initial conditions of the nutrient–plankton model are the two important factors on the dynamics of the system behavior. But here, our results obtained in this paper indicate that except for concentration of the input nutrient and initial values of the system (1.2), the intra-specific parameter and toxin liberation parameter also affect the dynamical properties of the model (1.2). Comparing with paper (Sharma et al. 2014), a Hopf-bifurcation arises also at the interior equilibrium, but the conditions for its occurrence here at are more restrictive, involving also the intra-specific competition, recycling of the zooplankton and toxin liberation. Besides, differently from literatures (Fan et al. 2013; Sharma et al. 2014), our numerical investigations show that the nutrient recycling delays can induce stability switches, such that the interior equilibrium switches from stable coexistence equilibrium to stale periodic orbits, to stable coexistence equilibrium again and so on(see Figs. 3, 4). This phenomenon is ecologically important and especially can lead to potentially dramatic shifts to the system dynamics. In biological terms, our finding has ecological significance in the estuarine system. There are jungles and forest adjacent to the estuary, which are the main source of productivity. The nutrients come from the litterfall which can be decomposed after a period of time. The tide not only collects the nutrient from the litters but also mixes them into the estuarine water (Wang and Wang (2007)). Our research indicates that delays in the decomposition of litterfall cause destabilization of this system.

Unfortunately, we cannot give a complete mathematical analysis of the asymptotic stability of the positive equilibrium for model (1) with different delay, i.e., and . We shall leave the problems as future work.