Study of the chemostat model with non-monotonic growth under random disturbances on the removal rate

Abstract: We revisit the chemostat model with Haldane growth function, here subject to bounded random disturbances on the input flow rate, as often met in biotechnological or waste-water industry. We prove existence and uniqueness of global positive solution of the random dynamics and existence of absorbing and attracting sets that are independent of the realizations of the noise. We study the longtime behavior of the random dynamics in terms of attracting sets, and provide first conditions under which biomass extinction cannot be avoided. We prove conditions for weak and strong persistence of the microbial species and provide lower bounds for the biomass concentration, as a relevant information for practitioners. The theoretical results are illustrated with numerical simulations.


Introduction
The chemostat refers to a laboratory device used for the growth of micro-organisms in a culture environment [1,2], that has been regarded as an idealization of the nature to study microbial ecosystems in stationary stage [3]. It turned out to be an important investigation field due to a large number of applications, especially in waste water treatment [4,5] but also in ecological and environmental sciences (see [6][7][8][9][10]).
It is worth mentioning that the chemostat has been subject to a large number of scientific publications and books, not only in biology and ecology but also in mathematics. Indeed there exists a specific research area about the so-called "theory of the chemostat" [11,12] where many researchers have been involved in the last years. This interest has been strengthened by the fact that the chemostat device can be mathematically modeled in a simple way which reproduces quite faithfully real bio-processes.
Let us recall quickly in what consists the original chemostat device. It is composed of three tanks, the feed bottle, the culture vessel and the collection vessel, which are interconnected by pumps (see Figure 1). The substrate is stored in the feed bottle and provided with a given flow rate to the culture vessel, where interactions between nutrient and microbial biomass take place. The media from the culture vessel is also withdrawn towards the collection vessel with the same flow rate, to keep the volume in the culture vessel constant. where s = s(t) and x = x(t) denote the concentration of substrate and species, respectively, s in is the input concentration of nutrient, D is the removal rate (also called input flow rate), and µ(·) is the specific growth function describing the kinetics of the nutrient consumption by the bacterial species. Here, we assume that the yield coefficient of the conversion of the substrate into biomass is equal to 1 (that is always possible to impose by a change of the unit of the biomass concentration). More specifically, we consider in the present work the Haldane growth function where k s is the affinity constant and k i a parameter modeling the growth inhibition under large concentrations of substrate. Since it will be useful later, let us define the number s m := arg max µ(s) = k i k s . (1.4) Many works have been dedicated to this classical deterministic model of the chemostat (see for instance [11,12]) but most of time the removal rate D is kept constant, although it is well-known that in practice it is frequently subject to disturbances (see for instance [13] where some chronicles of time varying removal rates are depicted).
Motivated by this fact, it has been proposed to model the perturbations on the input flow rate in the chemostat model (1.1)-(1.2) by D + Φ(z * (θ t ω)), where z * (θ t ω) denotes the Ornstein-Uhlenbeck process (introduced in more detail in Section 2) and Φ is a bounded function defined as where d > 0 (see, for instance, [14,15] for other possible functions). Indeed, any odd function Φ which is increasing and such that lim z→+∞ Φ(z) = d can be considered. In this way, once practitioners provide an interval [D l , D r ] ⊂ R, 0 < D l < D < D r < ∞, (typically obtained from observations) we can define d = D r − D = D − D l and then the perturbed input flow is bounded for every time and any realization of the noise, i.e., (1.6) Then, the random chemostat model becomes One may wonder the reason why we consider this way of modeling bounded random fluctuations to perturb the input flow D in the model (1.1)-(1.2), instead of considering other stochastic process, such as for instance the well-known standard Wiener process. Indeed, this way has been typically used to model real perturbations and provides several advantages from both mathematical and biological points of view, see for instance [15].
On the one hand, this way of modeling noise fits in a loyal way the bounded variations of the input flow rate observed in real life. On the contrary, the Wiener process is unbounded with probability one which leads to arbitrary large (possibly negative) values of the corresponding perturbed input flow rate, which is not realistic at all from the biological point of view. We refer readers to [15,16] where the authors explain the relevant drawbacks found when perturbing the input flow in the classical deterministic chemostat (1.1)-(1.2) with a Wiener process (where µ is a Monod growth function).
Moreover, the approach proposed in this paper allows to prove the persistence of the bacterial species (under some conditions on the growth function), as it is observed by practitioners on very long time periods despite variations of the input flow rate. This is not the case when considering the Wiener process where persistence cannot be ensured (see [16] and [17] where the Wiener process is used to model disturbances on the input flow and environmental perturbations in the classical deterministic chemostat).
In industrial setup, large concentrations of the input substrate s in can be observed and it is also wellknown that bacterial species may suffer from growth inhibition under very large concentrations s. The non-monotonic growth function (1.3) precisely models this fact (see [18]). Differently to the classical case, where the growth function µ is assumed to be increasing (as this is the case for the Monod function), the dynamics of the deterministic chemostat model with such non-monotonic growth function may exhibit a bi-stability for certain values of the dilution rate D (see for instance [11]). Depending on the initial condition, the state of the system converges asymptotically to the wash-out of the biomass (which is not a desirable state) or to a positive equilibrium. This kind of instability is observed in practice and present an issue in industrial applications because it requires a good monitoring of the system to detect if the state belongs the attraction basin of the wash-out equilibrium [19][20][21][22]. Most of the time, practitioners prefer to size the process to avoid such a behavior, i.e. such that the system admits an unique globally stable equilibrium (see also [23,24]). This mathematically amounts to have the following condition D < min(µ(s in ), µ(s m )) (1.9) (this result is recalled later). The purpose of the present work is to study the behavior of the random dynamics when the constant removal rate is replaced by D + Φ(z * (θ t ω)) where D satisfies condition (1.9).
The realizations of this variable may satisfy or not condition (1.9) at some times t. If not, one may wonder if this could lead the biomass to extinction. This question is of primer importance for the practitioners for the good health of the bio-process. In other words, for a nominal removal rate that satisfies condition (1.9), is the persistence of the biomass always guaranteed, even when the realizations of the noise provide effective values of the removal rate that do not satisfy this condition? Precise definitions of persistence in the framework of the chemostat will be given later.
Let us underline that stochastic modeling of the chemostat has received a great attention in the literature, considering different kinds of demographic noise [25][26][27][28][29][30], but few works have dealt with noise on the input, and much less in the case of considering Haldane consumption function, whereas this is quite natural for an open system as the chemostat, which is often the main source of fluctuations.
The paper is organized in the following way: in Section 2 we provide preliminaries and classical results about the deterministic chemostat model (1.1)-(1.2). In Section 3 we study the properties of the solutions of the random chemostat model with Haldane consumption kinetics (1.7)-(1.8). Then, in Section 4 we give conditions to ensure both uniform weak and strong persistence of the species. In Section 5, several simulations are presented to support the theoretical study. Finally, we present some conclusions in Section 6.

