On the Stability of Periodic Solutions in the Perturbed Chemostat

We study the chemostat model for one species competing for one nutrient using a Lyapunov-type analysis. We design the dilution rate function so that all solutions of the chemostat converge to a prescribed periodic solution. In terms of chemostat biology, this means that no matter what positive initial levels for the species concentration and nutrient are selected, the long term species concentration and substrate levels closely approximate a prescribed oscillatory behavior. This is significant because it reproduces the realistic ecological situation where the species and substrate concentrations oscillate. We show that the stability is maintained when the model is augmented by additional species that are being driven to extinction. We also give an input-to-state stability result for the chemostat tracking equations for cases where there are small perturbations acting on the dilution rate and initial concentration. This means that the long term species concentration and substrate behavior enjoys a highly desirable robustness property, since it continues to approximate the prescribed oscillation up to a small error when there are small unexpected changes in the dilution rate function.


1.
Introduction. The chemostat model provides the foundation for much of current research in bio-engineering, ecology, and population biology [6,8,10,11,21]. In the engineering literature, the chemostat is known as the continuously stirred tank reactor. It has been used for modeling the dynamics of interacting organisms in waste water treatment plants, lakes and oceans. In its basic setting, it describes the dynamics of species competing for one or more limiting nutrients. If there are n species with concentrations x i for i = 1, . . . , n and just one limiting nutrient with concentration S and dilution rate D > 0, then the model takes the form where µ i denotes the per capita growth rate of species i andṗ is the time derivative of any variable p. (In much of the paper, we simplify our notation by omitting the arguments of the functions. For instance, when no confusion can arise from the context, we denote S(t) simply by S.) The functions µ i depend only on the nutrient concentration, and are zero at zero, continuously differentiable and strictly increasing, although non-monotone functions have been the subject of research as well. The conversion of nutrient into new biomass for each species i happens with a certain yield γ i ∈ (0, 1) and the natural control variables in this model are the input nutrient concentration S in and the dilution rate D. The latter variable is defined as the ratio of the volumetric flow rate F (with units of volume over time) and the reactor volume V r which is kept constant. Therefore it is proportional to the speed of the pump that supplies the reactor with fresh medium containing the nutrient. The equations (1) are then straightforwardly obtained from writing the mass-balance equations for the total amounts of the nutrient and each of the species, assuming the reactor content is well-mixed. The full model (1) is illustrated in Figure 1. In the present work, we consider the case where there is just one species with concentration x, in which case the equations (1) take the form (but see Theorem 3 below for results on chemostats with disturbances, and Section 8 for models involving several species). We assume S in is a given positive constant, while the per capita growth rate µ is a Monod function (which is also known as a Michaelis-Menten function) taking the form for certain positive constants m and a that we specify later. The dilution rate is an appropriate continuous positive periodic function we also specify below. SinceṠ ≥ 0 when S > 0 is near zero, one can readily check that (2) leaves the domain of interest positively invariant (i.e., trajectories for (2) starting in X remain in X for all future times); see Theorem 3 for a more general invariance result for perturbed chemostats.
Since we are taking S in to be a fixed positive constant, we rescale the variables to reduce the number of parameters. Using the change of variables and dropping bars, we eliminate S in and γ and so obtain the new dynamics again evolving on the state space X = (0, ∞) 2 . Motivated by the realistic ecological situation where the species concentrations are oscillating, we solve the following biological problem: Biological Problem B1: For a prescribed oscillatory behavior for the species concentration and substrate level given by a time-periodic pair (S r (t), x r (t)), design a dilution rate function D(t) such that if this choice of D(t) is used in the chemostat (4), then all solution pairs (S(t), x(t)) for the substrate levels and corresponding species levels obtained from solving (4) (i.e. for all possible initial values) closely approximate (S r (t), x r (t)) for large times t. See also Problem (SP) in Section 4 below for a precise mathematical statement of the preceding problem. In the language of control theory, solving Biological Problem B1 means we will prove the stability of a suitable periodic reference signal for the species concentration t → x r (t) in (4) using an appropriate time-periodic dilution rate D(t); see [15] for the fundamental ideas from control theory we need in the sequel.
Since D(t) is proportional to the speed of the pump which supplies the chemostat with medium containing nutrient, implementation of the prescribed oscillatory behavior requires that we control the pump in a very precise way. In practice this control process is prone to errors, and the actual pump speed will be subject to small fluctuations which we will model by replacing D(t) by D(t) + u 1 (t) in the chemostat equations, where u 1 (t) models the error. It is therefore of interest to study the effect of these small fluctuations on the periodic behavior. Preferably this effect will be small, and the resulting behavior is not too different from the prescribed periodic behavior. We will show that this is indeed the case by actually quantifying how small the deviations are, relying on the well-known control-theoretic notion of Input-to-State Stability or ISS; see [23,24] and Remark 2 for details about ISS. Summarizing this a bit more formally, we will solve the following biological problem (which we state in a more precise mathematical way in Section 5): Biological Problem B2: For the prescribed oscillatory behavior (S r (t), x r (t)) and dilution rate D(t) obtained in Biological Problem B1, quantify how the substrate and species levels (S(t), x(t)) in the chemostat model (4) are affected by unexpected changes in the dilution rate, and show that the convergence of (S(t), x(t)) to the oscillatory behavior (S r (t), x r (t)) is robust to small changes in D(t).
Our solution to Biological Problem B2 will be a special case of our more general input-to-state stability result for the chemostat tracking equations, assuming the dilution rate and initial concentration are both perturbed by noise terms of small enough magnitude.
In the next section, we briefly review the literature focusing on what makes our approach different. In Section 3, we fix the reference signal we wish to track. In Section 4, we precisely formulate the definitions and the stability problem we are solving. We state our main stability theorem in Section 5 and we discuss the significance of our theorem in Section 6. We prove our stability result in Section 7. In Section 8, we show that the stability is maintained when there are additional species that are being driven to extinction. We validate our results in Section 9 using a numerical example. We conclude in Section 10 by summarizing our findings.
2. Review of the Literature and Comparison with Our Results. The behavior of the system (1) is well understood when S in and D are positive constants, as well as cases where n = 2 and either of these control variables is held fixed while the other is periodically time-varying. See [13,20] for periodic variation of S in and [4] for periodic variation of D and the general reference [21] on chemostats. When both S in and D are constants, the so-called "competitive exclusion principle" holds, meaning that at most one species survives. Mathematically this translates into the statement that system (1) has a steady state with at most one nonzero species concentration, which attracts almost all solutions; see [21]. This result has triggered much research to explain the discrepancy between the (theoretical) competitive exclusion principle and the observation that in real ecological systems, many species coexist.
The results on the periodically-varying chemostat mentioned above should be seen as attempts to explain this paradox. They involve chemostats with n = 2 species, and their purpose is to show that an appropriate periodic forcing for either S in (t) or D(t) can make the species coexist, usually in the form of a (positive) periodic solution. Few results on coexistence of n > 2 species are available. An exception is [19], where a periodic function S in (t) is designed (with D kept fixed) so that the resulting system has a (positive) periodic solution with an arbitrary number of coexisting periodically varying species. The stability properties of this solution are not known.
More recent work has explored the use of state-dependent but time invariant feedback control of the dilution rate D to generate coexistence; see [6,8] for monotone growth rate functions in the n = 2 species case, and [7] for the n = 3 species case. The paper [11] considers feedback control when the growth rate functions are non-monotone. In [12], [16], and [18], coexistence is proved for models taking into account intra-specific competition. In these models, the usual growth functions µ i (S) are replaced by functions µ i (S, x i ) which are decreasing with respect to the variable x i . All the results discussed so far apply to a more general model than (4) involving n > 1 species. This is because the main purpose of these papers is to investigate environmental conditions under which the competitive exclusion principle fails and several species can coexist.
Here we will not consider any coexistence problems. Our main objective is to provide a proof of stability of a periodic solution based on a Lyapunovtype analysis and to investigate the robustness properties of the periodic solution with respect to perturbations. As an illustration we show that the stability of the periodic solution is robust with respect to additional species that are being driven to extinction, or to small disturbances on the initial nutrient concentration or dilution rate. These features set our work apart from the known results on periodically forced chemostat models which do not rely on the construction of a Lyapunov function. Proving stability in the chemostat usually relies on reduction and monotonicity arguments, and not so often on Lyapunov functions (but see for instance Theorem 4.1 in [21] which uses a Lyapunov function introduced in [14] and more recently [12]).
Finally we point out that closely related to our results is [10] where a single-species chemostat with a continuous and bounded (but otherwise arbitrary) function S in (t) and constant dilution rate is investigated; there it is shown that two positive solutions converge to each other. However, the proof is not based on a Lyapunov function. The advantage of having a Lyapunov function is that it can be used to quantify the effect of additional noise terms on the stability of the unperturbed dynamics. In fact, to our knowledge, our work provides the first input-to-state stability analysis of chemostats whose dilution rates and initial concentrations are perturbed by small noise; see Remark 2 for a discussion on the importance of input-to-state stability in control theory and engineering applications.
3. Choosing a Reference Trajectory. We first choose the dilution rate D = D(t) that will give a reference trajectory (S r (t), x r (t)) for (4) which we show to be stable. We assume a growth rate with constants m, a > 0 as follows, in which the constants a and m and the variable S are all dimensionless by the change of coordinates used to obtain the normalized equations (4) so the units do not matter: For the sake of computational simplicity, we choose a sinusoidal reference trajectory but the extension to more general reference trajectories can be handled by similar methods; see Remark 5 for details. Simple calculations show that (4) admits the trajectory which we refer to as a reference trajectory when we choose Condition (5) then provides constantsD, D o > 0 such that for all t ≥ 0:D See Figure 2 for the graph of D(t) for m = 10 and a = 1 2 .  (7) 4. Definitions and Statement of Stability Problem. We wish to solve the following stability problem which is merely a restatement of Biological Problem B1 above in precise control theoretic terms: (SP) Given any trajectory (S, x) : [0, ∞) → (0, ∞) 2 for (4) corresponding to the dilution rate D(t) from (7) and µ as in (5) (i.e. for any initial value for (S, x)), show that the corresponding deviation (S(t),x(t)) := (S(t)−S r (t), x(t)−x r (t)) of (S, x) from the reference trajectory (6) asymptotically approaches (0, 0) as t → +∞.
We will solve (SP) by proving a far more general tracking result for a single species chemostat acted on by a disturbance vector u = (u 1 , u 2 ) : [0, ∞) → R 2 as follows: .
We will quantify the extent to which the reference trajectory (6) tracks the trajectories of (9) which will solve Biological Problem B2 from the introduction. To this end, we need to introduce a priori bounds on u 1 and u 2 ; see Remark 8. Our main theoretical tool will be the input-to-state stability (ISS) property [22] which is one of the central paradigms of current research in nonlinear stability analysis; see Remark 2. The relevant definitions are as follows.
We let K ∞ denote the set of all continuous functions γ : [0, ∞) → [0, ∞) for which (i) γ(0) = 0 and (ii) γ is strictly increasing and unbounded. We let KL denote the class of all continuous functions β : for which (I) β(·, t) ∈ K ∞ for each t ≥ 0, (II) β(s, ·) is non-increasing for each s ≥ 0, and (III) β(s, t) → 0 as t → +∞ for each s ≥ 0. Consider a general control-affine dynamiċ We always assume that such solutions are uniquely defined on all of [t o , ∞) (i.e., (10) is forward complete and O is positively invariant for this system) and that there exists Θ ∈ K ∞ such that |F (y, t)| + |G(y, t)| ≤ Θ(|y|) everywhere, where | · | is the usual Euclidean norm. For example, is the solution of (9) for the disturbance u = (u 1 , u 2 ) = α(t) satisfying the initial condition ( Definition 1. We call (10) input-to-state stable (ISS) provided there exist β ∈ KL and γ ∈ K ∞ such that Here |α| ∞ denotes the essential supremum of α ∈ U . By causality, the ISS condition (11) is unchanged if |α| ∞ is replaced by the essential supremum |α| [to,t] of α restricted to [t o , t]. In particular, (11) says y(t; t o , y o , Z) → 0 as t → +∞ for all initial values y o and initial times t o , where Z is the zero disturbance α(t) ≡ 0.
Remark 2. The theory of ISS systems originated in [22]. ISS theory provides the foundation for much current research in robustness analysis and controller design for nonlinear systems, and has also been used extensively in engineering and other applications [1,2,5,17,22,24]. The ISS approach can be viewed as a unification of the operator approach of Zames (e.g. [25,26]) and the Lyapunov state space approach. The operator approach involves studying the mapping (t o , y o , α) → y(·; t o , y o , α) of initial data and control functions into appropriate spaces of trajectories, and it has the advantages that it allows the use of Hilbert or Banach space techniques to generalize many properties of linear systems to nonlinear dynamics. By contrast, the state space approach is well suited to nonlinear dynamics and lends itself to the use of topological or geometric ideas. The ISS framework has the advantages of both of these approaches including an equivalent characterization in terms of the existence of suitable Lyapunov-like functions; see Remark 7 below. For a comprehensive survey on many recent advances in ISS theory including its extension to systems with outputs, see [24].
To specify the boundū on our disturbances u = (u 1 , u 2 ), we use the following constants whose formulas will be justified by the proof of our main stability result: 5. Statement of Theorem. From now on, we assume the disturbance vector u = (u 1 , u 2 ) in (9) takes all of its values in a fixed square control set of the form where c, κ, and C 1 are in (12) (but see Remark 8 for related results under less stringent conditions on the disturbance values). We will prove the following robustness result: Theorem 3. Choose D(t), µ, and (S r , x r ) as in (5)-(7). Then the corresponding solutions of (9) satisfy Moreover, there exist β ∈ KL and γ ∈ K ∞ such that the corresponding transformed error vector satisfies the ISS estimate (11) for all α ∈ U , t 0 ≥ 0, t ≥ t 0 , S 0 > 0, and 6. Discussion on Theorem 3. Before proving the theorem, we discuss the motivations for its assumptions, and we interpret its conclusions from both the control theoretic and biological viewpoints.
Remark 5. Theorem 3 says that in terms of the error signals y, any componentwise positive trajectory of the unperturbed chemostat dynamics (9) converges to the nominal trajectory (6), uniformly with respect to initial conditions. This corresponds to putting α ≡ 0 in (11). It also provides the additional desirable robustness property that for an arbitrary U-valued control function α ∈ U , the trajectories of the perturbed chemostat dynamics (9) are "not far" from (6) for large values of time. In other words, they "almost" track (6) with a small overflow γ(|α| ∞ ) from the ISS inequality (11). Similar results can be shown for general choices of x r and D(t). For example, we can choose any x r (t) that admits a constant ℓ > 0 such that for all t ≥ 0 and S r = 1 − x r . In this case, we take the dilution rate which is again uniformly bounded above and below by positive constants. The proof of this more general result is similar to the proof of Theorem 3 we give below except with different choices of the constants c and κ.

Remark 6. The robustness result
of Theorem 3 differs from the classical ISS condition in the following ways: 1. For biological reasons, negative values of the nutrient level S and the species level x do not make physical sense. Hence, only componentwise positive solutions are of interest. Therefore, (15) is not valid for all (S 0 , x 0 ) ∈ R 2 but rather only for (S 0 , x 0 ) ∈ (0, ∞) × (0, ∞). 2. Our condition (15) provides an estimate on the transformed error component ln(x(t; t 0 , (S 0 , x 0 ), α)) − ln(x r (t)) instead of the more standard error x(t; t 0 , (S 0 , x 0 ), α) − x r (t). Our reasons for using the transformed form of the error are as follows. The function ln(x) goes to −∞ when x goes to zero. This property is relevant from a biological point of view. Indeed, in the study of biological systems, it is important to know if the concentration of the species is above a strictly positive constant when the time is sufficiently large or if the concentration admits zero in its omega limit set. In the first case, the species is called persistent. The persistency property is frequently desirable, and it is essential to know whether it is satisfied. Hence, the function ln(x(t; t 0 , (S 0 , x 0 ), α)) − ln(x r (t)) has the desirable properties that (a) it goes to +∞ if x(t; t 0 , (S 0 , x 0 ), α) does, (b) it is equal to zero when x(t; t 0 , (S 0 , x 0 ), u) is equal at time t to the value of x r , and (c) it goes to −∞ if x(t; t 0 , (S 0 , x 0 ), u) goes to zero. Therefore, roughly speaking, if the species faces extinction, then it warns us.

Remark 7. Our proof of Theorem 3 is based on a Lyapunov type analysis.
Recall that a C 1 function V : R n ×[0, ∞) → [0, ∞) is called an ISS Lyapunov function (ISS-LF) for (10) provided there exist γ 1 , γ 2 , γ 3 , γ 4 ∈ K ∞ such that 1. γ 1 (|y|) ≤ V (y, t) ≤ γ 2 (|y|) and 2. V t (y, t) + V y (y, t)[F (y, t) + G(y, t)u] ≤ −γ 3 (|y|) + γ 4 (|u|) hold for all y ∈ O, t ≥ 0, and u ∈ U. The function V we will construct in the proof of Theorem 3 is not an ISS-LF for the chemostat error dynamics because of the specificities of the state space, which preclude the existence of the necessary functions γ 1 , γ 2 ∈ K ∞ in Condition 1 above. Hence, we cannot directly apply the result that the existence of an ISS Lyapunov function implies that the system is ISS e.g. [9, Theorem 1] to prove our theorem. Instead, we prove our Theorem 3 directly from the decay inequality satisfied by the time derivative of V along the trajectories. The proof that the decay inequality implies ISS is very similar to that part of the proof of [9, Theorem 1], so we only sketch that part of our proof in the appendix.
Remark 9. If our control set U is replaced by the larger control set U ♯ := [−ū, +ū] 2 for any fixed constantū ∈ (0, min{1, D o }), then the error dynamics (18) instead satisfies the less stringent integral ISS property. The relevant definitions are as follows. We say that (10) is integral input-to-state stable (iISS) provided there exist δ 1 , δ 2 ∈ K ∞ and β ∈ KL such that everywhere for all measurable essentially bounded functions α : [0, ∞) → U ♯ . This condition is less restrictive than ISS since e.g.ẏ = − arctan(y) + u is iISS but not ISS [3]. An iISS-LF for (10) (with controls in U ♯ ) is then defined to be a C 1 function V : R n × [0, ∞) → [0, ∞) for which there are γ 1 , γ 2 , γ 4 ∈ K ∞ and a positive definite function γ 3 (i.e., γ 3 : [0, ∞) → [0, ∞) is continuous and zero only at zero) such that Conditions 1-2 in Remark 7 hold everywhere. This is less restrictive than the ISS-LF condition since γ 3 need not be of class K ∞ . Arguing as in the proof of Theorem 3 up through (28) and solving the appropriate constrained minimum problem to get γ 3 shows that V = L 3 satisfies the iISS Lyapunov function decay condition (namely, Condition 2 from Remark 7) for the error dynamics (18) and the control set U ♯ using γ 3 (s) = C 1 (e −s − 1) 2 /2. Therefore, this system is in fact iISS, by a slight variant of the proof of the iISS estimate in [3,Theorem 1]. We leave the details to the reader.
8. Stability in the Presence of Several Species. Theorem 3 shows that the stability of the reference trajectory (6) is robust with respect to small perturbations of the dilution rate and initial concentration. To further demonstrate the robustness of our results, we next show that the stability of (6) is also maintained when the model (4) is augmented to include additional species that are being driven to extinction, in the following sense.
We assume for simplicity that u 1 ≡ u 2 ≡ 0. Consider the augmented system where µ is as in (5) and ν i is continuous and increasing and satisfies ν i (0) = 0 for i = 1, 2, . . . , n. The variables y i represent the levels of n additional species. We choose D and D o as in (7) and (8), and we assume ν i (1) < D o for i = 1, 2, . . . , n. (This assumption is, in a sense, natural because one can easily check that it ensures that each species concentration y i converges to zero. Indeed, the fact that, for all t ≥ 0, D(t) > D o > 0 and µ(S(t))x(t) + n i=1 ν i (S(t))y i (t) ≥ 0 ensures, in combination with the inequalities ν i (1) < D o , that there exists an instant T > 0 and a constant c > 0 such that for all t ≥ T and for i = 1, 2, . . . , n,ẏ i (t) ≤ −cy i (t). This implies that y i 's converge to 0 exponentially.) We show that the transformed error (z,ξ,ỹ) : between any componentwise positive solution (S, x, y) of (33) and the reference trajectory (S r , x r , 0, . . . , 0) = 1 2 − 1 4 cos(t), 1 2 + 1 4 cos(t), 0, . . . , 0 converges exponentially to the zero vector as t → +∞.
To this end, notice that in the coordinates (34), the system (33) becomes by the same calculations that led to (18). Set L 1 (ξ) = exp(ξ) −ξ − 1 where exp(r) := e r . Since e ξr ≥ 1/4, z r ≡ 1, 0 ≤ z r −e ξr = S r = 1/2−cos(t)/4 ≤ 1, and 0 ≤ z − e ξ =z + 1 − e ξ ≤ 2 +z 2 , we deduce that the derivative of L 3 (ξ,z) := L 1 (ξ) + 4m aD oz 2 along the trajectories of (35) satisfieṡ where the second inequality is by (16) with p = |z|, q = |eξ − 1|, and d = 4 and the last inequality used the relation J 2 + K 2 ≥ −2JK for real values J and K. On the other hand, since ν i (1) < D o for each i, the form of the dynamics for S and the nonnegativity of µ and the ν i 's along our componentwise positive trajectories imply that there exist ε > 0 and T ≥ 0 such that (i) S(t) ≤ 1 + ε for all t ≥ T and (ii) ν i (1 + ε) < D o for all i = 1, 2, . . . , n. We deduce that, for all i = 1, 2, . . . , n and for all t ≥ T , where δ = D o − max{ν i (1 + ε) : i = 1, 2, . . . , n} > 0. Hence, each y i (t) converges exponentially to zero. Next notice that along each pair (ξ(t),z(t)), the function ∆(ξ,z) := ma(eξ − 1) 2 16(a + 1)(a + 2 +z 2 ) is positive if and only ifξ = 0. By (36) and (37), the time derivative of along the trajectories of (35) satisfieṡ provided t > T where T is chosen as above. (The second inequality in (39) follows because for any nonnegative a k , we get a k ≤ ( n i=1 a 2 i ) 1/2 which we sum and then square to get ( n i=1 a i ) 2 ≤ n 2 n i=1 a 2 i . The last inequality used ν i (S(t)) ≤ ν i (1 + ε) < D o for all t ≥ T and our choice of A. ) It is tempting to surmise from (39) and the structure of L 4 that L 4 is a Lyapunov function for (35) since then we could use standard Lyapunov function theory to conclude that (ξ(t),z(t), y(t)) asymptotically converges to zero. However, such an argument would not be technically correct, since the state space of (35) is not R n+2 (because the original augmented chemostat model (33) is only defined for componentwise nonnegative values of the state). Instead, we argue as follows (in which we may assume for simplicity that the initial time for the augmented error dynamics is zero).

Remark 10.
Notice that −L 4 is bounded below by a quadratic of the form c|(ξ,z, y)| 2 along the trajectories of (35), and that L 4 is bounded above and below by such quadratics along the trajectories as well, since the trajectories are bounded. From this fact and (39), one can deduce that the trajectories (ξ(t),z(t), y(t)) converge exponentially to zero.

9.
Simulation. To validate our convergence result, we simulated the dynamics (9) with the initial values x(0) = 2 and S(0) = 1 and the reference trajectory x r (t), using the parameters m = 10 and a = 1 2 and t o = 0. In this case, the lower bound on D(t) provided by (8) is D o = 7/3. It follows from Remark 9 that the convergence of x(t) to x r (t) is robust to disturbances that are valued in [−ū,ū] 2 for any positive constantū < min{1, D o } = 1, in the sense of integral input-to-state stability. Moreover, using the estimate (28), one easily checks that in this case, estimate (iISS) on p.14 above holds with δ 2 (r) = 2C 2 r; cf. the proof of [3, Theorem 1]. For our simulation, we took the disturbance u 1 (t) = 0.5e −t on the dilution rate and u 2 (t) ≡ 0. This gave the plot of x(t) and x r (t) against time in Figure 3. Our simulation shows that the state trajectory x(t) closely tracks the reference trajectory x r (t) even in the presence of small disturbances and so validates our findings.

10.
Conclusions. The chemostat model is a useful framework for modeling species competing for nutrients. For the case of one species competing for one nutrient and a suitable time-varying dilution rate, we proved stability of an appropriate reference trajectory. Moreover, we found that the stability was maintained even if the model is augmented with other species that are being driven to extinction, or if there are disturbances of appropriately small magnitude acting on the dilution rate and input nutrient concentration. Appendix. For completeness, we provide the slight variant of the classical ISS arguments needed to finish the proof of Theorem 3. Multiplying through (32) by e C 5 l and applying the standard "variation of parameters" formula to [t o , t] ∋ l → V (ξ(l),z(l)) (by integrating between t 0 ≥ 0 and t ≥ t o ) gives V (ξ(t),z(t)) ≤ e (t 0 −t)C 5 V (ξ(t 0 ),z(t 0 )) + C 2 |u| [to,t] , where we enlarged C 2 without relabeling. We deduce that L 1 (ξ(t)) + κ D o −ūz 2 (t) ≤ ln 1 + e (t 0 −t)C 5 V (ξ(t 0 ),z(t 0 )) + C 2 |u| [to,t] where L 1 is defined in (19). Since e r − 1 − r ≥ 1 2 r 2 and ln(1 + r) ≤ r for all r ≥ 0, we deduce from the formula for V that 1 2ξ 2 (t) + κ D o −ūz 2 (t) ≤ e (t 0 −t)C 5 Ω |(ξ(t 0 ),z(t 0 ))| + C 2 |u| [to,t] where Ω(r) = e e r −1−r+ κ Do−ū r 2 − 1 .
The desired ISS estimate (11)  This completes the proof of Theorem 3.