Mathematics in Applied Sciences and Engineering https://ojs.lib.uwo.ca/mase Volume 1, Number 4, December 2020, pp.333-354 https://doi.org/10.5206/mase/10842 A PREDATOR-PREY MODEL IN THE CHEMOSTAT WITH HOLLING TYPE II RESPONSE FUNCTION TEDRA BOLGER, BRYDON EASTMAN, MADELEINE HILL, AND GAIL S. K. WOLKOWICZ Abstract. A model of predator-prey interaction in a chemostat with Holling Type II predator func- tional response of the Monod or Michaelis-Menten form, is considered. It is proved that local as- ymptotic 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 re- sponses is also given. 1. 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 Received by the editors 30 June 2020; accepted 1 December 2020; published online 22 December 2020. 2010 Mathematics Subject Classification. Primary 34C25, 34C60, 34D05, 34D23 Secondary 92B05, 34C23, 92D40. Key words and phrases. Hopf bifurcation, global dynamics, predator-prey system, chemostat, Monod form. G.S.K.W was supported in part by a Natural Science and Engineering Council of Canada (NSERC) Discovery Grant with Accelerator supplement. 333 https://ojs.lib.uwo.ca/mase https://dx.doi.org/https://doi.org/10.5206/mase/10842 334 T. BOLGER, B. EASTMAN, M. HILL, AND G.S.K. WOLKOWICZ 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 de- scribed 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 PREDATOR-PREY CHEMOSTAT MODEL: HOLLING TYPE II RESPONSE 335 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]. 2. 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: S′(t) = (S0 −S)D −xp(S) x′(t) = x(−D + cp(S)) −yq(x) y′(t) = y ( −D + γq(x) ) S(0) ≥ 0, x(0) ≥ 0, y(0) ≥ 0. (2.1) where S0 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 S0,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. (2.2) A function of the form p(S) = mS where m is a positive constant, satisfies all of the above properties. 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, q(0) = 0, q′(x) > 0 for x ≥ 0 and q′′(x) < 0 for x ≥ 0. (2.3) 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: q(x) −xq′(x) > 0 for all x > 0. (2.4) To simplify (2.1), the yield constants c and γ can be scaled out by performing a change of variables. Let Ŝ = σS, x̂ = ξx, ŷ = ηy and t̂ = τt. Then, taking c = σ ξ , γ = ξ η , Ŝ0 = σS, D̂ = D τ , m̂ = m ξτ , â = aξ τη , and b̂ = ξb, and removing hats to simplify notation, reduces system (2.1) to: 336 T. BOLGER, B. EASTMAN, M. HILL, AND G.S.K. WOLKOWICZ S′(t) = (S0 −S)D −mSx x′(t) = x(−D + mS) −y ax b + x y′(t) = y ( −D + ax b + x ) S(0) ≥ 0, x(0) ≥ 0, y(0) ≥ 0. (2.5) 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 C1, 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: S′ + x′ + y′ = (S + x + y)′ = D ( S0 − (S + x + y) ) (2.6) is a first order ODE that has the solution S(t) + x(t) + y(t) = e−Dt (( S(0) + x(0) + y(0) ) −S0 ) + S0. (2.7) Thus S(t) + x(t) + y(t) ≤ max{S(0) + x(0) + y(0),S0}. 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 S0 exponentially as t tends to infinity. This tells us that for any point (S,x,y) in the omega limit set, ω = { (S,x,y) ∈ R3+ : ∃{tn} ∞ n=1 , tn →∞ as n →∞, ( S(tn),x(tn),y(tn) ) → (S,x,y) as n →∞ } , it follows that S(tn) + x(tn) + y(tn) → S + x + y = S0. Hence, ω is restricted to the simplex{ (S,x,y) : S + x + y = S0 } , a two dimensional set. We can then obtain three equivalent limiting sys- tems by eliminating one variable in our system using the fact that S(t) + x(t) + y(t) = S0. 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 F(x) = −Dx + mx(S0 −x) mx + q(x) , we obtain the limiting system, x′(t) = x ( −D + m(S0 −x−y) ) −y ax b + x , ( mx + ax b + x )( F(x) −y ) y′(t) = y ( −D + ax b + x ) x(0) ≥ 0, y(0) ≥ 0. (2.8) PREDATOR-PREY CHEMOSTAT MODEL: HOLLING TYPE II RESPONSE 337 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. 3. The Prey Nullcline The prey nullcline is given by the continuously differentiable function F(x) = −Dx+mx(S0−x) mx+q(x) . 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 F(0) = lim x→0+ F(x) = (mS0 −D)b a + mb . (3.1) We will assume that mS0 > D, so that F(0) > 0. The slope of the prey nullcline is given by F ′(x) = −(b + x)2 m2 −a (b + 2 x− S0 ) m−aD (m (b + x) + a) 2 . (3.2) Define K such that F(K) = 0, namely K = mS 0−D m > 0. Then, F ′(K) = − S0m + bm−D S0m + bm + a−D < 0 (3.3) and F ′(0) = −b2m2 + S0am−abm−aD (bm + a)2 , (3.4) with the sign determined by the sign of the numerator. For x ≥ 0, the second derivative of the prey nullcline is: F ′′(x) = − 2ma(S0m + bm + a−D) (bm + mx + a)3 < 0, (3.5) and the third derivative is: F ′′′(x) = 6m2a(S0m + bm + a−D) (bm + mx + a)4 > 0. (3.6) Though not discussed any further, it is clear that the higher order derivatives will alternate in sign. Since F is concave down for x ≥ 0, it has at most one hump on [0,K]. Define M such that F(M) = max x∈[0,K] F(x). Remark 3.1. It follows immediately that (a) if F ′(0) ≤ 0, then M = 0, and (b) if F ′(0) > 0, then 0 < M < K and F ′(M) = 0. It is useful to know where M lies in this case. Lemma 3.1. Assume mS0 > D. Then, M < mS 0−D 2m . 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, F ′(M) = 0. By (3.2), the derivative of F at mS 0−D 2m is F ′ (mS0 −D 2m ) = − ( (mS0−D)2 4m )( m + q′ ( mS0−D 2m )) ( mS0−D 2 + q ( mS0−D 2m ))2 < 0 338 T. BOLGER, B. EASTMAN, M. HILL, AND G.S.K. WOLKOWICZ (a) F ′(0) < 0 (b) F ′(0) = 0 (c) F ′(0) > 0 Fig. 1. Permissible configurations of the prey nullcline. In plots (a) and (b) F(x) is strictly decreasing for all x ∈ [0,K], whereas in plot (c) F(x) has a single local extremum, a local maximum, at M ∈ (0,K). since q′(x) > 0 for x > 0. Thus, since F is concave down on [0,K] and F ′(M) = 0, then M < mS0−D 2m . � 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: M = −mb−a + √ S0am + abm + a2 −aD m . (3.7) Fig 1 illustrates how the shape of the graph of F depends on the sign of F ′(0). In Section 4, we select S0 as the bifurcation parameter and investigate how changes in S0 affect the shape of the graph of F . Since, ∂ ∂S0 F(x) = m(b + x) bm + mx + a > 0 for x > 0, (3.8) F is an increasing function of S0 for each fixed x. F ′ is also an increasing function of S0 for each fixed x > 0 since ∂ ∂S0 F ′(x) = ma (mb + mx + a)2 > 0. (3.9) Remark 3.2. By equations (3.8) and (3.9), as S0 increases, the local maximum of F moves up and to the right. It also is important to note that since F ′(x) increases with S0, one can find a fixed S0 value so that F ′(x) = 0, namely: S0 , S0(x) = (b + x) 2 m2 + a (b + 2 x) m + aD am , (3.10) for any fixed x > 0. If we take x = x∗ in S0(x), we get: S0crit , S 0(x∗) = (abm + aD −D2)(bm + a−D) m(a−D)2 . (3.11) Then, taking S0 = S0crit forces F ′(x∗) = 0. Taking S0 = S0crit will prove especially useful in our analysis. PREDATOR-PREY CHEMOSTAT MODEL: HOLLING TYPE II RESPONSE 339 4. 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). Non- negativity and boundedness of solutions is a prerequisite of any reasonable model of the chemostat. There are three possible critical points of (2.8). The first two, (0, 0) and (K, 0), always exist. Define x∗ = q−1(D) = bD a−D > 0 since we are assuming that a > D, and let y ∗ = F(x∗) > 0. Then the interior equilibrium, (x∗,y∗), exists only if x∗ < K. If x∗ = K, then (x∗,y∗) coalesces with (K, 0). The following is the stability analysis for these critical points. The Jacobian matrix of (2.8) is J(x,y) = ( F ′(x) ( q(x) + mx ) + ( q′(x) + m )( F(x) −y ) − ( mx + q(x) ) yq′(x) −D + q(x) ) . (4.1) Evaluating it at the zero equilibrium, we obtain: J(0, 0) = (( q′(0) + m ) F(0) 0 0 −D ) . (4.2) 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), J(K, 0) = ( F ′(K) ( q(K) + mK ) − ( mK + q(K) ) 0 −D + q(K) ) . (4.3) 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 A , J(x∗,y∗) = ( F ′(x∗)(D + mx∗) −(mx∗ + D) y∗q′(x∗) 0 ) , (4.4) which has characteristic equation: λ2 −λF ′(x∗)(D + mx∗) + y∗q′(x∗)(D + mx∗) = 0. (4.5) 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. Proposition 4.1. System (2.8) has a supercritical Hopf bifurcation occurring when F ′(x∗) = 0. 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 λ(S0) = F ′(x∗)(D + mx∗) ± √( F ′(x∗)(D + mx∗) )2 − 4y∗q′(x∗)(D + mx∗) 2 (4.6) 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 S0 = S0crit from (3.11). Also, the imaginary part of 340 T. BOLGER, B. EASTMAN, M. HILL, AND G.S.K. WOLKOWICZ Fig. 2. Bifurcation diagram of system (2.8) with parameters a = 2, b = 0.5, m = 2, and D = 1.3788. 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 S0 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 S0 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 S0crit and cross the imaginary axis at S 0 = S0crit, 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: w = − 4am ( (b + x∗)2(b + 1 2 x∗)m2 + 2(a− 1 4 D)(b + x∗)2m + a ( (a + D)b + Dx∗ ))( (b + S0)m + a−D ) ( (b + x∗)m + a )4 (b + x∗)2 , (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. 5. 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 C , ∮ Γ div(x′,y′) dt = ∮ Γ ( mx + q(x) ) F ′(x) dt. PREDATOR-PREY CHEMOSTAT MODEL: HOLLING TYPE II RESPONSE 341 Proof. Let g(x) = −D+m(S0−x). Then F(x) = xg(x) mx+q(x) , and thus, F ′(x) = g(x)+xg′(x) mx+q(x) − xg(x) ( m+q′(x) )( mx+q(x) )2 . This means that C = ∮ Γ div(x′,y′) dt = ∮ Γ [( xg′(x) + g(x) −y ( m + q′(x) )) + ( −D + q(x) )] dt = ∮ Γ [( xg′(x) + g(x) + (x′ −xg(x) mx + q(x) )( m + q′(x) )) + y′ y ] dt = ∮ Γ [( mx + q(x) ) F ′(x) + x′ ( m + q′(x) mx + q(x) ) + y′ y ] dt = ∮ Γ [( mx + q(x) ) F ′(x) + d dt ln ( mx + q(x) ) + d dt ln(y) ] dt = ∮ Γ ( mx + q(x) ) F ′(x) dt. � 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, ∆ ,∇· (hf) = yβ−1H(x) mx + q(x) (5.1) 342 T. BOLGER, B. EASTMAN, M. HILL, AND G.S.K. WOLKOWICZ where H(x) = ( mx + q(x) ) F ′(x) + β ( −D + q(x) ) . (5.2) 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 −D + q(x) < 0 for x ∈ [0,M], while F ′(x) < 0 and −D + q(x) ≥ 0 for x ∈ [x∗,K]. 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. Proposition 5.3. Let β(x) = − ( mx+q(x) ) F′(x) −D+q(x) . If there exits β > 0 such that max x∈[0,M] β(x) ≤ β ≤ min x∈[x∗,K] β(x), then H(x) ≤ 0, for all x ∈ [0,K]. Remark 5.1. Any positive β would work for x ∈ [M,x∗], since both F ′(x) ≤ 0 and −D + q(x) ≤ 0. First consider when x∗ ∈ (M,K]. For system (2.8) with a Monod functional response, the function β(x) is given by: β(x; S0) = ( (b + x)2m2 + a(b + 2x−S0)m + aD ) x ((a−D)x− bD)(m(b + x) + a) . (5.3) Note that β(x) is parameterized by S0. We next look to see how changes in S0 affect β(x). Since, ∂ ∂S0 β(x; S0) = − amx (a−D)(x−x∗)(bm + mx + a) , (5.4) β(x) is increasing with respect to S0 if x < x∗ (and hence, if x < M) and decreasing with respect to S0, if x > x∗. Using (3.10), we can find the value of S0 that forces x∗ = M, namely S0crit from (3.11). Since x∗ > M with the original S0, S0crit > S 0 by Remark 3.2. Thus, β(x; S0crit) > β(x; S 0) when x < M and β(x; S0crit) < β(x; S 0) when x > x∗. The function β(x; S0crit) is given by: β(x; S0crit) = mx(2abm + amx−Dbm−Dmx + 2a2 − 2aD) (bm + mx + a)(a−D)2 . (5.5) Note that the function β(x; S0crit) in (5.5) is a continuously differentiable function. Let βcrit be the value of β(x) with S0 = S0crit so that x = x ∗ = M, namely: βcrit = β(x ∗; S0crit) = 2Dbm (a−D)2 > 0. (5.6) Lemma 5.4. If x∗ > M, then max x∈[0,M] β(x; S0) < βcrit < min x∈[x∗,K] β(x; S0). Proof. Consider the derivative with respect to x of β(x; S0crit): β′(x; S0crit) = (2a−D)b2m3 + 2(a−D)bm3x + (a−D)m3x2 (bm + a)(a−D)2 + abm2(4a− 3D) + 2am2x(a−D) + 2a2m(a−D) (bm + a)(a−D)2 (5.7) > 0 for x ≥ 0, hence, β(x; S0crit) 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; S0) = 0 < βcrit, and lim x→x∗+ β(x; S0) = ∞ > βcrit. PREDATOR-PREY CHEMOSTAT MODEL: HOLLING TYPE II RESPONSE 343 (a) S0 < S0crit (b) S 0 > S0crit Fig. 3. Graphs of β(x;S0) (solid), β(x;S0crit) (dashed), βcrit (dash-dotted) and the vertical asymptote x = x∗ (dotted) with parameters a = 2, b = 0.5, m = 2, and D = 1.3788. Plot (a) illustrates when S0 < S0crit, i.e., x ∗ > M, and hence, the coexistence equilibrium is locally asymptotically stable, the function β(x;S0) can be used to show that if β = βcrit, that the function H(x) defined in (5.2) is of one sign (see Proposition (5.3), Remark 5.1, and Lemma 5.4)), and hence, can be used to prove global stability of the coexistence equilibrium using the Dulac criterion. Plot (b) illustrates that when S0 > S0crit, and hence, the coexistence equilibrium is unstable, as expected, the function H(x) changes sign for any choice of β. Thus, altogether, max x∈[0,M] β(x; S0) < βcrit < min x∈[x∗,K] β(x; S0). � Finally, we consider the case when x∗ = M. Then β(x; S0) is precisely β(x; S0crit), an increasing function with respect to x, and thus, max x∈[0,M] β(x; S0) = βcrit = min x∈[x∗,K] β(x; S0). Fig 3 demonstrates that the function β(x; S0) 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 S0 < S0crit, i.e., x ∗ < M, which is precisely when the coexistence equilibrium is locally asymptotically stable, but when S0 > S0crit, 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 sur- rounds 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. � 344 T. BOLGER, B. EASTMAN, M. HILL, AND G.S.K. WOLKOWICZ (a) S0 < x∗ + D m (b) x∗ + D m < S0 < S0crit (c) S0 > S0crit 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 coex- istence 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. 6. 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 (S0 −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) ∈ R3 : S,x,y > 0,S + x + y = S0}. 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 (S0, 0, 0), (S0 −K,K, 0) and (S0 −x∗ −y∗,x∗,y∗), respectively, for the 3D system. In terms of local PREDATOR-PREY CHEMOSTAT MODEL: HOLLING TYPE II RESPONSE 345 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 (S0 − 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 = S0 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) ∈ R3 : σS,ξx,ηy > 0,σS + ξx + ηy = σS0} simplex, and it is orbitally asymptotically stable. 7. 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 346 T. BOLGER, B. EASTMAN, M. HILL, AND G.S.K. WOLKOWICZ Fig. 5. Graphs of the three mathematical forms considered in [14]: the Monod form fM (x) = aM bM +x (dashed), the Ivlev form fI(x) = aI(1 − e−bI x) (solid) and Hyperbolic tangent form fT (x) = aT tanh(bT x) (dotted). Parameters were chosen by nonlinear least-squares fits to Ivlev’s response functions with aI = 1 and bI = 2 (i.e. aM = 1.14, bM = 0.37, aT = 0.99 and bT = 1.48). 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. Acknowledgement: This research is based on collaboration with T.B, M.H., and B.E. during their graduate studies at McMaster University under the supervision of G.S.K.W. Appendix Appendix A Modified Huang’s Theorem We use a slightly adapted version of Huang’s theorem that better applies to system (2.8) by account- ing for invariance of the first quadrant. Theorem A.1 (Huang, 1988). Consider a system of the form: dx dt = φ(x) (F(x) −π(y)) , (A.1) dy dt = ρ(y)ψ(x). PREDATOR-PREY CHEMOSTAT MODEL: HOLLING TYPE II RESPONSE 347 If (i) ∃ K > x∗ such that F(K) = 0 and (x − K)F(x) < 0 for x 6= K, ∃ 0 < x∗ < K such that ψ(x∗) = 0, i.e. there is a positive equilibrium point (x∗,y∗), (ii) all functions in (A.1) are C1 in the interior of R3+, and F ′(x) is continuous in the interior of R3+, (iii) φ(0) = π(0) = ρ(0) = 0, φ′(x) > 0 and ψ′(x) > 0 for x > 0, ρ′(y) > 0 and π′(y) > 0 for y > 0, and (iv) H(x) = −F ′(x)φ(x)/ψ(x) 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: x′(t) = ξ(x) ( F(x) −γ(y) ) y′(t) = η(y) ( −D + q(x) ) x ≥ 0, y ≥ 0. (B.1) (B.1) 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) = x(1− x K ) q(x) , and q(x) = a tanh(bx). Theorem B.1. Assume there exists an equilibrium (x∗,y∗) with positive components satisfying q(x∗) = D and η(y∗) = F(x∗), and that the following assumptions hold: (i) ∃K > 0 such that F(K) = 0, F(x) > 0 for 0 ≤ x ≤ K, F(x) < 0 for x > K, and F(x) ∈ C1 ( [0,K] ) . (ii) F has at most one local extremum, a local maximum, for x ∈ [0,K] assumed to occur at x = M. (iii) η(y) ∈ C1(R+), η(0) = 0 and η′(y) > 0 for y > 0. (iv) γ(y) ∈ C1(R+), γ(0) = 0 and γ′(y) > 0 for y > 0. (v) ξ(x) ∈ C1(R+), ξ(0) = 0 and ξ′(x) > 0 for x > 0. (vi) q(x) ∈ C1(R+), q(0) = 0, and q′(x) > 0 for x > 0. (vii) F ′(x∗) ≤ 0. Then (x∗,y∗) is globally stable provided there exists a β > 0 such that: max x∈[0,M] − ξ(x)F ′(x)( −D + q(x) ) inf y≥0 η′(y) ≤ β ≤ min x∈[x∗,K] − ξ(x)F ′(x)( −D + q(x) ) sup y≥0 η′(y) . (B.2) 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 348 T. BOLGER, B. EASTMAN, M. HILL, AND G.S.K. WOLKOWICZ function to be h(x,y) = ξ(x)−1η(y)β−1, where β > 0. The resulting divergence is then ∆ ,∇· (hf) = η(y)β−1H(x) ξ(x) (B.3) where f = ( x′ y′ ) and H(x) is given by: H(x) = ξ(x)F ′(x) + β ( −D + q(x) ) η′(y). (B.4) 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 β > 0. Observe that d dx γ−1 ( F(x) ) = F ′(x) γ′ ( γ−1 ( F(x) )). (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 H(x) ≤ 0 ∀x ∈ [0,M], it is then required that β ≥ max x∈[0,M] − ξ(x)F ′(x)( −D + q(x) ) inf y≥0 η′(y) . (B.6) However, β can not be too large either. Since x∗ ≥ M, −D + q(x) ≥ 0 and F ′(x) ≤ 0 for x ∈ [x∗,K], then β ≤ min x∈[x∗,K] − ξ(x)F ′(x)( −D + q(x) ) sup y≥0 η′(y) (B.7) must also hold to guarantee that H(x) ≤ 0 ∀x ∈ [x∗,K]. We do not need to worry about the region [M,x∗] since both F ′(x) ≤ 0 and −D + q(x) ≤ 0 there, so any positive β would work. Thus, if there exists β > 0 that satisfies both (B.6) and (B.7), 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 assump- tions. 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) = x(1− x K ) 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 7→ Dm̄ and b 7→ x∗(m̄−1), where m̄ > 1. Then q(x) = Dm̄x x∗(m̄−1)+x. We first show that these functions satisfy the assumptions of Theorem B.1. The prey nullcline is given by the function: F(x) = (1 − x K ) ( x∗(m̄− 1) + x ) Dm̄ (B.8) which has the second derivative F ′′(x) = − 2 Dm̄K < 0. (B.9) PREDATOR-PREY CHEMOSTAT MODEL: HOLLING TYPE II RESPONSE 349 Since F(x) is concave down for all x ∈ [0,K], it has at most one local extremum, a local maximum, in this region and thus, satisfies assumption (ii). We find the maximum to occur at: M = 1 2 ( K −x∗(m̄− 1) ) . (B.10) It is also clear that F(x) > 0 for all x ∈ [0,K], F(K) = 0, and F(x) < 0 x > K, satisfying assumption (i). Furthermore, ξ′(x) = q′(x) = Dm̄x∗(m̄−1)( x∗(m̄−1)+x )2 > 0 for x > 0, and ξ(0) = q(0) = 0. Thus both (v) and (vi) are also satisfied. Finally, γ(0) = η(0) = 0, and γ′(y) = η′(y) = 1 > 0 for y > 0, satisfying assumptions (iii) and (iv). 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 β(x) = − ξ(x)F ′(x)( −D + q(x) ) η′(y) = − x ( K − 2x−x∗(m̄− 1) ) DK(m̄− 1)(x−x∗) = − 2x(M −x) DK(m̄− 1)(x−x∗) . (B.11) We can then determine where the maximum and minimum of β(x) occur, by finding its critical values. Taking the first derivative, we get: β′(x) = 2(x2 − 2xx∗ + Mx∗) DK(m̄− 1)(x−x∗)2 . (B.12) Note that the sign of β′(x) depends on the numerator. The roots of x2 − 2xx∗ + Mx∗ = 0 are x+,− = x ∗ ± √ x∗(x∗ −M), (B.13) 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]. Lemma B.2. Let β− = β(x−) > 0. Then max x∈[0,M] β(x) = β−. Furthermore, β ′(x) < 0 for x∗ < x < x+. Proof. First consider x−. We proceed using proof by contradiction to show that x− ∈ [0,M]. Suppose that x− > M. Then x∗ − √ x∗(x∗ −M) > M − √ x∗(x∗ −M) > M −x∗ > 0, which is and x− ∈ [0,M], we can conclude that β− is a local maximum, and hence, that max x∈[0,M] β(x) = β−. (B.14) Now consider x+. Clearly x+ > x ∗. Moreover, β′(x) has a removable singularity at x = x∗ and since sign ( lim x→x∗ β′(x∗) ) = x∗(M −x∗) < 0, (B.15) β(x) is decreasing between x∗ and x+. � 350 T. BOLGER, B. EASTMAN, M. HILL, AND G.S.K. WOLKOWICZ Lemma B.3. Let β+ = β(x+) > 0 and βK = β(K) > 0. If x ∗ ≤ K 2 2K−M , then min x∈[x∗,K] β(x) = β+. Otherwise, if x∗ > K 2 2K−M , then min x∈[x∗,K] β(x) = βK. 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, β(K) = 2K2 − 2Kx∗ + Mx∗ = 2K2 + (x∗)2 − m̄(x∗)2 − 3Kx∗ = 2K(K −x∗) + ( −Kx∗ + (x∗)2(1 − m̄) ) = 2K(K −x∗) + x∗ ( (1 − m̄)x∗ −K ) = 2K(K −x∗)2 + x∗ ( − (K − 2M) −K ) = 2K(K −x∗) + 2x∗(M −K) = 2[K(K −x∗) + x∗(M −x∗ + x∗ −K)] = 2[(K −x∗)2 + x∗(M −x∗)] ≥ 0 This means that in this case, min x∈[x∗,K] β(x) = β+. (B.16) 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 x∈[x∗,K] β(x) = βK. (B.17) � Lemma B.4. For our positive quantities β−, β+ and βK, (a) β+ > β− when M ≤ x+ ≤ K (b) βK > β− when x+ > K. Proof. (a) β+ −β− = 8x∗(x∗ −M) DK(m̄− 1) √ x∗(x∗ −M) > 0 (B.18) (b) βK −β− = 2 ( (x∗ −K)2 + x∗(x∗ −M) )√ x∗(x∗ −M) + 4x∗(K −x∗)(x∗ −M) DK(m̄− 1)(K −x∗) √ x∗(x∗ −M) > 0 (B.19) In any case, there exists a positive β between the two quantities when x∗ ∈ (M,K]. � Finally consider when x∗ = M. Then β′(x) = 2(x−x∗)2 DK(m̄−1)(x−x∗)2 > 0 for x 6= x ∗ and furthermore lim x→x∗ β′(x) = 2 DK(m̄−1) > 0, hence, β(x) is increasing on [0,K]. Define β(x∗) = lim x→x∗ β(x) = lim x→x∗ 2x(x−x∗) DK(m̄− 1)(x−x∗) = 2x∗ DK(m̄− 1) . (B.20) Thus, max x∈[0,M] β(x) = β(x∗) = min x∈[x∗,K] β(x). PREDATOR-PREY CHEMOSTAT MODEL: HOLLING TYPE II RESPONSE 351 Remark B.2. Our above examples were both the case that η(y) = y. It is conceivable that other increasing functions will work, for instance η(y) = δy + αy γ+y , where α, δ and γ are all positive. Then sup y≥0 η′(y) = δ + α γ and inf y≥0 η′(y) = δ, both positive, finite values that would not blow up (B.2). 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, S0crit. 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 J = P−1AP = ( α β −β α ) where α± iβ are the eigenvalues of A. Lemma C.1. Let the eigenvector for α + iβ be PRe + iPIm. Then P = ( PRe PIm ) . Proof. We have: P−1P = P−1 ( PRe PIm ) = ( P−1PRe P −1PIm ) . By the identity P−1P = I, this means that P−1PRe = ( 1 0 ) and P−1PIm = ( 0 1 ) . Since PRe + iPIm is the eigenvector for α+iβ, we know A(PRe +iPIm) = (α+iβ)(PRe +iPIm), hence, APRe = αPRe−βPIm and APIm = βPRe + αPIm. Thus, J = P−1AP = ( P−1APRe P −1APIm ) = ( P−1(αPRe −βPIm) P−1(βPRe + αPIm) ) = ( αP−1PRe −βP−1PIm βP−1PRe + αP−1PIm ) = ( α ( 1 0 ) −β ( 0 1 ) β ( 1 0 ) + α ( 0 1 )) = ( α β −β α ) . � Near S0crit, the eigenvalues of (4.4) are: l+,− = (mx∗ + D)F ′(x∗) ± i √ − ( (mx∗ + D)F ′(x∗) )2 + 4q′(x∗)y∗(mx∗ + D) 2 . (C.1) Thus by Lemma C.1, we can construct the transformation matrix P by taking the real and imaginary parts of the eigenvector corresponding to l+, namely P =  F′(x∗)(mx∗+D) 2y∗q′(x∗) √ (mx∗+D) ( −F′(x∗)2(mx∗+D)+4y∗q′(x∗) ) 2y∗q′(x∗) 1 0   . (C.2) 352 T. BOLGER, B. EASTMAN, M. HILL, AND G.S.K. WOLKOWICZ With this P , we can rewrite system (2.8) so that the linear part is in real Jordan canonical form as follows: ( x y )′ = A ( x y ) + H.O.T.’s P−1 ( x y )′ = P−1A ( x y ) + P−1 H.O.T.’s P−1 ( x y )′ = P−1APP−1 ( x y ) + P−1 H.O.T.’s P−1 ( x y )′ = JP−1 ( x y ) + P−1 H.O.T.’s and thus, by letting ( u v ) = P−1 ( x y ) , we transform system (2.8) to: ( u v )′ = J ( u v ) + H.O.T.’s , ( f(u,v) g(u,v) ) . (C.3) In our case, P−1 ( x′ y′ ) =   y ( −D + q(x) ) 2y∗q′(x∗) ( q(x)+mx )( F(x)−y ) −F′(x∗)(mx∗+D)y ( −D+q(x) ) √ (mx∗+D) ( −F′(x∗)2mx∗−F′(x∗)2D+4y∗q′(x∗) )   = (f(x,y) g(x,y) ) , (C.4) where x(u,v), and y(u,v). One can find x and y in terms of u and v by inverting the change of variables, ( x y ) = P ( u v ) =  F′(x∗)(mx∗+D)u+ √ (mx∗+D) ( −F′(x∗)2mx∗−F′(x∗)2D+4y∗q′(x∗) ) v 2y∗q′(x∗) u   . (C.5) Thus, we finally compute the vague attractor condition (4.7) using the chain rule: W ′′′(x∗,y∗,S0crit) = 3π 4|β| ( fuuu + fuvv + guuv + gvvv ) + 3π 4|β|2 ( −fuv(fuu + fvv) + guv(guu + gvv) + fuuguu −fvvgvv ) = (mx∗ + D) ( 2F ′′(x∗)q′(x∗)2 + 2F ′′(x∗)q′(x∗)m + F ′′′(x∗)q′(x∗)D ) y∗q′(x∗)2 − (mx∗ + D) ( q′′(x∗)F ′′(x∗)mx∗ −F ′′′(x∗)q′(x∗)mx∗ + q′′(x∗)F ′′(x∗)D ) y∗q′(x∗)2 . (C.6) Since only the sign of W ′′′ matters, one can factor the positive constant mx ∗+D y∗q′(x∗) , obtaining w = (mx∗ + D)F ′′′(x∗) − (mx∗ + D)F ′′(x∗)q′′(x∗) q′(x∗) + (2q′(x∗) + 2m)F ′′(x∗), (C.7) which is consistent with the expression found in Wolkowicz [38] for, in their notation, p(x) = mx+q(x). At last, substituting the Monod response function yields (4.7). PREDATOR-PREY CHEMOSTAT MODEL: HOLLING TYPE II RESPONSE 353 References [1] C. Aldebert and D. B. Stouffer, Community dynamics and sensitivity to model structure: towards a probabilistic view of process-based model predictions, J. R. Soc, Interface 15 (2018), 20180741. [2] B. S. Attili and S. F. Mallak, Existence of limit cycles in a predator-prey system with a functional response of the form arctan(ax), Commun. Math. Anal. 1 (2006), no. 1, 33–40. [3] I. Bendixson, Sur les courbes définies par des équations différentielles, Acta Math. 24 (1901), no. 1, 1–88. [4] T. Bolger, Chemostat predator-prey model with Ivlev functional response, Master’s thesis, McMaster University, 2017. [5] G. J. Butler, H. I. Freedman, and P. Waltman, Uniformly persistent systems, Proc. Amer. Math. Soc. 96 (1986), no. 3, 425–430. [6] G. J. Butler and G. S. K. Wolkowicz, Predator-mediated competition in the chemostat, J. Math. Biol. 24 (1986), 167–191. [7] K. S. Cheng, Uniqueness of a limit cycle for a predator-prey system, SIAM J. Math. Anal. 12 (1981), 541–548. [8] K. S. Cheng, S.-B. Hsu, and S. S. Lin, Some results on global stability of a predator-prey system, J. Math. Biol. 12 (1981), 115–126. [9] W. A. Coppel, Stability and asymptotic behaviour of differential equations, D. C. Heath, New York, 1965. [10] S. Daoussis, Complicated dynamics of a three species chemostat food chain, Ph.D. thesis, McMaster University, 1997. [11] B. Eastman, Sensitivity to predator response functions in the chemostat, Master’s thesis, McMaster University, 2017. [12] B. Ermentrout, Simulating, analyzing, and animating dynamical systems: A guide to xppaut for researchers and students, SIAM, Philadelphia, 2002. [13] H. I. Freedman and G. S. K. Wolkowicz, Predator–prey systems with group defence: the paradox of enrichment revisited, Bull. Math. Biol. 48 (1986), 493–508. [14] G. F. Fussmann and B. Blasius, Community response to enrichment is highly sensitive to model structure, Biol. Lett. 1 (2005), 9–12. [15] S. R. Hansen and S. P. Hubbell, Single-nutrient microbial competition: Qualitative agreement between experimental and theoretical forecast outcomes., Science 213 (1981), 972–979. [16] G. W. Harrison, Global stability of predator–prey interactions, J. Math. Biol. 8 (1979), 159–171. [17] A. Hastings and Powell T., Chaos in a three species food chain, Ecology 72 (1991), 896–903. [18] M. Hill, A predator-prey model in the chemostat with Holling type II dynamics, Master’s thesis, McMaster University, 2016. [19] J. Hofbauer and W.-H. J. So, Multiple limit cycles for predator-prey models, Math. Biosci. 99 (1990), 71–75. [20] C. S. Holling, The components of predation as revealed by a study of small-mammal predation of the European pine sawfly, Can. Entomol. 91 (1959), 293–320. [21] S.-B. Hsu, On global stability of a predator–prey system, Math Biosci. 39 (1978), 1–10. [22] X.-C. Huang, Uniqueness of limit cycles of generalised Liénard systems and predator-prey systems, Journal of Physics A: Mathematical and General 21 (1988), no. 13, L685–L691. [23] V. S. Ivlev, Experimental ecology of the feeding of fishes, Yale University Press, New Haven, CT, 1961. [24] A. D. Jassby and T. Platt, Mathematical formulation of the relationship between photosynthesis and light for phy- toplankton, Limnol. Oceanogr. 21 (1976), 540–547. [25] R. E. Kooij and A. Zegeling, A predator–prey model with Ivlev’s functional response, J. Math. Anal. Appl. 198 (1996), 473–489. [26] M. Kot, Elements of mathematical ecology, Cambridge University Press, Cambridge, 2001. [27] Y. Kuang and H. I. Freedman, Uniqueness of limit cycles in Gause-type models of predator-prey systems, Math. Biosci. 88 (1988), no. 1, 67–84. [28] L. P. Liou and K. S. Cheng, On the uniqueness of a limit cycle for a predator-prey system, SIAM J. Math. Anal. 19 (1988), no. 4, 867–878. [29] Z. Lu and G. S. K. Wolkowicz, Global dynamics of a mathematical model of competition in the chemostat: General response functions and differential death rates, SIAM J. Appl. Math. 52 (1992), no. 1, 222–233. [30] Maple, (2018), Maplesoft, a division of Waterloo Maple Inc., Waterloo, Ontario. [31] J. E. Marsden and M. McCracken, The hopf bifurcation and its applications, Applied Mathematical Sciences, vol. 19, Springer-Verlag, New York, NY, 1976. [32] MATLAB, (r2019a), The MathWorks Inc., Natick, Massachusetts. [33] J. D. Murray, Mathematical biology: I. an introduction, third ed., Springer-Verlag, New York, 2002. [34] M. L. Rosenzweig and R.H. MacArthur, Graphical representation and stability conditions of predator–prey interac- tions, Am. Nat. 97 (1963), 209–223. 354 T. BOLGER, B. EASTMAN, M. HILL, AND G.S.K. WOLKOWICZ [35] G. Seo and G. S. K. Wolkowicz, Existence of multiple limit cycles in a predator-prey model with arctan(ax) as functional response, Commun. Math. Anal. 18 (2015), 64–68. [36] G. Seo and G. S. K. Wolkowicz, Sensitivity of the dynamics of the general Rosenzweig-MacArthur model to the mathematical form of the functional response: A bifurcation approach, J. Math. Biol. 76 (2018), 1873–1906. [37] H. L. Smith and P. Waltman, The theory of the chemostat: Dynamics of microbial competition, Cambridge University Press, New York, 1995. [38] G. S. K. Wolkowicz, Bifurcation analysis of a predator–prey system involving group defence, SIAM J. Appl. Math. 48 (1988), 592–606. [39] , Successful invasion of a food web in a chemostat, Math. Biosci. 93 (1989), 249–268. Department of Mathematics and Statistics, McMaster University, ON, Canada L8S 4K1 Current address: 425 Parkside Dr, Waterdown ON, Canada L8B0Y6 E-mail address: tedra.bolger@gmail.com Department of Mathematics and Statistics, McMaster University, ON, Canada L8S 4K1 Current address: Department of Applied Mathematics, University of Waterloo, ON, Canada N2L 3G1 E-mail address: b2eastma@uwaterloo.ca Department of Mathematics and Statistics, McMaster University, ON, Canada L8S 4K1 Current address: 404 King St. W Apt. 420 Kitchener, ON, Canada N2G 4Z9 E-mail address: Madeleine 0715@hotmail.com Corresponding author, Department of Mathematics and Statistics, McMaster University, ON, Canada L8S 4K1 E-mail address: wolkowic@mcmaster.ca 1. Introduction 2. The Model 3. The Prey Nullcline 4. Local Analysis 5. Global Analysis 6. Dynamics of the 3D System 7. Discussion Appendix A. Modified Huang's Theorem Appendix B. Extension of Hsu's theorem Appendix C. Analysis of the Hopf bifurcation References