Preliminaries
In this section we recall briefly some results that are useful in this paper. For the sake of clarity, we split this preliminary section in three different parts: the first one concerns classical results of the deterministic chemostat model (1.1)-(1.2), the second one is about the Ornstein-Uhlenbeck process and the third one recalls the definitions of persistence.

The deterministic chemostat model
The next proposition recalls the classical results about the chemostat model (1.1)-(1.2) when µ is a non-monotonic function. We refer readers to [11,12] for proofs and more details.
Proposition 2.1. Assume that there existsŝ ∈ (0, s in ) such that the function µ is increasing on (0,ŝ) and decreasing on (ŝ, s in ). Define the break-even concentrations λ − (D), λ + (D) as follows 2) presents a bi-stability between E − and E 0 . From any initial condition in R + × R + excepted on a set of null measure, the solution converges asymptotically to E − or E 0 .

Remark 1.
In practice, only the second case is desirable because it guarantees that in any situation the wash-out of the biomass is avoided.

Remark 2. For the Haldane expression (1.3), one has explicit expressions of the functions λ
Let us recall that the concept of break-even concentrations has been revisited in the context of stochastic models of the chemostat [31, 32] but we will not need it here. We keep the classical deterministic definition.
In the rest of the paper, we shall consider that we are in conditions of Proposition 2.1, that is µ non-monotonic on the interval [0, s in ] (otherwise the analysis is similar to monotonic growth function and cannot present bi-stability). Throughout the paper we shall consider the following hypothesis Assumption 1. There exists s m ∈ (0, s in ) such that the function µ is increasing on (0, s m ) and decreasing on (s m , s in ).

