A Predator-Prey model in the chemostat with Holling Type II response function

A model of predator-prey interaction in a chemostat with Holling Type II functional and numerical response functions of the Monod or Michaelis-Menten form is considered. It is proved that local asymptotic stability of the coexistence equilibrium implies that it is globally asymptotically stable. It is also shown that when the coexistence equilibrium exists but is unstable, solutions converge to a unique, orbitally asymptotically stable periodic orbit. Thus the range of the dynamics of the chemostat predator-prey model is the same as for the analogous classical Rosenzweig-MacArthur predator-prey model with Holling Type II functional response. An extension that applies to other functional rsponses is also given.


Introduction
In this paper, we analyze a predator-prey model in the chemostat with Holling Type II predator response function of Monod form and prey response function of mass action form. The chemostat is a widely-used apparatus used in the study of microbial biology. It is helpful for the study of microbial growth and interactions under nutrient limitation in a controlled environment. Chemostats can be used as a guide for identifying the dynamical nature of population interactions that may be present in a more complex system such as a lake. There are many articles related to the study of the chemostat, both from the experimental and the modelling point of views (see for example [6], [15], [29] and [37]). Here, we look at a model of a basic chemostat setup in which a single, essential, non-reproducing nutrient is supplied to a growth chamber from a nutrient reservoir at a constant rate. A population of microorganisms, designated as the prey, lives in the growth chamber and it feeds on the nutrient and a predator population predates on the prey population. The growth chamber is assumed to be well-stirred. It is assumed that the inflow rate of nutrient is the same as the outflow rate from the growth chamber to the waste reservoir so that the volume of the growth chamber remains constant and all of the contents of the growth chamber are removed in proportion to their amount in the growth chamber. How the amounts of the nutrient, prey population, and predator populations change as time changes are all modelled.
The more familiar predator-prey models, introduced in Rosenzweig and MacArthur [34], involve models of predator-prey interactions in which the prey population reproduces and involves only two equations, one for the prey and one for the predator. In [34], models for which the prey nullcline can have at most a single local extremum, a local maximum, were considered. Such a model is often referred to as the Rosenzweig-MacArthur model, and it has been well studied (see, for example, [26] and [33]). The case where the prey nullcline has both a local maximum and a local minimum is possible and has also been considered as a generalized version of the Rosenzweig-MacArthur model (see for example, [1], [13], [14], [36] and [38]). Collectively, the Rosenzweig-MacArthur model and its generalizations are known as the classical predator-prey models.
The mathematical expressions needed for the response functions in these models are not usually known in practice. Biologists typically collect sets of data that can be fitted by many functions. Following the work of Holling [20], modellers have used three main forms to describe response functions in predator-prey models. Holling Type I refers to a mass-action response function, that is linearly increasing. Holling Type II responses are increasing and concave down. Holling Type III are sigmoidal. The dynamics of predator-prey models are directly influenced by the types of response functions chosen. Recently, in Fussmann and Blasius [14], it was shown that the range of possible dynamics can be different for classical predator-prey models modelled by different forms of Holling Type II response functions.
Here, we determine the global dynamics of a predator-prey model in the chemostat with predator response function of Holling Type II (Monod form) and prey response function of mass action form. We then compare the dynamics of this predator-prey model in the chemostat with that of the analogous classical predator-prey model.
Harrison [16] considered a wide class of classical predator-prey models and obtains sufficient (but not necessary) conditions for the global stability of the coexistence equilibrium. He also proved that when the coexistence equilibrium is unstable, at least one periodic orbit exists. In the special case that the prey grow logistically in the absence of predators and the predator response function is of Holling Type I form, Hsu [21] proved that the coexistence equilibrium is globally asymptotically stable whenever it exists. He also proves that when the predator response function is of Monod form, then the coexistence equilibrium is globally asymptotically stable whenever it is locally asymptotically stable, and Liou and Cheng [28] and Kuang and Freedman [27] proved that when the coexistence equilibrium exists and is unstable, it is surrounded by a unique periodic orbit.
In this paper, we focus on the analogous model of predator-prey interaction in a chemostat. In Wolkowicz [39], a food web in a chemostat was considered. In the special case that the model studied in [39] is a food chain that includes one resource, and only a single prey, and a single predator population, the model is of the basic form studied here. A Lyapunov function was used to prove that the coexistence equilibrium is globally asymptotically stable whenever it exists, under the assumption that the predator response function is Holling Type I and the prey response function is either Holling Type I or Holling Type II (of Monod form).
In the model considered in this paper, we assume instead that the prey response function is Holling Type I and the predator response function is Holling Type II (of Monod form). In this case, we prove that the dynamics are more complicated. In particular, we prove that whenever the coexistence equilibrium is locally asymptotically stable, it is globally asymptotically stable, and that whenever the coexistence equilibrium is unstable, there is a unique orbitally asymptotically stable periodic orbit. We also prove that the change in stability occurs by means of a Hopf bifurcation that is always supercritical. We thus show the similarity between the dynamics of the classical predator-prey model in which the prey is assumed to grow logistically in the absence of the predator population and its analogous chemostat predator-prey counterpart model. The chemostat model is of one dimension higher, since the resource that the prey population consumes in order to grow is also modelled and the growth of the prey in the absence of the predator population depends instead on the abundance of the resource. This paper is organized in the following manner. The predator-prey model in a chemostat is described in Section 2, where three equivalent lower dimensional limiting systems are also derived. The limiting system that most resembles the classical predator-prey model is then chosen to be the focal system. Properties of the prey nullcline of this system are derived in Section 3. Preliminary analytic results appear in Section 4, where it is also determined that the system undergoes a supercritical Hopf bifurcation. In Section 5, we prove that the coexistence equilibrium is globally asymptotically stable whenever it is locally asymptotically stable, and also that when the coexistence equilibrium is unstable, there is a unique, orbitally asymptotically stable periodic orbit. Finally, a discussion of the similarities between chemostat models and their analogous, classical predator-prey models is given in Section 7. Appendices provide a modification of a theorem due to Huang [22], an extension of a theorem due to Hsu [21], and the Hopf bifurcation analysis. The bifurcation diagrams were done using XPPAUT [12] and the simulations were done using Matlab [32].

