Competition between a nonallelopathic phytoplankton and an allelopathic phytoplankton species under predation

We propose a model of two-species competition in the chemostat for a single growth-limiting, nonreproducing resource that extends that of Roy [38]. The response functions are specified to be Michaelis-Menten, and there is no predation in Roy's work. Our model generalizes Roy's model to general uptake functions. The competition is exploitative so that species compete by decreasing the common pool of resources. The model also allows allelopathic effects of one toxin-producing species, both on itself (autotoxicity) and on its nontoxic competitor (phytotoxicity). We show that a stable coexistence equilibrium exists as long as (a) there are allelopathic effects and (b) the input nutrient concentration is above a critical value. The model is reconsidered under instantaneous nutrient recycling. We further extend this work to include a zooplankton species as a fourth interacting component to study the impact of predation on the ecosystem. The zooplankton species is allowed to feed only on the two phytoplankton species which are its perfectly substitutable resources. Each of the models is analyzed for boundedness, equilibria, stability, and uniform persistence (or permanence). Each model structure fits very well with some harmful algal bloom observations where the phytoplankton assemblage can be envisioned in two compartments, toxin producing and non-toxic. The Prymnesium parvum literature, where the suppressing effects of allelochemicals are quite pronounced, is a classic example. This work advances knowledge in an area of research becoming ever more important, which is understanding the functioning of allelopathy in food webs.

1. Introduction.The term "plankton" is used to describe freely-floating and weakly-swimming marine and freshwater organisms.It was coined by the German scientist Victor Hensen in 1887 (Thurman [46]).Plankton are divided into broad functional groups, among them phytoplankton that live near the surface of the water where there is sufficient light to support photosynthesis, and zooplankton that feed on other plankton.The microscopic and unicellular plants, phytoplankton, are consumed by zooplankton, the animals, which in turn are eaten by larger organisms.Plankton are at the base of the food chain in the aquatic environment, and are responsible for much of the oxygen present in the Earth's atmosphere: half of the total amount produced by all plant life.They also absorb carbon dioxide from their surrounding environment.The highly diverse nature of phytoplankton 788 JEAN-JACQUES KENGWOUNG-KEUMO communities seems to contradict the competitive exclusion principle, which states that when two species compete for the same resource, only one will survive.The modeling of plankton populations is an essential tool in improving our understanding of the physical and biological processes that contribute to the complexity of these systems.
The term allelopathy was first coined by the Austrian plant physiologist Molisch [29] to explain the effect of ethylene on fruit ripening.Rice [36] defines allelopathy as the effects of one plant (including micro-organisms) on the growth of another plant through the release of chemical compounds (called allelochemicals by Whittaker and Feeny [47]) into the environment.At the 1996 meeting of the International Allelopathy Society (IAS), allelopathy was redefined as any process involving secondary metabolites produced by plants, algae, bacteria, and fungi that influences the growth and development of agricultural and biological systems.
Mathematical modeling of plankton populations goes back to Riley [37] and Hallam [14,15,16].The first mathematical model to represent allelopathic interactions between competing species was introduced by Maynard-Smith [27].Improvements and refinements of Maynard-Smith's model have followed several directions (see for examples Nakamaru and Iwasa [32], An et al. [1] and references therein).
The two functions of allelopathy are autotoxicity (where the toxin affects the growth of the toxic species itself) and phytotoxicity (the effect of toxin on the growth of another species).We assume that the phytotoxic term depends upon the product of square of concentration of non-toxic species with square of concentration of toxic species in accordance with Solé et al. [44].In addition, we suppose that the autotoxic term is a linear function of the concentration of the toxic species as per Sinkkonen [42].
To better understand the dynamics of plankton, one should not neglect the interactions of zooplankton.Many authors have worked in this direction.Edwards and Brindley [11] investigated the bifurcational structure of a simple plankton model with zooplankton mortality modeled by −cZ m , 1 ≤ m ≤ 2. They showed explicitly how cycles can persist for 1 < m < 2. In addition, m = 2 does not preclude the existence of cycles or chaos.Edwards [12] examined the behavior of two nutrient-phytoplankton-zooplankton-detritus models to help understand the factors that most influence the dynamics of such models.He further showed that the addition of a detritus compartment has little impact on the nature of the qualitative dynamics that were found for the corresponding nutrient-phytoplanktonzooplankton model.Mukhopadhyay and Bhattacharryya [31] examined a model of nutrient-phytoplankton-zooplankton interaction with spatial heterogeneity.They proved that phytoplankton species with low diffusivity and zooplankton functional response with half-saturation constant can control algal blooms.Ruan [41] studied plankton nutrient models with both instantaneous and delayed nutrient recycling.He successively chose the nutrient input concentration and the maximal zooplankton ingestion rate as bifurcation parameters to show that the positive equilibrium loses its stability via a Hopf bifurcation as these parameters are varied through respective critical values.Jang and Baglama [22] explored nutrient-phytoplankton interaction with both instantaneous and delayed nutrient recycling and zooplankton mortality modeled by −cZ 2 .Unlike other ecological models for which delays can destabilize the system (see Roy et al. [39] and Piotrowska et al. [34]), their numerical simulations suggested that delayed nutrient recycling can actually stabilize the nutrient-phytoplankton system.Chakraborty and Chattopadhyay [9] proposed and analyzed four models of nutrient-phytoplankton-zooplankton populations to observe the dynamics of such models in the presence of additional food.Here the phytoplankton are toxic to the zooplankton population.However, little has been done in modeling competition, allelopathy, predation, and instantaneous nutrient recycling in the same food chain.The model we propose incorporates competition, allelopathy, predation, and instantaneous nutrient recycling.
With the use of mathematical models, Roy [38] demonstrates theoretically that the stable coexistence of two species competing for a single nutrient in a homogeneous medium would be possible provided one of the two species has a sufficiently strong allelopathic effect on the other.The uptake functions are specified to be Michaelis-Menten in Roy's model.The contents of this paper are largely devoted to extending the results of Roy [38].More specifically, we extend the nutrient-nontoxic phytoplankton-toxic phytoplankton model of Roy [38] to general uptake functions and include a zooplankton species as a fourth interacting component.We carefully structure this work in a way that presents the incremental effects of adding the increased complexity and realism to the model.This paper is organized as follows.In the next section, we describe the main model.We then present some preliminary results, and study existence and local stability of steady states.Ecological interpretations of inequalities are followed by some global results.We further extend this work to include a zooplankton species to study the impact of predation on the competition-allelopathy model.The zooplankton species predates only on the two phytoplankton species which are its perfectly substitutable resources.We reconsider each of the models under the effects of instantaneous nutrient recycling.We utilized Matlab and Mathematica Version 10 to run extended simulations to support our analytical findings.