The Ornstein-Uhlenbeck process
We present here briefly the Ornstein-Uhlenbeck (OU) process. For more details we refer readers to [13,33,34].
The Ornstein-Uhlenbeck (OU) process is a stationary mean-reverting Gaussian stochastic process defined as where ω denotes a standard Wiener process in a probability space (Ω, F , P), β is the mean reversion constant representing the strength with which the process is attracted by the mean, ν > 0 is the volatility constant describing the variation or the size of the noise and θ t denotes the usual Wiener shift flow given by We note that the OU process (2.1) can be obtained as the stationary solution of the Langevin equation dz + βzdt = νdω.
Typically, the OU process (2.1) can model the position of a particle by taking into account its friction in a fluid (which is the main difference with the typical standard Wiener process). Indeed, it can be considered as a generalization of the standard Wiener process and provides a link between the standard Wiener process (β = 0, ν = 1) and no noise at all (β = 1, ν = 0).
From now on we consider β and ν fixed and z * (θ t ω) the OU process defined above. We recall in the next proposition some of its properties.

Persistence in the chemostat
We recall here the definitions of (uniform) persistence (see for instance [37]) that we consider in the present work .

Properties of the solutions of the random dynamics
In this section we study the random chemostat model (1.7)-(1.8) presented in the introduction. We prove the existence and uniqueness of a global positive solution and provide results about the existence of absorbing and attracting sets which, in addition, are deterministic (i.e. that do not depend on the realization of the noise). In addition, we derive first conditions under which extinction of species cannot be avoided.
Proof. Let us first write system (1.7)-(1.8) as and whence F is locally Lipschitz respect to u ∈ X. Therefore, for each realization of the noise, du/dt = L(θ t ω)u + F(u, θ t ω) is an non-autonomous differential equation with a right member Lipschitz with respect to u and continuous with respect to t (recall that z * (θ t ω) is continuous respect to t and Φ is also continuous). Therefore, the solution of the Cauchy problem admits an unique local solution of system (1.7)-(1.8) from the theory of ordinary differential equations.
Let us first show that both x and s remain in X for any u 0 ∈ X. To this end, we notice first that x(t) = 0 solves (1.7). By uniqueness of the Cauchy problem, we deduce that any other solution is such that x(t) 0 for any t. Notice also that one has, thanks to (1.6), which proves that the axis s = 0 is repulsive in X. This demonstrates the positiveness of the solution of system (1.7)-(1.8). Now we prove that both the substrate and the microorganisms concentrations remain bounded for every time. To this end, define v(t) := s(t) + x(t). Then the variable v is solution of the following differential equation and, thanks to (1.6), dv dt ≤ D r s in − D l v, By comparison of solutions of scalar ODEs (see [38]) we obtain Then, v is forward bounded and since v = s + x ≥ 0, we deduce that both the solution s and  Proof. Define the variable q(t) := s(t) − s in + x(t). Then q satisfies the following differential equation From (1.6), we have q(0)e −D r t ≤ q(t; 0, ω, q(0)) ≤ q(0)e −D l t and taking limit when t goes to infinity we have lim t→+∞ q(t; 0, ω, q(0)) = 0, which means that, for any ε > 0, there exists T (ω, ε) > 0 such that This proves that is a (forward) deterministic absorbing set for system (1.7)-(1.8). As ε > 0 is arbitrary, the proof is done.
Remark 3. Let us underline that the attracting set (3.1) obtained in Theorem 3.2 does not depend on the event ω ∈ Ω. This is another particularity that we obtain when considering this way of modeling random bounded realizations.
Our aim now is to provide first conditions on the parameters of system (1.7)-(1.8) under which extinction of the species occurs.

