Uniqueness of Limit Cycles for a Limiting Case of the Chemostat: Does It Justify the Use of Logistic Growth Rates?

On infinitesimally short time intervals various processes contributing to population change tend to operate independently so that we can simply add their contributions (Metz and Diekmann, The dynamics of physiologically structured populations, 1986, p. 3). This is one of the cornerstones for differential equations modeling in general. Complicated models for processes interacting in a complex manner may be built up, and not only in population dynamics. The principle holds as long as the various contributions are taken into account exactly. In this paper we discuss commonly used approximations that may lead to non-removable dependency terms potentially affecting the long run qualitative behavior of the involved equations. We prove that these terms do not produce such effects in the simplest and most interesting biological case, but the general case is left open. Our main result is a rather complete analysis of an important limiting case. Once complete knowledge of the qualitative properties of simple models is obtained, it greatly facilitates further studies of more complex models. A consequence of our analysis is that standard methods can be applied. However, the application of those methods is far from straightforward and require non-trivial estimates in order to make them valid for all values of the parameters. We focus on making these proofs as elementary as possible.


Introduction
A number of ecological phenomena can be studied under chemostat conditions, cf.Smith and Waltman [39, p. xi].For instance, a lake ecology may have rivers transferring a limiting resource to the lake and rivers diluting this resource.A phenomenological model containing (1.1) Here s > 0 is the substrate, x > 0 is the prey having the substrate s as its limiting resource and y > 0 is a predator feeding on the prey x.The parameters C > 0, D > 0, a > 0, b > 0, m > 0, A > 0, B > 0, and M > 0 stand for concentration, dilution rate, search rate for the prey, handling time for the prey (cf.Holling [18]), conversion factor for the prey, search rate for the predator, handling time for the predator and conversion factor for the predator, respectively.A logistic approximation is often made when studying (1.1) and related systems, cf.Kooi, Boer, and Kooijman [23] and references therein.More precisely, consider the function It satisfies Ḣ = −DH meaning that the surface H = 0 is asymptotically invariant for (1.1).A study of the system on this surface allows for reducing it to a planar predator-prey system as follows ẋ = ax(mC Such reductions can be made rigorously under certain conditions, see Smith and Waltman [39,Appendix F].See also Thieme [41] for a nice example describing how such reductions may break down.If b = 0, then the growth function of (1.2) is non-logistic and given by h(x) = ax(mC − x) in the absence of predators.The second equality of (1.3) explicitly clarifies its vertical and oblique asymptotes, respectively.The growth function is concave down at the left-hand side of its vertical asymptote.Therefore, it is unimodal (has one local maximum) on the interval The last inequality is identically true.Thus, no qualitative differences can occur in the single species case when doing logistic approximations.The above expressions should be compared to the corresponding expressions in the logistic case.When b → 0 the function is unimodal on [0, (amC − D)/a] when D < amC.The growth rate and carrying capacity parameters are thus, given by r = amC − D and K = (amC − D)/a, respectively, cf.Fitzsimmons, Schoustra, Kerr, and Kassen [8], and are strongly related.Equivalents of the predator-prey system (1.2) have been studied in Smith and Waltman [39, Sections 3.2-3.6]and Kuang [25].The results were that local stability implies global stability and that uniqueness of limit cycles was proved for a certain range of parameters.However, it is still not known whether the limit cycle is unique for all parameters of the system.
If the handling time for the prey is small b ≈ 0, the denominator in the first term of the first equation in (1.2) vanishes and a logistic approximation is often made for the growth rate of the prey feeding on the substrate.We get (1.4) In this context, the system (1.5) has often been used as an approximation of the system (1.4).Note that the coupling term axy M is set to zero in this procedure, too.This is entirely in concordance with well-established principles in differential equation modeling.As soon as logistic growth is accepted for the prey-population x, then one is expected to add the contributions of the predation process to the resulting equations.The equations (1.5) contain the result of such a procedure.The longrun qualitative dynamical properties including possibilities for limit cycle bifurcations of (1.5) are well-known trough a number of well-established results (cf.Cheng, Hsu, and Lin [3], Liou and Cheng [32] and Kuang and Freedman [26]).The objective of this paper is to elucidate the influence of the above coupling term on those results.The numerical part of the work by Kooi, Boer, and Kooijman [23] reported differences in fixed point bifurcations for still longer food-chains as consequences of similar simplifications.The methods used in this paper are the Lyapunov functions introduced by Ardito and Ricciardi [1] and LaSalle's [28] theorem, cf.Hale [14, p. 317], Hirsch, Smale, and Devaney [17, p. 199], and Wiggins [42,Section 8.3].We use the Kuang and Freedman [26] version of Zhang's [46] theorem (see also Cherkas and Zhilevich [4]) for proving uniqueness of limit cycles.The result is that the qualitative properties of the bifurcation diagram are robust with respect to the coupling term mentioned above.
However, it is not straightforward to do the estimates required for application of these methods in particular cases.We have therefore limited the presentation here to the simplest and most interesting case (logistic growth combined with Holling [18] specialist predator functional response).In the presentation of the proofs we have focused on finishing the relevant inequalities to the best forms whereas chains of equalities that involve extensive but elementary algebra have been given low priority.Our impression is that equalities are more straightforward to check than inequalities.

Reparametrization
We begin our study of the qualitative differences between the systems (1.5) and (1.4) by reparametrizing the systems in order to eliminate some of the parameters involved (cf.Keener [22]).We put and the system under consideration takes the form (2.1) The known system (1.5) corresponds to (2.1) with ǫ = 0 and seven parameters have been reduced to four.We assume in the rest of the paper that all parameters and variables are non-negative.Non-negative variables can be assumed, since by uniqueness of solutions, no solutions can intersect the solutions at ξ = 0 and η = 0.The parameters C and D were those used as bifurcation parameters in Kooi, Boer, and Kooijman [23].

Isocline form and equilibria
Most of the theorems and results regarding the system (2.1) are formulated for its isocline form ξ, η ≥ 0. The involved functions are required to meet some conditions.We specify them as follows: (A-I) f , ψ, and Condition (A-I) grants the local existence and uniqueness of solutions.In fact, only Lipschitz conditions are needed for this conclusion.However, we may need continuous derivatives of higher orders for some results below.These conditions are valid anyway for the specific system (2.1).Indeed, for (2.1), we have Proof.By uniqueness of solutions, the solutions of (3.1) cannot intersect the solutions at η = 0 and ξ = 0. Therefore, the solutions remain positive.For boundedness, note first that all solutions enter the region 0 < ξ < 1 by (A-II)-(A-III) and that we have η < 0 for 0 < ξ < λ by (A-IV).Next consider the function (cf.Lindström) [30]) The last integral diverges as η → ∞.Therefore, all level curves of this function intersect the line ξ = λ.The derivative along the solutions of (3.1) for the above function is Since F(ξ) + ψ(ξ) is bounded on λ < ξ < 1, this quantity is negative for η large enough.That is, select κ such that V < 0 for V(ξ, η) > κ and we have a trapping region ensuring bounded solutions according to Figure 5.1(a).
Remark 3.2.Because all positive solutions enter the region 0 < ξ < 1 by the above argument, we usually check the conditions of the theorems alluded to below in that interval.
We make a summary of the local bifurcation results.
Proof.Existence and number of equilibria follow from (A-II)-(A-IV).The Jacobian matrix of (3.1) takes the form meaning that the Jacobian matrices evaluated at the first two equilibria take diagonal or triangular forms.The stability properties of the last equilibrium are usually referred to as the Rosenzweig-MacArthur [38] criterion.It shows that the properties of F ′ are important.Returning to (3.2), this function can be written in two forms where the latter one makes use of the zeroes of the numerator.We note that ξ − < 0.Moreover, ξ + is a strictly increasing function of β that takes the value zero when β = 1 + ǫ.Thus, the numerator of (3.3) has a positive zero if β > 1 + ǫ, otherwise not.Theorem 3.3 implies that λ can be used as a bifurcation parameter and λ ≥ 1 implies (1, 0) locally stable, max(0, ξ + ) < λ ≤ 1 implies (λ, F(λ)) to be locally stable, and implies that all equilibria are unstable.Standard analysis gives now that This relates our results to (2.1) with ǫ = 0.The coupling term in the system (2.1) therefore stabilizes the system locally.If the term containing ǫ > 0 is present, then λ must be smaller in order to ensure possibilities for all non-negative equilibria to be unstable.The local bifurcation routes are qualitatively the same (as λ decreases, possible destabilization occurs), only quantitative differences occur (the parameter λ must be still smaller to ensure destabilization to take place as ǫ > 0).This phenomenon is usually called the "Paradox of enrichment" in the biological literature, cf.Rosenzweig [37].We conclude this section by proving the following lemma.
Proof.By Theorem 3.1 all solutions remain bounded and positive.The unstable manifold of (1,0) can therefore not intersect the stable manifold of itself or the origin.The Butler-McGehee theorem (Smith and Waltman [39, Section 1.3]) implies that no ω-limit set can contain (1,0) or the origin.By the Poincaré-Bendixson theorem the ω-limit set is either the fixed point (λ, F(λ)) or a limit cycle that must surround (λ, F(λ)) by index theory, see e.g.Jordan and Smith [21, Chapter 3] or Dumortier, Llibre, and Artés [7, Chapter 6].In both cases, trajectories starting in the positive cone have some part of the ray ξ = λ, 0 < η ≤ F(λ) in their limit set.