2.
The Competition-Allelopathy model.We will consider population interactions in a chemostat environment.For a detailed description of the chemostat and its application in biology and ecology, the reader is referred to Hsu et al. [20], Smith and Waltman [43].
Two-species compete exploitatively for a single nonreproducing resource.Our model also incorporates allelopathic effects, and can be written In these equations P i (t) is the biomass of the i th population of phytoplankton in the culture vessel at time t, i = 1, 2. Population P 1 is assumed to be nontoxic, while population P 2 is assumed to be toxic.The concentration of the nonreproducing resource in the culture vessel at time t is denoted by N (t), while N 0 is the concentration of resource N in the feed vessel.
The removal rate m 1 of nontoxic phytoplankton P 1 is the sum of washout rate D and the specific death rate 1 , so that, m 1 = D + 1 .The removal rate m 2 of toxic phytoplankton P 2 is the sum of washout rate D, the specific death rate 2 , and the autotoxic coefficient a 2 , so that, m 2 = D + 2 + a 2 as per Sinkkonen [42].
As in Wolkowicz and Lu [48], it is interesting to note that the analysis of the model requires no assumptions on the signs of the i 's and a 2 , as long as the m i 's all remain positive.This leaves the i 's and a 2 open to other interpretations.For instance, a negative i describes an additional food source for the i th population of phytoplankton while a positive i accounts for further deleterious effects (such as sinking and mixing) on the i th population of phytoplankton.Finally, a zero i means that there is no intrinsic death of the i th population of phytoplankton.A negative, zero, and positive a 2 indicate respectively stimulatory effects, no effects, and inhibitory effects of toxins produced by P 2 on the growth of conspecies.
As in Solé et al. [44], we express the phytotoxic interactions as γP 2 1 P 2 2 , where γ denotes the phytotoxic coefficient.In system (1) the response functions f i (N ) represent the per capita rate of conversion of nutrient to biomass of population P i per unit of population P i as a function of the concentration of nutrient N .We assume that the rate of conversion of nutrient to biomass is proportional to the amount of nutrient consumed, so that the consumption rate of resource N per unit of population P i is of the form 1 γi f i (N ), where γ i is the growth yield constant (number of phytoplankton per unit of nutrient).We make the following assumptions concerning the response functions f i : The break-even concentration for population P i on nutrient N is obtained by setting dPi dt = 0 = f i (N ) − m i and solving for N .By the monotonicity assumptions, the solution λ i is a uniquely defined positive extended real number provided we assume λ i = ∞ if f i (N ) < m i for all N ≥ 0. Finally, let µ i denote the maximal growth rate of population P i on resource N , so that lim Lotka-Volterra kinetics (or Holling type I), Michaelis-Menten kinetics (or Holling type II), and sigmoidal kinetics (Holling type III or multiple saturation dynamics) are prototypes of response functions f i found in the literature (Aris and Humphrey [3], Boon and Laudelout [6], Edwards [12], Jost et al. [23], Monod [30], Ruan [40], Wolkowicz and Lu [48], Yang and Humphrey [49]).The half-saturation constant K i of the i th phytoplankton species for nutrient is given by f i (K i ) = µ i /2 and so represents the resource concentration supporting growth at half the maximal growth rate.Half-saturation constants and maximal growth rates can be measured experimentally (Hansen and Hubbell [17]).System (1) was considered by Roy [38] under the assumption that the response functions f i are Michaelis-Menten and the yield constants γ i equal 1.
2.1.Preliminary results.The first lemma is a statement that solutions of (1) are positive and bounded.These are minimal requirements for a reasonable model of the chemostat.(b) Given any δ > 0, for all solutions N (t) of ( 1) N (t) ≤ N 0 +δ, for all sufficiently large t.(c) If there exists a t 0 ≥ 0 such that N (t 0 ) ≤ N 0 , then N (t) < N 0 for all t > t 0 .
Proof of (b).Let δ > 0 be given.From the first equation of (1) we have Hence, approaches 0 as t tends to infinity, N (t) ≤ N 0 + δ for all sufficiently large t.
Proof of (c).Suppose there exists a first time t > t 0 such that N ( t) = N 0 and N (t) < N 0 for all t 0 ≤ t < t.Then dN dt ( t) ≥ 0. However, by (1), γi f i (N ( t)) < 0, a contradiction.The Fundamental Existence-Uniqueness Theorem (see, for example, Perko [33]) and Lemma 2.1(a) ensure that solutions of (1) exist uniquely for all time.
In the absence of allelopathic effects (γ = 0), system (1) reduces to the model studied by Hsu [21] in the case of two competing species.As such, the system exhibits the competitive exclusion principle, which Hardin [18] states as "complete competitors cannot coexist".That is, two species cannot coexist if they compete for a single resource available in growth-limiting amounts.
Of note, if γ = 0 then only one species survives, the one with the lower breakeven concentration.

