STABILITY AND ROBUSTNESS ANALYSIS FOR A MULTISPECIES CHEMOSTAT MODEL WITH DELAYS IN THE GROWTH RATES AND UNCERTAINTIES

. We study a chemostat model with an arbitrary number of competing species, one substrate, and constant dilution rates. We allow delays in the growth rates and additive uncertainties. Using constant inputs of certain species, we derive bounds on the sizes of the delays that ensure asymptotic stability of an equilibrium when the uncertainties are zero, which can allow persistence of multiple species. Under delays and uncertainties, we provide bounds on the delays and on the uncertainties that ensure input-to-state stability with respect to uncertainties.


Introduction.
1.1. Preliminaries. This paper continues our work (which we began in [8,9,22,25,26,27,28,30,39]) on control and other methods to ensure desired asymptotic behaviors in chemostat models, such as the coexistence of multiple competing species, convergence to equilibria, or input delay compensation. Our work is strongly motivated by the ubiquity of chemostats in biological and engineering settings that are of compelling ongoing interest. The chemostat is a laboratory device and a mathematical model for the continuous culture of microorganisms. It was introduced primarily in the works [34] of Monod and [37] of Novick and Szilard from 1950. In the past few decades, chemostat models have been studied extensively, because of their role in biotechnology, ecology, and microbiology as ideal representations for modeling cell or microorganism growth, natural environments such as lakes, and wastewater treatment processes [5]; see [42] for an overview of the chemostat literature.
The classical model of competition in the chemostat is described by the system   where n ≥ 2 microbial species (whose concentrations are denoted by x 1 , . . . , x n ) are in competition for a nutrient with concentration s [42]. The positive constants s in and D are called the input nutrient concentration and the dilution rate, respectively, and Y i is a positive yield constant related to the conversion rate of the substrate into new biomass for each i. The continuously differentiable functions µ i for i = 1, . . . , n are strictly increasing, satisfy µ i (0) = 0, and describe the consumption of the nutrient by species i. The model also assumes that the growth of the ith species is proportional to the consumption of the nutrient.
It is well known (e.g., from [12,42]) that if the preceding conditions hold, and if which is the competitive exclusion principle [12,15,33,35]. The constant µ −1 i (D) is called the break-even concentration for the ith species, namely the minimal nutrient concentration that ensures a positive growth for the ith species. Condition (2) says that the species can be ordered by their competitive ability, which is determined by their break-even concentrations. Condition (3) implies that only the nth species persists, because it only requires the lowest concentration of the nutrient to have positive growth.
However, it is commonly observed in experiments that multiple competing species can persist in chemostats with one limiting substrate. Numerous methods and theories have been developed to generate or explain coexistence in chemostats, such as crowding effects, feedback controls in which the inputs are functions of the state variables instead of being constant (e.g., in [8,27]), flocculation [11], heterogeneity properties of the medium (as noted in [10,16,36,38]), impulsive use of substrates (as explained in [31,32,47]), intra-species competition [22], multiple substrates [10,20], and deterministic or stochastic time-varying inputs (as explained, e.g., in [2,3,7,14,19,29,30,41,45]). In this work, we study another approach, based on an alternative model that we describe next.
in order to promote the coexistence of all species and provides sufficient conditions on the nonnegative constant species inputs x 0 i for i = 1, . . . , n that can ensure the coexistence of multiple species, which puts [39] outside the scope of the previously cited works. One approach that was used in [39] involved polytopic Lyapunov functions [9], which are related to other types of piecewise continuous or piecewise differentiable Lyapunov functions such as those in [46].
The system (4) can be seen as a limiting dynamics for a chain of interconnected chemostats, by the following argument. If we consider the system (1) with n = 2 under condition (2), then the competitive exclusion principle will imply that the concentration of the first species will converge to 0 as t → +∞. To promote the coexistence of the two species, the first species (which is the less advantaged competitor) is cultivated in a first chemostat, whose dynamics is described by the system whose asymptotic behavior is described by (3) with n = 1, that is If the output of (5) becomes the input of (1) with n = 2, then we obtain the coupled system where x ij is the concentration of the ith species in the jth chemostat, while s i is the concentration of the same nutrient in the ith chemostat. The last three equations of (6) can be written as the perturbed system   where the perturbations q 0 (t) = D[s 1 (t) − s * ] and q 1 (t) = D[x 11 (t) − x * 1 ] converge to zero as t → +∞, so (7) is asymptotically autonomous to   which is a particular case of (4) with n = 2, s in = s * , x 0 1 = x * 1 and x 0 2 = 0. A generalization for a chain of n chemostats can be done in a recursive way, by considering the competitive abilities described in (2).
In this article, we propose a generalization of the model (4) by taking into in account (i) the presence of additive perturbations and (ii) the effect of delays in the growth responses [31,46]. Therefore, the system (4) is changed to The functions δ i : [0, +∞) → [d i ,d i ] are measurable and essentially bounded (for i = 0, 1, . . . , n) and have known constant lower and upper bounds d i andd i , respectively. The system (7) is a particular case of (9) but the functions δ i can describe uncertainties such as unmodeled features (for example, disturbances in the chain of chemostats mentioned above), external perturbations, or uncertainties in the input concentrations which commonly occur in applications. For example, a perturbed biomass input x 0 i + ∆ i (t) can be represented by an additive perturbation δ i (t) = D∆ i (t). Also, a disturbance ∆ D added to the dilution rate can be captured by setting δ . . , n, and uncertainties in the uptake functions µ i or in s in can be modeled in an analogous way using suitable choices of the δ i 's.
Our assumptions will imply that d 0 > −Ds in and d i ≥ −Dx 0 i for i = 1, 2, . . . , n; see Section 2 for our assumptions, in which the lower bounds d i are allowed to be negative for each i ≥ 1 such that x 0 i > 0 (so the disturbances need not be positive valued). Therefore, all solutions of (9) with positive valued initial conditions are in the positive orthant for all positive instants, so (9) has the state space X = (0, +∞) n+1 . A key component of our analysis is that we will prove a uniform persistence condition, where we compute positive lower bounds on the x i (t) values for all i such that x 0 i > 0 and for sufficiently large t; see Lemma 3.3. The model (9) assumes the existence of a time interval [0, τ i ] necessary for the ith microbial species metabolize the nutrient. The delays τ i have been reported experimentally in several works, including [1,4,40]. Other delays can arise in our model but will not be considered in this note (but their effects can be incorporated into our δ i 's). For example, in [13,17,23] it is assumed that there is mortality of microbial biomass and a fraction is recycled into nutrient with some time delay. The mortality is also considered in [6,21,44], where the growth of biomass term x i (t)µ i (s(t − τ i )) of (9) is replaced by e −diτi x i (t − τ i )µ i (s(t − τ i )). Through a change of coordinates (based on a scaling the x i 's, x 0 i 's, and δ i 's for i ≥ 1), the model (9) becomes     ṡ where for simplicity we kept the same notation, and (10) will be the subject of this paper.
1.3. Structure of article. In the next section, we provide our theorem for (10), which we prove in Section 3 and illustrate in an example in Section 4. Our work is novel in its use of the more general model (10) (which we believe had not been studied in the presence of nonzero delays or uncertainties) and our use of a new Lyapunov-Krasovskii functional method that is beyond the scope of our prior Lyapunov function designs (such as those in [24,39], which were confined to undelayed systems). Also, the equilibrium that we stabilize is in the boundary of the state space X when at least one of the x 0 i 's is 0. Moreover, while our previous work [39] required that the nth species input be x 0 n = 0, here we allow a range of possible x 0 n values, so we cover a much broader class of equilibria than [39].
2. Definitions, assumptions, and main result. Our main result provides sufficient conditions ensuring input-to-state stability (or ISS) properties for the dynamics for the error variable with respect to the disturbance vector δ = (δ 0 , δ 1 , . . . , δ n ), for a large class of possible equilibrium points E * = (s * , x * ), where x * = (x 1 * , . . . , x n * ) and x = (x 1 , . . . , x n ). The ISS framework is used extensively across engineering; see [18] for the definition of ISS for undelayed systems without state constraints. To allow delays and state constraints, we use a variant of the usual ISS property. To explain this variant, we first need several definitions. Let K ∞ be the set of all continuous strictly increasing unbounded functions γ : [0, +∞) → [0, +∞) such that γ(0) = 0. Also, KL is the set of all continuous functionsβ : [0, +∞) × [0, +∞) → [0, +∞) such that (i) for each t ≥ 0, the function f (s) =β(s, t) is of class K ∞ and (ii) for each s ≥ 0, the function g(t) =β(s, t) is nonincreasing and satisfies lim t→+∞ g(t) = 0. We also set q t ( ) = q(t + ) for all ≤ 0, t ≥ 0, and functions q for which the equality is defined. By ISS of a delay system of the formq(t) = F(q t , δ(t)) with respect to a pair (D, S), we mean that there exist functions β ∈ KL and γ ∈ K ∞ such that holds for all t ≥ 0, all solutions q(t) of the system that have initial states q(0) ∈ S, and all measurable essentially bounded functions δ : [0, +∞) → D. Here and in the sequel, we assume that the initial conditions of our systems are constant, | · | is the usual Euclidean norm, | · | [0,t] is the essential supremum over [0, t], and | · | ∞ is the essential supremum. Later we specialize the preceding definitions to cases where q = E, and where q is an error vector in a different set of variables that we introduce later. Our first assumption is: and the constants x 0 i are all nonnegative. Also, the constant D satisfies 0 < D < µ n (s in ).
Assumption 1 can always be satisfied for all choices of the constant D ∈ (0, sup s≥0 µ n (s)), by first fixing s * such that µ i (s * ) < D for all i, and then choosing the x 0 i 's such that D < µ n (s in ) and s in > 0, i.e., we view the x 0 i 's and s in as constant controls. Although this differs from the usual treatment of controls where the controls depend on time (i.e., open loop controls) or on the state (as in feedback controls), our use of constant controls is sufficient for our delays and robustness analysis (and is included in the usual framework of open loop or feedback controls, by specializing the usual framework to cases where the controls are constant ones). By the symmetry of the dynamics for the x i 's, we can replace the condition D < µ n (s in ) by the requirement that D < µ i (s in ) for some i, by renumbering the species. Assumption 1 implies that when the δ i 's and τ i 's are 0, the system (10) admits the equilibrium E * = (s * , x * ), where Then E * ∈ [0, +∞) n+1 , and when the x 0 i 's are all positive, we have E * ∈ (0, +∞) n+1 . From Assumption 1, we have µ n (s * ) < µ n (s in ), so since µ n is strictly increasing, we have s * ∈ (0, s in ). Our assumption on the measurable essentially bounded uncertainties δ i (t) is as follows, where P = {i ∈ {1, . . . , n} : x 0 i > 0}: . . , n} \ P. Assumption 2 allows d i =d i = 0 for i = 0, 1, . . . , n which corresponds to cases where the δ i 's are all zero, but is far more general, e.g., because thed i 's for i ≥ 1 can be as large as we want. Notice that in important cases where the x 0 i 's are all positive, there is no systematic positive bias in the disturbance values, since for instance, we can allow δ to be take its values in sets such as that are symmetric hypercubes centered at the origin. To state our assumption on the constant delays τ i , we assume that the µ i 's have the Monod's form with m i > 0 and a i > 0 being respectively the maximal growth rate and the half saturation constants for all i ∈ {1, 2, . . . , n}, and we set where the constants ≥ s in will be specified later. Finally, we use the positive constants The motivation for the constant Γ 0 is that it is a lower bound for the function for all s ∈ [0,s ]. The function (18) will play an important role in the proof of our theorem, but is not needed to state our theorem. Our main result is: Theorem 2.1. If Assumptions 1-2 hold, and if there exists a constants ≥ s in such that each of the constant delays τ i satisfies where M is defined by then for all constants x > 0, the dynamics for the error (11) are ISS with respect Our proof of Theorem 2.1 will show how the constants s 1 and α i described in (17) are related to the lower bounds of the solutions of (10) for arbitrarily large values of t. Note also that M depends on the τ i 's for all i = 1, . . . , n and is equal to zero when these delays are 0, which implies that there are positive valuesτ i such that (19) is satisfied when the delays τ i all satisfy 0 ≤ τ i ≤τ i . Since However, sinces ≥ s in and x > 0 are arbitrary, it follows that when the δ i 's are zero, Theorem 2.1 implies that all solutions (s(t), x(t)) of (10) starting in X = (0, +∞) n+1 remain in X at all positive times and satisfy lim t→+∞ (s(t), x(t)) = (s * , x * ). This ensures uniform persistence of the ith species for all i ∈ P, i.e., x i has a positive lower bound (and lim t→+∞ x i (t) = 0 for all i ∈ {1, . . . , n} \ P). Our results are new, even in the special cases where the delays τ i or the uncertainties δ i are all zero, because [39] did not include delays or uncertainties. Since we do not restrict the values ofd i ≥ 0 for i ≥ 1, we obtain ISS under arbitrarily large upper bounds on the disturbances δ i in the x subdynamics. The functions β ∈ KL and γ ∈ K ∞ in the final ISS estimate will depend ons, x, and the disturbance boundsd i and d i from Assumption 2, and our proof of Theorem 2.1 can be used to provide an algorithm for constructing β and γ. Moreover, as noted in Section 1.2, several types of uncertainties (including in the dilution rate or in the concentration of species input) can be captured by our additive uncertainties δ i , so our ISS results can estimate or measure the effects of these types of uncertainties.

