Global dynamics of the chemostat with different removal rates and variable yields.

. In this paper, we consider a competition model between n species in a chemostat including both monotone and non-monotone growth functions, distinct removal rates and variable yields. We show that only the species with the lowest break-even concentration survives, provided that additional technical conditions on the growth functions and yields are satisﬁed. We construct a Lyapunov function which reduces to the Lyapunov function used by S. B. Hsu [SIAM J. Appl. Math., 34 (1978), pp. 760-763] in the Monod case when the growth functions are of Michaelis-Menten type and the yields are constant. Various applications are given including linear, quadratic and cubic yields.

1. Introduction. In this paper we study the global dynamics of the following model of the chemostat in which n populations of microorganisms compete for a single growth-limiting substrate: where S(0) ≥ 0 and x i (0) > 0, i = 1 · · · n and S 0 , D and D i are positive constants. In these equations, S(t) denotes the concentration of the substrate at time t; x i (t) denotes the concentration of the ith population of microorganisms at time t; f i (S) represents the uptake rate of substrate of the ith population; p i (S) represents the per-capita growth rate of the ith population and so the function y i (S), defined by y i (S) = pi(S) fi(S) is the growth yield; S 0 and D denote, respectively, the concentration of substrate in the feed bottle and the flow rate of the chemostat; each D i represents the removal rate of the ith population. For general background on model (1.1), in the constant yield case y i (S) = Y i , the reader is referred to the monograph of Smith and Waltman [19].
The global analysis of this model was considered by Hsu, Hubbell and Waltman [6], in the Monod case [15] when the growth functions are of Michaelis-Menten form, (1.2) and the yields are constant y i (S) = Y i , and D i = D for i = 1 · · · n. The authors showed that only the species with the lowest break-even concentration survives. Thus the competitive exclusion principle (CEP) holds: only one species survives, namely the species which makes optimal use of the resources. Hsu [5] applied a Lyapunov argument to give a simple and elegant proof of the result in [6] for the case of different removal rates D i . The Lyapunov function V H used by Hsu is Yi ai ai−Di , i = 1 · · · n, x * 1 = DY 1 and λ 1 = b1D1 a1−D1 is the lowest break-even concentration of the species.
Wolkowicz and Lu [21] extended the results of [5] by allowing more general growth functions. These authors used the Lyapunov function where α i , i = 2 · · · n are positive constants to be determined. They identified a large class of growth functions, including many prototypes of growth functions often found in the literature, where the constant α i in (1.4) can always be found. Despite the fact the α i cannot be found for all growth functions, the work of Wolkowicz and Lu [21] represents a major step in the extension of the result of Hsu [5] to general growth functions.
The CEP has also been proved under a variety of hypotheses by Armstrong and McGehee [2], Butler and Wolkowicz [3], Wolkowicz and Xia [22] and Li [9]. The hypotheses used in [2,3,6,5,9,21,22] are summarized in Table 1 of [8]. However, the problem is not yet completely solved: the CEP holds for a large class of growth functions but an important open question remains: is the CEP true assuming only that the f i are monotone with no restriction on the D i ? This major open problem remains unresolved, see in particular [8,20]. For other studies and complements on the use of Lyapunov techniques in the chemostat, see [4,10,12,13,14].
The variable yield case was considered, for n = 1, 2 by Pilyugin and Waltman [16], with a particular interest to linear and quadratic yields, and by Huang, Zhu and Chang [7]. The model (1.1), with variable yields, was considered by Arino, Pilyugin and Wolkowicz [1]. For biological motivations concerning the dependence of the yields on the substrate, see [1,16] and the references therein.
Notice that, in the case when the growth functions are of Michaelis-Menten form (1.2), the Lyapunov function (1.4) does not reduce to the Lyapunov function (1.3). Our aim in this paper is to extend the Lyapunov function (1.3) of Hsu [5] to the chemostat with a more general class of growth functions and variable yields. Our Lyapunov function is given by where α i , i = 2 · · · n are positive constants to be determined. This Lyapunov function is just a multiple of the Lyapunov function (1.3) that Hsu used in [5] in the Monod case, see Section 3.1. It is also a multiple of the one used in [22] Page 1039 or [20] Section 3.3, in the case of one species, general growth function and constant yield, see Section 3.2.
The paper is organized as follows. In section 2 we prove our main result (see Theorem 2.2) and we compare it with the main result in [21] (see Theorem 2.3), where the yields are assumed to be constant. It should be noticed that, in the case hal-00418676, version 2 -11 May 2010 when the yields are constant, our result follows from the result in [21]. Actually, both theorems 2.2 and 2.3 are corollaries of a more general result, which is valid in the case when the yields are variable [18]. In Section 3.1 we consider the Monod model with constant yields. In Section 3.2 we consider the one species case and we show that our Lyapunov function can be used to obtain the same result as in [1]. In Section 3.3 we show that for the Monod model with constant yields replaced by either linear or quadratic functions of S, under certain additional technical assumptions, the CEP still holds (see Corollary 3.1). In Section 3.4 we consider the model of Pilyugin and Waltman [16] which was used to demonstrate that a periodic orbit was possible in the case of variable yield model. In this model, with two species, where one yield is constant and the other is cubic in S, we show that our Lyapunov function can be used to prove that for some values of the parameters the CEP holds (see Corollary 3.2). In Section 3.5 we identify a class of growth functions, including Lotka-Volterra and Michaelis-Menten growth functions where our Lyapunov function works. Concluding remarks are given in Section 4.
2. Global asymptotic stability. We make the following assumptions on the functions p i and f i : • p i (0) = f i (0) = 0 and for all S > 0, p i (S) > 0 and f i (S) > 0. Following Butler and Wolkowicz [3], we make the following assumptions on the form of the growth functions p i : there exist positive extended real numbers λ i and µ i with It is known (see Theorem 4.1 [1]) that the non-negative cone is invariant under the flow of (1.1) and all solutions are defined and remain bounded for all t ≥ 0. System (1.1) can have many equilibria: the washout equilibrium E 0 = (S 0 , 0, · · · , 0), which is locally exponentially stable if and only if for all i = 1 · · · n, S 0 / ∈ [λ i , µ i ] and the equilibria E * i and E * * i where all components of E * i and E * * i vanish except for the first and the (i + 1)th, which are The equilibrium E * i lies in the non-negative cone if and only if λ i ≤ S 0 . If λ i < λ j for all i = j and F ′ i (λ i ) < 0 then it is locally exponentially stable. It coalesces with E 0 when λ i = S 0 . The equilibrium E * * i lies in the non-negative cone if and only if µ i ≤ S 0 and is locally exponentially unstable if it exists. Its coalesces with E 0 when µ i = S 0 . Besides these equilibria, the system (1.1) can have a continuous set of non-isolated equilibria in the non-generic cases where two or more of the break-even concentrations are equal. In what follows we assume, that λ 1 < λ 2 ≤ · · · ≤ λ n , and λ 1 < S 0 < µ 1 . Hence E 0 is locally exponentially unstable and the equilibrium , lies in the non-negative cone. It is locally exponentially stable if and only if F ′ 1 (λ 1 ) < 0. We consider the global asymptotic stability of E * 1 .