The Model
Let S(t) denote the concentration of the nutrient in the growth chamber at time t, and let x(t) and y(t) denote the density of prey (that feeds off this nutrient) and predator populations, respectively. We consider the following system of autonomous ordinary differential equations as a model of predator-prey interaction in a well-stirred chemostat: where S 0 denotes the concentration of the nutrient in the nutrient reservoir, and c and γ are yield constants (respectively, to the prey's consumption of the nutrient and the predator's consumption of prey). To ensure that the volume of this vessel remains constant, D denotes both the rate of inflow from the nutrient reservoir to the growth chamber, as well as the rate of outflow from the growth chamber. We assume the species specific death rates are insignificant with respect to the flow rates and they are ignored. It is also assumed that the functions p and q are continuously differentiable (additional smoothness assumptions are given below) and that S 0 , D, c, and γ are all positive constants.
The rate of conversion of nutrient to biomass is given by the function p(S), and is assumed to satisfy p(0) = 0, p(S) > 0 for S > 0 and p (S) > 0 for S ≥ 0. For the remaining of this paper, we consider p to be of this form. The function q(x) denotes the predator response function. It is assumed that q(x) has properties similar to p(S). In particular, The Monod functional response q(x) = ax b+x satisfies these properties and will be the focus in the remainder of this paper. Note that the following inequality holds for that form: To simplify (2.1), the yield constants c and γ can be scaled out by performing a change of variables. LetŜ = σS,x = ξx,ŷ = ηy andt = τ t. Then, taking c = σ ξ , γ = ξ η ,Ŝ 0 = σS,D = D τ ,m = m ξτ ,â = aξ τ η , andb = ξb, and removing hats to simplify notation, reduces system (2.1) to: Remark 2.1. Since the system has four variables (S, x, y and t), we are able to scale out two more parameters for a total of four. However, this is not needed to complete our analysis.
Lemma 2.1. The solutions of S(t), x(t) and y(t) of (2.5) are non-negative and bounded.
Proof. It is important to first note that since the vector field is C 1 , existence and uniqueness of solutions hold. S(t) is positive for all t > 0 since S(τ ) = 0 for any τ ≥ 0 implies that S (τ ) > 0. Furthermore, x(t) > 0 and y(t) > 0 for all t > 0 since the (S, 0, y)-plane and (S, x, 0)-plane are invariant with respect to solutions of (2.5). Hence, by the uniqueness of solutions, they cannot be reached in finite time by trajectories for which x(0) > 0 or y(0) > 0, respectively.
Next consider the sum of our differential equations: is a first order ODE that has the solution Thus S(t) + x(t) + y(t) ≤ max{S(0) + x(0) + y(0), S 0 }. Since all three components are non-negative, this implies each component is bounded above.
Nonetheless, (2.5) is still a 3D system of equations, and so is harder to analyze than the classical predator-prey model. From the above proof, we can see that the sum of solutions S(t) + x(t) + y(t) converges to S 0 exponentially as t tends to infinity. This tells us that for any point (S, x, y) in the omega limit set, , t n → ∞ as n → ∞, S(t n ), x(t n ), y(t n ) → (S, x, y) as n → ∞ , it follows that S(t n ) + x(t n ) + y(t n ) → S + x + y = S 0 . Hence, ω is restricted to the simplex (S, x, y) : S + x + y = S 0 , a two dimensional set. We can then obtain three equivalent limiting systems by eliminating one variable in our system using the fact that S(t) + x(t) + y(t) = S 0 . In this case, it is most useful to eliminate S, as this yields a 2D predator-prey system with nullclines that resemble those of the classical predator-prey model. In doing so, if we define we obtain the limiting system, (2.8) Since (2.8) closely resembles the classical predator-prey model, it will be the main system we analyze throughout this paper. That it has the same dynamics as the 3D systems (2.5) and therefore (2.1), is justified in Section 6.
The predator nullcline is the vertical line x = q −1 (D) = bD a−D and is in the first quadrant if and only if a > D. We will assume, unless otherwise stated, that a > D so that the coexistence equilibrium point exists in the interior of the first quadrant. The prey nullcline, F (x), is studied in detail in the next section.