Steady states: Existence and local stability. Steady states of model (1) are solutions of:
In what follows, λ i is the unique positive solution of the equation Three of the steady states are readily identified and are given by: E 0 = (N 0 , 0, 0), E λ1 = (λ 1 , P1 , 0), and We say that a steady state does not exist if any one of its components is negative.E 0 always exists, whereas for each i ∈ {1, 2} a necessary and sufficient condition on the parameters for feasibility of E λi is N 0 > λ i .Note that when N 0 = λ i , E 0 and E λi coalesce.
If any other equilibria exist, they must be interior equilibria (for which N , P 1 , and P 2 are all positive).By the third equation of (1), we must have N = λ 2 , so that an interior equilibrium Ê = (λ 2 , P1 , P2 ) of (1) corresponds to solutions Note that there is no interior equilibrium when λ 1 ≥ λ 2 .Using a method similar to that of Roy [38] we establish the existence of zero, one or two interior equilibria as a function of the value of the input nutrient concentration N 0 .The second equation of system (5) 2 in the first quadrant convex to the origin and satisfying lim P2→∞ F 2 (P 2 ) = 0 and lim The first equation of system (5) gives a straight line in the first quadrant with slope −γ 2 f 1 (λ 2 )/γ 1 m 2 and P 2 -intercept γ 2 (N 0 − λ 2 )D/m 2 .These curves may or may not intersect in the first quadrant.If N 0 is too small, the two curves do not intersect and model ( 1) does not have an interior equilibrium point.As N 0 is increased, there is a critical value N c 0 for which the straight line is tangent to the curve, and model ( 1) has precisely one equilibrium given by Ẽ = (λ 2 , P1 , P2 ), where The critical value is computed using system (5), and is given by As N 0 is increased beyond N c 0 , the straight line intersects the curve in two points and (1) has two interior equilibria Ê = (λ 2 , P 1 , P 2 ) and Ê = (λ 2 , P 1 , P 2 ), where 0 < P 1 < P1 , P2 < P 2 < 3 P2 /2, P1 < P 1 < 3 P1 , and 0 < P 2 < P2 (see Figure 1).Recall P1 = γ 1 (N 0 − λ 2 )D/3f 1 (λ 2 ), and P2 = 2γ 2 (N 0 − λ 2 )D/3m 2 are respectively the P 1 − and P 2 −coordinates of the unique interior equilibrium Ẽ when N 0 = N c 0 .We now investigate the local stability properties of (1) through an examination of the linearized system about each of the equilibria.The variational matrix of (1), denoted by V (N, P 1 , P 2 ), is given by At E 0 = (N 0 , 0, 0) we have  1) in the (P 1 , P 2 )-plane.For N 0 < N c 0 there is no interior equilibrium (the graph of the first equation of system (5) does not intersect the curve).For N 0 = N c 0 there is a unique interior equilibrium (the graph of the first equation of system ( 5) is tangent to the curve).For N 0 > N c 0 there are precisely two interior equilibria (the graph of the first equation of system (5) intersects the curve in two distinct points).
In the remaining cases, we use the following notations for the interior equilibria when N 0 > N c 0 : Ê = (λ ) has at least one eigenvalue with positive real, so that Ê is unstable (saddle point).
Roy [38] asserts that , when N 0 > N c 0 and the functions f i are Michaelis-Menten one of these two coexistence equilibria is locally stable whereas the other is unstable.Our results are consistent with those of Roy [38] in that Ê = (λ 2 , P 1 , P 2 ) is unstable and we are able to give conditions under which Ê = (λ 2 , P 1 , P 2 ) is locally asymptotically stable.
We summarize the results of this subsection in the following theorem.
Theorem 2.1.1. E 0 always exists.It is locally asymptotically stable for (1), if and only if N 0 < λ i for i = 1, 2.
(8) 2.3.Ecological interpretations of inequalities.This subsection gives ecological interpretations of inequalities resulting from the local stability results for our model (1) (see Theorem 2.1).By statement 1 of Theorem 2.1, the species-free steady state is locally asymptotically stable if and only if f i (N 0 ) < m i , for i = 1, 2. That is, the growth rate of species P i is strictly less than its removal rate m i , even when the culture vessel is held at the input nutrient concentration N 0 .Thus neither P 1 nor P 2 can survive at this input level of nutrient.
Statement 2 means that the growth rate of species P 2 is strictly less than its removal rate when the nutrient level in the culture vessel is held at λ 1 .Thus, species P 2 cannot compensate for the rate at which it is being removed from competition.It makes biological sense that only P 1 avoids extinction for initial conditions in a neighborhood of E λ1 .
The biological interpretation of the local stability conditions for E λ2 is symmetrical to that of E λ1 .