Global stability
We now continue with global stability.Our results follow from the results in Smith and Waltman [39, Section 3.5] in a limiting case but our method is different (and simpler) anyway.
For θ ≥ 0, we introduce the range of Lyapunov functions see Ardito and Ricciardi [1] and Lindström [31].We introduce and deduce that the total derivative of W θ along the solutions of (3.1) is given by Therefore, the sign condition implies that LaSalle's [28] invariance principle can be applied.We also note that W θ is a first integral of 3) are compared to each other.For θ = 0, all solutions (4.3) in the positive quadrant are closed orbits, and hence the level curves of (4.1), see Figure 5.1(c).For θ > 0, (4.3) has non-closed solutions in the positive cone.However, all solutions intersecting the ray ξ = λ, 0 < η ≤ F(λ) are closed orbits.See Figure 5.1(d) and cf.Lemma 3.4.We point out that this was not mentioned in the original paper by Ardito and Ricciardi [1].
If F(ξ) − Fθ (ξ) decreases, then the sign-condition (4.2) is clearly satisfied.We differentiate and obtain the requirement (cmp.e.g.(3.2) and (3.3)) we can choose θµ = 2ξ + β > 0 and the above inequality is equivalent to It is straightforward to check the equality above and it can be derived using a Taylor expansion around ξ = ξ + .The inequality holds since λ ≥ ξ + and We can thus, formulate the following theorem.