3.1.
Changes of variables and preliminary state bounds. We setτ = max i τ i . Since (10) is forward complete on X = (0, +∞) n+1 , we can first fix any solution (s(t), x(t)) of the system all of whose components are positive for all t ≥ 0, with x n (0) > x. Then Assumption 2 implies that s(t) ≤s for all t ≥ 0, since we would haveṡ(t) < 0 for all t for which s(t) >s and s(0) ≤s. We use the new variables for all i ∈ {1, 2, . . . , n}. Then we obtain the new system in our new variables, where we also set ∆ 0 = δ 0 for consistency. We next prove three lemmas which produce three functions T i ∈ K ∞ for i = 1, 2, 3, whose class K ∞ properties will be used later to build an ISS estimate that is valid for all times Our first lemma is: Proof. Assumptions 1-2 imply that for all t ≥ 0, we have s(t) > 0 and Hence, to prove the lemma, it suffices to consider the case where s(0) > s in +d 0 /D and consequentlys(0) = 0. In this case, it follows that for all t ≥ 0 such that s( ) > s in +d 0 /D for all ∈ [0, t], we also havė s( ) < 0. Therefore, for any such t, we can use the facts that D < µ n (s in ) and µ n is nondecreasing, combined with the inequalities Dx 0 n + d n ≥ 0 and s * < s in ≤ s( ) for all ∈ [0, t], to deduce that for all ∈ [0, t], we have α n ( ) ≥ α n (0) ≥ xe −µn(s * )τn and D(s in − s( )) + ∆ 0 ( ) ≤ 0, so (by integratingṡ( ) on [0, t]) which implies that t ≤ |s(0)|e 2µn(s )τn /(µ n (s in )x), where we used the fact that −µ n (s(p)) ≥ −µ n (s ) for all p ∈ [0, t]). Therefore, there is a time t * ∈ 0, 2|s(0)|e 2µn(s )τn /(µ n (s in )x) (24) such that s(t * ) ≤ s in +d 0 /D (since by the preceding argument, the largest possible time t such that s(t) > s in +d 0 /D is |s(0)|exp(2µ n (s )τ n )/(µ n (s in )x), so at any later times t, it must be the case that s(t) ≤ s in +d 0 /D), which by the first part of the proof, implies that s(t) ≤ s in +d 0 /D for all t ≥ t * . Hence, we can choose T 1 (r) = 2e 2µn(s )τn r/(µ n (s in )x).
Consider the operator Then the constants L i satisfy L i ≤L i for all i, where the constantL i 's were defined in (17); this follows from our assumptions that the initial conditions are constant combined with the facts that µ i 's are nondecreasing and s( ) ≤s for all ≥ 0, which gives µ i (s( )) ≤ µ i (s ) for all ≥ −τ i and all i. Pick any constant λ 1 > 1 such that our delay upper bounds from (19) are satisfied with C replaced by λ 1 C in the formulas from (17). Such a λ 1 exists, because of the strictness of the inequalities in (19), by picking λ 1 > 1 close enough to 1. We next prove: Lemma 3.2. If Assumptions 1-2 hold, then there is a T 2 ∈ K ∞ such that σ(t) ≤ λ 1 C for all t ≥ T 2 (|E(0)|).
We now fix a constant λ 2 ∈ (0, 1) such that our delay bounds from (19) hold with C replaced by λ 1 C in (17) and with s 1 and the α i 's replaced by respectively for all i in (17). As in the λ 1 case, the existence of such a λ 2 follows from the strictness of the inequalities in our conditions on the delays, by choosing λ 2 ∈ (0, 1) and λ 1 > 1 both to be close enough to 1. The following lemma provides useful positive lower bounds on s(t), and on the α i (t)'s for all i ∈ P: Proof. For all t ≥ T 2 (|E(0)|), Lemma 3.2 and our choice of σ in (25) give α i (t) ≤ λ 1Li C for all i ∈ {1, 2, . . . , n}. Hence, for all t ≥ T 2 (|E(0)|), the (s, α) dynamics described by (22), and the Monod's description (15) for the growth functions, givė where we also used the fact that Dx 0 i + d i > 0 for all i ∈ P to get (32). The right side of (31) is bounded below by the positive value (Ds in + d 0 )(1 − λ 2 ) if t is such that s(t) ≤ s λ . Also, for each i ∈ P, the right side of (32) is bounded below by the positive value (1 − λ 2 )(Dx 0 i + d i )e −µi(s * )τi if t is such that α i (t) ≤ α iλ . Hence, for any t 0 ≥ 0 such that s(t 0 ) ≥ s λ , we have s(t) ≥ s λ for all t ≥ t 0 . Also, for each i ∈ P, and for any t 0 ≥ 0 such that α i (t 0 ) ≥ α iλ , we have α i (t) ≥ α iλ for all t ≥ t 0 . Hence, to construct T 3 , it suffices to choose T 3 ∈ K ∞ such that the following conditions hold: (i) If s(0) < s λ , then s(t) ≥ s λ for some t ∈ (0, T 3 (|E(0)|) −τ ] (which implies that s(t) ≥ s λ for all t ≥ T 3 (|E(0)|) −τ , by the preceding argument) and (ii) for each i ∈ P such that α i (0) < α iλ , we have α i (t) ≥ α iλ for some t ∈ (0, T 3 (|E(0)|) −τ ] (which implies that α i (t) ≥ α iλ for all t ≥ T 3 (|E(0)|) −τ , also by the preceding argument).
To satisfy the preceding conditions (i)-(ii), first note that if we choose any constant T L > 0 such that  (21) gives s(0) ≥ λ 2 s * ≥ s λ and α i (0) ≥ x i (0)e −τiµi(s * ) ≥ λ 2 x i * e −τiµi(s * ) ≥ α iλ for all i ∈ P, which imply that s(t) ≥ s λ and α i (t) ≥ α iλ all hold for all t ≥ 0 and i ∈ P. Hence, we can choose where the formula for the restriction of T 3 to [0, T M ] was chosen to ensure that T 3 is of class K ∞ .