Global results.
In this subsection, we establish the global asymptotical stability of boundary equilibria due to the inadequacy of the resource supply.We first establish the competition-independent extinction of P i .The proof uses the following result due to Miller [28].
Proof of Lemma 2.3.Choose δ > 0 so that N 0 + δ < λ i .By Lemma 2.1(b), N (t) < N 0 + δ for all sufficiently large t.From the second and third equations of system (1), and by monotonicity properties of uptake functions f i , we have for all sufficiently large t.Hence by the definition of λ i (that is, f i (λ i ) = m i ), dP i (t) dt < 0 for all sufficiently large t.Also, P i (t) is bounded below.It follows from Lemma 2.2 (b) that P i (t) → 0 as t → 0. However, lim sup t→∞ f i (N (t)) < f i (N 0 + δ) < m i so that the only possibility is that P i (t) → 0 as t → ∞.
We are now in a position to prove that E 0 is a global attractor when it is the only steady state.
The next theorem gives conditions under which E λi is globally asymptotically stable.
Proof.Here we need to prove only (a) as the proof of (b) is similar.Take Let Ω(Q) denote the omega limit set of the orbit through Q.By the hypothesis and Lemma 2.3, any P = (N, P 1 , P 2 ) ∈ Ω(Q) satisfies P 2 = 0. On (N, P 1 , 0) ∈ R 3 + the system reduces to By an argument comparable to that given in Hsu [21], N (t) → λ 1 and P 1 (t) → P1 = γ1(N0−λ1)D m1 . Therefore, {E λ1 } ∈ Ω(Q).Since (10) has no periodic orbits and the boundary is acyclic, it follows from Lemma 4.3 in Thieme [45] that E λ1 is globally asymptotically stable for (1).