Uniqueness of limit cycles
If F ′ (λ) > 0, then there exists at least one limit cycle for (3.1) with (A-I)-(A-IV) due to the Poincaré-Bendixson theorem and Theorem 3.1.If ǫ = 0, then this limit cycle is unique by Theorem 4.3 in Kuang and Freedman [26] for (3.1) with (3.2).The theorem basically shows that Zhang's [46] theorem for Liénard equations can be applied to the predator-prey system (3.1) after a sequence of coordinate transformations.The conditions are stronger than (A-I)-(A-IV) in order to validate the coordinate transformations and uniqueness.Alternative methods for obtaining reminiscent conditions are given by Cheng [2], Liou and Cheng [33], and Ding [6].
Other methods for proving uniqueness of limit cycles for Liénard systems have since then been made available through various coordinate transformations, including Sansone's method (Xiao and Zhang [43]).A nice survey on how to choose between various methods was given by Hasík [15].However, the lack of examples carried out carefully in this paper demonstrates a clear difficulty.Although a large variety of methods exist, all of them tend to lead to formidable algebraic estimations when checking the essential criteria as they are applied to important special cases of biological relevance.This holds especially as one tries to verify the conditions for all possible parameter values.In our case the essential condition is verifying constant sign of a quartic polynomial for all possible values of the three involved parameters.We encountered a similar problem already when proving global stability, but then the polynomial was cubic.
We shall use Kuang and Freedman's [26] theorem to prove that the limit cycle is unique for ∀ǫ > 0. It is difficult to verify its crucial condition for all relevant parameter values.Using (3.3) it takes the form for ξ = λ in 0 ≤ ξ ≤ 1.The criterion states that the divergence integrated around an assumed outer limit cycle is smaller than if it was integrated around an assumed inner limit cycle.The remaining conditions are satisfied for our selection of involved functions (3.2).It is sufficient to check that the numerator of (5.1) is negative.In our case, this condition is equivalent to that after Taylor expansion around ξ = λ is equivalent to (5.3) Straightforward computations verify the expressions above.The last term of (5.3) is negative definite since We remark that the last inequality holds because It remains to prove that the three first terms of (5.3) form a negative contribution.This quadratic expression has negative leading term and may be estimated from above by lines tangent to it at ξ = 0 and ξ = λ.These two lines are given by 2) and (3.3)) and (from (5.3)), respectively.We evaluate L 0 and L λ at points allowing for simple computations as follows Now, since λ < ξ + , the three first terms in (5.3) give a negative total contribution.Therefore, the limit cycle existing for 0 < λ < ξ + is unique.We have thus, proved the following theorem.(c) For θ = 0, all solutions of (4.3) are closed orbits. (d) For θ > 0 there may be non-closed solutions of (4.3), but all solutions passing the cross-section ξ = λ, 0 < η ≤ F(λ) are closed.The ×-marks denote two saddles in a heteroclinic loop.