hal-00418676, version 2 -11 May 2010
Before presenting the results, we need the following lemma, Lemma 2.1. The solutions S(t), x i (t), i = 1 · · · n of (1.1) with positive initial conditions are positive and bounded, and if λ i < S 0 < µ i for some i = 1 · · · n, then S(t) < S 0 for all sufficiently large t.
Proof. The proof is similar to the proof of Lemma 2.1 in [21] obtained for the model (1.1) in the case where the yields are constant. We have the following result.
Then the equilibrium E * 1 is globally asymptotically stable for system (1.1) with respect to the interior of the positive cone.
Proof. From Lemma 2.1 it follows that there is no loss of generality in restricting our attention to 0 ≤ S < S 0 . Consider the function V = V (S, x 1 , · · · , x n ) given by (1.5), where α i are the positive constants satisfying (2.1). The function V is continuously differentiable in the positive cone and positive except at the point E * 1 , where it is equal to 0. The derivative of V along the trajectories of (1.1) is given by . First, note that, using hypotheses 1 and 3, the first term of the above sum is always non-positive for 0 < S < S 0 and equals 0 for Thus θ i (S) < 0 for every S ∈]0, S 0 [, provided that the numbers α i satisfy (2.1). Hence V ′ ≤ 0 and V ′ = 0 if and only if S = λ 1 and x i = 0 for i = 2 · · · n. By the Krasovskii-LaSalle extension Theorem, the ω-limit set of the trajectory is E * 1 . In the case when the yields are constant, y i (S) = Y i , (1.1) takes the form Using the Lyapunov function (1.4), Wolkowicz and Lu (see Theorem 2.3 in [21]) proved the following result Theorem 2.3. Assume that 1. λ 1 < λ 2 ≤ · · · ≤ λ n , and Then the equilibrium E * 1 is globally asymptotically stable for system (2.2) with respect to the interior of the positive cone. It should be noticed that, in the case when the yields are constant, conditions (2.3) are consequences of conditions (2.1) in Theorem 2.2. Indeed, we have Thus, hypotheses 2 and 3 of Theorem 2.2 imply hypothesis 2 of Theorem 2.3. Hence, in the case when the yields are constant, Theorem 2.2 follows from Theorem 2.3. It is of interest to identify classes of growth functions where conditions (2.1) are satisfied, and hence Theorem 2.2 can be applied. We give below a result which will be used in the following section to verify easily that conditions (2.1) are satisfied. This proposition is similar to Corollary 2.4 in [21].

hal-00418676, version 2 -11 May 2010
We consider the case where, for all i = 1 · · · n, a i > D i . In that case: By Proposition 2.4, the conditions (2.1) are satisfied. Since the first derivative of the function F 1 (S) is negative. Hence, hypothesis 3 in Theorem 2.2 is satisfied. The global stability of the equilibrium E * 1 of (3.1) follows from Theorem 2.2. This result was obtained by Hsu [5], using the Lyapunov function (1.3). Notice that, in this case, the Lyapunov function (1.5) is simply V = a1−D1 a1Y1 V H where V H is the Lyapunov function (1.3) used by Hsu [5].
If λ 1 < S 0 < µ 1 and hypothesis 3 in Theorem 2.2 is satisfied then the equilibrium is globally asymptotically stable with respect to the interior of the positive quadrant. This result follows from Theorem 2.2 since in the case where n = 1 the condition (2.1) is obviously satisfied. The global asymptotic stability of E * 1 was obtained previously by Arino, Pilyugin and Wolkowicz [1]. These authors used the following Lyapunov function They proved (see [1], Theorem 2.11) that if 1 − f1(S)(S 0 −λ1) f1(λ1)(S 0 −S) has exactly one sign change for S ∈ (0, S 0 ) then E * 1 is globally asymptotically stable. The condition on the change of sign is equivalent to hypothesis 3 in Theorem 2.2. Notice that the Lyapunov function we obtain is not proportional to the Lyapunov function V AP W considered in [1]. However, in the case when the yields is constant the global asymptotic stability of the equilibrium E * 1 of the system

hal-00418676, version 2 -11 May 2010
can be obtained, using the following Lyapunov function (see [22] page 1039 or [20], Section 3.3) In this case we simply have V = Y 1 V W BL , where V is our Lyapunov function (1.5).
3.3. Michaelis-Menten growth functions and linear or quadratic yields. Consider the particular case of (1.1), where the growth functions p i (S) are given by (1.2), and the yields y i (S) = p i (S)/f i (S) are linear where Y i > 0 and c i ≥ 0. System (1.1) takes the form Corollary 3.1. Consider system (3.5) where the yields are given by (3.3) or (3.4). Assume that 1. λ 1 < λ 2 ≤ · · · ≤ λ n and Then the equilibrium E * 1 is globally asymptotically stable for (3.5) with respect to the interior of the positive cone.

hal-00418676, version 2 -11 May 2010
Under hypothesis 2 there exists α i satisfying (2.1). The result follows by Theorem 2.2. For quadratic yields (3.4) we have Thus Next, the proof is mutatis mutandis the same as the proof given above for the case of linear yields (3.3). This result contains as a particular case the result of Hsu [5] which corresponds to the case where the yields are constant. Indeed, for constant yields c i = 0, so that hypotheses 2 and 3 in Corollary 3.1 are satisfied. Remark. For linear or quadratic yields the function F 1 (S) is not monotone in general on the interval ]0, S 0 [, and it is not easy to give a condition on the parameters for which hypothesis 3 in Corollary 3.1 holds. However, in each example, the graphical depiction of this hypothesis is very simple as shown in Fig. 3.2. 3.4. Pilyugin-Waltman's example. This system was given in [16] as a model of the competition in the chemostat exhibiting limit cycles. The existence of the limit cycles is a consequence of the variable yield in the model. The model takes the form x2 120 (3.6) In their study Pilyugin and Waltman [16] fixed c = 50 and considered m 2 as a bifurcation parameter. They showed that for m 2 ≥ 9.85 the system exhibits sustained oscillations. In this section we fix m 2 = 10 and we consider c ≥ 0 as a bifurcation parameter. In this case we have

hal-00418676, version 2 -11 May 2010
Straightforward computations lead to the formula    Fig. 3.4, left). The function g 2 (S) is defined by For c ≥ 0, the function w 2 (S) is non-decreasing. By Proposition 2.4, the condition (2.1) with i = 2 holds (see Fig. 3.4, right), and the result follows from Theorem 2.2.
Pilyugin and Waltman showed by numerical simulations that their system exhibits limit cycles in the case where c = 50 and m 2 ≥ 9.85 (see Fig. 4 in [16]). The example was revisited by Huang, Zhu and Chang [7] who claimed that the limit cycle of the system should remain only on the face x 2 = 0 (see [7], Remark 2). We do not agree with this claim. We performed ourselves numerical simulations and actually the limit cycle is contained within the positive cone as shown in Fig. 4 in [16] and not in the face x 2 = 0 as claimed in [7]. Huang, Zhu and Chang [7] made a simple modification by replacing 2S/(0.7 + S) with 2S/(0.71 + S) in (3.6) and obtained an example exhibiting competitive exclusion. The model takes the form x2 120 It is claimed, without proof, in [7] that the equilibrium E * 1 is globally asymptotically stable. Hypothesis 3 in Theorem 2.2 is not satisfied (see Fig. 3.5, left) and we cannot prove the global asymptotic stability of E * 1 . However an explanation of the high sensitivity when 0.7 is replaced by 0.71 is easy to find. Actually the plots of the function F 1 (S) in the case of (3.6), where c = 50 and (3.7) are very similar (see Fig. 3.5, left), but a magnification of the neighborhood of the value S = λ 1 shows the differences (see Fig. 3.5, center and right). In (3.6), F ′ 1 (λ 1 ) > 0. Hence the equilibrium E * 1 is locally exponentially unstable. In (3.7), F ′ 1 (λ 1 ) < 0. Hence the equilibrium E * 1 is locally exponentially stable.