The Prey Nullcline
The prey nullcline is given by the continuously differentiable function . The properties of this function play a key role in our analysis. In this section, we determine the properties of F (x) when the predator response function is of Monod form, q(x) = ax b+x . Define We will assume that mS 0 > D, so that F (0) > 0. The slope of the prey nullcline is given by and with the sign determined by the sign of the numerator. For x ≥ 0, the second derivative of the prey nullcline is: (3.5) and the third derivative is: Though not discussed any further, it is clear that the higher order derivatives will alternate in sign.  Proof. If we assume that M = 0, the result follows, since then, mS 0 −D 2m > 0 under our assumptions. Next assume that M > 0. Then, Lemma 3.1 tells us that M is always less than half the distance to K. Furthermore, with a Monod response function, one can actually find an explicit expression for M , namely: Fig 1 illustrates how the shape of the graph of F depends on the sign of F (0). In Section 4, we select S 0 as the bifurcation parameter and investigate how changes in S 0 affect the shape of the graph of F . Since, F is an increasing function of S 0 for each fixed x. F is also an increasing function of S 0 for each fixed Remark 3.2. By equations (3.8) and (3.9), as S 0 increases, the local maximum of F moves up and to the right.
It also is important to note that since F (x) increases with S 0 , one can find a fixed S 0 value so that F (x) = 0, namely: for any fixed x > 0. If we take x = x * in S 0 (x), we get: Then, taking S 0 = S 0 crit forces F (x * ) = 0. Taking S 0 = S 0 crit will prove especially useful in our analysis.

