Bounded random fluctuations on the input flow in chemostat models with wall growth and non-monotonic kinetics

http://www.aimspress.com/journal/Math AIMS Mathematics, 6(4): 4025–4052. DOI: 10.3934/math.2021239 Received: 17 November 2020 Accepted: 19 January 2021 Published: 04 February 2021 Research article Bounded random fluctuations on the input flow in chemostat models with wall growth and non-monotonic kinetics Tom´ as Caraballo1and Javier L´ opez-de-la-Cruz2,* 1Departamento de Ecuaciones Diferenciales y An´ alisis Num´ erico, C/Tarfia s/n, Facultad de Matem´ aticas, Universidad de Sevilla, 41012 Sevilla, Spain 2Departamento de Matem´ atica Aplicada a las TIC, Escuela T´ ecnica Superior de Ingenieros Inform´ aticos, Campus de Montegancedo, Universidad Polit´ ecnica de Madrid, 28660 Boadilla del Monte, Madrid, Spain *Correspondence: Email: javier[email protected]. Dedicated to the memory of Mar´ıa Jos´e Garrido Atienza Abstract: This paper investigates a chemostat model with wall growth and Haldane consumption kinetics. In addition, bounded random fluctuations on the input flow, which are modeled by means of the well-known Ornstein-Uhlenbeck process, are considered to obtain a much more realistic model fitting in a better way the phenomena observed by practitioners in real devices. Once the existence and uniqueness of global positive solution has been proved, as well as the existence of deterministic absorbing and attracting sets, the random dynamics inside the attracting set is studied in detail to provide conditions under which persistence of species is ensured, the main goal pursued from the practical point of view. Finally, we support the theoretical results with several numerical simulations. Keywords: chemostat; wall growth; non-monotonic kinetics; absorbing set; real noise Mathematics Subject Classification: 34F05, 34K99 1. Introduction Chemostat refers to a laboratory device, invented at the same time in the 1950s by Monod (see [1]) and Novick and Szilard (see [2]), very used in practice for growing microorganisms in a cultured environment. The chemostat has been regarded as an idealization of nature to study microbial ecosystems at steady state, which is an important and interesting problem for researchers from many areas due to its large number of applications, for instance, it can be used to study genetically altered microorganisms (see e.g. [3, 4]), waste-water treatment process (see e.g. [5, 6]), models of mammalian large intestine 4026 (see e.g. [7,8]) and plays an important role in theoretical ecology (see e.g. [9–20]). It is worth mentioning that the chemostat has been subject to a large number of scientific publications, not only in applied sciences but also in Mathematics. It is due to the fact that it is very interesting also as a mathematical object, in fact, it constitutes a very active branch of applied mathematics which, in addition, proposes a recent formal framework called the theory of the chemostat centered around a well-identified community of mathematicians. The reputation of the chemostat is mainly due to the fact that it is a very simple device that allows to develop many different works and reproduces the real experiments in a very loyal way. Those are some of the reasons that encouraged us to consider this model as part of our research. Concerning the biological aspects, the chemostat device consists of three tanks, the feed bottle, the culture vessel and the collection vessel, that are interconnected by means of pumps. The substrate stored in the feed bottle is pumped to the culture vessel, where the interactions between the substrate and the species take place. In addition, another flow is pumped from the culture vessel to the collection vessel in order to keep the volume of the culture vessel constant. The classical mathematical model used to describe the dynamics into the culture vessel in the chemostat device consists of a system of two ordinary differential equations (see e.g. [21,22]) given by ds dt =D(sin −s)−µ(s)x,(1.1) dx dt =−Dx +µ(s)x,(1.2) where s=s(t) and x=x(t) denote the concentration of substrate (or nutrient) and species (or microorganisms), respectively, sin is the input concentration of substrate, Dis the removal rate, also called input flow, and µis a function describing the consumption kinetics of the consumer species. Among the different assumptions made when dealing with the deterministic model (1.1)–(1.2), one is to assume that the input flow Dis fast enough such that species are removed to the collection vessel before having the opportunity of sticking to the wall of the culture vessel. However, this phenomenon, introduced by Pilyugin and Waltman (see [23]) and known as wall growth, is quite usual and it can be easily observed in practice. In fact, these observations motivate us to go further from the classical deterministic chemostat model (1.1)–(1.2) and consider the corresponding chemostat model with wall growth given by the following differential system (see [22]) ds dt =D(sin −s)−µ(s) yx1−µ(s) yx2+bνx1,(1.3) dx1 dt =−(ν+D)x1+µ(s)x1−r1x1+r2x2,(1.4) dx2 dt =−νx2+µ(s)x2+r1x1−r2x2,(1.5) where s=s(t) denotes again the concentration of the nutrient and we divide in this case the total concentration of species x=x(t) as follows: on the one hand, x1=x1(t) refers to the concentration of the species in the liquid media (the also called planktonic biomass) and x2=x2(t) is the concentration of the species sticked on the wall of the culture vessel (the attached biomass). AIMS Mathematics Volume 6, Issue 4, 4025–4052. 4027 In addition, the rate b∈(0,1) describes the fraction of dead biomass which is recycled and ν > 0 is the collective death rate coefficient of the microbial biomass representing factors such as diseases, aging, etc. Apart from that, r1>0 and r2>0 represent the rates at which the species stick on to and shear offfrom the walls of the culture vessel, respectively. Finally µis, as in the classical chemostat, the consumption function of the consumer species and y>0 denotes the yield coefficient of transformation of substrate into biomass. As explained before, we consider Haldane kinetics in this paper, a non-monotonic function introduced by Haldane (see [24]) and defined by µ(s)=¯µ0s s+ks+1 ki s2 ,(1.6) where ¯µ0corresponds to the growth rate coefficient of the species, ksdenotes the so-called Monod constant and kiis the inhibition constant modeling the inhibition of species to take substrate when it is at high concentration. We would also like to highlight that, in industrial setup, large input concentrations sin can be observed and it is also very well-known that microorganisms may suffer growth inhibition under high concentrations of substrate. This is precisely what the non-monotonic consumption function (1.6) models (see [25]). There are many papers in the literature dedicated to investigate different versions of the deterministic classical chemostat model, either (1.1)–(1.2) or (1.3)–(1.5), even though a few of them take into account the wall growth, see e.g. [14, 21, 22, 26–29]. Nevertheless, most of them need to assume restrictions that are quite strong from the biological point of view. For instance, a very representative example of this is the case of the input flow D, that is typically considered to be constant in spite of the fact that it is often subject to suffer bounded random fluctuations in real devices (see Figure 1 were we display the evolution of the input flow in a real experiment). 0 10 20 30 40 50 60 70 80 90 100 0.1 0.15 0.2 0.25 0.3 Figure 1. Evolution of input flow in real chemostat. Motivated by the previous facts, our aim in this work consists in modeling and investigating the dynamics of the deterministic classical chemostat with wall growth (1.3)–(1.5) and Haldane consumption kinetics where the input flow is subject to random disturbances. To this end, since the perturbed input flow should be bounded (as observed in real devices) we propose to model the random fluctuations on AIMS Mathematics Volume 6, Issue 4, 4025–4052. 4028 the input flow by means of the Ornstein-Uhlenbeck process, which has proved to be a useful tool when modeling bounded noises (see [30–34] for more details about modeling bounded noises). To this end, we replace the input flow Din the deterministic system (1.3)–(1.5) by D+ Φ(z∗(θtω)), where z∗(θtω) denotes the Ornstein-Uhlenbeck process (which will be introduced in the next section) and Φ:R→[−d,d], where d>0, is a function given by Φ(z)=2d πarctan(z).(1.7) Hence, the resulting random chemostat is given by the following random differential system ds dt =(D+ Φ(z∗(θtω)))(sin −s)−µ0s s+ks+1 ki s2 x1−µ0s s+ks+1 ki s2 x2+bνx1,(1.8) dx1 dt =−(ν+D+ Φ(z∗(θtω)))x1+¯µ0s s+ks+1 ki s2 x1−r1x1+r2x2,(1.9) dx2 dt =−νx2+¯µ0s s+ks+1 ki s2 x2+r1x1−r2x2,(1.10) where µ0:=¯µ0/yand µ0≥¯µ0. We recall that modeling bounded perturbations is a very interesting topic for most of researchers in applied sciences when trying to find a way to model real noises, since they are random but bounded. Apart from that, as it will be observed in this paper, such bounded perturbations will allow us to ensure the persistence of the species, one of the most important goals pursued by practitioners. In fact, persistence is the object of study in many papers (see [30–36]). The rest of the paper is organized in the following way: in Section 2 we present some brief preliminaries about the classical deterministic chemostat with Haldane kinetics, the Ornstein-Uhlenbeck process, whose properties will be continually used in the paper, and the way in which bounded perturbations are modeled in this work; in Section 3 we investigate the existence and uniqueness of global positive solutions of the random chemostat (1.8)–(1.10), the existence of deterministic absorbing and attracting sets and we provide conditions under which both extinction of species and persistence can be ensured; in Section 4 we display several numerical simulations to support the theoretical work; finally, in Section 5 we present some final comments. 2. Preliminaries In this section we recall briefly some classical results concerning the deterministic chemostat model (1.1)–(1.2) with non-monotonic growth and present the Ornstein-Uhlenbeck (OU) process, providing some results about relevant properties that will be used further and will help us to set up the way to model bounded random perturbations. AIMS Mathematics Volume 6, Issue 4, 4025–4052. 4029 2.1. The deterministic chemostat model with non-monotonic growth The next proposition presents the classical results about the classical chemostat model (1.1)-(1.2) with non-monotonic consumption function µ. We refer every interested reader to [21, 37] for proofs and more details. Proposition 2.1. Assume that there exists ˆs∈(0,sin)such that the function µis increasing on (0,ˆs) and decreasing on (ˆs,sin). Define the break-even concentrations λ−(D),λ+(D)as follows λ−(D)=min{s∈[0,ˆs] ; µ(s)≥D},D∈[0, µ(ˆs)] λ+(D)=max{s∈[ˆs,sin] ; µ(s)≥D},D∈[µ(sin), µ(ˆs)] 1. If D > µ(ˆs), the system (1.1)-(1.2) possesses a unique equilibrium E0:=(0,sin), which is globally asymptotically stable on R2 +. 2. If D < µ(sin), the system (1.1)-(1.2) admits an unique positive equilibrium E−:=(sin − λ−(D), λ−(D)) which is globally asymptotically stable on R? +×R+. 3. If D ∈[µ(sin), µ(ˆs)], the system (1.1)-(1.2) presents a bi-stability between E−and E0. From any initial condition in R? +×R+excepted on a set of null measure, the solution converges asymptotically to E−or E0. 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.6), one has explicit expressions of the functions λ± λ±(D)= (µ0−D)ki±q(µ0−D)2k2 i−4D2kski 2D. Let us note that, even though the concept of break-even concentrations has been revisited in the context of stochastic models of the chemostat [38, 39], 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,sin] (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 sm∈(0,sin) such that the function µis increasing on (0,sm) and decreasing on (sm,sin). 2.2. The Ornstein-Uhlenbeck process The OU process is a stationary mean-reverting Gaussian stochastic process defined as the following random variable z(t, ω) :=z∗(θtω)=−βγ 0 Z −∞ eβsθtω(s)ds,t∈R, ω ∈Ω, β, γ > 0,(2.1) AIMS Mathematics Volume 6, Issue 4, 4025–4052. 4030 where ωdenotes a standard Wiener process in a probability space, namely (Ω,F,P), β > 0 is the mean reverting constant describing how our system reacts under perturbations, γ > 0 is the volatility constant representing the variation of the noise and θtdenotes the usual Wiener shift flow given by θtω(·)=ω(·+t)−ω(t),t∈R. We notice that the OU process (2.1) can be obtained as the stationary solution of the Langevin equation dz +βzdt =γdω. (2.2) From now on, we consider βand γfixed. Once presented the OU process, we state in the following proposition some essential properties which will be crucial to study this paper. Proposition 2.2 (See [40]).There exists a θt-invariant subset e Ω∈ F of Ωof full P−measure such that for ω∈e Ωand β, γ > 0, we have 1. the random variable |z∗(ω)|is tempered with respect to {θt}t∈R, i.e., for a.e. ω∈e Ω, lim t→∞ e−ηtsup t∈R|z∗(θ−tω)|=0,for all η > 0. . 2. the mapping (t, ω)→z(t, ω) :=z∗(θtω)=−βγ 0 Z −∞ eβsω(t+s)ds+ω(t) is a stationary solution of (2.2) with continuous trajectories; 3. for any ω∈˜ Ωone has lim t→±∞ |z∗(θtω)| t=0; lim t→±∞ 1 tZt 0 z∗(θsω)ds =0; lim t→±∞ 1 tZt 0|z∗(θsω)|ds =E[|z∗|]<∞; From now on, we will denote e Ωby Ω. 2.3. Modeling bounded random perturbations Once introduced the OU process, our aim in this section is to present a way to model mathematically bounded random perturbations that fit the real disturbances observed in labs. To this end, assume we are given an interval [Dl,Dr] (this interval is typically provided by practitioners from observation), 0 <Dl<D<Dr<∞and let us define d=Dr−D=D−Dl>0. Consider now the function Φ:R→[−d,d] defined by Φ(z)=2d πarctan(z), AIMS Mathematics Volume 6, Issue 4, 4025–4052. 4031 as in (1.7). Thus, Dl<D+ Φ(z∗(θtω)) <Drfor all t∈R,(2.3) whence the random perturbations in system (1.8)–(1.10) remain bounded (as the real ones) in the desired interval. In addition, the random perturbations in system (1.8)–(1.10) have the following ergodic property. Proposition 2.3. Let Φ(z)=2d πarctan(z)be defined as in (1.7). Then lim t→+∞ 1 tZt 0 Φ(z∗(θsω))ds =0,a.s. in Ω.(2.4) Proof. From the definition of Φit yields that ZΩ|Φ(z∗(ω))|dP(ω)=ZΩ 2d πarctan(z∗(ω))dP(ω)≤d|Ω|=d whence Φ(z∗(·)) ∈L1(Ω,F,P). Therefore, since Pis invariant by θ(see [26, 41]), from the Birkhoff ergodic theorem we obtain lim t→+∞ 1 tZt 0 Φ(z∗(θsω))ds =E[Φ(z∗(ω))],a.s. in Ω. Then, it remains to prove that E[Φ(z∗(ω))] =0. In fact, E[Φ(z∗(ω))] =Z+∞ −∞ Φ(x)fOU(x)dx =0, where fOU denotes the density function of the random variable z∗(·), which is an even function since it is gaussian, and Φis an odd function.  Remark 3. From the proof of Proposition 2.3 we notice that we could consider any odd measurable function Φsuch that lim z→+∞Φ(z)=d>0. Remark 4. We would like to point out that the resulting random system (1.8)–(1.10) generates a random dynamical system and then we could use the theory of random dynamical systems and pullback attractors (see [26]). Nevertheless, we will study the solutions of the random system (1.8)–(1.10) for fixed events ω∈Ωsince, thanks to this, every result in this work can be proved in forward sense, which is more natural from the biological point of view than the pullback one. We refer readers to [30, 31, 33, 34] for more details about modeling bounded perturbations with the OU process not only in chemostat models but also when dealing with other models in population dynamics. Moreover, readers can compare the results in those works with the ones in [42,43] where the usual Wiener process was used to model stochastic perturbations in different versions of the chemostat model. AIMS Mathematics Volume 6, Issue 4, 4025–4052. 4032 3. Random chemostat with wall growth and Haldane consumption kinetics In this section we investigate the random chemostat model with wall growth (1.8)–(1.10) presented in the introductory section. We first prove that there exists a unique global positive solution of system (1.8)–(1.10). After that, by studying the asymptotic behavior of system (1.8)–(1.10), we prove the existence of absorbing and attracting sets (forwards in time) that do not depend on the noise. Then, in order to have more detailed information about the random dynamics inside the attracting set, we go further and investigate the long-time behavior of the state variables involved in the model. Thanks to this fact, we provide conditions under which species go to extinction and also conditions to ensure the persistence of species. 3.1. Existence and uniqueness of global positive solution In this section we prove that system (1.8)-(1.10) has a unique global positive solution that remains positive for every positive initial condition. From now on X=n(s,x1,x2)∈R3:s≥0,x1≥0,x2≥0odenotes the positive cone in the three dimensional space. Theorem 3.1. For every initial condition u0:=(s0,x10,x20)∈ X, random system (1.8)-(1.10) possesses a unique global solution u(·; 0, ω, u0)=(s(·; 0, ω, s0),x1(·; 0, ω, x10),x2(·; 0, ω, x20)) ∈ C1([0,+∞); X) which remains in Xfor every t >0. Proof. Let us write our random system (1.8)-(1.10) as du dt =L(θtω)u+F(u, θtω),(3.1) where*u=(s,x1,x2)t, L(θtω)= −(D+ Φ(z∗(θtω))) bν0 0−(ν+D+ Φ(z∗(θtω)) −r1)r2 0r1−ν−r2  and F(u, ω)=(sin(D+ Φ(z∗(θtω))) −(µ(s)/y)(x1+x2), µ(s)x1, µ(s)x2)t. On the one hand, from Proposition 2.2 (ii) one has that t7→ z(t, ω) :=z∗(θtω) is continuous. Apart from that, it is easy to check that F(·, ω)∈ C1(X) (recall that µis given by (1.6)) whence we deduce that Fis locally Lipschitz with respect to the state variables (s,x1,x2). Hence, the right side of (3.1) is continuous respect to the time and locally Lipschitz respect to uand, from classical arguments from the theory of ODEs, it yields that the random system (1.8)–(1.10) possesses a unique local solution. Let us now define the new state variable y=s+x1+x2as the sum of the substrate and both species concentrations. By differentiation, one has that ysatisfies the following random differential inequality dy dt ≤Drsin −my,(3.2) *We note that tdenotes the transposed of the vector. AIMS Mathematics Volume 6, Issue 4, 4025–4052. 4033 where m=min{Dl, ν}and we used (2.3) and the fact that b<1. Then, by solving (3.2), we obtain y(t; 0, ω, y0)≤y0e−mt +Drsin m1−e−mt(3.3) for all t>0, any ω∈Ωand every initial condition y0=s0+x10 +x20 >0. Thus, y=s+x1+x2is bounded from above (and from below by zero, as we will prove next). As a consequence, one can define the unique local solution of the random system (1.8)-(1.10) on the whole interval t∈(0,+∞) and then it is in fact a global one. It remains to prove that every state variable remains in Xfor any initial condition u0=(s0,x10,x20)∈ X. To this end, let us consider in the rest of the proof e s≥0, e x1≥0 and e x2≥0. Assume that we reach the plane x2=0 at some point (e s,e x1,0). Since our initial condition is in the positive cone Xthen the rest of the state variables remain positive. By computing the derivative of x2 at (e s,e x1,0) we obtain dx2 dt (s,x1,x2)=(e s,e x1,0) ≥0 (3.4) whence we deduce that the solution will not cross the plane x2=0 and then x2will always remain non-negative. By similar arguments one also has dx1 dt (s,x1,x2)=(e s,0,e x2) =r2x2≥0,(3.5) when assuming that the solution reaches some point in the plane x1=0 and ds dt (s,x1,x2)=(0,e x1,e x2) =(D+ Φ(z∗(θtω)))sin +bνe x1≥0 (3.6) when reaching s=0. Summing up, s,x1and x2remain non-negative and then the unique global solution of the random system (1.8)–(1.10) is in Xfor every time.  3.2. Existence of absorbing and attracting sets In this section, we study the long-time behavior of the state variables involved in system (1.8)–(1.10) to prove the existence of absorbing and attracting sets (forwards in time). Theorem 3.2. For any ε > 0, the random system (1.8)–(1.10) possesses an absorbing set (forwards in time) given by Bε:=((s,x1,x2)∈ X :s+x1+x2≤Drsin m+ε),(3.7) where m =min{Dl, ν}. As a consequence, A:=((s,x1,x2)∈ X :s+x1+x2≤Drsin m)(3.8) AIMS Mathematics Volume 6, Issue 4, 4025–4052. 4040 3.4. Persistence of species In this section we provide conditions on the parameters involved in system (1.8)–(1.10) to ensure the persistence of species. To this end, we firstly prove that the total concentration of species persists and then we use this fact and some properties proved about the dynamics of the proportion to prove coexistence of each group of species, the planktonic biomass and the attached biomass. Theorem 3.6. Assume ¯µ0 ν+Drprzl>zr+¯µ0ks+(zr)2 ki¯µ0 ,(3.46) where every constant was already defined. Then, the attracting set (3.8) is reduced to A:=((s,x1,x2)∈ X :s+x1+x2≤Drsin m,s≥¯s,x≥¯x),(3.47) where ¯s:=−a1−qa2 1−4a2a0 2a2 and ¯x:=−b1−qb2 1−4b2b0 2b2 ,(3.48) with a2:=−Dr−2¯µ0 µ0 bνpl,a1:=Dlsin −Drks−zr+zlbνpl µ0−2¯µ0 µ0 ksbνpl, a0:=Dlsinks+zlksbνpl µ0 ,b2:=−µ2 0 ki¯µ0 , b1:=−¯µ0µ0 ν+Drpr+µ0+2zrµ0 ki¯µ0 ,b0:=¯µ0zl ν+Drpr−zr−¯µ0ks−(zr)2 ki¯µ0 . Proof. We start proving the lower bound for the concentration of the species. To this end, from (3.11), we have ds dt =(D+ Φ(z∗(θtω)))(sin −s)−µ0s s+ks+1 ki s2 x+bνpx ≥Dlsin −Drs−µ0s s+ks+1 ki s2 1 µ0 (zr+ε−¯µ0s)+bνpl1 µ0 (zl−ε−¯µ0s) ≥Dlsin −Drs−s(zr+ε−¯µ0s) s+ks +1 µ0 bνpl(zl−ε−¯µ0s) (3.49) where we used (2.3) and (3.33). Let us consider ε < ¯µ0¯s, where ¯s>0 such that Dlsin −Dr¯s−¯szr ¯s+ks +1 µ0 bνpl(zl−2¯µ0¯s)≥0.(3.50) AIMS Mathematics Volume 6, Issue 4, 4025–4052. 4041 Hence, from (3.49) one has ds dt ≥Dlsin −Drs−s(zr+¯µ0¯s−¯µ0s) s+ks +1 µ0 bνpl(zl−¯µ0¯s−¯µ0s) (3.51) and, from (3.51), we obtain ds dt s=¯s≥0 (3.52) Then, our goal now is to prove the existence of some ¯s>0 such that (3.50) holds. To do this, by doing some calculations, it is easy to verify that (3.50) is equivalent to a2¯s2+a1¯s+a0≥0,(3.53) where a2,a1and a0are defined as in the statement of Theorem 3.6. Define now the map F:s∈[0,+∞)7→ F(s)∈Rwhere F(s)=a2s2+a1s+a0. Since a2<0, then Fdefines a concave parabolic function. In addition, a0>0 then it follows that there exists s∗>0 such that F(s)≥0 for every s∈[0,s∗] and F(s)<0 when s∈(s∗,+∞). Hence, (3.50) fulfills by considering ¯s=s∗>0, where the value s∗can be easily calculated by solving F(s∗)=0. As a consequence, we find that ¯s=−a1−qa2 1−4a2a0 2a2 ,(3.54) where a0,a1and a2are defined as in the statement of Theorem 3.6. The second part of the proof consists of proving the existence of a lower bound for the total concentration of species. Thanks to this fact, we will be able to ensure that, at least, some of both species will persist. To this end, from (3.12) one has dx dt =−(ν+(D+ Φ(z∗(θtω)))p)x+¯µ0s s+ks+1 ki s2 x ≥−ν−Drpr+¯µ0(zl−ε−µ0x) zr+ε−µ0x+¯µ0ks+(zr)2+ε2+µ2 0x2−2zrµ0x+2zrε ki¯µ0  x,(3.55) where we used (2.3) and (3.46). Let us consider ε<δ¯x, where ¯x>0 such that ¯µ0(zl−δ¯x−µ0¯x) zr+δ¯x−µ0¯x+¯µ0ks+(zr)2+δ2¯x2+µ2 0¯x2−2zrµ0¯x+2zrδ¯x ki¯µ0 ≥ν+Drpr.(3.56) Then, the right hand side of (3.55) can be bounded from below by −ν−Drpr+¯µ0(zl−δ¯x−µ0x) zr+δ¯x−µ0x+¯µ0ks+(zr)2+δ2¯x2+µ2 0x2−2zrµ0x+2zrδ¯x ki¯µ0  x. AIMS Mathematics Volume 6, Issue 4, 4025–4052. 4042 Then, we have dx dt x=¯x≥0, whence ¯x>0 provides a lower bound for the total concentration of species. Hence, our goal now is to prove that there exists ¯x>0 such that (3.56) is fulfilled. To this end, by making some calculations, from (3.56) one obtains b2(δ) ¯x2+b1(δ) ¯x+b0≥0,(3.57) where b2(δ)=− 1 ki¯µ0 (δ2+µ2 0)!,(3.58) b1(δ)=−(µ0+δ)¯µ0 ν+Drpr−δ+µ0+2zrµ0 ki¯µ0−2δzr,(3.59) b0=¯µ0 ν+Drprzl−zr−¯µ0ks−(zr)2 ki¯µ0 .(3.60) Define now the map G:x∈[0,+∞)7→ G(x)∈Rwhere G(x)=b2(δ)x2+b1(δ)x+b0. Since b2(δ)<0 for every δ > 0, then Gdefines a concave parabolic function for every δ > 0. In addition, b0>0 as long as (3.46) is satisfied. Then, for every δ > 0, there exists x∗ δ>0 such that G(x)≥0 for every x∈[0,x∗ δ] and G(x)<0 when (x∗ δ,+∞). Hence, the value x∗ δcan be easily calculated by solving G(x)=0 whence x∗ δ=−b1(δ)−qb2 1(δ)−4b2(δ)b0 2b2(δ),(3.61) where b2(δ), b1(δ) and b0were defined above. Hence, (3.50) fulfills when considering ¯x=lim δ→0x∗ δ>0.  Remark 5. Recall that the chemostat model (1.8)-(1.10) studied in this work does not only incorporate wall growth but also nutrient recycling (coefficient b). Then, one could wonder how this feature is favorable or not to species persistence, i.e., the difference with case b=0. To answer this question, since zr(defined as in (3.26)) is the only parameter in (3.46) (condition to ensure persistence of species) that involves b, and the value of zrdecreases if b=0, we conclude that b=0 would help condition (3.46) to be true, i.e., b=0 would favor the persistence of the species. Finally, we prove that not only the total concentration of species remains positive as long as (3.46) holds true but also both the planktonic biomass and the attached biomass persist. To this end, we use Theorem 3.6 and the dynamics of the proportion. AIMS Mathematics Volume 6, Issue 4, 4025–4052. 4043 Theorem 3.7. Assume that (3.46) holds true. Then, for every initial condition x10 >0and x20 >0, ω∈Ωand t >0large enough, the following lower bounds hold true x1(t; 0, ω, x10)≥¯x1and x2(t; 0, ω, x20)≥¯x2,(3.62) where ¯x1:=pl¯x and ¯x2:=(1 −pr) ¯x,(3.63) i.e., both the planktonic biomass and the attached biomass persist. Proof. From the definition of the proportion (3.10), we have x1(t; 0, ω, x10)=p(t; 0, ω, p0))x(t; 0, ω, x0),(3.64) for any initial condition x10 >0, p0∈(0,1) and x0=x10 +x20 >0. Then, from Theorem 3.3 and assuming (3.46), one can apply Theorem 3.6 and obtain x1(t; 0, ω, x10)≥pl¯x>0 (3.65) for every tlarge enough, any ω∈Ωand every initial condition x10 >0, where pland ¯xare defined in (3.16) and (3.48). Similarly, from the definition of the proportion (3.10) one has x2(t; 0, ω, x20)=(1 −p(t; 0, ω, p0)x(t; 0, ω, x0) (3.66) for every tlarge enough, any ω∈Ωand every initial condition x20 >0, p0∈(0,1) and x0=x10 +x20 > 0. Then, from Theorem 3.3 and assuming (3.46), we can use Theorem 3.6 and one has x2(t; 0, ω, x20)≥(1 −pr) ¯x>0, for every t, any ω∈Ωand every initial condition x20 >0, where prand ¯xare defined in (3.22) and (3.48).  4. Numerical simulations In this section we depict several numerical simulations to support the theoretical study presented in the previous section and help every reader to understand better the dynamics of the random chemostat model with wall growth (1.8)–(1.10). Every figure in this section shows three different panels: the biggest one displays the three dimensional phase plane of system (1.8)–(1.10); the one placed on the right at the top shows the dynamics of the species in the liquid media versus the substrate; finally, the panel on the right at the bottom displays the dynamics of the species sticked on the wall of the culture vessel versus the substrate. In such a way, it will be easy to observe in the little panels when the species go to extinction or persist. Apart from that, in every panel a label next to an arrow is depicted to point the initial condition. In every figure we plot five realizations of the solution of the random chemostat (1.8)–(1.10) with s0=2, x10 =5 and x20 =2 as initial condition. In addition, we set β=1 and γ=1 to simulate the realizations of the OU process. The rest of parameters will be specified further in every case. AIMS Mathematics Volume 6, Issue 4, 4025–4052. 4044 In Figure 2 we consider sin =4, ks=1.7, µ0=8, ki=7, ¯µ0=7.5, b=0.1, ν=0.4, r1=0.4, r2=0.8. Concerning the removal rate, we set D=0.7, Dr=0.8, Dl=0.7 and d=0.1. One can check that condition (3.46) is fulfilled (128.3377 >114.5576) and then persistence of both groups of species is ensured, as observed in the figure. 1.4 6 1.6 1.8 52 2 4 2.2 1.5 2.4 31 2.6 20.5 10 0 0.5 1 1.5 2 1 2 3 4 5 6 0 0.5 1 1.5 2 1.4 1.6 1.8 2 2.2 2.4 2.6 Figure 2. Persistence in random chemostat (1.8)-(1.10) with sin =4, ks=1.7, µ0=8, ki=7, ¯µ0=7.5, b=0.1, ν=0.4, r1=0.4, r2=0.8, D=0.7, Dr=0.8, Dl=0.6, d=0.1, β=1, γ=1 and initial condition (s0,x10,x20)=(2,5,2). Figure 3 displays a zoom of the different panels of Figure 2 to let the reader observe the long-time random dynamics of system (1.8)–(1.10) more clearly, in fact, we can notice how the solutions of the system reach some place in which the attracting set Agiven by (3.47) is. 1.4 1.45 2.4 1.5 1.55 1 2.2 1.6 1.65 0.8 2 1.7 0.6 1.75 1.8 1.8 0.4 1.6 0.2 1.4 0 0.1 0.15 0.2 0.25 0.3 1.8 2 2.2 2.4 2.6 2.8 0 0.1 0.2 0.3 0.4 0.5 1.45 1.5 1.55 1.6 1.65 1.7 1.75 1.8 1.85 Figure 3. Persistence in random chemostat (1.8)-(1.10) with sin =4, ks=1.7, µ0=8, ki=7, ¯µ0=7.5, b=0.1, ν=0.4, r1=0.4, r2=0.8, D=0.7, Dr=0.8, Dl=0.6, d=0.1, β=1, γ=1 and initial condition (s0,x10,x20)=(2,5,2). AIMS Mathematics Volume 6, Issue 4, 4025–4052. 4045 In this case we can observe how the substrate concentration sdecreases at the beginning while the concentration of both species x1and x2increases until both species reach their maximum values. Then, both species start decreasing and, as a consequence, the concentration of substrate increases slightly up to reaching the attracting set . In Figure 4 we change a little bit the value of the input flow Dto observe its influence on the dynamics of the system. Then we set again sin =4, ks=1.7, µ0=8, ki=7, ¯µ0=7.5, b=0.1, ν=0.4, r1=0.4, r2=0.8. Concerning the removal rate, in this case we increase the input flow to D=0.7, and set Dr=1.9, Dl=1.5 and d=0.2. It is not difficult to check that condition (3.46) in Theorem 3.5 is verified (122.7993 >95.2622), and then persistence of species is ensured, as observed in the figure. 2 5.5 2.5 5 3 2 4.5 1.8 3.5 1.6 4 1.4 4 3.5 1.2 4.5 1 30.8 0.6 2.5 0.4 20.2 0 0.5 1 1.5 2 2 2.5 3 3.5 4 4.5 5 5.5 0 0.5 1 1.5 2 2 2.5 3 3.5 4 4.5 Figure 4. Persistence in random chemostat (1.8)-(1.10) with sin =4, ks=1.7, µ0=8, ki=7, ¯µ0=7.5, b=0.1, ν=0.4, r1=0.4, r2=0.8, D=1.7, Dr=1.9, Dl=1.5, d=0.2, β=1, γ=1 and initial condition (s0,x10,x20)=(2,5,2). In order to have a better panoramic of the dynamics of the system near the attracting set, Figure 7 shows a zoom of every panel of Figure 4. 2.4 4 2.6 2.8 3 3.2 0.5 3.5 3.4 3.6 0.45 3.8 4 30.4 4.2 4.4 0.35 2.5 0.3 0.25 20.2 0.2 0.22 0.24 0.26 0.28 0.3 0.32 0.34 2 2.1 2.2 2.3 2.4 2.5 0.2 0.22 0.24 0.26 0.28 0.3 0.32 0.34 2.5 3 3.5 4 4.5 Figure 5. Persistence in random chemostat (1.8)-(1.10) with sin =4, ks=1.7, µ0=8, ki=7, ¯µ0=7.5, b=0.1, ν=0.4, r1=0.4, r2=0.8, D=1.7, Dr=1.9, Dl=1.5, d=0.2, β=1, γ=1 and initial condition (s0,x10,x20)=(2,5,2). AIMS Mathematics Volume 6, Issue 4, 4025–4052. 4046 The main difference between the behavior of our system in this last case, respect to the first one, is due to the change in the removal rate D, which is one of the most important parameters when modeling chemostats, due to its influence on the dynamics of the device. To be more precise, recall that species sticked on to the wall of the culture vessel cannot be removed by the flow passing through the pumps, differently to the species in the medium that can be easily removed from their tank (specially when increasing the removal rate as in this case). This is the reason why we observe how fast the first species decreases whereas the second species behaves in a total different way. Now, in Figure 6 we present the case of extinction of both species studied in Theorem 3.5. To this end, we do not change the value of the input flow respect to the value in the last figure but we consider sin =4, ks=0.4, µ0=4, ki=0.5, ¯µ0=1.4, b=0.1, ν=1.4, r1=0.4, r2=0.8. It is easy to notice that we decreased considerably the the Monod constant ks, the inhibition kiand the consumption of both species whereas we increased the death rate ν. Then it is not a surprise that condition (3.38) in Theorem (3.5) is fulfilled (3.1>2.7548) and condition (3.46) in Theorem 3.5 is not true (1.6261 ≯41.9151), whence we have extinction of both species, as one can observe. 0 5 0.2 0.4 0.6 4 0.8 4 1 3.5 1.2 3 1.4 3 1.6 2.5 1.8 22 2 1.5 11 0.5 00 01234 0 1 2 3 4 5 01234 0 0.5 1 1.5 2 Figure 6. Extinction in random chemostat (1.8)-(1.10) with sin =4, ks=0.4, µ0=4, ki=0.5, ¯µ0=1.4, b=0.1, ν=1.4, r1=0.4, r2=0.8, D=1.7, Dr=1.9, Dl=1.5, d=0.2, β=1, γ=1 and initial condition (s0,x10,x20)=(2,5,2). Finally, we show that bi-stability also happens when considering wall growth, as we do in this work. To this end, we set sin =10, ks=1.8, µ0=3.4, ki=7.7, ¯µ0=2.3, b=0.1, ν=0.4, r1=0.4, r2=0.8. Concerning the removal rate, in this case we take D=1.7, Dr=1.3, Dl=2.1 and d=0.4. It is easy to check that µ(sin)≤D≤µ(sm) (1.3717 <1.7<1.7285), the same condition under which bi-stability is observed in the deterministic framework (see Proposition 2.1 (3)). The initial condition will be specified later. AIMS Mathematics Volume 6, Issue 4, 4025–4052. 4047 In Figure 7 we can notice how both the planktonic biomass and the attached biomass go to extinction when considering their initial concentrations x10 and x10 closed to the wash-out. 0 0.01 0.002 0.009 0.004 0.008 10 0.006 0.007 9.5 0.006 9 0.008 8.5 0.005 0.01 8 0.004 7.5 0.012 0.003 7 6.5 0.002 6 0.001 5.5 05 5678910 0 0.001 0.002 0.003 0.004 0.005 0.006 0.007 0.008 0.009 0.01 5678910 0 0.002 0.004 0.006 0.008 0.01 0.012 Figure 7. Extinction in random chemostat (1.8)-(1.10) with sin =10, ks=1.8, µ0=3.4, ki=7.7, ¯µ0=2.3, b=0.1, ν=0.4, r1=0.4, r2=0.8, D=1.7, Dr=1.3, Dl=2.1, d=0.4, β=1, γ=1 and initial condition (s0,x10,x20)=(0.01,5,0.01). On the other hand, in Figure 8 we consider the initial condition far enough from the wash-out. Thus, we observe persistence of both the planktonic biomass and the attached biomass. 1 4 2 3 3.5 4 8 3 5 7 6 6 2.5 7 5 8 24 3 1.5 2 11 12345678 1 1.5 2 2.5 3 3.5 4 12345678 1 2 3 4 5 6 7 8 Figure 8. Persistence in random chemostat (1.8)-(1.10) with sin =10, ks=1.8, µ0=3.4, ki=7.7, ¯µ0=2.3, b=0.1, ν=0.4, r1=0.4, r2=0.8, D=1.7, Dr=1.3, Dl=2.1, d=0.4, β=1, γ=1 and initial condition (s0,x10,x20)=(4,5,1). AIMS Mathematics Volume 6, Issue 4, 4025–4052. 4048 Eventually, Figure 9 presents a zoom of the previous case in order to help the reader to see how the dynamics inside the attracting set of our model is. 4 4 4.5 3.9 5 3.8 4 5.5 3.7 6 3.5 3.6 6.5 3.5 3 7 3.4 2.5 3.3 2 3.2 3.1 1.5 3 1.2 1.4 1.6 1.8 2 2.2 2.4 3 3.1 3.2 3.3 3.4 3.5 3.6 3.7 1.2 1.4 1.6 1.8 2 2.2 2.4 5.4 5.6 5.8 6 6.2 6.4 6.6 6.8 7 7.2 Figure 9. Persistence in random chemostat (1.8)-(1.10) with sin =10, ks=1.8, µ0=3.4, ki=7.7, ¯µ0=2.3, b=0.1, ν=0.4, r1=0.4, r2=0.8, D=1.7, Dr=1.3, Dl=2.1, d=0.4, β=1, γ=1 and initial condition (s0,x10,x20)=(4,5,1). 5. Conclusions In this final section of our paper we would like to summarize the results obtained in Section 3 and observed in the numerical simulations displayed in Section 4. We have considered the chemostat model (1.8)–(1.10) with wall growth and Haldane consumption kinetics, under bounded perturbations on the removal rate, motivated by observations in real devices in industrial setup and biotechnology. To this end, we make use of a function and the Ornstein-Uhlenbeck process to ensure the random fluctuations be bounded in some interval (that is typically observed in practice). Thanks to standard techniques from the classical theory of ordinary differential equations and, of course, due to the properties satisfied by the bounded noise Φ(z∗(θtω)), we proved in Theorem 3.1 that the random chemostat (1.8)–(1.10) possesses a unique global solution which, in addition, remains positive as long as we consider positive initial conditions. Then, we proved in Theorem 3.2 the existence of an absorbing set for the solutions of system (1.8)– (1.10). This allowed us to ensure the existence of an attracting set (forwards in time) for the solutions of system (1.8)–(1.10). In addition, we proved that both absorbing and attracting set enjoys a very nice property: They are deterministic, i.e., they do not depend on the realization of the noise. After that, we went further and we investigated the internal structure of the deterministic attracting set obtained in Theorem 3.2 in order to have more detailed information about the long-time dynamics of our system. To this end, we proved two theorems that can be interesting for both mathematicians AIMS Mathematics Volume 6, Issue 4, 4025–4052. 4049 and biologists: On the one hand, Theorem 3.5 provides conditions under which extinction of both species cannot be avoided as well as some conditions leading in upper bounds for both species; on the other hand, Theorem 3.6 and 3.7 states a condition under which not only the persistence of the total concentration of species but also the persistence of both the planktonic biomass and the attached biomass is ensured. This last point is undoubtedly the most meaningful since it is the main goal pursued by practitioners. We would also like to highlight again the importance of considering the non-monotonic growth rate (1.6), since it allows us to study mathematical models that really fit better the devices displayed in real life, as underlined in the introduction. Concerning the non-monotonic growth of species, it is worth to mention some aspects observed in this work compared to the classical deterministic model (1.1)–(1.2) with Haldane consumption function. More precisely, the dynamics of the deterministic chemostat (1.1)–(1.2) with non-monotonic consumption function may exhibit a bi-stability for certain values of the removal rate D(see [21]), as explained in Section 2.1. In such a case, the concentration of the biomass converges asymptotically either to the wash-out equilibrium (which is not desirable since it means the extinction of the species) or to a positive steady state. This type of instability is observed in practice and entails an issue in industrial applications since it requires a good monitoring of the system to prevent it from going to the wash-out equilibrium (see [44–47]). Then, practitioners normally decide on avoiding this behavior by sizing the process such that the system has an unique globally stable equilibrium (see [48, 49]) which mathematically means to have the following condition D<min{µ(sin), µ(sm)}.(5.1) In addition, we supported every theoretical result with several numerical simulations depicting the random dynamics of the system when considering different values of the parameters involved in the model. Moreover, this also allowed us to illustrate the bi-stability for some values of the parameters involved in the model and the huge importance that the removal rate Dhas on the dynamics of the chemostat giving another justification of our interest in introducing the random perturbations on it. Finally, we would like to highlight that the Ornstein-Uhlenbeck process has proved again to be a very powerful tool when modeling real bounded fluctuations, not only when studying chemostats but also when dealing with any model involving bounded perturbations, as the ones in real life. Acknowledgments This work has been partially supported by Proyecto de Excelencia P12-FQM-1492 from Junta de Andaluc´ ıa, Project PGC2018-096540-B-I00 from Ministerio de Ciencia, Innovaci´ on y Universidades (Spanish government), Project I+D+i FEDER US-1254251 from Consejer´ ıa de Econom´ ıa y Conocimiento (Junta de Andaluc´ ıa) and Project P18-FR-4509 from Consejer´ ıa de Econom´ ıa y Conocimiento (Junta de Andaluc´ ıa). We thank the referees for their helpful comments which allowed us to add some details that made this work more elegant and interesting, not only for mathematicians, but also for people from other areas who are interested in this research. AIMS Mathematics Volume 6, Issue 4, 4025–4052.