A Nutrient-Prey-Predator Model: Stability and Bifurcations
In this paper we consider a model of a nutrient-prey-predator system in a chemostat with general functional responses, using the input concentration of nutrient as the bifurcation parameter. We study the changes in the existence of isolated equilibria and in their stability, as well as the global dynamics, as the nutrient concentration varies. The bifurcations of the system are analytically verified and we identify conditions under which an equilibrium undergoes a Hopf bifurcation and a limit cycle appears. Numerical simulations for specific functional responses illustrate the general results.
We consider a mathematical model of two-species predator-prey interaction in the chemostat under nutrient limitation. With the exception of one nutrient, all nutrients that the prey species requires are supplied to the growth vessel from the feed vessel in ample supply. The predator species grows exclusively on the prey. With the concentration of the limiting nutrient, the concentration of prey (say, phytoplankton), and the concentration of predator (say, zooplankton), we consider the following model:
for initial conditions , , and .
The concentration of the growth-limiting nutrient in the feed vessel is denoted , and will be the bifurcation parameter in our analysis. is the input rate from the feed vessel to the growth vessel as well as the washout rate from the growth vessel to the receptacle, so that constant volume is maintained. The parameters and are the removal rates of and , respectively, from the growth vessel, incorporating the washout rate and the intrinsic death rates of and . Our analysis does not necessarily require that and are positive; however, and should be positive. The yield coefficient gives the amount of prey biomass produced per unit of nutrient consumed, while gives the amount of predator biomass produced per unit of prey biomass consumed.
The function represents the per capita consumption rate of nutrient by the prey populations as a function of the concentration of available nutrient; similarly, the function represents the per capita consumption rate of the prey by the predator as a function of available prey. These functions are assumed to satisfy , , 2, for all , and for all . We further assume that and are bounded. To avoid the case of washout due to an inadequate resource, we assume that , and to avoid the case of an inadequate prey, we assume that . Define and to be the unique numbers satisfying
The number is the break-even concentration of nutrient, at which the growth and removal of phytoplankton balance; is similarly interpreted. The number plays a central role in our investigation, so we assume that and are also defined. From the perspective of , due to the boundedness assumptions on and , and are perturbations of .
From a functional point of view, and are right inverses of and , respectively. Accordingly, on their respective domains and are as differentiable as and . We note for later use
Kuang and Li  studied this system with general functional responses and distinct values of , , and . However, they fixed the input nutrient concentration, whereas we have this as a parameter. With the hypothesis that , they provide numerical criteria for the stability of a coexistence equilibrium, and prove that a cycle exists when the equilibrium is unstable [7, Theorem 3.2]. When the hypothesis does not hold, they provide numerical evidence that stability of the coexistence equilibrium breaks down and a cycle appears. A similar model was studied in  with functional response of Holling type I and of Holling type II, demonstrating the existence of a Hopf bifurcation in response to varying nutrient concentration. These results inspire our work. Our goal is to prove analytically that the system undergoes a Hopf bifurcation without restricting the forms of the uptake functions or the values of the removal rates and .
This paper is organized as follows. In section 2, conditions for existence and local stability of predator-free equilibria are obtained for general functional responses and . In section 3, we study stability of a coexistence equilibrium which appears as the parameter increases. We quote a version of the Hopf bifurcation theorem, stating the result in a limited form most useful for us. Accordingly, application of this theorem requires us to develop specific information about the behavior of eigenvalues of the linearizations as the paramter increases. In section 4, we use a smooth change of variables to reach a situation where the Hopf bifurcation theorem in three dimensions applies, enabling us to conclude the existence of cycles. In section 5, our results are illustrated using simulations arising by choosing rate functions of Holling type II and Holling type III forms. In section 6, we explain in detail the approximations and estimates that support our work in the latter part of section 3.
2 Steady States and Their Stability
To begin, we establish that the solutions of (1) are nonnegative and bounded. These are minimum requirements for a reasonable model of the chemostat. We then develop conditions for the existence and stability of equilibria. We conclude the section by proving uniform persistence in the sense of  when is sufficiently large and the initial values and are positive.
All solutions , , and of (1) are nonnegative and bounded.
The plane is invariant for system (1). Therefore, by existence and uniqueness, if then for all . Similarly, since , the plane is invariant, so implies for all .
Suppose . If there exists a with , then there is a least such number, say . Then since . Consequently, there is such that , which is a contradiction to the choice of .
For the boundedness of solutions, set . From system (1), it follows that
where . Then
Since , , and are nonnegative, the boundedness of implies the boundedness of , , and . ∎
There are at most three biologically relevant equilibria of system (1) depending on the value of . The equilibria and the conditions of their existence are summarized in the following theorem.
Thus, for only the equilibrium exists; for , there are two equilibria ; and, when , there are three equilibria , , and .
From the -equation, either or (so that ). If , the -equation yields
which implies either or (so that ).
If and , then the -equation gives , so that . Thus, the washout equilibrium is given by and exists for . Note that as increases, moves along the -axis in . This is the proof of the first part of the theorem.
If and , then in the -equation gives
Note that for all , and so the single-species equilibrium exists in the positive cone for all . This is the proof of the second part of the theorem.
We now determine the critical value of at which first appears in the positive cone. Equation (7) implies
Since is strictly increasing, if and only if , and for . Substituting in (6), we obtain the equation
Thus, is positive when , and the coexistence equilibrium exists in the positive cone for all . This proves part three of the theorem. ∎
In Theorems 2.3, 2.4, and 2.5 we investigate the stability of the equilibria of system (1) by finding the eigenvalues of the associated Jacobian matrices. The Jacobian matrix of the system (1) takes the form
We first summarize the stability of in the following theorem. Here, the breakeven concentration of nutrient given in (2) plays a critical role.
The equilibrium point is locally asymptotically stable if and unstable if . When , is globally asymptotically stable with respect to solutions initiating in . That is, the plane is , the stable manifold of .
The Jacobian at is
so that the eigenvalues of are , and . The stability of now follows from (2) and the fact that is strictly increasing: when and when . To see that when , consider the Lyapunov function
Clearly, and if . The time derivative of at a point on a trajectory of system (1) is
for . Thus is globally asymptotically stable in the plane . ∎
For , , so that and coalesce (see equation (4)). When , enters the positive cone. We summarize the stability of in the following theorem. Note that the critical value of given in (5) now plays a central role.
The equilibrium point is locally stable if and unstable if . When , is globally asymptotically stable with respect to solutions initiating in . That is, the plane is , the stable manifold of .
The Jacobian matrix at is
The determinant of the upper lefthand -by- submatrix is positive and its trace is negative, so its eigenvalues have negative real parts. The third eigenvalue is . If , so that , then , and is locally stable. Similarly, if , so that , then and is unstable with one dimension of instability.
To see that when , consider the Lyapunov function introduced by Hsu in 
Notice that ,
precisely when . Moreover,
Therefore, is the only critical point of and it is a local minimum, so that for all .
If , then and , so
Thus, when .
If , then . When ,
while for , then . In either case, , and so when .
Said another way, as increases through , enters the positive cone by passing through .
With , the Jacobian matrix at takes the form
The eigenvalues of satisfy
The coexistence equilibrium point
is asymptotically stable if and only if and .
Since is positive, this follows from the Routh-Hurwitz criterion. ∎
Assume that . Then system (1) is uniformly persistent with respect to all solutions satisfying and .
We first show that . If and , then . But this is impossible, for then it follows from the -equation that as and this, in turn, contradicts the fact that is bounded.
Now, suppose while . Then there exists a sequence of local minima of satisfying as . Thus,
, since is a local minimum, and
as , since .
From the -equation we have
Rearranging and using the facts that as , is continuous and , we get
a contradiction. Hence, .
Choose . Then is a nonempty, compact invariant set with respect to system (1). We claim and are not in . Suppose . Since , is an unstable hyperbolic equilibrium point. By theorem 2.3, is globally asymptotically stable with respect to solutions initiating in the plane . Since , . By the Butler-McGehee lemma (see lemma A1, ), there exists , so that . For such an initial condition, the governing system is , and . But then becomes unbounded as . This is a contradiction to the compactness of , and so .
Suppose is in . Since , is an unstable hyperbolic equilibrium point. By theorem 2.4 is globally asymptotically stable with respect to solutions initiating in the plane . Since , . By the Butler-McGehee lemma, there exists , so that . If , then , a contradiction. Thus , and this implies is unbounded as , contradicting the compactness of . Therefore, .
3 Stability of the coexistence equilibrium
In this section we study the coexistence equilibrium point
first when and varies and then after relaxing the assumption . In particular, we study the eigenvalues of the Jacobians at the coexistence equilibria as increases. To prepare for the study of the evolution of the equilibrium , we observe the following consequence of the implicit function theorem [5, p.122].
For , , and there is a local parameterization of the locus of coexistence equilibria
defined on an interval containing and a disc around that is smooth to the smaller of the degree of smoothness of and the degree of smoothness of .