Local Analysis
In this section, we consider some of the properties of the different invariant sets associated with (2.8). It is first important to note that non-negativity and boundedness of solutions for system (2.8) follows immediately from the non-negativity and boundedness of solutions of the 3D system (2.5). Nonnegativity and boundedness of solutions is a prerequisite of any reasonable model of the chemostat.
The Jacobian matrix of (2.8) is Evaluating it at the zero equilibrium, we obtain: Since q(x) is an increasing function and F (0) > 0, the diagonal entries have opposite signs. Hence, (0, 0) is always a saddle. Similarly for (K, 0), One can see that if x * > K, then q(K) − D = q(K) − q(x * ) < 0 and hence, (K, 0) would be a global attractor since F (K) < 0 as well. However, there is no coexistence equilibrium in this case. Therefore we assume that x * < K, which means the diagonal entries of (4.3) have opposite signs, and hence, (K, 0) is also a saddle. Finally, we consider the coexistence equilibrium point. Note once again, the coexistence equilibrium exists if and only if x * < K. In this case the Jacobian is given by which has characteristic equation: Since D + mx * > 0, the constant term is positive and the roots of (4.5) have negative real part if and only if F (x * ) < 0. Thus, (x * , y * ) is locally asymptotically stable when F (x * ) < 0 and unstable when F (x * ) > 0. Note that unless x * = K or F (x * ) = 0 all critical points are hyperbolic. When x * = K, we determine that the equilibrium point is asymptotically stable using standard phase plane analysis. Proof. Since system (2.8) is planar, any periodic orbit must surround an equilibrium point by the Poincaré-Bendixson Theorem [3]. By the non-negativity of this system, the only equilibrium a periodic orbit can surround is (x * , y * ). From phase plane analysis, any periodic orbit would lie in the set {(x, y) : 0 < x < K and y > 0}. The eigenvalues of the variational matrix about (x * , y * ) are When the coexistence equilibrium exists, these eigenvalues are ply imaginary if and only if F (x * ) = 0. The condition F (x * ) = 0 can be achieved by fixing S 0 = S 0 crit from (3.11). Also, the imaginary part of Solid lines correspond to stable equilibria, dashed lines correspond to unstable equilibria, filled circles correspond to stable periodic orbits and empty circles correspond to unstable periodic orbits. As S 0 increases there is a transfer of stability from (0, 0) to (K, 0) to (x * , y * ) by means of transcritical bifurcations, and finally a transfer of stability from (x * , y * ) to a stable periodic orbit by means of a Hopf bifurcation.
the eigenvalues is non-zero. Furthermore, the transversality condition holds, since the derivative with respect to S 0 of the real part of the eigenvalue at the Hopf bifurcation is positive by (3.9). Thus, the eigenvalues are complex in a neighbourhood of S 0 crit and cross the imaginary axis at S 0 = S 0 crit , implying that a Hopf bifurcation occurs there. The direction and stability of the bifurcating periodic orbit is determined by the sign of the following quantity, called the vague attractor condition: (4.7) which was determined using the algorithm in Marsden and McCracken [31] as outlined in Appendix C. Under our assumptions on the parameters, w < 0, hence, the Hopf bifurcation is always supercritical.
The bifurcation diagram in Fig 2 summarizes each of these local stability results.