Conditions for persistence
In this section our aim is to provide conditions to ensure persistence of the species in the sense defined in Section 2.3.
Hence, it is enough to show that E[Φ(z * (ω))] = 0. In fact, we have where f OU denotes the density function of the random variable z * (·), which is an even function since it is gaussian, and Φ is an odd function.
Remark 4. From Proposition 4.1 we notice that we could model random perturbations as in this paper, i.e., by means of Φ(z * (θ t ω)), and the ergodic property (4.1) remains true as long as Φ is an odd measurable function. is fulfilled. Then, the random chemostat model (1.7)-(1.8) is weakly persistent, that is there exits ε > 0 such that for any initial condition x(0) ∈ X, any realization satisfies lim sup t→+∞ x(t; 0, ω, x(0)) ≥ ε.
Then, from (1.8), one can write the following inequalitẏ and by integration between T (ω, ε) and t > T (ω, ε) one obtain Remark 5. Theorem 4.1 proves the weak (uniform) persistence of species as long as D < µ(s in ) is fulfilled, the same condition that guarantees the persistence in the deterministic case (see Proposition 2.1). However, let us underline that we do not impose the upper bound D r of the variations of the removal rate to fulfill this inequality. This means that one could have realizations of the disturbances such that the effective value of the removal rate is above µ(s in ) on large periods of time making the species to be arbitrary closed to the extinction but it will always persist.
In the next result we consider a stronger condition on the upper bound D r of the removal rate which ensures the strong persistence of the species. is fulfilled. Then, the random chemostat model (1.7)-(1.8) is strongly persistent and the set Proof. Take ε > 0 such that µ(s) > D r for any s ∈ [s in − ε, s in + ε] and posit From now on, we will s(t), x(t) and q(t) instead of s(t; 0, ω, s(0)), x(t; 0, ω, x(0)) and q(t; 0, ω, q(0)) to make the readability easier even though we recall that every state variable depends on the noise. From Theorem 3.2, we know that the variable q(t) = s(t) − s in + x(t) converges asymptotically to zero for any realization, and thus any solution s(t) of the random model (1.7)-(1.8) converges to the set [0, s in ]. Consequently, there exists T (ω, ε) > 0 such that If s(t) ∈ [s in − ε, s in + ε] for any t > T (ω, ε), then one has, from (1.8),ẋ(t) > ηx(t) for any t > T (ω, ε), which implies that x is unbounded, thus a contradiction. We deduce that there exists a finite timeT (ω) ≥ T (ω, ε) such that s(T (ω)) ≤ s in − ε/2. On another hand, from (1.7), s(t) can be written as the solution of the non-autonomous dynamics Note that one has and then (from the definition of κ). We deduce that the set [0, s in − ε/2] is forward invariant for the semi-flow {σ = F(t, σ), t > T (ω, ε)}. Therefore, one has Then, once identified σ with s, from (4.5) and from the comparison of solutions of scalar ordinary differential equations [38], one has the inequality s(t) ≤ s + (t) for any t >T (ω), where s + (t) is solution of the Cauchy problem Note that the solution s + (t) belongs to the interval [0, s in ] for any t ≥T (ω) (and is thus bounded) and that its dynamics is asymptotic autonomous with limiting dynamics Under assumption D r < µ(s in ), one has necessarily λ − (D r ) < s in and the property is fulfilled. One finally obtains that any solution of (4.6) in [0, s in ] is such that s † (t) → λ − (D r ) when t → +∞. From the theory of asymptotically autonomous dynamical systems [39], one concludes that s + (t) converges also to λ − (D r ) when t → +∞, which proves that one has Now we consider the last situation when µ(s in ) < D < µ(s m ) which corresponds to the bi-stability in the deterministic case. In the random framework, one cannot guarantee persistence nor wash-out of the biomass. As we shall see later on simulations, the asymptotic behavior of the solutions depends on the initial condition and the realization of the noise. However we show that an upper bound on the biomass can be provided.  Proof. Let us remark that in this proof we will write x(t), s(t) and q(t), or simply x, s and q, instead of x(t; 0, ω, x(0)), s(t; 0, ω, s(0)) and q(t; 0, ω, q(0)) for the sake of simplicity.

