A NUTRIENT-PREY-PREDATOR MODEL: STABILITY AND BIFURCATIONS

. We model a nutrient-prey-predator system in a chemostat with general functional responses, using the input concentration of nutrient as the bifurcation parameter. We study changes in the existence and the stability of isolated equilibria, as well as changes in the global dynamics, as the nutrient concentration varies. The bifurcations of the system are analytically veriﬁed and we identify conditions under which an equilibrium undergoes a Hopf bifur- cation and a limit cycle appears. Numerical simulations for speciﬁc functional responses illustrate the general results.


Introduction.
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 N the concentration of the limiting nutrient, P the concentration of prey (say, phytoplankton), and Z the concentration of predator (say, zooplankton), we consider the following model: dN/dt = (µ − N )D − P f 1 (N ) dP/dt = γ 1 P f 1 (N ) − D 1 P − Zf 2 (P ) (1) dZ/dt = γ 2 Zf 2 (P ) − D 2 Z for initial conditions N (0) = N 0 > 0, P (0) = P 0 ≥ 0, and Z(0) = Z 0 ≥ 0. The concentration of the growth-limiting nutrient in the feed vessel is denoted µ, and will be the bifurcation parameter in our analysis. D 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 D 1 = D + 1 and D 2 = D + 2 are the removal rates of P and Z, respectively, from the growth vessel, incorporating the washout rate D and the intrinsic death rates of prey biomass produced per unit of nutrient consumed, while γ 2 gives the amount of predator biomass produced per unit of prey biomass consumed.
The function f 1 (N ) represents the per capita consumption rate of nutrient by the prey populations as a function of the concentration of available nutrient; similarly, the function f 2 (P ) represents the per capita consumption rate of the prey by the predator as a function of available prey. These functions are assumed to satisfy f i (0) = 0, i = 1, 2, f 1 (N ) > 0 for all N ≥ 0, and f 2 (P ) > 0 for all P ≥ 0. We further assume that f 1 (N ) and f 2 (P ) are bounded. To avoid the case of washout due to an inadequate resource, we assume that lim N →∞ f 1 (N ) > D 1 /γ 1 , and to avoid the case of an inadequate prey, we assume that lim P →∞ f 2 (P ) > D 2 /γ 2 . Define λ P (D 1 ) and λ Z (D 2 ) to be the unique numbers satisfying f 1 λ P (D 1 ) = D 1 /γ 1 and f 2 λ Z (D 2 ) = D 2 /γ 2 . (2) The number λ P (D 1 ) is the break-even concentration of nutrient, at which the growth and removal of phytoplankton balance; λ Z (D 2 ) is similarly interpreted. The number D plays a central role in our investigation, so we assume that λ P (D) and λ Z (D) are also defined. From the perspective of D, due to the boundedness assumptions on f 1 and f 2 , D 1 and D 2 are perturbations of D.
Lemma 1.1. From a functional point of view, λ P and λ Z are right inverses of γ 1 f 1 and γ 2 f 2 , respectively. Accordingly, on their respective domains λ P and λ Z are as differentiable as f 1 and f 2 . We note for later use Variations on this system have been investigated by many authors over a long period of time. For example, Gard [4] studied a version with a general function representing the rate of nutrient removal and developed a numerical condition for the presence of a cycle. Li and Kuang [9] studied this system with general functional responses and distinct values of D, D 1 , and D 2 . With the hypothesis that D=D 1 =D 2 , they provide numerical criteria for the stability of a coexistence equilibrium, and prove that a cycle exists when the equilibrium is unstable [9,Theorem 3.2]. When the hypothesis D=D 1 =D 2 does not hold, they provide numerical evidence that stability of the coexistence equilibrium breaks down and a cycle appears. We note that the input nutrient concentration is fixed in both [4] and [9], whereas we have this as a parameter µ. A similar model was studied in [12] with functional response f 1 (N ) of Holling type I and f 2 (P ) 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 system (1) undergoes a Hopf bifurcation without restricting a priori either the forms of the uptake functions or the values of the removal rates D 1 and D 2 .
This paper is organized as follows. In Section 2, conditions for the existence of equilibria are obtained for general functional responses f 1 (N ) and f 2 (P ). We discuss the stability of predator-free equilibria and begin the discussion of the stability of a coexistence equilibrium which appears as the parameter µ increases. In Section 3, we study stability of the coexistence equilibrium and prove the existence of a Hopf bifurcation to an orbitally asymptotically stable cycle. This material is quite technical and comprises our main contribution to understanding the system. The analysis depends on Theorem 3.4, a paraphrase of the Hopf bifurcation theorem as developed in [5], stated in a form most useful for us. Accordingly, application of this theorem requires us to develop specific information about the behavior of the eigenvalues of the linearizations at the coexistence equilibrium as the parameter µ increases. We first make the assumption D=D 1 =D 2 , which allows us to factor the characteristic polynomial of the linearization at the coexistence equilibrium. We give conditions on the uptake functions that enable us to complete the verification of the hypotheses of the Hopf bifurcation theorem for this case in Theorem 3.7. Under additional conditions on the uptake functions, we can prove that there is a unique value of the bifurcation parameter at which a Hopf bifurcation to an orbitally asymptotically stable cycle occurs. Relaxing the condition D=D 1 =D 2 , we use the inverse function theorem to factor the characteristic polynomials of the linearizations at coexistence equilibria E 2 (µ, D 1 , D 2 ) for (D 1 , D 2 ) near (D, D). By continuity of the factorization, the eigenvalue conditions of the Hopf theorem are fulfilled. In Lemma 3.10 we state a uniform approximation result that implies the transversality and stability conditions of the Hopf theorem are also satisfied. Although the proof of Lemma 3.10 is straightforward, it is technical, so we postpone explanation of the proof to Section 5. In Theorem 3.11 we complete the verification of the hypotheses of Theorem 3.4 for the general case, establishing that the behavior in the case D=D 1 =D 2 propagates to the general situation. In Section 4 our results are illustrated using simulations arising by choosing rate functions of Holling type II and Holling type III forms.
2. Steady states and their stability. To begin, we establish that the solutions of system (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 [2] when µ is sufficiently large and the initial values P 0 and Z 0 are positive. Lemma 2.1. All solutions N (t), P (t), and Z(t) of system (1) are nonnegative and bounded.
Suppose N 0 > 0. If there exists a t > 0 with N (t) = 0, then there is a least such number, say t 0 . Then N (t 0 ) = µD > 0 since f 1 (0) = 0. Consequently, there is t < t 0 such that N (t) < 0, which is a contradiction to the choice of t 0 .
For the boundedness of solutions, set Since N (t), P (t), and Z(t) are nonnegative, the boundedness of U (t) implies the boundedness of N (t), P (t), and Z(t).
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. Theorem 2.2. Let λ P (D 1 ) and λ Z (D 2 ) be as in (2). The equilibria of the system (1) satisfy the following conditions: 1. The washout equilibrium E 0 = (µ, 0, 0) exists for all µ > 0.
If Z = 0 and P = 0, then the N -equation gives 0 = (µ − N )D, so that N = µ. Thus, the washout equilibrium is given by E 0 = (µ, 0, 0) and exists for µ > 0. Note that as µ increases, E 0 moves along the N -axis in R 3 ≥0 . This is the proof of the first part of the theorem.
If Z = 0, then P = λ Z (D 2 ) in the N -and P -equations, and equations (6) and (7) define N and Z as implicit functions of µ, D 1 , and D 2 .
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 E 0 in the following theorem. Here, the breakeven concentration of nutrient given in (2) plays a critical role. Theorem 2.3. The equilibrium point E 0 is locally asymptotically stable if µ < λ P (D 1 ) and unstable if µ > λ P (D 1 ). When µ > λ P (D 1 ), E 0 is globally asymptotically stable with respect to solutions initiating in {(N, P, Z) ∈ R 3 + | P = 0}. That is, the plane P = 0 is m + (E 0 ), the stable manifold of E 0 .
With f 2 λ Z (D 2 ) = D 2 /γ 2 , refer to (9) to see that the Jacobian matrix at E 2 (µ, D 1 , D 2 ) takes the form The eigenvalues of J E 2 (µ, D 1 , D 2 ) satisfy and Theorem 2.5. The coexistence equilibrium point is asymptotically stable if and only if a 1 > 0 and a 1 a 2 > a 3 .
Proof. Since a 3 is positive, this follows from the Routh-Hurwitz criterion.
We conclude this section with a significant strengthening of Lemma 2.1. We use the concept of uniform persistence introduced in [2]. Theorem 2.6. Assume that µ > µ c1 (D 1 , D 2 ). Then system (1) is uniformly persistent with respect to all solutions satisfying P 0 > 0 and Z 0 > 0.
Proof. Recall from Lemma 2.1 that all solutions of system (1) for which P 0 > 0 and Z 0 > 0 are positive and bounded. We first show that lim inf t→∞ N (t) > 0.
If lim inf t→∞ N (t) = 0 and lim sup t→∞ N (t) = 0, then lim t→∞ N (t) = 0. But this is impossible, for then it follows from the N -equation that N (t) → µD > 0 as t → ∞ and this, in turn, contradicts the fact that N (t) is bounded.
Rearranging and using the facts that N (τ n ) → 0 as n → ∞, f 1 is continuous and Then ω X(0) is a nonempty, compact invariant set with respect to system (1). We claim E 0 = (µ, 0, 0) and ). For such an initial condition, the governing system is This is a contradiction to the compactness of ω X(0) , and so is an unstable hyperbolic equilibrium point. By Theorem 2.4 E 1 (µ, D 1 ) is globally asymptotically stable with respect to solutions initiating in the plane , neither of which can be true. Thus, lim inf t→∞ P (t) > 0 and lim inf t→∞ Z(t) > 0, and it follows from the main result of [2] that system (1) is uniformly persistent.
3. Stability of the coexistence equilibrium. In this section we study the coexistence equilibrium point as µ varies, first when D=D 1 =D 2 and then after relaxing this 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 E 2 (µ, D 1 , D 2 ), we observe the following consequence of the implicit function theorem [7, p.122].
there is a local parameterization of the locus of coexistence equilibria defined on an interval containing µ and a disc around (D 1 , D 2 ) that is smooth to the smaller of the degree of smoothness of f 1 (N ) and the degree of smoothness of f 2 (P ).
Proof. Set P = λ Z (D 2 ) in the N -and P -equations of system (1). Let Then we want to parametrize the set G −1 (0, 0). The derivative of G is represented by the matrix Observe that, for fixed µ 0 > µ c1 (D 1 , D 2 ), the first two columns of DG at the point Since f 1 (N ) > 0, the first two columns are linearly independent and so the implicit function theorem applies. There exists a ball B 1 around (µ 0 , D 1 , D 2 ) and a ball B 2 Moreover, the functions N (µ, D 1 , D 2 ) and Z(µ, D 1 , D 2 ) have the same degree of differentiability as does G, which is the minimum of the degrees of differentiability of f 1 (N ) and f 2 (P ). To explain how the differentiability of f 2 enters, note that the computation of ∂N/∂D 2 and ∂Z/∂D (3). This all gives us a smooth parameterization of the equilibrium locus, as desired.
Remark 3.2. From the point of view of calculus, each of the independent variables µ, D 1 , D 2 is on an equal footing with the others. However, viewed through the lens of the structure of system (1), we can describe the variable µ as a control parameter and the variables D 1 and D 2 as experimental parameters fixed at some earlier time. This distinction informs our analysis where we first consider the model when D=D 1 =D 2 as µ varies, and subsequently vary D 1 and D 2 in a neighborhood of D.
Proof. Since D 1 and D 2 are fixed throughout this proof, we drop these symbols and write simply N (µ) and Z(µ). To show N (µ) > 0, recall equation (6) from Theorem 2.2: Differentiating (6) with respect to µ and rearranging, we get (7), obtaining Differentiating with respect to µ and rearranging, we get Since N (µ) and f 1 (N ) are positive, it follows that Z (µ) is positive. This proves part one. To prove part two, note that equation (6) implies Since f 1 (N ) is bounded, lim µ→∞ N (µ) = ∞. From (7) we have Before we turn to a study of the stability properties of the coexistence equilibrium, we include a partial paraphrase of the C L Hopf bifurcation theorem as stated in [5, p.16]. Since our goal is to make an application of this result to the coexistence equilibrium E 2 , and because the verification of the hypotheses is lengthy, we refer to our paraphrase to keep track of progress.
Theorem 3.4. Consider a system dX/dt = F (X, µ) with X ∈ R n and µ a real parameter. If, 1. for µ in an open interval containing µ c (characterized in 3 below), F (0, µ) = 0 and 0 ∈ R n is an isolated equilibrium point of dX/dt=F (X, µ); 2. all partial derivatives of the components F of the vector F of orders ≤ L+2, (L ≥ 2) exist and are continuous in X and µ in a neighborhood of (0, µ c ) in where α(µ c ) = 0 and α (µ c ) = 0; 4. the remaining n−2 eigenvalues of J(0, µ c ) have strictly negative real parts, then the system dX/dt = F (X, µ) has a family of periodic solutions.
Moreover, if α (µ c ) > 0, then the family of periodic solutions arises for µ > µ c and for µ sufficiently close to µ c these periodic solutions are orbitally asymptotically stable.
Remark 3.5. For the purposes of the proof in [5] the authors assume the critical value of the bifurcation parameter is µ c = 0 and that the isolated critical point is at the origin. A straightforward linear adjustment of the bifurcation variable moves the critical value of the bifurcation parameter to 0, but a more serious issue with hypothesis 1 is that, for us, the coordinates of the equilibrium (1) are changing with the parameter µ. Figure 1 shows an example. In [5, p.14] a linear change of variables is recommended to translate the equilibrium smoothly varying with the bifurcation parameter back to the origin. The reader interested in the details of this adjustment is referred to the preprint [1], where we use the inverse function theorem to overcome this difficulty. After the change of variables, hypothesis 1 is Figure 1. A curve of coexistence equilibria satisfied and the Jacobian at the origin is simply a conjugate of the Jacobian at Hypothesis 2 is fulfilled by imposing more differentiability conditions on the functions f 1 (N ) and f 2 (P ) at an appropriate point in the exposition. Verification of hypotheses 3 and 4 in the statement of Theorem 3.4 is the most involved part of our process and occupies the rest of this section.
Hypotheses 3 and 4 concern the eigenvalues of the Jacobian J(E 2 ), but the characteristic polynomial of J(E 2 ) given in Equation (10) is not immediately amenable to factorization, a standard approach to eigenvalue study. Motivated by the choice of the auxiliary variable U (t) defined in Lemma 2.1, we conjugate the Jacobian J(E 2 ) in (10) by to obtain a more convenient form, as follows. Using where First we make the assumption that D=D 1 =D 2 , so that 1 = 2 =0; that is, we assume the death rates of P and Z are negligible with respect to washout rate D. The results of [9] suggest this is a useful initial assumption. Since we regard D as fixed for the discussion, we will abbreviate N (µ, D, D) by N (µ) and Z(µ, D, D) by Z(µ).
Similarly, we will abbreviate E 2 (µ, D, D) by E 2 (µ). From (14) we can explicitly compute the characteristic polynomial and eigenvalues of J E 2 (µ) , since conjugation does not change them. The characteristic polynomial of J E 2 (µ) is, therefore, where and The next result amplifies Theorem 2.5 in the case D=D 1 =D 2 .
In [9,Theorem 3.2] it is shown that, under the assumption D=D 1 =D 2 , cycles exist for certain values of the parameters. To prove this, Li and Kuang observe that there is a limiting plane for solutions of the system, given by in terms of our variables. Using their form of (20), they eliminate explicit consideration of the nutrient equation and reduce the number of variables to two. The Poincaré-Bendixson theorem is then applied to the resulting two-dimensional limit system, providing a cycle for the (N, P, Z)-system in the plane (20). Now we amplify their result by proving cycles arise from Hopf bifurcations. From the factorization of the characteristic polynomial of J E 2 (µ) given in (15) it is immediate that the Jacobian at E 2 (µ) has one negative eigenvalue. Thus, hypothesis 4 of Theorem 3.4 for E 2 (µ) is satisfied in the case D=D 1 =D 2 . Now we verify hypothesis 3 for this situation.
1. If f 1 (N ) is continuous and D/ γ 2 λ Z (D) > f 2 λ Z (D) , then there exists a value µ c2 > µ c1 (D, D) for which A(µ c2 ) = 0. Consequently, when µ = µ c2 , the Jacobian has a conjugate pair of imaginary eigenvalues. 2. If, in addition, f 1 is twice differentiable with respect to N and f (2) 1 N (µ c2 ) < 0, then A (µ c2 ) > 0. Combining this with part 1, we have that hypothesis 3 of Theorem 3.4 is satisfied for E 2 (µ). 3. If, in addition, f 1 and f 2 are four times continuously differentiable, then all parts of Theorem 3.4 are satisfied. Moreover, the last sentence in the statement also applies, and we conclude there is a Hopf bifurcation at µ c2 to an orbitally asymptotically stable cycle. 4. If, in addition to the hypotheses in parts 1, 2, and 3, f 1 (N ) < 0 for all N , then there is a unique Hopf bifurcation at µ c2 to a stable cycle. 2 (P ) < 0, i.e., that f 2 is is concave down, then the slope of the secant that passes through the points (0, 0) and λ Z (D), (D/γ 2 ) is greater than the slope of the tangent line to the graph of f 2 at λ Z (D); that is, D/ γ 2 λ Z (D) > f 2 λ Z (D) . This will be the case, for example, when f 2 (P ) is Holling Type II. However, the one-point condition D/ γ 2 λ Z (D) − f 2 λ Z (D) > 0 may hold even if f 2 is not concave down. We will give such an example in Section 4.
Concerning part 2 of the theorem, our argument for the existence of a parameter value µ c2 such that A(µ c2 ) = 0 is based on the intermediate value theorem, so we cannot a priori exclude the possibility that there is more than one parameter value meeting this condition. If f > 0 at some zero of A(µ), the sign of A at that zero cannot be determined. Consequently, it may be possible that more than one Hopf bifurcation occurs and that some may result in bifurcations to unstable cycles.
Concerning part 4, the assumption that f 1 (N ) < 0 for all N is fufilled if f 1 is a function of Holling type II and the possibility of a multiplicity of zeroes of A(µ) is eliminated. As such, the assumption is natural in some respects, although one can conceive of examples not fulfilling the assumption.