Global Analysis
To establish the global dynamics of system (2.8), first examine periodic orbits. The following lemma together with the Poincaré criterion will be used to show that when periodic orbits exist, they must surround the local maximum, M, F (M ) .
Lemma 5.1. Let Γ be any periodic orbit of (2.8). Then This means that Proposition 5.2. Any periodic orbit of (2.8) must surround the local maximum of F (x).
Proof. Assume F (x) > 0 for the entire portion of F inside of Γ. Then C > 0 by Lemma 5.1, implying that Γ is an unstable periodic orbit by the Poincaré criterion [9]. Since any periodic orbit must surround (x * , y * ), it follows that F (x * ) > 0 and so (x * , y * ) would also be unstable, which is impossible. Using a similar argument, it also follows that F (x) < 0 for the entire portion of F inside of Γ is impossible. Thus, the slope of the portion of the prey nullcline inside any periodic orbit cannot be entirely of the same sign, i.e., it must change sign, and therefore any periodic orbit must surround the local maximum of F (x), M, F (M ) .
Global stability of the coexistence equilibrium point for the classical predator-prey model was studied by Harrison [16]. If one denotes the coexistence equilibrium point as (x * , y * ), he proved that this equilibrium is globally asymptotically stable, if it is locally asymptotically stable, i.e., in phase space it lies on prey nullcline, F (x), where its slope is negative (F (x * ) < 0), and as well, it lies below the vertical line y = F (0), i.e., y * < F (0). Hence, he only obtained a sufficient condition for the global stability.
Hsu [21] theorized more generally that if the prey nullcline is concave down, then wherever the interior equilibrium is locally stable, it is also globally stable. Although this is not true for all systems (some counterexamples were found in [19]), we prove for (2.8), that local asymptotic stability implies global stability by using the Dulac criterion and the Poincaré-Bendixson theorem, that the equilibrium destabilizes via a supercritical Hopf bifurcation, and that if a periodic orbit exists it is unique.
Let (x * , y * ) be locally asymptotically stable, i.e. x * ∈ [M, K]. Since the Hopf bifurcation at x * = M is supercritical, (x * , y * ) is also locally asymptotically stable there. We start using a similar argument to the proof given for Theorem 3.3 in Hsu [21]. However, it was pointed out in [8], that the proof in Hsu [21] is "not rigorously correct" and so we make modifications. Choose h(x, y) = mx+q(x) −1 y β−1 as the auxiliary function to use with the Dulac criterion, where the value β > 0 will be determined.
Here, the function h(x, y) is defined in the interior of the first quadrant. Let f = x y . Then, In the interior of the first quadrant, y β−1 > 0 and mx + q(x) > 0, hence, ∆ changes sign if and only if H(x) changes sign. As opposed to the argument in [21], it is important to note that F (x) ≥ 0 and . Therefore, the choice of β is critical in order to ensure that H(x) does not change sign. As a consequence, we have the following proposition.
. If there exits β > 0 such that max Note that β(x) is parameterized by S 0 . We next look to see how changes in S 0 affect β(x). Since, β(x; S 0 ).
Proof. Consider the derivative with respect to x of β(x; S 0 crit ): > 0 for x ≥ 0, hence, β(x; S 0 crit ) is an increasing function on [0, K]. This, together with (5.4) and Remark 3.2, tells us that β crit > β(x; S 0 crit ) > β(x; S 0 ) for x < M , and β crit < β(x; S 0 crit ) < β(x; S 0 ) for x > x * . Finally at the boundaries, β(M ; S 0 ) = 0 < β crit , and lim x→x * + β(x; S 0 ) = ∞ > β crit .  can be used to determine when the Dulac criterion can be used to obtain global stability of the coexistence equilibrium using the function H(x) defined in (5.2). In particular, a value β crit can be chosen so that the Dulac criterion can be used to prove global stability of the coexistence equilibrium if S 0 < S 0 crit , i.e., x * < M , which is precisely when the coexistence equilibrium is locally asymptotically stable, but when S 0 > S 0 crit , i.e., x * < M , and the coexistence equilibrium is unstable, no value β crit can be chosen.
Theorem 5.5. Consider system (2.8). If (x * , y * ) is locally asymptotically stable, then it is globally asymptotically stable. If (x * , y * ) is unstable, then there exists a unique, stable periodic orbit that surrounds the point (M, F (M )).
Proof. Since β crit satisfies Proposition 5.3, H(x) ≤ 0 and consequently ∆ ≤ 0 for 0 ≤ x ≤ K. Therefore, since ∆ does not change sign, by the Dulac criterion it follows that system (2.8) has no nontrivial closed orbits lying entirely in the first quadrant. Thus, since we have proved that all orbits are bounded (see Lemma 2.1), by the Poincaré-Bendixson theorem, (x * , y * ) is globally asymptotically stable.
To determine uniqueness of the periodic orbit when x * < M , apply a modified version of Huang's theorem from [22], stated as Theorem A.1 in Appendix A by taking φ(x) = mx + q(x), ψ(x) = q(x) − D, π(y) = y and ρ(y) = y. Then conditions (i)-(iii) of Theorem A.1 are satisfied. For condition (iv) of Theorem A.1, notice that the function H(x) is identical to β(x) from (5.3). We had β (x) > 0 from (5.7), (5.4), and Remark 3.2, so we also have H (x) > 0. Thus, condition (iv) of Theorem A.1 is satisfied, and so when x * < M , there is a unique periodic orbit and it is orbitally asymptotically stable. Fig. 4. Sample trajectories of (2.8) with parameters a = 2, b = 0.5, m = 2, and D = 1.3788. Plot (a) depicts globally asymptotically stable convergence to (K, 0), when x * > K. Plot (b) depicts convergence to the globally asymptotically stable interior equilibrium, when it is to the right of the local maximum. Plot (c) depicts an orbitally asymptotically stable periodic orbit surrounding the unstable interior equilibrium, when it is to the left of the local maximum.

