Modeling and analysis of random and stochastic input flows in the chemostat model

Manuscript submitted to doi:10.3934/xx.xx.xx.xx AIMS’ Journals Volume X, Number 0X, XX 200X pp. X–XX MODELING AND ANALYSIS OF RANDOM AND STOCHASTIC1 INPUT FLOWS IN THE CHEMOSTAT MODEL2 Tom´ as Caraballo, Maria-Jos´ e Garrido-Atienza, Javier L´ opez-de-la-Cruz Dpto. de Ecuaciones Diferenciales y An´alisis Num´erico Facultad de Matem´aticas, Universidad de Sevilla C/ Tarfia s/n, Sevilla 41012, Spain and Alain Rapaport∗ MISTEA, Univ. Montpellier, Inra, Montpellier SupAgro 2, place Pierre Viala 34060 Montpellier, France (Communicated by the associate editor name) Abstract. In this paper we study a new way to model noisy input flows in the chemostat model, based on the Ornstein-Uhlenbeck process. We introduce a parameter βas drift in the Langevin equation, that allows to bridge a gap between a pure Wiener process, which is a common way to model random disturbances, and no noise at all. The value of the parameter βis related to the amplitude of the deviations observed on the realizations. We show that this modeling approach is well suited to represent noise on an input variable that has to take non-negative values for almost any time. 1 Introduction. Chemostat refers to a laboratory device used for growing mi-3 croorganisms in a cultured environment and has been regarded as an idealization4 of nature to study microbial ecosystems at steady state, which is a really impor-5 tant and interesting problem since they can be used to study genetically altered6 microorganisms, waste water treatment (see e.g. [16,25]) and play an important7 role in theoretical ecology (see e.g. [3,15,23,24,26]). The simplest chemostat8 device consists of three interconnected tanks called feed bottle,culture vessel and9 collection vessel. The nutrient is pumped from the first tank to the culture vessel,10 where the interactions between the species and the nutrients take place, and there11 is also another flow being pumped from the culture vessel to the third one such12 that the volume of the culture vessel remains constant. Derivation and analysis of13 chemostat models are well documented in [19,21] and references therein.14 Some standard assumptions for simple chemostat models are usually imposed,15 for instance, it is common to suppose that the availability of the nutrient and its16 supply rate are fixed. However, they are very strong restrictions since the real world17 2010 Mathematics Subject Classification. Primary: 37N25, 92D25, 34F05, 37H10; Secondary: 37B55, 60H10. Key words and phrases. chemostat model, Ornstein-Uhlenbeck process, random dynamical system, random attractor. Partially supported by FEDER and Ministerio de Econom´ıa y Competitividad under grant MTM2015-63723-P, Junta de Andaluc´ıa under the Proyecto de Excelencia P12-FQM-1492 and VI Plan Propio de Investigaci´on y Transferencia de la Universidad de Sevilla. ∗Corresponding author: Alain Rapaport. 1 2 T. CARABALLO, M.-J. GARRIDO-ATIENZA, J. L ´ OPEZ-DE-LA-CRUZ, A. RAPAPORT is non-autonomous and stochastic. This is one of the reasons which encourage us to1 go further away from deterministic systems and study stochastic chemostat models.2 Let us consider the classical chemostat model3 ds dt = (sin −s)Q V−µ(s)x, (1.1) dx dt =−Q Vx+µ(s)x, (1.2) where s(t) and x(t) denote concentrations of the nutrient and the microbial biomass,4 respectively; sin denotes the input concentration, Qis the input flow, Vis the5 volume of liquid media inside the culture vessel and the ratio Q/V , which is also6 denoted by D, is called dilution rate. We notice that all parameters are supposed to7 be positive and a function Holling type-II, µ(s) = ms/(a+s), is used as functional8 response of the microorganisms describing how the nutrient is consumed by the9 species, where ais the half-saturation constant and mdenotes the maximal specific10 growth rate of microorganisms. Sometimes we will refer to this function as the11 consumption function or growth function as well.12 There exist many different ways of modeling stochasticity and/or randomness in13 some deterministic model, see e.g. [4,5,6,18,22,28,29,30,32,33]. Considerations14 of stochastic processes in the chemostat model have already been tackled in the15 literature, but mainly on the growth function (see, for instance, [7]). This appears16 particularly relevant when the number of individual bacteria could be small, with17 a risk of extinction of the biomass population in finite time. Nevertheless, sudden18 extinctions in continuous cultures that are well supervised about a nominal regime19 are quite rare in practice. On another hand, fluctuations on the input flow that20 brings permanently resources to the bacterial population in continuous cultures are21 much likely to be observed. In the present work, we focus on the way to model22 these random fluctuations, taking into consideration that the effective flow rate has23 to stay non-negative. From the biological point of view, the fact of introducing a24 noisy term in the input flow of our chemostat model is a really interesting problem25 found in the laboratory since, for instance, it reflects the presence of particles of dirt26 inside the pumps or temporary clogs at the input or output of the chemostat. Then,27 it is well known that continuous flow is often subjected to random fluctuations with28 time.29 In this paper we shall follow two different approaches to perturb the input flow30 in the chemostat model. As the volume Vis constant, it is equivalent to have31 disturbances on the dilution rate D=Q/V instead of considering them on the input32 flow. On the one hand, we will consider a perturbation by making use of the well-33 known standard Wiener process such that we replace Dby the stochastic term D+34 α˙ω(t) in the deterministic system (1.1)-(1.2), where ωdenotes the standard Wiener35 process and α > 0 represents the intensity of the noise. Thus, the corresponding36 stochastic system is given by37 ds =(sin −s)D−msx a+sdt +α(sin −s)dω(t), dx =−Dx +msx a+sdt −αxdω(t). As will be discussed in Section 3, this common approach will lead into some38 drawbacks, for instance, the persistence of the microorganisms can never be ensured39 RANDOM AND STOCHATIC CHEMOSTAT MODELS 3 and some state variables can take negative values. Therefore, this stochastic process1 does not seem to be the best way to represent disturbances on the input flow in2 the chemostat model, due to the fact that it provides unrealistic situations from3 the biological point of view since it would mean that the pumps are revering flow,4 which is not reasonable at all.5 On the other hand, we will consider a suitable Ornstein-Uhlenbeck (O-U) process6 to perturb the dilution rate. Particularly, we are interested in replace Dby the7 random term D+αz∗ β,ν (θtω), where z∗ β,ν(θtω) denotes some suitable O-U process8 which will be carefully introduced later and α > 0 represents again the intensity of9 noise. In such a way, the resulting random model is given by the following system10 of random differential equations11 ds dt = (sin −s)D+αz∗ β,ν (θtω)−µ(s)x, (1.3) dx dt =−D+αz∗ β,ν (θtω)x+µ(s)x. (1.4) Concerning the O-U process, some essential properties will be provided which12 will allow us to set up a new framework and, moreover, to make calculations in13 the next section. To sum up some of the main ingredients to be used, for every14 fixed event ω, it will be possible to choose βω∈Rsuch that the corresponding15 realizations of the perturbed input flow, D+αz∗ βω,ν (θtω), will remain for every t∈R16 inside some strictly positive band which should be previously fixed, for instance, by17 practitioners. In such a way, for every fixed event ω, the resulting random chemostat18 model will be given by19 ds dt = (sin −s)D+αz∗ βω,ν (θtω)−µ(s)x, (1.5) dx dt =−D+αz∗ βω,ν (θtω)x+µ(s)x. (1.6) As a consequence, since βω∈Rdepends on the event ωpreviously fixed, the20 solutions of system (1.5)-(1.6) may not generate a random dynamical system. Nev-21 ertheless, this does not represent any inconvenient for the analysis of the long time22 behavior of the random differential system (1.5)-(1.6), since it can be investigated23 for every fixed event ω. In fact, we will be able to obtain some results on forwards24 convergence (in time) of solutions, instead of the pullback convergence obtained25 within the framework of random dynamical systems.26 This new approach, which arises from the nature of the particular noise (the27 suitable O-U process), leads into another unusual technique which seems to be really28 interesting since, for instance, it allows us to guarantee the existence of compact29 and attracting sets (forwards in time) which are strictly inside the positive cone,30 whence we will ensure the persistence of the species in the sense that there exists31 a number η > 0 such that for any non null initial biomass x(0) each realization32 satisfies 1 33 lim inf t→+∞x(t)≥η > 0.(1.7) Needless to say that this is the principal goal pursued by biologists for modeling34 bounded disturbances in a biological framework, differently to other several previ-35 ous works as [22] where the authors consider disturbances in the chemostat model36 1We will refer to this sense throughout the whole paper when using the term “persistence”. In case of referring to another sense, we will specify the new definition. 4 T. CARABALLO, M.-J. GARRIDO-ATIENZA, J. L ´ OPEZ-DE-LA-CRUZ, A. RAPAPORT by means of the standard Wiener process and prove some results concerning the1 persistence of the biomass in the sense lim inft→+∞x(t)>0,which is weaker than2 (1.7).3 We also achieve some improvements comparing our results throughout this paper4 with the ones by Xu et al in [31] since, even though they consider stochastic noise on5 the dilution rate in the chemostat model, they need a condition on the parameter of6 the amplitude of the noise to ensure the persistence (see, for instance, Theorem 1.27 and Section 4 in [31] where the authors ensure the necessity of a smallness condition8 on the amount of noise α > 0) whereas, in our case, modeling the disturbances with9 the Ornstein-Uhlenbeck process, there is no discussion needed on the amplitude10 of the noise to ensure the persistence (which is, in addition, in the stronger sense11 (1.7)). Moreover, the authors in [31] prove the results in probability while we will12 prove all the results almost surely, i.e., for every realization in a set of events of full13 measure.14 The previous reasons constitute a few representative examples which support that15 this way of perturbing the dilution rate by using the Ornstein-Uhlenbeck process16 fits much better the real situations we wish to model. Apart from the advantages17 described above, we will also obtain some improvements with respect to the results18 obtained when analyzing the deterministic chemostat model (1.1)-(1.2), as we will19 explain in more detail in Section 2. To be more precise, in the deterministic setting20 the washout equilibrium (sin,0) is attractive if D=µ(sin) (see [19]) whereas, in our21 case by using the O-U process, it is possible to prove that there exists an attracting22 (forwards in time) set for the solutions of our system which has several points (in23 fact, all of them except the washout) inside the positive cone.24 The paper is organized as follows. In Sections 2and 3we will analyze the25 chemostat models perturbed by the Ornstein-Uhlenbeck process and the white noise,26 respectively. We will prove the existence and uniqueness of a global solution and27 we will also state some results regarding the existence of a compact absorbing set28 as well as an attracting one. Finally, in Section 4some numerical simulations which29 will support the results previously proved will be also presented. In addition, we30 present briefly some basic concepts and results concerning the theory of random31 dynamical systems (RDSs) in Appendix.32 2 The chemostat model with random input flow. In this section, we are33 interested in investigating system (1.3)-(1.4), where a random perturbation on the34 dilution rate has been introduced by means of a suitable O-U process, as explained35 in the introductory section. To this end, we will present a new framework which will36 be a bit different from the ones used in previous papers by Caraballo and several37 coauthors (see e.g. [7,8,9,11]), where the theory of random dynamical systems is38 used.39 2.1 The Ornstein-Uhlenbeck process. Let Wbe a two sided Wiener process. Kolmogorov’s theorem ensures that Whas a continuous version, that we will denote by ω, whose canonical interpretation is as follows: let Ω be defined by Ω = {ω∈ C(R,R) : ω(0) = 0}, Fthe Borel σ−algebra on Ω generated by the compact open topology (see [2] for details) and Pthe corresponding Wiener measure on F. We consider the Wiener shift flow given by θtω(·) = ω(·+t)−ω(t), t ∈R. RANDOM AND STOCHATIC CHEMOSTAT MODELS 5 Then, (Ω,F,P,{θt}t∈R) is a metric dynamical system (see Appendix for details).1 Now, let us introduce the following Ornstein-Uhlenbeck (O-U) process defined2 on (Ω,F,P,{θt}t∈R) as the random variable given by3 z∗ β,ν(θtω) = −βν 0 Z −∞ eβsθtω(s)ds, t ∈R, ω ∈Ω, β, ν > 0,(2.1) which solves the following Langevin equation (see [2,12,13])4 dz +βzdt =νdω(t), t ∈R.(2.2) The O-U process given by (2.1) is a stationary mean-reverting Gaussian stochastic5 process where β > 0 is a mean reversion constant that represents the strength with6 which the process is attracted by the mean or, in other words, how strongly our7 system reacts under some perturbation, and ν > 0 is a volatility constant which8 represents the variation or the size of the noise independently of the amount of the9 noise α > 0. In fact, the O-U process can describe the position of some particle10 by taking into account the friction, which is the main difference with the standard11 Wiener process and makes our model to be a better approach to the real ones,12 specially when modeling processes in microbiology as in our case. In addition, the13 O-U process can be understood as a kind of generalization of the standard Wiener14 process, which would correspond to take β= 0 and ν= 1 in (2.1).15 By taking into account the definition of both parameters βand νinvolved in the16 Langevin equation (2.2), we highlight the following relevant observations concerning17 the effect caused by each of them on the evolution of the process.18 2.1.1 Fixed β > 0.Then, the volatility of the process is larger if we consider a19 larger ν. However, the evolution of the process is smoother when we take a smaller20 value of ν. This is reasonable since νdecides the amount of noise introduced to dz,21 which measures the variation of the process, hence the process will be subjected to22 suffer many more changes when choosing a larger value of ν. We can observe this23 behavior in Figure 1where we simulate two realizations of the perturbed dilution24 rate D+αz∗ β,ν (θtω) with D= 2, α= 0.8, β= 2 and we consider ν= 0.1 (blue) and25 ν= 0.5 (orange).26 Figure 1. Realizations of the perturbed dilution rate with D= 2, α= 0.8 and β= 2 6 T. CARABALLO, M.-J. GARRIDO-ATIENZA, J. L ´ OPEZ-DE-LA-CRUZ, A. RAPAPORT 2.1.2 Fixed ν > 0.Then, the process tends to go further away from the mean1 value if we consider a smaller value of β. However, the attraction of the mean2 increases when taking a larger β, which is quite logical since βhas a huge influence3 on the drift of the Langevin equation (2.2). For instance, we can observe this4 behavior in Figure 2, where we simulate two realizations of the perturbed dilution5 rate D+αz∗ β,ν (θtω) with D= 2, α= 0.8, ν= 0.5 and we take β= 2 (blue) and6 β= 10 (orange).7 Figure 2. Realizations of the perturbed dilution rate with D= 2, α= 0.8 and ν= 0.5 Now, we establish the following result involving important ergodic properties8 held by the O-U process which will be used at several places in the sequel.9 Proposition 2.1. There exists a θt-invariant set e Ω∈ F of Ωof full P−measure10 such that for ω∈e Ωand β, ν > 0, we have11 (i) the random variable |z∗ β,ν (ω)|is tempered (see Definition 4.3).12 (ii) the mapping13 (t, ω)→z∗ β,ν(θtω) = −βν 0 Z −∞ eβsω(t+s)ds+ω(t) is a stationary solution of (2.2)with continuous trajectories;14 15 (iii) for any ω∈˜ Ωone has:16 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∗ β,ν ]<∞; (iv) finally, for any ω∈e Ω, lim β→∞ z∗ β,ν (θtω) = 0,for all t∈R. Remark 2.1. We note that the proof of (iv) can be found in [1](see Lemma 4.1)17 and we refer the readers to [2,12]for the proof of (i)-(iii).18 RANDOM AND STOCHATIC CHEMOSTAT MODELS 7 Then, if we restrict the metric dynamical system to e Ω, we obtain a new metric1 dynamical system, see [10]. For simplicity, we will denote this new metric dynamical2 system by the old symbols, namely (Ω,F,P,{θt}t∈R).3 Our aim in this section is to analyze the long-time behavior of system (1.3)-4 (1.4). To this end, let us first fix a strictly positive interval, namely (b,¯ b)⊂R,5 where ¯ b > b > 0. Thanks to the last item in Proposition 2.1, for each ω∈Ω, it is6 possible to choose β∈Rlarge enough such that the corresponding realization of7 the perturbed input flow, D+αz∗ β,ν(θtω), remains inside the interval (b,¯ b) for every8 t∈R. Nevertheless, it is not possible to ensure, from a theoretical point of view,9 that there exists some β∈Rsuch that almost all realizations of the perturbed input10 flow remains in (b,¯ b). Because of this reason, we will analyze our system (1.3)-(1.4)11 for every fixed ω∈Ω. As stated above, we know that it is possible to find βω∈R12 such that D+αz∗ βω(θtω)∈(b,¯ b) for every t∈R, then we need to analyze the13 following random system14 ds dt = (sin −s)D+αz∗ βω,ν (θtω)−µ(s)x, dx dt =−D+αz∗ βω,ν (θtω)x+µ(s)x. We would like to remark that the choice of βdepends on ω∈Ω, this is the reason to15 write βωin the previous system. Then, βωacts in practice as a control parameter.16 However, once an event ω∈Ω is fixed, we have that βωis also a fixed real number,17 thus we will rewrite βω=βand we will focus on analyzing the system (1.3)-(1.4) for18 every fixed ω∈Ω. The interesting fact is that the attracting sets for the solutions19 will not depend on ωeven though βdepends on ω, so we will obtain a non random20 set where all solutions for all realizations will approach to.21 In the rest of the section, we will prove the existence and uniqueness of a global22 solution of system (1.3)-(1.4) as well as the existence of a strictly positive forward23 attracting set for its solutions, whence we will ensure that the microorganism con-24 centration will converge asymptotically to a strictly positive interval or, in other25 words, we will be able to guarantee the persistence of the species.26 Before starting with the analysis previously motivated, let us recall the classical Monod expression µ(s) = ms a+s,for all s≥0, denoting the consumption for the specific growth rate function of the species. Then,27 we define the following constants28 s:= µ−1(b) and ¯s:= µ−1(¯ b) (2.3) which will be essential henceforth. In Figure 3where we plot the mapping s7→ µ(s)29 and overlap a realization of the perturbed input flow as well without taking into30 account the dependency of time.31 8 T. CARABALLO, M.-J. GARRIDO-ATIENZA, J. L ´ OPEZ-DE-LA-CRUZ, A. RAPAPORT D ¯ b b s¯s m s Figure 3. Realizations of the perturbed dilution rate, sand ¯s 2.2 The random chemostat model. We are interested in analyzing the fol-1 lowing random chemostat model2 ds dt =−D+αz∗ β,ν (θtω)s−µ(s)x+sin D+αz∗ β,ν (θtω),(2.4) dx dt =−D+αz∗ β,ν (θtω)x+µ(s)x, (2.5) where µ(s) = ms/(a+s) denotes the Holling type-II consumption function as the3 functional response of the microorganisms. Henceforth, ω∈Ω is fixed and β∈R4 is also a parameter which has been fixed such that D+αz∗ β,ν (θtω)∈(b,¯ b) for all5 t∈R.6 In this section, X={(x, y)∈R2:x≥0, y ≥0}denotes the positive cone.7 Theorem 2.1. For any initial pair v0:= (s0, x0)∈ X, system (2.4)-(2.5)possesses a unique global solution v(·; 0, ω, v0) := (s(·; 0, ω, v0), x(·; 0, ω, v0)) ∈ C1([0,+∞),X) with v(0; 0, ω, v0) = v0, where s0:= s(0; 0, ω, v0)and x0:= x(0; 0, ω, v0).8 Proof. We set v(·; 0, ω, v0) := (s(·; 0, ω, v0), x(·; 0, ω, v0)) such that system (2.4)-(2.5) can be rewritten as dv dt =L(θtω)·v+F(v, θtω), where9 L(θtω) = −(D+αz∗ β,ν(θtω)) −m 0−(D+αz∗ β,ν(θtω)) + m and F:X × [0,+∞)−→ R2is given by10 F(ξ, θtω) =    ma a+ξ1 ξ2+sinD+αsinz∗ β,ν(θtω) −ma a+ξ1 ξ2  , where ξ= (ξ1, ξ2)∈ X.11 Since z∗(θtω) is continuous, Lgenerates an evolution system on R2. Moreover,12 we notice that F(·, θtω)∈ C1(X × [0,+∞); R2) which implies that it is locally13 RANDOM AND STOCHATIC CHEMOSTAT MODELS 9 Lipschitz with respect to (ξ1, ξ2)∈ X. Therefore, system (2.4)-(2.5) possesses a1 unique local solution. Now, we prove that the unique local solution of system (2.4)-2 (2.5) is defined for any forward time and is, in fact, a unique global one. To this3 end, we define the new state variable q(t) := s(t) + x(t)−sin. Then, qsatisfies the4 following differential equation5 dq dt =−D+αz∗ β,ν (θtω)q, (2.6) whose solution is given by6 q(t; 0, ω, q(0)) = q(0)e−Dt−αRt 0z∗ β,ν (θsω)ds.(2.7) It is straightforward to check that qdoes not blow up at any finite time, thanks to the positiveness of the dilution rate and the ergodic properties of the O-U process (see Proposition 2.1, (iii)), moreover qis bounded. In addition, after solving (2.5) we have the following upper bound for the x-equation x(t; 0, ω, x(0)) ≤x(0)e−(D−m)t−αRt 0z∗ β,ν (θsω)ds, since µ(s)≤mfor any s≥0.7 In conclusion, xis also bounded by an expression which does not blow up at8 any finite time. Therefore, sdoes not blow up either and we can conclude that our9 chemostat model (2.4)-(2.5) possesses a unique global solution.10 Moreover, since x≡0 solves (2.5) and every realization of our noise remains in some strictly positive interval, we conclude that ds dt s=0 =D+αz∗ β,ν (θtω)sin >0. Thus, we can ensure the unique solution of system (2.4)-(2.5) to be in the positive11 cone Xfor every initial value v0∈ X . 12 Now, we are interested in proving the existence of an attracting set. From now13 on, F⊂ X denotes a bounded set in the positive cone.14 Theorem 2.2. For any ε > 0, there exists a deterministic compact absorbing set15 Bε⊂ X for the solution of our system (2.4)-(2.5), i.e., there exists TF(ω, ε)>0such16 that for every given initial pair v0∈F, the solution corresponding to v0remains17 inside Bεfor all t≥TF(ω, ε).18 Proof. Consider again q(t) = s(t) + x(t)−sin. Then, thanks to (2.7), we obtain19 lim t→+∞q(t; 0, ω, q(0)) = 0.(2.8) Thus, given v0∈Fand any ε > 0, there exists TF(ω, ε)>0 such that −ε≤q(t; 0, ω, q(0)) ≤ε for every t≥TF(ω, ε).20 Then,21 Bε:= {(s, x)∈ X :sin −ε≤s+x≤sin +ε}.(2.9) is a compact absorbing set in X. 22 Therefore, thanks to Theorem 2.2, we have that B0:= {(s, x)∈ X :s+x=sin} 16 T. CARABALLO, M.-J. GARRIDO-ATIENZA, J. L ´ OPEZ-DE-LA-CRUZ, A. RAPAPORT then we can find some ε(ω)>0 small enough such that σ(t) is strictly decreasing1 and2 −(¯ D+αz∗(θtω)) −m(sin +σ(t)e−αz∗(θtω)) a+sin +σ(t)e−αz∗(θtω)κ(t)>0 (3.7) holds for all t∈[¯ t−ε(ω),¯ t). Hence, from (3.7) we have dσ dt (¯ t−ε(ω)) >0. Consequently, there exists some δ(ω)>0, small enough, such that σ(t) is strictly3 increasing for all t∈[¯ t−ε(ω),¯ t−ε(ω)+δ(ω)), which clearly contradicts the unique-4 ness of solution. Hence, (3.6) holds true for all t∈Rand we can also ensure that5 σis bounded.6 7 Since σ+κand σare bounded in both cases, κis also bounded. Hence, the8 unique local solution of system (3.3)-(3.4) is a unique global one. Moreover, the9 unique global solution of system (3.3)-(3.4) remains in e Xfor every initial value in10 e Xsince κ≡0 solves the same system.11 Finally, the mapping ϕu:R+×Ω×e X → e Xgiven by12 ϕu(t, ω)u0:= u(t; 0, ω, u0),for all t≥0, u0∈e X, ω ∈Ω, defines an RDS generated by the solution of (3.3)-(3.4). The proof of this statement13 follows similarly to the one of Theorem 2.1, hence we omit it.14 Theorem 3.2. There exists a tempered compact random absorbing set B0(ω)∈15 E(e X)for the RDS {ϕu(t, ω)}t≥0, ω∈Ω.16 Proof. Thanks to (3.5), we have17 eq(t; 0, θ−tω, eq(0)) = eq(0)e−¯ Dt−αR0 −tz∗(θsω)ds t→+∞ −→ 0. Then, for all ε > 0 and u0∈E(θ−tω), there exists TE(ω)>0 such that, for all t≥TE(ω), −ε≤eq(t; 0, θ−tω, eq(0)) ≤ε. If we assume that σ(t)≥0 for all t≥0, which corresponds to Case 1 in the proof of Theorem 3.1, since κ(t)≥0 for all t≥0, we have that B1 ε(ω) := {(σ, κ)∈ X :σ≥0, σ +κ≤ε} is a tempered compact random absorbing set in e X.18 In the other case, i.e., if there exists some t∗>0 such that σ(t∗) = 0, which corresponds to Case 2 in the proof of Theorem 3.1, we proved that σ(t; 0, θ−tω, σ(0)) >−(a+sin)eαz∗(ω). Hence, we obtain that B2 ε(ω) := n(σ, κ)∈e X:−ε−(a+sin)eαz∗(ω)≤σ≤0,−ε≤σ+κ≤εo is a tempered compact random absorbing set in e X.19 In conclusion, defining20 Bε(ω) = B1 ε(ω)∪B2 ε(ω) =n(σ, κ)∈ X :−ε≤σ+κ≤ε, σ ≥ −(a+sin)eαz∗(ω)−εo, RANDOM AND STOCHATIC CHEMOSTAT MODELS 17 (see Figure 5) we obtain that Bε(ω) is a tempered compact random absorbing set1 in e Xfor every ε > 0.2 σ κ −(a+sin)eαz∗(ω)−ε ε −ε B2 ε(ω) B1 ε(ω) Figure 5. Absorbing set Bε(ω) σ κ −(a+sin)eαz∗(ω) B0(ω) Figure 6. Absorbing set B0(ω) Then, we have that B0(ω) := n(σ, κ)∈e X:σ+κ= 0, σ ≥ −(a+sin)eαz∗(ω)o is a tempered compact random absorbing set as we wanted to prove (see Figure3 6).4 Therefore, thanks to Proposition 4.1, it follows directly that system (3.3)-(3.4)5 possesses a unique pullback random attractor such that A(ω)⊂B0(ω).6 Proposition 3.1. The pullback random attractor of system (3.3)-(3.4)consists of7 a singleton component given by A(ω) = {(0,0)}as long as8 ¯ D > µ(sin) (3.8) holds true.9 Proof. We would like to note that the result in this proposition follows trivially if10 σremains always positive (Case 1 in the proof of Theorem 3.1) since in that case11 both σand κare positive and σ+κtends to zero when tgoes to infinity, thus the12 pullback random attractor is directly given by A(ω) = {(0,0)}.13 Due to the previous reason, we will only present the proof when there exists some14 t∗>0 such that σ(t∗) = 0, which implies that σ(t)<0 for all t > t∗. Hence, since15 µ(s) = ms/(a+s), from (3.4) we have16 dκ dt ≤ −(¯ D+αz∗)κ+msin a+sin κ, which allows us to state the following inequality17 κ(t;t∗, θ−tω, κ(t∗)) ≤κ(t∗)e−¯ D−msin a+sin (t−t∗)−αRt∗ −tz∗(θsω)ds,(3.9) where the right side of (3.9) tends to zero when tgoes to infinity, as long as (3.8)18 is fulfilled. Therefore A(ω) = {(0,0)}.19 Finally, defining the mapping T: Ω ×e X −→ e Xby T(ω, ζ) = (ζ1−sin)eαz∗(ω), ζ2eαz∗(ω), 18 T. CARABALLO, M.-J. GARRIDO-ATIENZA, J. L ´ OPEZ-DE-LA-CRUZ, A. RAPAPORT it is a homeomorphism whose inverse is given by the following expression T−1(ω, ζ) = sin +ζ1e−αz∗(ω), ζ2e−αz∗(ω). If we now define1 ϕv(t, ω)v0:= T−1(θtω, ϕu(t, ω)T(ω, v0)) =T−1(θtω, u(t; 0, ω, u0)) =v(t; 0, ω, v0) by Proposition 4.2 it turns out that ϕvis the RDS associated to the stochastic2 chemostat model (3.1)-(3.2) with pullback random attractor b A(ω)⊂b B0(ω), where3 b B0(ω) := n(s, x)∈e X:s+x=sin, s ≥ −ao.(3.10) In addition, under (3.8), the pullback random attractor for (3.1)-(3.2) reduces4 to a singleton set b A(ω) = {(sin,0)}, which means that the microorganisms become5 extinct.6 We remark that it is not possible to provide conditions which ensure the persis-7 tence of the microbial biomass even though our numerical simulations show that we8 obtain it for many different values of the parameters involved in the system, as we9 will present in Section 4. Moreover, we will also compare the numerical simulations10 concerning this section with those corresponding to the previous one.11 4 Numerical simulations and final comments. In this section we would like12 to highlight some comments about both ways of modeling stochasticity and ran-13 domness in a chemostat model and show some numerical simulations which will14 support the results proved throughout Sections 2and 3.15 We remark that we use the Euler-Maruyama method which is a simple general-16 ization of the Euler method for ODEs to stochastic differential equations. The main17 difference is the discretization of the term dW(t) = W(τj)−W(τj−1) for some par-18 tition {τj}of the time interval, where we make use of the fact that such a difference19 is a Gaussian variable with mean zero and variance τj−τj−1. For more detailed20 information about the definition of the numerical scheme necessary to obtain the21 numerical simulations we refer the readers to [7,8,9,20].22 In every simulation, the dashed lines represent the solution of the deterministic23 systems, i.e., the behavior of the stochastic/random system after taking α= 0,24 whereas the continuous lines correspond to different realizations of the solution of25 the corresponding stochastic/random system.26 Now, we will show some simulations concerning the random chemostat model27 studied in Section 2. In each of the following figures three panels are displayed: the28 left one shows the phase plane and the general dynamics of the chemostat model;29 the two panels on the right side help us to see two important zones in the phase30 plane.31 In Figure 7we set D= 2, sin = 4, a= 0.6, m= 5, α= 0.5, β= 1, ν= 0.732 and initial values s(0) = 2, x(0) = 5 for the nutrient and the microorganism,33 respectively. In this case (2.11) holds true and this is the reason why we can observe34 the persistence of the species. We can also see how the realizations are approaching35 the line s+x=sin, as proved in (2.8).36 RANDOM AND STOCHATIC CHEMOSTAT MODELS 19 Figure 7. Persistence of the species in the random chemostat model In contrast to the last case, in Figure 8we take D= 3.5, sin = 2, a= 0.8, m=1 0.5, α= 0.5, β= 1, ν= 0.7 and initial values s(0) = 2.5, x(0) = 5 for the nutrient2 and the microorganisms, respectively. Then we can see that the microorganisms3 extinguish what is not surprising due to the fact that condition (2.10) is fulfilled.4 We can also see here how the realizations are approaching the line s+x=sin, as5 proved in (2.8).6 Figure 8. Extinction of the species in the random chemostat model The following simulations represent the stochastic chemostat analyzed in Section7 3, where the dilution rate is perturbed by a white noise. On the one hand, in the8 left picture we take D= 3.5, sin = 1, a= 0.8, m= 1.5 and α= 0.5. In this9 case, condition (3.8) is fulfilled so we can observe that the microorganisms become10 extinct. On the other hand, in the right picture we choose D= 2, sin = 1, a= 0.6,11 m= 5 and α= 0.5. In this second case ¯ D < µ(sin) holds, then we can see12 that the species persist even though it is not possible to prove mathematically the13 persistence. In addition, we can also see that some realizations in the right picture14 take negative values which is another significant drawback.15 20 T. CARABALLO, M.-J. GARRIDO-ATIENZA, J. L ´ OPEZ-DE-LA-CRUZ, A. RAPAPORT Figure 9. Stochastic chemostat model. Extinction (left) and persistence (right) Next we present two figures where we overlap a typical realization of the solution1 of system (3.1)-(3.2) and another one of the solution of system (2.4)-(2.5) such that2 we can notice much more easily the differences between the simulations concerning3 Sections 2and 3. In each figure we plot a big panel where the general dynamics4 can be seen and four smaller panels which correspond to two different zooms of5 two interesting places of the realizations, the dynamics around (s, x) = (2,2) and6 (s, x) = (2,0) in Figure 10 and the dynamics about (s, x) = (0.2,5.5) and (s, x) =7 (0.4,3.75) in Figure 11.8 In Figure 10 we plot a typical realization when perturbing the dilution rate with9 the Wiener process (orange) and two different ones when perturbing with the O-U10 process for β= 2 (red) and β= 0.5 (green). In this case, we take sin = 2, D= 3.5,11 a= 0.8, m= 0.5, α= 0.8, σ= 0.8, x(0) = 5 and s(0) = 2.5. We can observe that12 (2.10) and (3.8) are both fulfilled then the microorganisms become extinct, as we13 already proved in previous sections.14 Figure 10. Comparison in case of extinction Eventually, in Figure 11 we plot again a typical realization when perturbing15 the dilution rate with the Wiener process (orange) and two different ones when16 RANDOM AND STOCHATIC CHEMOSTAT MODELS 21 perturbing with the O-U process for β= 2 (red) and β= 0.5 (green) but now1 we take sin = 4, D= 2, a= 0.6, m= 5, α= 0.15, σ= 0.8, x(0) = 5 and2 s(0) = 2. In this case (3.8) does not hold true, thus it is not possible to ensure the3 persistence of the species (in the chemostat model perturbed by using the standard4 Wiener process) although numerically it can be obtained for the previous values5 of the parameters. In addition, if D < µ(sin) and ¯s < sin hold true, then we6 can ensure the persistence of the microbial biomass when perturbing the chemostat7 model by means of the O-U process. Moreover, we can observe that every realization8 is approaching the line s+x=sin, as proved in Sections 2(see (2.16)) and 3(see9 (3.10)).10 Figure 11. Comparison in case of persistence In conclusion, we can observe that the Ornstein-Uhlenbeck process provides us a11 really useful tool when modeling stochasticity and randomness since it allows us to12 set up mathematical models which guarantee the positiveness of the variable and13 therefore it better suits to represent reality. This new framework could also allow us14 to revisit the persistence of species under input disturbances, in case of competition15 between several species, which will be the topic of a future work.16 Appendix. Although very good references (see e.g. [2]) providing very detailed17 information about random dynamical systems (RDSs) can be found in the literature,18 we recall briefly some useful definitions and results to make our presentation as much19 self-contained as possible.20 Let (X,k · kX) be a separable Banach space.21 Definition 4.1. An RDS on Xconsists of two ingredients: (a) a metric dynami-22 cal system (Ω,F,P,{θt}t∈R)where (Ω,F,P)is a probability space and a family of23 mappings θt: Ω →Ωsatisfying24 (1) θ0=IdΩ,25 (2) θs◦θt=θs+tfor all s, t ∈R,26 (3) the mapping (t, ω)7→ θtωis measurable,27 (4) the probability measure Pis preserved by θt, i.e., θtP=P28 22 T. CARABALLO, M.-J. GARRIDO-ATIENZA, J. L ´ OPEZ-DE-LA-CRUZ, A. RAPAPORT and (b) a mapping ϕ: [0,∞)×Ω× X → X which is (B[0,∞)× F × B(X),B(X))-1 measurable, such that for each ω∈Ω,2 (i) the mapping ϕ(t, ω) : X → X,x7→ ϕ(t, ω)xis continuous for every t≥0,3 (ii) ϕ(0, ω)is the identity operator on X,4 (iii) (cocycle property) ϕ(t+s, ω) = ϕ(t, θsω)ϕ(s, ω)for all s, t ≥0.5 Definition 4.2. Let (Ω,F,P)be a probability space. A random set Kis a measu-6 rable subset of X × Ωwith respect to the product σ−algebra B(X)× F. Moreover K7 will be said a closed or a compact random set if K(ω) = {x: (x, ω)∈K}, ω ∈Ω,8 is closed or compact for P−almost all ω∈Ω, respectively.9 Definition 4.3. A bounded random set K(ω)⊂ X is said to be tempered with respect to {θt}t∈Rif for a.e. ω∈Ω, lim t→∞ e−βt sup x∈K(θ−tω) kxkX= 0,for all β > 0; a random variable ω7→ r(ω)∈Ris said to be tempered with respect to {θt}t∈Rif for a.e. ω∈Ω, lim t→∞ e−βt sup t∈R |r(θ−tω)|= 0,for all β > 0. In what follows we use E(X) to denote the set of all tempered random sets of X.10 Definition 4.4. A random set B(ω)⊂ X is called a random absorbing set in E(X) if for any E∈ E(X)and a.e. ω∈Ω, there exists TE(ω)>0such that ϕ(t, θ−tω)E(θ−tω)⊂B(ω),for all t≥TE(ω). Definition 4.5. Let {ϕ(t, ω)}t≥0,ω∈Ωbe an RDS over (Ω,F,P,{θt}t∈R)with state11 space Xand let A(ω)(⊂ X)be a random set. Then A={A(ω)}ω∈Ωis called a12 global random E−attractor (or pullback E−attractor) for {ϕ(t, ω)}t≥0,ω∈Ωif13 (i) (compactness) A(ω)is a compact set of Xfor any ω∈Ω;14 (ii) (invariance) for any ω∈Ωand all t≥0, it holds ϕ(t, ω)A(ω) = A(θtω); (iii) (attracting property) for any E∈ E(X)and a.e. ω∈Ω, lim t→∞ distX(ϕ(t, θ−tω)E(θ−tω), A(ω)) = 0, where distX(G, H) = supg∈Ginfh∈Hkg−hkXis the Hausdorff semi-metric15 for G, H ⊆ X.16 Proposition 4.1. [See [14,17]] Let B∈ E(X)be a closed absorbing set for the continuous RDS {ϕ(t, ω)}t≥0,ω∈Ωthat satisfies the asymptotic compactness condition for a.e. ω ∈Ω, i.e., each sequence xnin ϕ(tn, θ−tnω)B(θ−tnω)has a convergent subsequence in Xwhen tn→ ∞. Then ϕhas a unique global random attractor A={A(ω)}ω∈Ωwith component subsets A(ω) = \ τ≥TB(ω)[ t≥τ ϕ(t, θ−tω)B(θ−tω). Remark 4.1. When the state space X=Rdas in this paper, the asymptotic com-17 pactness follows trivially.18 The next result ensures when two RDSs are conjugated (see [10,11,12]).19 RANDOM AND STOCHATIC CHEMOSTAT MODELS 23 Proposition 4.2. Let ϕube an RDS on X. Suppose that the mapping T: Ω×X →1 Xpossesses the following properties: for fixed ω∈Ω,T(ω, ·)is a homeomorphism2 on X, and for x∈ X , the mappings T(·, x),T−1(·, x)are measurable. Then the3 mapping4 (t, ω, x)→ϕv(t, ω)x:= T−1(θtω, ϕu(t, ω)T(ω, x)) is a (conjugated) RDS.5 REFERENCES6 [1] S. Al-azzawi, J. Liu and X. Liu, Convergence rate of synchronization of systems with additive7 noise, Discrete Contin. Dyn. Syst. Ser. B,22 (2017), 227–245.8 [2] L. Arnold, Random Dynamical Systems, Springer Berlin Heidelberg, 1998.9 [3] J. Barlow, W. Schaffner, F. de Noyelles, B. Peterson and J. Peterson, Continuous flow nutrient10 bioassays with natural phytoplankton populations, G. Glass (Editor): Bioassay Techniques11 and Environmental Chemistry, John Wiley & Sons Ltd., 1973.12 [4] F. Campillo, M. Joannides and I. Larramendy-Valverde, Stochastic modeling of the chemostat,13 Ecological Modelling,222 (2011), 2676–2689.14 [5] F. Campillo, M. Joannides and I. Larramendy-Valverde, Approximation of the Fokker–Planck15 equation of the stochastic chemostat, Mathematics and Computers in Simulation,99 (2014),16 37–53.17 [6] F. Campillo, M. Joannides and I. Larramendy-Valverde, Analysis and approximation of a18 stochastic growth model with extinction, Methodology and Computing in Applied Probability,19 18 (2016), 499–515.20 [7] T. Caraballo, M. J. Garrido-Atienza and J. L´opez-de-la-Cruz, Some Aspects Concerning21 the Dynamics of Stochastic Chemostats, vol. 69, 227–246, Springer International Publishing,22 Cham, 2016.23 [8] T. Caraballo, M. J. Garrido-Atienza and J. L´opez-de-la-Cruz, Dynamics of some stochastic24 chemostat models with multiplicative noise, Communications on Pure and Applied Analysis,25 16 (2017), 1893–1914.26 [9] T. Caraballo, M. J. Garrido-Atienza, J. L´opez-de-la-Cruz and A. Rapaport, Corrigendum to27 ”Some aspects concerning the dynamics of stochastic chemostats”, 2017, arXiv:1710.0077428 [math.DS].29 [10] T. Caraballo, M. J. Garrido-Atienza, B. Schmalfuss and J. Valero, Asymptotic behavior of a30 stochastic semilinear dissipative functional equation without uniqueness of solutions, Discrete31 and Continuous Dynamical Systems - Series B,14 (2010), 439–455.32 [11] T. Caraballo and X. Han, Applied Nonautonomous and Random Dynamical Systems, Applied33 Dynamical Systems, Springer International Publishing, 2016.34 [12] T. Caraballo, P. E. Kloeden and B. Schmalfuss, Exponentially stable stationary solutions for35 stochastic evolution equations and their perturbation, Applied Mathematics and Optimiza-36 tion,50 (2004), 183–207.37 [13] T. Caraballo and K. Lu, Attractors for stochastic lattice dynamical systems with a multi-38 plicative noise, Front. Math. China,3(2008), 317–335.39 [14] T. Caraballo, G. Lukaszewicz and J. Real, Pullback attractors for asymptotically compact40 non-autonomous dynamical systems, Nonlinear Analysis: Theory, Methods & Applications,41 64 (2006), 484–498.42 [15] I. F. Creed, D. M. McKnight, B. A. Pellerin, M. B. Green, B. A. Bergamaschi, G. R. Aiken,43 D. A. Burns, S. E. G. Findlay, J. B. Shanley, R. G. Striegl, B. T. Aulenbach, D. W. Clow,44 H. Laudon, B. L. McGlynn, K. J. McGuire, R. A. Smith and S. M. Stackpoole, The river as a45 chemostat: fresh perspectives on dissolved organic matter flowing down the river continuum,46 Canadian Journal of Fisheries and Aquatic Sciences,72 (2015), 1272–1285.47 [16] G. D’Ans, P. Kokotovic and D. Gottlieb, A nonlinear regulator problem for a model of bio-48 logical waste treatment, IEEE Transactions on Automatic Control,16 (1971), 341–347.49 [17] F. Flandoli and B. Schmalfuss, Random attractors for the 3D stochastic navier-stokes equation50 with multiplicative white noise, Stochastics and Stochastic Reports,59 (1996), 21–45.51 [18] J. Grasman, M. D. Gee and O. A. V. Herwaarden, Breakdown of a chemostat exposed to52 stochastic noise, Journal of Engineering Mathematics,53 (2005), 291–300.53 24 T. CARABALLO, M.-J. GARRIDO-ATIENZA, J. L ´ OPEZ-DE-LA-CRUZ, A. RAPAPORT [19] J. Harmand, C. Lobry, A. Rapaport and T. Sari, The Chemostat: Mathematical Theory of1 Micro-organisms Cultures, Wiley, Chemical Engineering Series, John Wiley & Sons, Inc.,2 2017.3 [20] D. J. Higham., An algorithmic introduction to numerical simulation of stochastic differential4 equations, SIAM Review,43 (2001), 525–546.5 [21] S. B. Hsu, S. Hubbell and P. Waltman, A mathematical theory for single-nutrient competition6 in continuous cultures of micro-organisms, SIAM Journal on Applied Mathematics,32 (1977),7 366–383.8 [22] L. Imhof and S. Walcher, Exclusion and persistence in deterministic and stochastic chemostat9 models, Journal of Differential Equations,217 (2005), 26–53.10 [23] H. W. Jannasch, Steady state and the chemostat in ecology, Limnology and Oceanography,11 19 (1974), 716–720.12 [24] J. Kalff and R. Knoechel, Phytoplankton and their dynamics in oligotrophic and eutrophic13 lakes, Annual Review of Ecology and Systematics,9(1978), 475–495.14 [25] J. W. M. La Rivi`ere, Microbial ecology of liquid waste treatment, in Advances in Microbial15 Ecology, vol. 1, Springer US, 1977, 215–259.16 [26] E. Rurangwa and M. C. J. Verdegem, Microorganisms in recirculating aquaculture systems17 and their management, Reviews in Aquaculture,7(2015), 117–130.18 [27] H. L. Smith and P. Waltman, The theory of the chemostat: dynamics of microbial competition,19 Cambridge University Press, 1995.20 [28] L. Wang and D. Jiang, Periodic solution for the stochastic chemostat with general response21 function, Physica A: Statistical Mechanics and its Applications,486 (2017), 378–385.22 [29] L. Wang, D. Jiang and D. O’Regan, The periodic solutions of a stochastic chemostat model23 with periodic washout rate, Communications in Nonlinear Science and Numerical Simula-24 tion,37 (2016), 1–13.25 [30] C. Xu and S. Yuan, An analogue of break-even concentration in a simple stochastic chemostat26 model, Applied Mathematics Letters,48 (2015), 62–68.27 [31] C. Xu, S. Yuan and T. Zhang, Asymptotic behavior of a chemostat model with stochastic28 perturbation on the dilution rate, Abstract and Applied Analysis,2013 (2013), 1–11.29 [32] D. Zhao and S. Yuan, Critical result on the break-even concentration in a single-species30 stochastic chemostat model, Journal of Mathematical Analysis and Applications,434 (2016),31 1336–1345.32 [33] D. Zhao and S. Yuan, Break-even concentration and periodic behavior of a stochastic chemo-33 stat model with seasonal fluctuation, Communications in Nonlinear Science and Numerical34 Simulation,46 (2017), 62–73.35 E-mail address:[email protected] E-mail address:[email protected] E-mail address:[email protected] E-mail address:[email protected]