Proof. Consider the expression
given in (19). Since , this implies that there exists an N * such that, if N > N * , then In addition, N (µ) is increasing without bound by Theorem 3.3, so there is an Since A µ c1 (D, D) < 0 and A(µ * ) > 0, there is a number µ c2 > µ c1 (D, D) such that A(µ c2 ) = 0. Note that when µ = µ c2 the discriminant of the quadratic factor of the characteristic polynomial in (15) is so its roots are indeed purely imaginary. This proves part one.
For part two, by continuity of the discriminant as a function of µ, there is a neighborhood of µ c2 on which the discriminant is negative. Continuing, differentiate A(µ) with respect to µ to obtain 1 (N (µ c2 )) < 0. Thus, A (µ c2 ) > 0. Combining parts one and two means that hypothesis 3 of Theorem 3.4 holds at µ c2 for the case D=D 1 =D 2 . This completes the proof of part 2.
With f 1 and f 2 four times continuously differentiable, all parts of Theorem 3.4 are satisfied. In particular, the last sentence of Theorem 3.4 also applies, and we conclude that at µ c2 there is a Hopf bifurcation to an orbitally asymptotically stable cycle.
For part four, if f Let us now discuss weakening the condition D=D 1 =D 2 . Intuitively, for (D 1 , D 2 ) sufficiently close to (D, D) the eigenvalues of J E 2 (µ, D 1 , D 2 ) should exhibit behavior similar to those of J E 2 (µ, D, D) . It is a well-known consequence of the inverse function theorem that, for a degree n polynomial with n distinct real roots, the roots and, hence, the linear factors, are smooth functions of the coefficients. In our situation we may specialize this fact in the form of the following lemma. Lemma 3.9. Let P 1 (x) − = {(α − x) | α < 0} be the space of polynomials of degree 1 in x, with leading coefficient −1 and negative constant term, let P 2 (x) − = {β −γx+x 2 | γ 2 −4β < 0} be the space of monic quadratic polynomials in x with real coefficients and having a complex conjugate pair of roots, and let P 3 (x) − = {p 0 + p 1 x + p 2 x 2 − x 3 } be the space of cubic polynomials in x with leading coefficient −1 and real coefficients.
Then the multiplication map M : P − 1 ×P − 2 → P − 3 is locally a diffeomorphism. To explain how Lemma 3.9 comes into play, consider the map χ : M 3,3 (R) → P − 3 which takes as input a real-valued 3 by 3 matrix R and produces its characteristic polynomial χ(R) = det(R − xI). The map χ has a coordinate expression by taking the coefficients in degrees 0, 1, and 2. These coefficients are polynomials in the matrix entries, so χ is smooth. Now look at Suppose we are in the situation of Theorem 3.7, where it is easy to factor the characteristic polynomial of J E 2 (µ, D, D) , and we have seen the Jacobian has a negative real eigenvalue and a complex conjugate pair of eigenvalues as µ varies near a potential bifurcation value µ c2 . The observed factorization of the characteristic polynomial of J E 2 (µ, D, D) explicitly inverts the polynomial multiplication M at particular points, and Lemma 3.9 implies the characteristic polynomial of J E 2 (µ, D 1 , D 2 ) has a similar factorization smoothly depending on the entries of J E 2 (µ, D 1 , D 2 ) , for (D 1 , D 2 ) sufficiently close (D, D).
More precisely, we have proved in Theorem 3.7 that there is an interval of parameters µ in which the characteristic polynomials of J E 2 (µ, D, D) have a complex conjugate pair of roots in addition to the eigenvalue −D<0, so they are in the image of the multiplication map M : . For (D 1 , D 2 ) close to (D, D) the entries in J E 2 (µ, D 1 , D 2 ) are close to the entries in J E 2 (µ, D, D) , so the characteristic polynomial of J E 2 (µ, D 1 , D 2 ) is close to the characteristic polynomial of J E 2 (µ, D, D) . To see this explicitly, refer to the formulas (11), (12), and (13); to account for the normalization of the characteristic polynomial to leading coefficient −1, multiply each expression by −1. Therefore, in view of Lemma 3.9, the characteristic polynomial of J E 2 (µ, D 1 , D 2 ) has a decomposition of the same form as that of the characteristic polynomial of J E 2 (µ, D, D) . Written formally, the decomposition is We remind the reader that we think of µ as a control parameter, adjustable by the experimenter, and D 1 and D 2 as experimental parameters, set at the beginning of an experiment. To bring out this distinction, we have written the components of the formal factorization of the characteristic polynomial of J E 2 (µ, D 1 , D 2 ) as α(D 1 , D 2 )(µ), β(D 1 , D 2 )(µ), and γ(D 1 , D 2 )(µ).
The map M is infinitely differentiable, so the local inverse M −1 is also. In particular, for a fixed value µ, the coefficients of the decomposition are smooth functions of (D 1 , D 2 ) defined in a neighborhood of (D, D).
We are primarily interested in the behavior of the eigenvalues of the Jacobian J E 2 (µ, D 1 , D 2 ) in the neighborhood of a value of µ where the Jacobian J E 2 (µ, D, D) has a purely imaginary pair of eigenvalues. Let µ c2 meet this condition, as provided by part 1 of Theorem 3.7. The discriminant of the quadratic factor of the characteristic polynomial of J E 2 (µ, D 1 , D 2 ) is At µ c2 and for (D 1 , D 2 ) sufficiently close to (D, D), this is close to the expression which is negative, because it's the discriminant of the quadratic factor of the characteristic polynomial of J E 2 (µ, D, D) given in (21). Therefore, at µ c2 , the dis- is also negative. By continuity of the discriminant as a function of µ, there is an interval [µ c2 − δ 0 , µ c2 + δ 0 ] through which it is negative. Therefore, the characteristic polynomial of J E 2 (µ, D 1 , D 2 ) has a complex conjugate pair of roots on this interval, satisfying part of hypothesis 3 of Theorem 3.4.
However, preservation of the transversality condition is more challenging and requires the uniform approximation of γ (D 1 , D 2 ) by γ (D, D) when (D 1 , D 2 ) is sufficiently close to (D, D). This is the content of Lemma 3.10, the proof of which is relegated to Section 5 so as not to disturb the flow of the exposition. Lemma 3.10. Assume f 1 is three times continuously differentiable and f 2 is two times continuously differentiable. Then there exists a µ-interval [µ−δ, µ+δ] on which the µ-derivative γ (D 1 , D 2 )(µ) is uniformly approximated by γ (D, D)(µ) = A (µ). In fact, there exists a constant C such that for any µ ∈ [µ c2 −δ, µ c2 +δ]  D). Consequently, the Jacobian J E 2 (µ, D 1 , D 2 ) has a conjugate pair of imaginary eigenvalues when µ=µ c2 (D 1 , D 2 ).  Remark 3.12. Wolkowicz [11] considers a model in which the uptake function of the predator is Lotka-Volterra, i.e., f 2 (P ) = mP . There it is shown that no Hopf bifurcation can occur, regardless of the functional response of the prey. Since f 2 (P ) so defined cannot satisfy the condition D/ γ 2 λ Z (D) > f 2 λ Z (D) , the result is consistent with our own.
Provided that (D 1 , D 2 ) is sufficiently close to (D, D), parts 3 and 4 of the theorem hold by the same reasoning as for the corresponding parts of Theorem 3.7.

