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.

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;. 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 phytoplanktonzooplankton 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 N 0 be the constant input of nutrient concentration and D be the washout rates for nutrient, phytoplankton and zooplankton, respectively. The constant delay parameter τ 1 , τ 2 and τ 3 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: α-nutrient uptake rate for the phytoplankton β-the maximum zooplankton ingestion rate d 1 -the natural death rate of phytoplankton d 2 -the natural death rate of zooplankton m 1 -the nutrient recycle rate after the death of phytoplankton population 0 < m 1 < 1 m 2 -the nutrient recycle rate after the death of zooplankton population 0 < m 2 < 1 k 1 -the conversion factor from nutrient to phytoplankton 0 < k 1 < 1 k 2 -the conversion factor from phytoplankton to zooplankton 0 < k 2 < 1 k 3 -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 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 t − τ 3 . So, the zooplankton mortality due to the toxic phytoplankton is described by the term p(t − τ 3 )z(t). In model (1), the term ρp(t−τ 3 )z(t) a+p(t−τ 3 ) 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). υp 2 is the reduction term for the phytoplankton population. (1) Here we observe that, if there is no delay (i.e., τ i = 0 ) and k 2 < k 3 , then ż < 0. If k 2 > k 3 and β(k 2 − k 3 ) − (D + d 2 ) < 0, then we also get ż < 0. Hence, throughout our analysis, we assume that From the standpoint of biology, we are only interested in the dynamics of model (1) in the closed first octant R 3 + . In accordance with the biological meaning, the initial conditions of the system (1) are taken as follows where τ = max{τ 1 , τ 2 , τ 3 }.
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, X < N 0 + ǫ 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.
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.

Then it is locally asymptotically stable if the following inequality hold
Proof The characteristic equation about E 0 (N 0 , 0, 0) is given by Therefore, all roots of (9) have negative real parts. By the Routh-Hurwitz criterion we obtain that the coexistence equilibrium E * (N * , p * , z * ) 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.
Proof Define a positive definite function Calculating the derivative of V 0 along the positive solution of system (1) we have = 0 if and only if (N , p, z) = (N 0 , 0, 0). Thus E 0 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.

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. τ = τ 1 = τ 2 = τ 3 , and m 2 = 0, 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).
From "Model (1) without delay" section 3.2 we know that the coexistence equilibrium E * (N * , p * , z * ) is locally asymptotically stable for τ = 0 if (7) holds. For τ � = 0, the linearization of system (1) at E * (N * , p * , z * ) is Then the associated characteristic equation of (16) is In the following, we study the Hopf bifurcation of the coexistence equilibrium. Now for τ � = 0, if = iω(ω > 0) is a root of G( , τ ) = 0, then we have Separating the real and imaginary parts, we have Adding up the squares of both equations, we obtain Denot r = ω 2 , then (19) becomes 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 g(0) = b 3 = a 2 3 − a 2 5 > 0 we can easily get the following lemma. (20)
Obviously, if � > 0, then is the local minimum of g(r). Thus, we get the following result. (20) has positive roots if and only if r 1 > 0 and g(r 1 ) ≤ 0.

Lemma 4 Equation
Proof The sufficiency is obvious. We only need to prove the necessity. Otherwise, we assume that either r 1 ≤ 0 or r 1 > 0 and g(r 1 ) > 0. Since g(r) is increasing for r ≥ r 1 and g(0) = b 3 > 0, g(r) has no positive real zeros for r 1 ≤ 0. If r 1 > 0 and g(r 1 ) > 0 , since is the local maximum value, it gives that g(r 1 ) < g(r 2 ). Hence, by g(0) = b 3 > 0 we obtain that g(r) has no positive real zeros. This completes the proof.
From above discussion, we get the following lammas.  Assume that Eq. (20) has positive roots. Without loss of generality, we suppose that it has two positive roots, denoted by u 1 , u 2 , respectively. Then (19) has two positive roots, say By (18) we have Let where k = 1, 2; j = 0, 1, 2, ....
Then ±iω k is a pair of purely imaginary roots of (17), τ = τ j k . define Therefore, applying Lemmas 1 and 5 to (17), we obtain the following lemma.
Lemma 7 Suppose g ′ (r 0 ) � = 0. If the conditions of Lemma 6 (b) are satisfied, then dτ and g ′ (r 0 ) have the same sign.
This completes the proof.

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