Fig 4 illustrates
convergence to the globally asymptotically stable equilibrium (K, 0) when the coexistence equilibrium does not exist and the existence of the two types of dynamics that are possible when a coexistence equilibrium point does exist: convergence to a globally asymptotically stable coexistence equilibrium point or convergence to an orbitally asymptotically stable periodic orbit that attracts all solutions with positive initial conditions except the unstable coexistence equilibrium point.
In Appendix B, we extend Theorem 5.5 to a more general class of predator-prey models, and give examples of systems that satisfy the hypotheses so that a β can be found in order to prove global stability of the coexistence equilibrium.

Dynamics of the 3D System
The dynamics of the 2D system (2.8) have been analyzed in great detail. We next justify why studying this system is equivalent to studying the 3D system (2.5), and hence, the original system (2.1).
Any point (x, y) in the 2D system (2.8) corresponds to a point (S 0 − x − y, x, y) in the 3D system (2.5), and so solutions of the 2D system correspond to solutions of the 3D system that lie on the 2D simplex S = {(S, x, y) ∈ R 3 : S, x, y > 0, S + x + y = S 0 }. Thus, the two systems share some of the same properties. The equilibrium points (0, 0), (K, 0) and (x * , y * ) of the 2D system correspond to (S 0 , 0, 0), (S 0 − K, K, 0) and (S 0 − x * − y * , x * , y * ), respectively, for the 3D system. In terms of local stability, the additional eigenvalue of the Jacobian matrix is negative, since solutions of the 3D system converge exponentially to S. Thus, the local stability of the 3D equilibrium points is the same as for the corresponding 2D equilibrium point. Since all periodic orbits of the 3D system must lie on S, the existence and number of periodic orbits of the two systems is the same.
From Theorem 5.5 and standard methods for asymptotically autonomous systems (see Smith and Waltman [37], or using the Butler-McGehee Lemma [5] directly), we obtain the following theorem.
Theorem 6.1. Consider system (2.5). If (x * , y * ) is locally stable for the 2D system (2.8), then it is globally stable, and (S 0 − x * − y * , x * , y * ) is globally asymptotically stable for the 3D system (2.5). If (x * , y * ) is unstable in 2D system (2.8), then there is a unique periodic orbit that lies on the S+x+y = S 0 simplex and it is orbitally asymptotically stable. Corollary 6.2. Consider the original system (2.1). If the coexistence equilibrium exists and is stable, then it is globally asymptotically stable, and if it is unstable, then there is a unique periodic orbit that lies on the {(S, x, y) ∈ R 3 : σS, ξx, ηy > 0, σS + ξx + ηy = σS 0 } simplex, and it is orbitally asymptotically stable.