Further applications.
In this section we describe a class of growth functions p i (S) and yields y i (S) for which constants α i satisfying (2.1) exist and hence Theorem 2.2 can be applied. It is convenient to use the notation Remark. We can take any functions P i (S) that are positive for 0 < S ≤ S 0 and satisfy P i (0) = 0 and use the righthand side of formulas (3.8) to define the functions p i (S). The function P i (S) must satisfy the condition P i (S) > S − λ i , so that le denominator in p i (S) remains positive. If we find a class of yield functions y i (S) such that the conditions (2.4) hold, where the functions w i (S), considered in Proposition 2.4, are given by then we can use Proposition 2.4 to obtain the global asymptotic stability of the equilibrium E * 1 . The Holling type II (Michaelis-Menten or Monod) growth functions correspond to the choice P i (S) = m i and m i > 1. The Holling type III (or sigmoidal) growth functions correspond to the choice P i (S) = (ai+λi)(bi+λi)S 2 (ai+bi)λiS+aibi(S+λi) . Here µ i = +∞. The prototype for a non-monotone growth function corresponds to the choice P i (S) = (ai+λi)(bi+λi) aibi−λiS . Here µ i = aibi λi . The growth functions (3.9-3.11) were considered by Wolkowicz and Lu [21] who indicated for each combination of them that it is always possible to find appropriate constants α W L i satisfying the criterion (2.3).
Hereafter we define two new classes of functions, which are not considered in the literature, for which our results apply. A class of monotone growth functions of the form (3.8) is obtained with P i (S) = α i S + αS 1+βS , where α > 0, β > 0 and α i ≥ 1. In this case we have A class of non-monotone growth functions of the form (3.8) is obtained with P i (S) = α i S 2 and α i > 1 4λi . In this case we have 4. Discussion. In this paper we considered a mathematical model (1.1) of n species of microorganisms in competition in a chemostat for a single resource. The model incorporates both monotone and non-monotone growth functions, distinct removal rates and variable yields. We demonstrated that the CEP holds for a large class of growth functions and yields.
In the case where the yields are constant, it is known [3] that the CEP holds provided that D i = D for all i, the set Q = i∈N ]λ i , µ i [ is connected, and S 0 ∈ Q, where N = {i : λ i < S 0 }. Wolkowicz and Lu [21] conjectured that this result can be extended to the case of different removal rates. Under hypothesis 1 in Theorem 2.2, it is clear that the set Q is connected, and S 0 ∈ Q. The condition λ 1 < λ i for i = 1 in hypothesis 1 can be stated without loss of generality, by labelling the populations such that the index i = 1 corresponds to the lowest break-even concentration, but the condition λ 1 < S 0 < µ 1 in hypothesis 1 cannot be stated without loss of generality. If µ 1 < S 0 , it is not possible to show the CEP by the methods that we used. To the best of our knowledge, in the case of different removal rates and non-monotone growth functions, the CEP has been proved only under the assumption S 0 < µ 1 [9,21,22]. However, Rapaport and Harmand [17] considered the case of two populations and proposed conditions on the growth functions such that the CEP holds under the condition µ 1 < S 0 . It should be interesting to extend their methods to more general cases. We leave this problem for future investigations.
In the case of constant yields, numerical simulations of model (1.1) have only displayed competitive exclusion. Our results concern also the case of variable yields, for which it is known [1,7,16] that more exotic dynamical behaviors, including limit cycles and chaos, are possible. Thus in the case of variable yields, it is of great importance to have criteria ensuring the global convergence to an equilibrium with at most one surviving species. Under certain technical restrictions, we extended the result of Hsu [5] to the case of linear or quadratic yields.
Our proof relies on the construction of non-strict Lyapunov functions, i.e. Lyapunov functions whose derivatives along the trajectories are non-positive. We conjecture that the strictification techniques of Chapter 5 of [11] can be used to construct strict Lyapunov functions, i.e. Lyapunov functions whose derivative along the trajectories are definite negative, which next can be used to establish some robustness properties. This can be the subject of further research.