Numerical simulations
In this section we present several numerical simulations to support our theoretical results in three different cases: extinction, strong persistence and weak persistence of species.
Two different figures are presented in each case. The first figure (with two panels) concerns the dynamics of the substrate (top) and the species (bottom) in the random chemostat model (1.7)-(1.8) for different realizations of the noise (in colored continuous lines) along with the case without noise i.e. the deterministic chemostat model (in blue dashed lines). In addition, in each panel we display a little box with a zoom to illustrate better the dynamics around the attracting set. The second figure offers information about the consumption function of the consumer species and those parameters involved in conditions to have extinction, weak persistence and strong persistence of species.
We first present the case of extinction in Figure 2  This is indeed quite intuitive as condition (3.3) basically means that the nominal value D of the removal rate is too large large compared to the growth kinetics (recall that D l > µ(s m ) implies D > µ(s m )), even though we start far away from the washout (s in , 0).   In this case D r < µ(s in ) holds true (see Figure 10) and then strong persistence of species is obtained from Theorem 4.1.    Moreover, we would like to notice that OU process proves again to be a powerful tool when modeling real (bounded) noises. In addition, it allows us to observe clearly the difference between weak and strong persistence.
Finally we present some numerical simulations to observe that both extinction and persistence can be obtained when D l < µ(s in ) < D is fulfilled, as in the deterministic case. However, the issue here depends also on the realization of the noise.

Conclusion
We have considered the chemostat model (1.1)-(1.2) with Haldane consumption kinetics, under bounded perturbations on the input flow rate, motivated by real cases in industrial setup and biotechnology. To this end, we use a bounded saturated function of the Ornstein-Uhlenbeck process that allows us to ensure the noise to be bounded in realistic intervals.
We prove in Theorem 3.1 existence and uniqueness of global positive solution of the corresponding random chemostat (1.7)-(1.8) by means of standard results from the theory of ODEs and thanks to properties of the Ornstein-Uhlenbeck process. Then, in Theorem 3.2 we establish the existence of an absorbing and attracting set which has the nice property to be deterministic, i.e. that does not depend on the realization of the noise.
We then focused on the long-time behavior of the random dynamics inside this attracting set. To this end, we first proved in Theorem 3.3 that extinction of species cannot be avoided as long as D l > µ(s m ) whatever is the input concentration s in . On the opposite, we proved in Theorem 4.1 the weakly uniform persistence of species when D < µ(s in ), which means that species can be temporarily arbitrary closed to the extinction but still persist even when having random disturbances in the input flow. The condition D < µ(s in ) ensures persistence in the deterministic case but the effective removal rate does not necessarily fulfill this condition depending on the realizations. Finally, we prove the strong persistence of the species in Theorem 4.2 under the stronger condition D r < µ(s in ). In this case, we provide an explicit lower bound for the asymptotic concentration of the species, a useful information for practitioners.
In addition, we support the theoretical results with several numerical simulations which depict the possible behaviors of the random dynamics. Moreover, this allows us to illustrate the difference between weak and strong persistence: once fixed every parameter we change the removal rate D to let conditions µ(s in ) > D and µ(s in ) > D r be true or not. Condition to have weak persistence (µ(s in ) > D) does not involve the noise and then the random realizations of the species can fluctuate closed to zero and we still ensure persistence whereas condition to have strong persistence (µ(s in ) > D r ) only fulfills for realizations of the perturbed removal rate inside a small enough interval.
To conclude, we would like to stress that the Ornstein-Uhlenbeck process has proved once again to be relevant tool to model random disturbances in a biological framework. This stochastic process fits in a quite loyal way the bounded fluctuations that are observed in practice, and justify the (weak or strong) persistence of the biomass despite the possible realizations of noise, as also observed in practice.
c 2020 the Author(s), licensee AIMS Press. This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0)