Discussion
A system of ODEs modeling predator-prey interactions in a chemostat was analyzed assuming a predator response function of Monod form. It was shown that whenever the coexistence equilibrium is locally asymptotically stable, it is also globally asymptotically stable, and whenever the coexistence equilibrium is unstable, there is a unique, orbitally asymptotically stable periodic orbit.
These results are consistent with the dynamics of the analogous classical predator-prey model with Monod predator response function, in whichthe resource and how the prey grows based on the amount of resource available is not modelled, but instead the prey is assumed to grow logistically in the absence of the predator population. This classical model has been studied extensively. See for example, [8,27,28,36] the resource is not modelled, but instead the prey is assumed to grow logistically in the absence of the predator. These authors found, just as for the predator-prey model in the chemostat studied in this paper, that whenever the coexistence equilibrium is locally asymptotically stable, it is globally asymptotically stable. It was also shown in Cheng [7] that when a periodic orbit exists around an unstable coexistence equilibrium in the classical predator-prey model with a Monod functional response, it is unique, and hence any nontrivial periodic orbit is asymptotically stable. However, there was a gap in the proof, that was later corrected by Liou and Cheng [28]. Another proof of the uniqueness of the limit cycle for the classical model in the case of the Monod functional response was also given in Kuang and Freedman [27].
That the range of dynamics of the model studied here, and of the analogous classical predator-prey model is basically the same, is not entirely surprising. For example, after Hastings and Powell [17] showed that a three-species food chain with Monod response functions in which the population at the lowest tropic level grows logistically in the absence of a predator population could have chaotic dynamics, Daoussis [10] showed this was also the case for the analogous chemostat food-chain model.
Classical predator-prey models have been shown to be sensitive to the mathematical form used to model the predator response function, even when the forms have the same qualitative shape by Fussmann and Blasius [14]. They considered three mathematical forms: the Monod form, the Ivlev form [23], and the Hyperbolic tangent form [24], all of which are monotone increasing and concave down and are nearly indistinguishable when appropriate parameters are chosen (see Fig 5). Fussmann and Blasius provided a similar figure and demonstrated that the qualitative and quantitative dynamics predicted by models with these response functions can be quite different. The model with the three . Functions with such a shape are said to be of Holling Type II form [20]. different response function forms was studied in more detail using a bifurcation theory approach in Seo and Wolkowicz [36].
In the case of the classical model with the Hyperbolic tangent response function, Seo and Wolkowicz [36] proved that the Hopf bifurcation is always supercritical when it occurs at the local maximum of the prey nullcline, but can be either super or subcritical when it occurs at the local minimum. They also studied the dynamics in more detail for all three forms of the response functions using a one and two parameter bifurcation approach, and found that in the Hyperbolic tangent case, two limit cycles surrounding a stable coexistence equilibrium can arise through a saddle-node bifurcation of limit cycles when the Hopf bifurcation at the local minimum is subcritical. Seo and Wolkowicz [35] also considered the classical predatory-prey model, with a functional response of arctan form and proved that when the coexistence equilibrium is locally asymptotically stable, more than one limit cycle is possible, providing a counterexample to a result in Attili and Mallak [2]. The classical predator-prey model with Ivlev response functions was analyzed by Kooij and Zegeling [25]. They proved that model has a similar range of dynamics as the model with Monod response function. The analogous chemostat models in the cases of the hyperbolic tangent and the arctan as the response function was studied in Eastman [11] and the analogous model with the Ivlev response function was studied in Bolger [4]. It was also shown that in these cases the analogous chemostat models have a similar range of dynamics when compared with their analogous classical predator-prey models.
there is a positive equilibrium point (x * , y * ), (ii) all functions in (A.1) are C 1 in the interior of R 3 + , and F (x) is continuous in the interior of φ (x) > 0 and ψ (x) > 0 for x > 0, ρ (y) > 0 and π (y) > 0 for y > 0, and is non-decreasing for 0 < x < x * and x * < x < K.
Then, system (A.1) has at most one limit cycle in the first quadrant, and, if it exists it is stable.

Appendix Appendix B Extension of Hsu's theorem
Consider the following predator-prey system: is a generalized version of the system that Hsu studied in [21]. Here, the prey nullcline is given by the function γ −1 F (x) , and the interior equilibrium is the unique point (x * , y * ) satisfying q(x * ) = D and γ(y * ) = F (x * ). Hsu [21,Theorem 3.3] conjectured that if (x * , y * ) is stable and the prey nullcline is concave down, then (x * , y * ) is globally stable. Since this conjecture is not true, as demonstrated by the counter example given by Hofbauer and So [19], we state the following theorem that was shown to be satisfied in Seo and Wolkowicz [36] for the classical predator-prey model with a Hyperbolic tangent response function. The following theorem can also be shown to be satisfied if, for example, ξ(x) = q(x), γ(y) = η(y) = y, F (x) = q(x) , and q(x) = a tanh(bx).
Then (x * , y * ) is globally stable provided there exists a β > 0 such that: Remark B.1. If M = 0 then the Lyapunov function in Harrison [16] can be used to show global stability.
Proof. To show that (x * , y * ) is globally stable, we will use a similar argument as Hsu. That is, we will show that there are no closed orbits in the first quadrant using the Dulac Criterion. Define the auxiliary function to be h(x, y) = ξ(x) −1 η(y) β−1 , where β > 0. The resulting divergence is then where f = x y and H(x) is given by: In the first quadrant, η(y) β−1 > 0 and ξ(x) > 0 by assumptions (iii) and (v), hence, ∆ will only change sign if H(x) changes sign. In Hsu's proof, he tried to show that H(x) ≤ 0 by choosing an appropriate .