Constant input of nutrient concentration
Nutrient uptake rate for the phytoplankton α 0.7 (1 day −1 ) Maximum zooplankton ingestion rate β 0.6 (1 day −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 (mg dm −1 ) Intra-specific competition coefficient υ 0.1 and g ′ (r 0 ) = −0.0728333887 � = 0 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 τ 0 = 1.7657, such that E * (5.4818, 1.9173, 14.7487) is locally stable for τ ∈ [0, 1.7657) and is unstable for τ > τ 0 . From our analytical findings we have seen that E * is locally asymptotically stable for τ < τ 0 . Figure 2 shows the simulation result for system (1) with τ = 1 < τ 0 . Interior equilibrium point looses its stability as τ passes through its critical value τ = τ 0 and a Hopf bifurcation occurs. A periodic solution is depicted in Fig. 2d, e. Next, We present some numerical results on the case of system (1) that τ 1 � = τ 2 � = τ 3 and m 2 � = 0. Take D = 0.381 day −1 , N 0 = 15 mg dm −1 , k 3 = 0.1 day −1 , a = 1 mg dm −1 , υ = 0.009 and other parameters the same as that in Table 1. With the help of this parameter set we obtain the interior equilibrium as E * (3.7040, 1.8108, 6.8715) . Let us fix τ 2 = 1, τ 3 = 2 and gradually increase the value of τ 1 . After some calculations one can find the minimum value of the delay parameter "τ 1 " for the model system (1) for which the stability behaviour changes and the first critical values are given by τ 0 1 = 3.9465, τ 1 1 = 6.7592 , such that E * (3.7040, 1.8108, 6.8715) is stable for τ 1 ∈ [0, 3.9465) and unstable for τ 1 ∈ [3.9465, 6.759). Figure 3 shows the simulation result for the model system (1) with τ 1 = 1 < τ 0 1 . Interior equilibrium point looses its stability as τ 1 passes through its critical value τ 1 = τ 0 1 and a Hopf bifurcation occurs, a stable Hopf-bifurcating periodic solution is depicted in Fig. 3d, e. The equilibrium point E * (3.7040, 1.8108, 6.8715) remains locally asymptotically stable whenever the delay parameter lies in the range (6.759, 10.6156). E * (3.7040, 1.8108, 6.8715) again switches from stability to instability as τ 1 passes through τ 1 = 10.6156 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 τ 0 1 = 3.9465, τ 1 1 = 6.759 and τ 2 1 = 10.6156, 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 τ 1 , τ 3 and varying the value of τ 2 or fixing τ 1 , τ 2 and varying τ 3 , the dynamical behavior of the system (1) explored by numerical simulation are the same as above, so we omit it here.  Fig. 1 a, b The asymptotical stability of the coexistence equilibrium E * (3.696, 5.6566, 19.3071) with τ = 1.

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  Fig. 2 a, b The asymptotical stability of the coexistence equilibrium E * (5.48181.917314.7487) with τ = 1 < τ 0 . c-e Coexistence equilibrium E * loses its stability when τ = 1.9 > τ 0 . Stable periodic solution arising from Hopf bifurcation at τ = τ 0 . Here τ 0 = 1.7657 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 τ i , τ j and gradually increase the value of τ l (i � = j � = l, i.j.l = 1, 2, 3) , the numerical simulations which we have performed show that there are threshold limit τ k l (l = 1, 2, 3, k = 1, 2, . . .) such that as the time delay crosses the threshold  Fig. 3 a, b The asymptotical stability of the coexistence equilibrium E * (5.48181.917314.7487) for (τ 1 , τ 2 , τ 3 ) with τ 1 = 1 < τ 0 1 . c Coexistence equilibrium E * loses its stability at (τ 1 , τ 2 , τ 3 ) with τ 1 = 4 > τ 0 1 . Stable periodic solution arising from Hopf bifurcation at (τ 1 , τ 2 , τ 3 ) with τ 1 = τ 0 1 . d, e Stable limit cycle is observed at (τ 1 , τ 2 , τ 3 ) with τ 1 = 4.2. Here τ 2 = 1, τ 3 = 2 and τ 0 1 = 3.9465 value τ k l , the delayed nutrient-plankton system enters into a Hopf bifurcation and we have a periodic orbit around the coexisting equilibrium point E * . The interior equilibrium point E * is stable whenever τ ∈ [0, τ 1 l ) ∪ [τ 2 l , τ 3 l ) ∪ · · · and unstable for τ ∈ [τ 1 l , τ 2 l ) ∪ [τ 3 l , τ 4 l , ) ∪ · · · , l = 1, 2, 3. 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 Hopfbifurcating 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 E * . 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  Fig. 4 a, b The coexistence equilibrium E * (5.48181.917314.7487) becomes stable for (τ 1 , τ 2 , τ 3 ) with τ 1 = 8 < τ 1 1 . c Coexistence equilibrium E * loses its stability at (τ 1 , τ 2 , τ 3 ) with τ 1 = 11.56 > τ 1 1 . Stable periodic solution arising from Hopf bifurcation at τ 1 = τ 1 1 . d, e Stable limit cycle is observed at (τ 1 , τ 2 , τ 3 ) with τ 1 = 12. Here τ 1 1 = 10.6156 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 E 2 under the condition N 0 > N 1 +N , which also ensures local asymptotic stability of the interior equilibrium E 2 . Our results shows that conditions ak 1 αm 2 d 2 β < D < aα and 0 < p * < p * 1 ensure the existence of the interior equilibrium E * . But for locally asymptotically stability of the E * , 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 k 3 , 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 E 2 have been proved by Fan et al. (2013) only for the special case of model (1) (with m i = 0, i = 1, 2) . 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 N 0 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 E * 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 E * for model (1) with different delay, i.e., τ 1 � = τ 2 � = τ 3 , and m 2 � = 0. We shall leave the problems as future work.