Construction of a Lyapunov-like functional.
We define the transformed error vectorÊ(t) = (s(t),α 1 (t), ...,α n (t)), and we choose the C 1 function for all i ∈ P, and Ψ i (α i ) = α i for all i ∈ {1, . . . , n} \ P. Since ν (s) =s/s and Ψ i (α i ) =α i /α i hold for all values ofs, theα i 's, and i ∈ P, it follows from the chain rule that along all solutions of (41) and for all i ∈ {1, 2, . . . , n}, we havė It follows that for all t ≥ T 3 (|E(0)|), the time derivative of V 1 along the solutions of (41) satisfieṡ where and where the positive constants q i = p i /c i = (D − µ i (s * ))/c i are defined in (17), and where the formulā follows from Lemmas 3.2-3.3 and the relationship between the ∆ i 's and the δ i 's in (21), and was used to bound the coefficients of the ∆ i 's in the curly braces in (42). Our choice of V 1 is motivated by the fact that if the delays and uncertainties are all 0, then V 1 is a Lyapunov function for the system (41). By Lemma 3.2, our decay estimate (43) on V 1 giveṡ where N (Ê(t)) is the quantity in curly braces in (43). By using the Mean Value Theorem, we can deduce the inequality for all X ∈ R. We apply (46) twice, by choosing X to be the integrals in (45). We deduce that when t ≥ T 3 (|E(0)|), theṅ and therefore alsȯ where the B i 's are defined in (17). Using the second equality in (39), and the fact that the function Γ from (35) satisfies Γ(s) ≥ Γ 0 for all s ∈ (0,s ], it follows that for all t ≥ T 3 (|E(0)|), we havė where the r i 's are defined in (17) and b i = De τiBi x 0 i m i (and where we used the formula µ i (s) − µ i (s * ) = c i m is /(a i + s)). We next use the positive lower bounds s λ and α iλ from Lemma 3.3, to upper bound the nonnegative valued terms contained between the brackets in (48), and then we will convert (48) into a decay condition on a suitable Lyapunov-like function.
3.4. Converting (48) into a decay estimate. We can use Jensen's inequality and Lemma 3.3 to check that for all t ≥ T 3 (|E(0)|), we have for all i ∈ P, where and where we also applied the inequality c 1 c 2 ≤ 1 2 c 2 1 + 1 2 c 2 2 and then Young's inequality c 1 c 2 ≤ 1 4 c 2 1 + c 2 2 , with c 1 and c 2 being the corresponding terms in curly braces in (49) in both applications, and we used the fact that s( ) ≤s for all ≥ 0 to get where M λ (τ ) Our bound (19) from Theorem 2.1 and our choices of the λ i 's imply thatτ M λ (τ )λ 1 < Γ 0 /2, by choosing λ 2 ∈ (0, 1) and λ 1 > 1 close enough to 1. Therefore, there is a constant v 0 such that along all solutions of (35) in our state space X = (0, ∞) n+1 and for all times t ≥ T 3 (|E(0)|), the time derivative of the function where we used the relations to find v 0 . The remainder of the proof consists of converting (53) into an ISS estimate in the transformed error variableÊ(t) = (s(t),α 1 (t), ...,α n (t)) that is valid for all t ≥ 0, which we then convert into an ISS estimate in the original error variable E = (s, x) − E * ; for details for the completion of the proof, see the appendix below.
We apply our theorem with the choices = s in . Then Assumption 1 is satisfied with s in = 1.34412, and the corresponding equilibria (14) are x 1 * = 1.29412 and x 2 * = 1.1. In this illustration, we will use the vectors δ(t) to model uncertainties in applying the constant input concentrations (which may occur in applications, because it may be difficult to maintain the inputs x 0 i at constant levels), so we can set δ 0 (t) = 0 and therefore choose d 0 =d 0 = 0. To isolate the effects of delays in one of the species, we assume that τ 1 is positive but that τ 2 = 0 (but analogous results can be obtained if the delay is in the dynamics for the second species, or in the dynamics for both species). Then Assumption 2 is satisfied with d 1 = −d 1 =d 2 = −d 2 = 0.1, and our requirements (19) on the delays τ i are satisfied with the choices = s in for all values τ 1 ∈ [0, 0.145).
We simulated the dynamics (10) using the command NDSolve in Mathematica [43], with the preceding choices of the parameters, the delays τ 1 = 0.14 and τ 2 = 0, and the disturbance vector δ(t) = (δ 0 (t), δ 1 (t), δ 2 (t)) = (0, −0.1 sin(t), 0.1 cos(t)), to isolate the effects of uncertainties in applying the constant input concentrations x 0 i in the x subdynamics. We report our results in the figures below, with the initial state (s(0), x 1 (0), x 2 (0)) = (0.2, 0.1, 1), and then with the initial state (s(0), x 1 (0), x 2 (0)) = (1.3, 0.2, 0.1). The figures show rapid convergence towards an oscillatory steady state, with a deviation from the equilibrium point (s * , x 1 * , x 2 * ) = (0.5, 1.29412, 1.1) that can be explained by the presence of the uncertainties δ 1 and δ 2 , and therefore help illustrate our theory. Proof. Choose any constant s d ≥ s in + (d 0 /D) such that µ i (s d ) > D for all i ∈ {1, 2, . . . , n}, and fix any solution (s(t), x(t)) of (10) whose initial value is componentwise positive. Then Lemma 3.1 provides a constant t a ≥ 0 such that s(t) ≤ s d for all t ≥ t a , and our assumptions ensure that all components of (s(t), x(t)) are nonnegative valued for all t ≥ 0. Hence, s(t) is bounded. Set D * = Ds in +d 0 . Then our assumptions giveṡ(t) ≤ D * for all t ≥ 0, which gives for all values m 1 ≥ 0 and m 2 ≥ m 1 . Also, for each i ∈ {1, 2, . . . , n} and all t ≥ t a +τ , we haveẋ i (t) ≤ x i (t)[µ i (s d ) − D] + Dx 0 i +d i for all t ≥ 0, which we can integrate to obtain Since lim k→∞ x j (t k ) = +∞, it follows from (65) that there is k p > 0 such that for all k ≥ k p , we have for all m ∈ I k . Combining (66) with (63), we obtain Since lim k→∞ x j (t k ) = +∞, it follows that there is a k ≥ k p such that s(t k −τ ) < 0. This yields a contradiction with the fact that s(t) ≥ 0 for all t ≥ 0.
6. Conclusions. We solved an input-to-state stabilization problem for a chemostat model with one limiting substrate, an arbitrary number of competing species, a constant dilution rate, delays in the uptake functions, and uncertainties. We used the constant species inputs and the input nutrient concentration as constant controls, and these controls can be chosen to input-to-state stabilize a large class of possible equilibria. In the special case where all of the constant inputs x 0 i are positive and the disturbances are zero, we first proved (in Lemma 3.3) that all solutions whose initial states are in the positive orthant are uniformly persistent, meaning, there is a positive lower bound on the species levels. Then by using a Lyapunov functional, we proved that the solutions asymptotically converge towards a positive equilibrium, which generalizes [39]. We use the x 0 i 's and the input nutrient concentration s in as constant controls.
To cope with delays in uptake functions or uncertainties, we used a new Lyapunov functional approach. The decay estimate for our Lyapunov-like function made it possible to prove robustness to uncertainties. We hope to generalize our work to larger classes of models with multiple species and multiple limiting substrates. Our choices of the uptake functions (15) made it possible to obtain lower and upper bounds in several intermediate steps of the construction of the Lyapunov functional. The Michaelis-Menten growth functions that we used provide a wide range of generality for our results, since the growth rates of many species are described by such functions. A more general result is still unsolved for more general uptake functions, due to the technical difficulties induced by the lack of monotonicity and more complex nonlinearities.
To extend (A.7) to obtain an ISS estimate that is valid for all t ≥ 0, first note that the structure (38) of theÊ dynamics and the Mean Value Theorem estimate (46), combined with our bounds on the µ i 's, s(t), and Γ and the global Lipschitzness of the µ i 's, provide a constantL (that is independent of the choice of the solution) such that where the maxima are over all i ∈ {1, 2, . . . , n}; this can be checked by comparing the formulas for the components of E andÊ. This proves the theorem.