4.
Simulations and examples. We present some simulations illustrating concretely our results. First we consider system (1) and take f 1 and f 2 to be Holling type II. We choose notations as follows. With these choices explicit formulas can be given for many quantities studied in earlier sections. For example, we obtain formulas for the numbers λ P (D 1 ) and λ Z (D 2 ) defined in Equation (2). For λ P (D 1 ) we have For rate functions of this type, Equation (6) determining N (µ, D 1 , D 2 ) reduces to a quadratic equation. In principle, one obtains explicit solutions for N (µ, D 1 , D 2 ) and the nonnegative solution is easily identified. Then Equation (7) for Z(µ, D 1 , D 2 ) is linear and easily solved for Z(µ, D 1 , D 2 ).
With these remarks we can turn to an illustration of Theorems 3.7 and 3.11 in a case using Holling type II rate functions. With these rate functions, the concavity of the function f 2 guarantees that the condition D/ γ 2 λ Z (D) > f 2 λ Z (D) is satisfied. First, we set D=D 1 =D 2 =1.0 and the remaining parameters to the following values.
With all quantities involved explicitly computed, the formula given in Equation (19) for the real part of the complex conjugate pair of eigenvalues of the linearization of the system at the coexistence equilibrium can be made explicit, though messy, and is easily plotted by a computer algebra system. This is the solid curve in Figure 2, which shows that a Hopf bifurcation occurs in the vicinity of µ = 0.6.  exhibits solutions of the system for µ slightly smaller and slightly larger than the bifurcation value µ c2 (D, D) ≈ 0.6. Figure 3. Before and after a Hopf bifurcation: D=D 1 =D 2 =1 and Holling type II rate functions Next, keep D=1 and set D 1 =1.2 and D 2 =1.3. The proof of Theorem 3.11 says that the graph of the function defined by taking the real part of the complex conjugate pair of eigenvalues associated with the coexistence equilibrium is an increasing function whose graph lies in a neighborhood of the curve we discussed in the preceding paragraph. We do not have an explicit formula for this function with the new values for D 1 and D 2 , but we can estimate its values by numerically computing the eigenvalues of the linearization along a sequence of µ-values. Figure 2 also shows a sequence of points derived from eigenvalue approximations when D 1 =1.2 and D 2 =1.3. Interpolating a curve through the plotted points, we see a Hopf bifurcation occurs in the vicinity of µ = 0.9. Figure 4 shows trajectories of this system for values of µ slightly smaller and slightly larger than the bifurcation value. Then the condition D/ γ 2 λ Z (D) > f 2 λ Z (D) of Theorem 3.11 is satisfied, but this is not a consequence of the concavity of the graph of f 2 . Determining λ P (D 1 ) and λ Z (D 2 ) in this example requires solving quadratic equations, so obtaining exact values is quite easy. Consequently, one can compute explicitly from (5) the value µ c1 (D 1 , D 2 ) beyond which the coexistence equilibrium exists. Locating a coexistence equilibrium E 2 (µ, D 1 , D 2 ) requires solving a cubic equation for N (µ, D 1 , D 2 ), for which it is more appropriate to use numerical methods. We choose a sequence of µ values starting beyond µ c1 (D 1 , D 2 ), approximate the coexistence equilibria and their Jacobian matrices J E 2 (µ, D 1 , D 2 ) , and numerically compute the real part of the complex pair of eigenvalues of the Jacobian for each µ value. Plotting the real part against µ produces Figure 5, which exhibits the expected change of sign and shows that a Hopf bifurcation occurs in the vicinity of µ = 7.25; trajectories for parameters slightly smaller and slightly larger than the bifurcation value are shown in Figure 6.
The conservative nature of system (1) for the case D=D 1 =D 2 implies that the orbitally asymptotically stable limit cycle lies in the limiting plane given by Equation (20). In the general case, there is no expectation that the limit cycle lies in a plane. That the cycles in Figures 4 and 6 are not planar is hard to see in the presented figures, because the scales of the axes are not uniform. However, calculations of estimates for the torsion of the cycles in these figures yield non-zero values, so the cycles are not planar curves, by a basic result of differential geometry.
Backtracking farther, we recapitulate Theorem 3.7, where the working assumptions are that a value D is fixed, a coexistence equilibrium E 2 (µ, D, D) exists, and that there is a range of parameters µ for which the linearizations at the coexistence equilibrium have a complex conjugate pair of eigenvalues. Moreover, at the parameter value µ c2 , the linearization of the system has a purely imaginary pair of eigenvalues. For any µ slightly smaller than µ c2 , the pair of complex eigenvalues has a negative real part and for any µ slightly larger than µ c2 the pair of complex eigenvalues has a positive real part. By Lemma 3.1 there is an interval I containing µ c2 and a disc ∆ centered at (D, D) such that the functions N (µ, D 1 , D 2 ) and Z(µ, D 1 , D 2 ) defined on I×∆ smoothly parametrize the locus of coexistence equilibria near E 2 (µ, D, D). Lemma 3.10). Assume f 1 is three times continuously differentiable and f 2 is two times continuously differentiable. Then there exists a µ-interval [µ−δ, µ+δ] on which the µ-derivative γ (D 1 , D 2 )(µ) is uniformly approximated by γ (D, D)(µ) = A (µ). In fact, there exists a constant C such that The proposition follows from a sequence of lemmas and estimates, given below. In the course of proving these results, we find it necessary to impose the differentiability conditions on f 1 and f 2 .
An essential ingredient in the process is to obtain bounds on magnitudes of the differences We can write where the projection π 3 to the third coordinate γ in P 1 (x) − ×P 2 (x) − is represented by 0, 0, 1 . The gradients ∇α, ∇β, and ∇γ are evaluated at where explicit expressions for a i (µ, D 1 , D 2 ) are given in (11), (12), and (13). The derivatives p 0 , p 1 , and p 2 are evaluated at (µ, D 1 , D 2 ).
Applying the convention of (24), we have Before we go farther, we compress taking the derivative D(M −1 ) at (p 0 , p 1 , p 2 ) and evaluating on the vector (p 0 , p 1 , p 2 ), writing Then the previous equation becomes We can estimate using operator norms computed in terms of the Euclidean metrics.
since the norm of a projection is 1. Now we use the triangle inequality to bound the last expression. Now we bound the two summands on the right side of the inequality in (26). The argument for the second summand is more complicated than that for the first, so we give a detailed explanation.
for a constant C D 2 (M −1 ) independent of (D 1 , D 2 , µ). We compute This takes care of the first factor on the righthand side of (27). For the remaining factor in the bounding term in (27)  Proof. The proof is similar to the proof of Lemma 5.3 but slightly simpler, because the mean value theorem is not needed. For details, see [1].
We can now prove the uniform convergence result.
Proof of Proposition 5.1. Combining the inequalities (25) and (26), we have Using the observations detailed in Lemma 5.4 and Lemma 5.3, Thus, there is a constant C such that We have already used a technique of obtaining bounds by splitting quantities. As we will continue to exploit the technique in the following results, we formulate the Lemma 5.5 for reference, leaving the elementary proof to the reader.
1. There is a constant c 1 such that 2. There is a constant c 2 such that 3. There are constants c 3 and c 4 such that |Q 1 (µ, D 1 , D 2 )| ≤ c 3 for (µ, D 1 , D 2 ) ∈ I×∆ and for dist (D 1 , D 2 ), (D, D) sufficiently small, and |Q 2 (µ, D, D)| ≤ c 4 for µ ∈ I. Then there is a constant c 5 such that, for dist (D 1 , D 2 ), (D, D) sufficiently small, As has been seen, the proof of 5.1 depends on the following three propositions.
Proposition 5.6. There are constants C 0 and C 0 such that Proposition 5.7. There are constants C 1 and C 1 such that Proposition 5.8. There are constants C 2 and C 2 such that The proofs of Propositions 5.6, 5.7, and 5.8 depend in turn on a number of elementary bounds and estimates, given below in Lemma 5.9 and Propositions 5.10, 5.11, and 5.12. We give the quick proofs of Lemma 5.9 and Propositions 5.10 and 5.11, because they are quite short, postponing the discussion of the many parts of 5.12 to the end of the section. After we state these results, we prove Proposition 5.6, leaving the proofs of Propositions 5.7 and 5.8 to the reader. Lemma 5.9. Given D, there is an interval J containing λ Z (D) such that for all P ∈ J. Thus, for all P ∈ J, f 2 (P ) is bounded away from zero, and, for all is bounded away from zero. Proof. By assumption on f 2 , f 2 λ Z (D) > 0. By continuity of f 2 , there is an interval J containing λ Z (D) such that, for all P ∈ J, is bounded by some constant on the interval I.
Proof. Each of the listed functions is continuous on the closed interval I, so each one is bounded.
Proposition 5.11. Each of the quantities is bounded by some constant on the domain I×∆.
Proof. Each of the listed functions is continuous on the compact set I×∆, so each one is bounded.
The proof of the next result depends on many more details of system (1) and consequences drawn from them. We present the details for items selected to be typical, referring the interested reader to the preprint [1] for details on other items.  D, D) .
. Proof of Proposition 5.6. After some reorganization, we have from (13) Proof of Propositions 5.7 and 5.8. The arguments are essentially the same as those given for Proposition 5.6; we refer the reader to the preprint [1] for details. Now we embark on the proof of 5.12. Part of this work is made easier by the fact that the rate functions in system (1) are explicitly linear in D, D 1 , and D 2 . On the other hand, because other quantities such as N , P , and Z depend implicitly on D, D 1 , and D 2 , the details of the analyses are somewhat lengthy.
Proof of bound 1. By definition and the mean value theorem applied to f 2 , we have for some number P 1 between λ Z (D 2 ) and λ Z (D). Consequently, Since P 1 is also close to λ Z (D), Lemma 5.9 gives f 2 (P 1 ) > f 2 λ Z (D) /2. Thus, for D 2 sufficiently close to D Proof of bound 2. .
We also observe that f 1 N (µ, D 1 , D 2 ) is bounded by a constant depending only on I×∆, because of the convexity of the closed disc ∆ centered at (D, D) from which we choose (D 1 , D 2 ).

Proofs of bounds 3 and 4.
The ingredients are the same as in the proof of bound 2, but the manipulations differ somewhat.
Proof of bound 5. This is very similar to the proof of bound 1.
Proof of bound 6. At this point the continuity of f .