(B.5)
Since the denominator is positive by assumption (iv), F (x) has the same sign as d dx γ −1 F (x) . The problem with the β Hsu chose is that it neglects the fact that F (x) ≥ 0 for x ∈ [0, M ] by assumption (vii). By assumptions (ii) and (vii), we can conclude that x * ≥ M . Despite that −D + q(x) ≤ 0 in [0, M ] by assumption (vi), if β is too small (as the one he chose), H(x) > 0. To guarantee that However, β can not be too large either. Since , then H(x) ≤ 0. Since H does not change sign, ∆ does not change sign. Hence, by the Dulac criterion, it follows that system (B.1) has no nontrivial orbits lying entirely in the first quadrant. Thus by the Poincaré-Bendixson theorem, (x * , y * ) is globally stable. We next include another example of a system where such a β is obtainable under the above assumptions. Consider the classical predator-prey model with a Monod response function. In this notation, this would be equivalent to letting ξ(x) = q(x), γ(y) = η(y) = y, F (x) = q(x) , and q(x) to be of Monod form.
Note that since γ(y) = y, γ −1 F (x) = F (x) and hence, that the interior equilibrium is given by x * , F (x * ) . We can then equivalently transform the parameters of the Monod response function by a → Dm and b → x * (m − 1), wherem > 1. Then q(x) = Dmx x * (m−1)+x . We first show that these functions satisfy the assumptions of Theorem B.1.
It is important to note that under assumption (i), the interior equilibrium (x * , y * ) lies in the positive quadrant. Moreover since η(0) = 0 and ξ(0) = 0, the point (0, 0) is also an equilibrium point. Thus, the x and y axes are both nullclines.
To satisfy assumption (vii), we require x * ∈ [M, K]. First assume x * ∈ (M, K]. Since all hypotheses of Theorem B.1 are met, we proceed to the find the β as described by (B.2). Let We can then determine where the maximum and minimum of β(x) occur, by finding its critical values.
Taking the first derivative, we get: Note that the sign of β (x) depends on the numerator. The roots of x 2 − 2xx * + M x * = 0 are and are both positive. The following two lemmas will be useful in determining which of these values is the local maximum of β (the other being the local minimum), and whether these values lie in our desired regions [0, M ] and [x * , K].
Proof. First consider x − . We proceed using proof by contradiction to show that x − ∈ [0, M ]. Suppose that x − > M . Then which is and x − ∈ [0, M ], we can conclude that β − is a local maximum, and hence, that Now consider x + . Clearly x + > x * . Moreover, β (x) has a removable singularity at x = x * and since sign lim β(x) is decreasing between x * and x + .
Proof. First assume that x * ≤ K 2 2K−M . Since β(x) is decreasing between x * and x + , if we can show that β(K) ≥ 0, then a local minimum occurs in [x * , K]. Furthermore, since x − is not in this region, it would mean that this local minimum must be β(x + ), and hence, x + ≤ K. Indeed, This means that in this case, min Now assume that x * > K 2 2K−M . Then, β(K) < 0, and hence, x + > K. Furthermore, since β(x) is decreasing between x * and x + , it sly is decreasing on [x * , K]. Thus the minimum value it attains on that region is at K, namely min Lemma B.4. For our positive quantities β − , β + and β K , In any case, there exists a positive β between the two quantities when x * ∈ (M, K]. β(x).

Appendix Appendix C Analysis of the Hopf bifurcation
Computation of (4.7) was done using the computer algebra system Maple [30], as provided in the supplementary material [18]. The following is a summary of the algorithm used, highlighting the main results.
The formula in Marsden and McCracken [31] is localized to where the Hopf bifurcation occurs, and thus, we assume that we are near the critical value of our bifurcating parameter, S 0 crit . To use this formula, we first need matrix (4.4) in real Jordan canonical form. That is, we need to find an invertible matrix P so that where α ± iβ are the eigenvalues of A.
Lemma C.1. Let the eigenvector for α + iβ be P Re + iP Im . Then P = P Re P Im .
Proof. We have: P −1 P = P −1 P Re P Im = P −1 P Re P −1 P Im .