Transfer of local stability and one-parameter bifurcation.
In this subsection we describe the evolution of equilibria into the nonnegative cone of R 3 and the consequent transfer of stability as N 0 is increased.In particular, whenever a new boundary steady state coalesces with an existing one, a transcritical bifurcation occurs.We have two scenarios to consider.Scenario 1: λ 1 < λ 2 .When N 0 < λ 1 only the washout equilibrium E 0 exists and is globally asymptotically stable.As N 0 is increased so that N 0 = λ 1 , E 0 and E λ1 coalesce.As N 0 is increased still further, so that λ 1 < N 0 < λ 2 , E λ1 bifurcates into the nonnegative cone of R 3 through E 0 .The washout equilibrium loses a degree of stability, and E λ1 is globally asymptotically stable for (1).As N 0 is increased so that N 0 = λ 2 , E 0 and E λ2 coalesce.As N 0 is increased still further so that λ 1 < λ 2 < N 0 < N c 0 , E λ2 bifurcates into the nonnegative cone of R 3 through E 0 .Here, E 0 loses another degree of stability, E λ1 loses a degree of stability, and E λ2 is locally asymptotically stable for (1).When N 0 = N c 0 the interior equilibrium Ẽ becomes feasible in the positive cone of R 3 .For N 0 > N c 0 , Ẽ undergoes a saddlenode bifurcation to generate Ê and Ê .Ê is a saddle point for model (1) while Ê is locally asymptotically stable for (1) under conditions given in Theorem 2.1.In other words, we have an unstable interior equilibrium and an asymptotically stable interior equilibrium.As such, persistence is not possible.Unlike Roy [38], our goal is not to identify the basin of attraction of the locally stable interior equilibrium.We are more interested in global stability and uniform persistence.Scenario 2: λ 2 < λ 1 .When N 0 < λ 2 only the washout equilibrium E 0 exists and is globally asymptotically stable.As N 0 is increased so that N 0 = λ 2 , E 0 and E λ2 coalesce.As N 0 is increased still further, so that λ 2 < N 0 < λ 1 , E λ2 bifurcates into the nonnegative cone of R 3 through E 0 .The washout equilibrium loses a degree of stability, and E λ2 is globally asymptotically stable.As N 0 is increased so that N 0 = λ 1 , E 0 and E λ1 coalesce.As N 0 is increased still further so that λ 2 < λ 1 < N 0 , E λ1 bifurcates into the nonnegative cone of R 3 through E 0 .Here, E 0 loses another degree of stability, E λ2 loses a degree of stability, and E λ1 is locally asymptotically stable for (1), attracting solutions from the interior.
Hence, system (1) cannot be uniformly persistent.
2.6.Effect of nutrient recycling.The model considered so far lacks the effects of nutrient recycling.The regeneration of nutrient due to bacterial decomposition of the dead biomass must be considered.We modify model (1) to incorporate the effect of nutrient recycled from each phytoplankton cell on its death.The mortality of nontoxic phytoplankton P 1 is the sum of its intrinsic death and the death due to phytotoxic effects.The mortality of toxic phytoplankton P 2 is due to its intrinsic death and the death due autotoxic effects.For simplicity, we assume limiting nutrient recycling is instantaneous.We denote by η i (assumed to be less than 1 and constant over time) the nutrient contents of a single cell of phytoplankton P i .Under these considerations, model (1) can be extended as follows: We assume that all model parameters are nonnegative.Adapting the method used for the previous model, we can show that model (11) has the following steady states: , and E * λ2 = (λ 2 , 0, γ2(N0−λ2)D m2−γ2η2( 2+a2 ) with N 0 > λ 2 and m 2 > γ 2 η 2 ( 2 + a 2 ).As before, there are zero, one or two interior equilibria depending on whether or not the magnitude of the input nutrient concentration equals or exceeds a modified critical value N c * 0 given by, The qualitative dynamics remain unchanged from those exhibited by system (1) when the effect of nutrient recycling is incorporated into the model.Only the critical values of the parameters at which transitions take place are affected.For instance, the eigenvalues of the Jacobian matrix at Thus the local stability of E 0 remains the same under nutrient recycling.The characteristic polynomial of the Jacobian matrix at where . The existence conditions of E λ1 and the Routh-Hurwitz criterion imply that the roots of this polynomial have negative real parts.Hence the local stability criteria of E λ1 and E λ1 are similar.Observe that along a given P i − direction, each equilibrium of model ( 11) appears below the corresponding equilibrium of model (1).
Our analytical findings are consistent with those of Roy [38].In the next section, we further extend this work to include predation of phytoplankton species by a zooplankton population.
3. The competition-allelopathy-predation model.We reconsider model ( 1) and add a zooplankton species Z as a fourth interacting component.The two phytoplankton populations are the only prey of zooplankton.We neglect the potential negative effect of toxic phytoplankton P 2 on the growth of zooplankton Z.With the above assumptions our model may be formulated as follows: In these equations, N (t), P 1 (t), P 2 (t), D, γ i , and m i have the same meanings as in model ( 1) and the f i 's satisfy assumptions of model (1).Z(t) is the biomass of zooplankton species at time t.Here, the removal rate c of zooplankton Z is the sum of washout rate D and the specific death rate ξ, so that c = D + ξ.The response function g i (P i ) represents the per capita rate of conversion of phytoplankton P i to biomass of population Z per unit of population Z as a function of the biomass of phytoplankton P i .We assume that the rate of conversion of biomass of phytoplankton P i to biomass of zooplankton Z is proportional to the amount of phytoplankton consumed, so that the consumption rate of phytoplankton P i per unit of population Z is of the form 1 ηi g i (P i ), where η i is a growth yield constant (number of zooplankton per unit of phytoplankton).Since P 1 and P 2 are perfectly substitutable resources for Z (see for examples, Butler and Wolkowicz [8], León and Tumpson [25], Rapport [35], and Ballyk and Wolkowicz [4]), the per capita growth rate of zooplankton as a function of P 1 and P 2 takes the form G(P 1 , P 2 ) = g 1 (P 1 ) + g 2 (P 2 ) for all P 1 ≥ 0 and P 2 ≥ 0. Following Li and Kuang [26], Ruan [40], Wolkowicz and Lu [48] and others, we make the following assumptions concerning the response functions g i : where ω i denotes the maximal growth rate of zooplankton Z on phytoplankton P i .
It will also be convenient to express g i (P i ) as where h i (P i ) is some positive and differentiable function.Since g i is continuously differentiable it follows that lim and so we define The breakeven concentration for population Z on phytoplankton P i is obtained by setting dZ dt = 0 = g i (P i )−c and solving for P i .By the monotonicity assumptions, the solution Λ i is a uniquely defined positive extended real number as long as we assume Λ i = ∞ if g i (P i ) < c for all P i ≥ 0.
The half-saturation constant L i of zooplankton for the i th phytoplankton is given by g i (L i ) = ωi 2 and represents the phytoplankton biomass P i supporting growth at half the maximal growth rate.System (13) was considered by Butler and Wolkowicz [7] under the assumptions that specific death rates are insignificant compared to the washout rate D ( i = 0, i = 1, 2, ξ = 0), Z feeds only on one phytoplankton P i and there is no allelopathic effect.Holt et al. [19] studied model (13) under linearity of response functions f i and g i , and in the absence of allelopathic effects.Grover and Holt [13] relaxed the linearity assumptions on the responses functions in Holt et al. [19] and included the Holling types I and II response functions.In addition, system (13) was considered by Li and Kuang [26] under the assumption that one of the P i 's is absent and there is no allelopathic interaction.As such, model ( 13) is a significant generalization of those previously considered.
3.1.Preliminary results.We first establish that solutions of (13) are positive and bounded.Proof of (a).The proof is similar to that of Lemma 2.1 (a) except that we have to replace the A i 's and t i 's in the proofs of P i (t) > 0 for all t, by the following Âi 's and ti 's, respectively: For i ∈ 1, 2, let ti = min {t > 0 : P i (t) = 0} and define Suppose now that Z(0) > 0, then dZ dt = (g Proof of (b).To prove boundedness of solutions of model ( 13), define T (t) = N (t)+ . From ( 13) we have ) exp(−D 0 t).Since lim t→∞ (T (0) − dN0 D0 ) exp(−D 0 t) = 0, it follows that for all > 0, the solutions (N (t), P 1 (t), P 2 (t), Z(t)) of ( 13) satisfy for sufficiently large t.
The Fundamental Existence-Uniqueness Theorem (see, for example, Perko [33]) and Lemma 3.1(b) ensure that solutions of (13) exist uniquely for all time.
3.2.Steady states: Existence and local stability.It is straightforward to prove that all the boundary steady states of model (1) along with their existence conditions are transferred to model (13) except that we add a zero Z component.In the remainder of this work, we denote by λ i and Λ i the unique positive solutions of f i (N ) = m i , and g i (P i ) = c, respectively.The steady states of model ( 13) along with their existence conditions are listed below.
No equilibrium of the form (N, 0, 0, Z) (for which N and Z are both positive) exists.
We will just investigate the local stability properties of ( 13) through an examination of the linearized system about the equilibria E Λ1 , and E Λ2 .
We assume that N 0 − Λ2m2 γ2D > N 2 > λ 2 and f 2 (N 2 ) > m 2 , so that E Λ2 exists.We examine the local stability properties of E Λ2 .The variational matrix of (13) evaluated at E Λ2 , is given by The corresponding characteristic polynomial is given by where The monotonicity of f 2 (N ) and g 2 (P 2 ), the positivity of N 2 , Λ 2 , and Z 2 , together with the Routh-Hurwitz criterion, ensure that the roots of the cubic factor have negative real parts if and only if Ã > 0 and Ã B > C. Hence, E Λ2 is locally asymptotically stable for (13) if and only if f 1 (N 2 ) < Z2 η1 g 1 (0) + m 1 , Ã > 0 and Ã B > C. The local stability analysis of E Λ1 = (N 1 , Λ 1 , 0, Z 1 ) is symmetrical to the analysis for E Λ2 .It is straightforward to show that the coefficients of the cubic factor of the corresponding characteristic polynomial are given by We summarize the results of this subsection in the following theorem which extends Theorem 2.1.

Global results.
In this subsection we list conditions under which E 0 , E λi , and E Λi are globally asymptotically stable for system (13) with respect to all solutions initiating in the positive cone of R 4 .
recycling takes the form: In these equations θ 1 , θ 2 , and θ 3 (assumed to be less than 1) are the nutrient contents of a single cell of species P 1 , P 2 , and Z, respectively.The variables and remaining parameters of model (31) are the same as in model (13).It is straightforward to verify that instantaneous nutrient recycling does not impact the qualitative behavior of model (13).Only the values at which transitions take place are altered.The explanations are similar to those for model (11) and therefore are omitted.
We first take N 0 = 0.05, so that N 0 < λ 1 < λ 2 .By Theorem 2.2, the species-free steady state E 0 is globally asymptotically stable for (1): all solutions of (1) tend to E 0 regardless of initial condition.One such solution is depicted in Figure 2.
We then increase N 0 to 0.1, so that λ 1 < N 0 < λ 2 .By Theorem 2.3, E λ1 is globally asymptotically stable for (1): all solutions of (1) tend to E λ1 regardless of initial condition.One such solution is depicted in Figure 3.
We then increase N 0 to 0.12, so that λ 1 < N 0 < λ 2 , and we pick ω 1 = 1, so that ω 1 < c.By Theorem 3.3, E λ1 is globally asymptotically stable for (13): all solutions of (13) tend to E λ1 regardless of initial data.One such solution is depicted in Figure 5.
We finally increase N 0 to 0.70. Figure 7 shows that there is a unique interior equilibrium for model (13).
. Graphs of N (t), P 1 (t), P 2 (t), and Z(t).The input nutrient concentration N 0 = 0.70 satisfies the assumptions on the existence of an interior equilibrium of model ( 13).
5. Concluding remarks.In this paper, we have extended the work of Roy [38] to general uptake functions, demonstrating again that allelopathy due to one toxinproducing species can lead to stable coexistence (while competitive exclusion of the weaker competitor would otherwise occur).Roy used the phytotoxic coefficient γ for his bifurcation analysis.We have used the input nutrient concentration N 0 as our bifurcation parameter because it can be controlled by the experimenter.As N 0 is increased, the model exhibits a biologically revelant evolution of equilibria into the nonnegative cone of R 3 and the consequent transfer of stability.We have studied local and global stability of each of the boundary equilibria E 0 , E λ1 , and E λ2 of model (1).When N 0 reaches at least the threshold N c 0 , we are able to study local stability of interior equilibria.Model (1) cannot exhibit uniform persistence because when N 0 > N c 0 , we have two distinct interior equilibria: one is a saddle and the other one is locally asymptotically stable.This result remains true when the toxic phytoplankton P 2 is the stronger competitor.We recall that when λ 2 < λ 1 (P 2 is stronger in competition than P 1 ) there is no interior equilibrium.
The incorporation of nutrient recycling into the model (1) modifies only the parameter bounds at which transitions occur, leaving the qualitative dynamics of the model (1) unchanged.
Predation can be responsible for diversity in ecosystems.Predation may promote, hinder or have no effect on interspecific competitive interactions (Chesson et al. [10]).We introduced predation of phytoplankton species by zooplankton in the model ( 1) and analyzed the modified system.Model (13) illustrates how predation increases the number of boundary equilibria of model (1) and brings diversity to interspecific competitive interactions.When the input nutrient concentration N 0 reaches a critical value, model (1) has one or two interior equilibria while our simulations indicate that model (13) has a unique interior equilibrium (see Figure 7).The system (31) extends model (13) to incorporate instantaneous nutrient recycling.
In a recent work (Kengwoung-Keumo [24]), we showed that model (13), in the absence of phytotoxic effects and nutrient recycling, exhibits uniform persistence.For future work, we could investigate if the persistence result still holds in the presence of phytotoxic effects and nutrient recycling.We could also incorporate excretion and/or the phytotoxic effects of toxins on the growth of zooplankton and study the impact(s) on the model (13).

Figure 1 .
Figure 1.Existence and non-existence of interior equilibria of model (1) in the (P 1 , P 2 )-plane.For N 0 < N c 0 there is no interior equilibrium (the graph of the first equation of system(5) does not intersect the curve).For N 0 = N c 0 there is a unique interior equilibrium (the graph of the first equation of system (5) is tangent to the curve).For N 0 > N c 0 there are precisely two interior equilibria (the graph of the first equation of system (5) intersects the curve in two distinct points).