Discussion
We have considered the system (2.1) in this paper.This system was derived from a threetrophic level (resource-prey-predator) chemostat model possessing an asymptotically invariant plane that allows reduction to two trophic levels (prey and predator).If the prey is unsaturated, logistic growth follows but a coupling term affecting the predators functional response cannot be removed.We prove that λ can be used as a bifurcation parameter and that (1, 0) is globally stable for λ ≥ 1, (λ, F(λ)) is globally stable for ξ + ≤ λ ≤ 1, and that a unique stable limit cycle exists for 0 < λ < ξ + .The bifurcations occurring at λ = The analysis shows it to be far from trivial in general to deduce what kind of effects such coupling terms might have on the existence of possible limit cycles.Obvious cases that are left open are saturated prey (giving rise to a non-logistic growth function that depends on the predator), generalist predator functional responses (Hassell [16,p. 39], Britton [45, p. 61] or the discontinuous variants in Krebs and Davies [24, p. 61]), and group defence cases cf.Geritz and Gyllenberg [9,10].In several cases, results granting global stability or uniqueness of limit cycles exist if logistic growth is assumed and coupling terms are neglected, see e.g.Lindström [29], Hwang [19,20], and Qiu and Xiao [36].Some authors are, however, critical to the foundation for some of the investigated functional responses, see e.g.Oksanen, Moen, and Lundberg [35], Diehl, Lundberg, Gardfjell, Oksanen, and Persson [5], and Gleeson [11].
We remark that examples demonstrating complicated limit cycles bifurcations have been created without considering the dependency terms considered here, see e.g.González-Olivarez and Rojas-Palma [12].It is therefore not evident that uniqueness of limit cycles is granted in predator-prey systems when they exist or that local stability implies global stability.Since such results are not granted even when the dependency terms are neglected, makes it even more complicated to elucidate the long-term dynamical effects of these terms in the general case.
Our analysis shows that a detailed investigation of many of these cases require further development of the methods, but this is usually difficult to do as long as a reasonable set of well-known special cases does not exist.Despite that many of the available methods provide possibilities for analyzing the different cases by elementary calculus methods, the arising expressions usually give rise to extensive algebraic manipulations when doing the required estimates.In this case we solved many of these problems through a careful selection of the quantities computed.
such a situation is given by ṡ = CD − Ds − axs
Figure 5.1(b) demonstrates the techniques used for proving global asymptotic stability in the case λ ≥ 1.The stability properties of the last equilibrium follows from the trace-determinant criterion, see e.g.Hirsch, Smale, and Devaney [17, Section 4.1] or Strang [40, p. 270-271].
[13,nd at λ = ξ + are a transcritical bifurcation, and a supercritical Hopf bifurcation, respectively, see e.g.Guckenheimer and Holmes[13, Section 3.4]orKuznetsov [27, Section 3.4].The result was that this destabilizing bifurcation pattern remains identical to the well-known case without coupling term, but that smaller values of λ are needed for destabilization.