# Study of the dynamics of two chemostats connected by Fickian diffusion with bounded random fluctuations

Tomás Caraballo, Javier López-De-La-Cruz, Alain Rapaport

► To cite this version:

Tomás Caraballo, Javier López-De-La-Cruz, Alain Rapaport. Study of the dynamics of two chemostats connected by Fickian diffusion with bounded random fluctuations. Stochastics and Dynamics, 2022, 22 (3), 10.1142/S0219493722400020 . hal-03359279

HAL Id: hal-03359279

https://hal.inrae.fr/hal-03359279

Submitted on 30 Sep 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
STUDY OF THE DYNAMICS OF TWO CHEMOSTATS
CONNECTED BY FICKIAN DIFFUSION WITH BOUNDED
RANDOM FLUCTUATIONS

TOMÁS CARABALLO

Dpto. de Ecuaciones Diferenciales y Análisis Numérico
C/ Tarfía s/n, Facultad de Matemáticas, Universidad de Sevilla
41012 Sevilla, Spain
caraball@us.es

JAVIER LÓPEZ-DE-LA-CRUZ*

Dpto. de Matemática Aplicada a las TIC
Escuela Técnica Superior de Ingenieros Informáticos
Campus de Montegancedo, Universidad Politécnica de Madrid
28660 Boadilla del Monte, Madrid, Spain
javier.lopez.delacruz@upm.es

ALAIN RAPAPORT

MISTEA, Univ. Montpellier, Inrae
Montpellier SupAgro 2, Place Pierre Viala
34060 Montpellier, France
alain.rapaport@inrae.fr

Received (Day Month Year)

Revised (Day Month Year)

Accepted (Day Month Year)

Communicated by (xxxxxxxxxx)

Dedicated to the memory of María J. Garrido Atienza**

This paper investigates the dynamics of a model of two chemostats connected by Fickian
diffusion with bounded random fluctuations. We prove the existence and uniqueness of

<footnote>*Corresponding author: javier.lopez.delacruz@upm.es.</footnote>

<footnote>**January 24, 2021 was a very sad day for us: our good friend and collaborator passed away on the day of her birthday. She not only was an excellent teacher and researcher, but a great person full of vitality and happiness. We will always remember her with all our deep appreciation and love.</footnote>
non-negative global solution as well as the existence of compact absorbing and attracting sets for the solutions of the corresponding random system. After that, we study the internal structure of the attracting set to obtain more detailed information about the long-time behavior of the state variables. In such a way, we provide conditions under which the extinction of the species cannot be avoided and conditions to ensure the persistence of the species, which is often the main goal pursued by practitioners. In addition, we illustrate the theoretical results with several numerical simulations.

Keywords: chemostat; random; diffusion; Ornstein-Uhlenbeck; interconnection.

AMS Subject Classification 2020: 34A12, 34F05, 34D99.

# 1. Introduction

The chemostat is a laboratory device used for growing microorganisms in a culture environment. It provides a very useful tool to investigate microbial ecosystems at steady state and possesses many applications in real life, for instance, it can be used to study genetically altered organisms (see [40, 42]), models of mammalian large intestine (see [21, 22]), waste-water treatment processes (see [15, 30]) and models in ecology regarding populations of microorganisms in lakes, wetlands, rivers or aquaculture systems (see [2, 13, 26, 27, 28, 37]).

Since its invention, thanks to Monod (see [34]) and Novick and Szilard (see [36]) simultaneously in 1950, the chemostat has been subject to a huge number of scientific publications, not only in applied mathematics but also in applied sciences, being the main focus of attention of many books and research articles.

The success enjoyed by the chemostat is mainly due to the fact that it is a quite simple device that allows scientists to develop several different works that reproduce the real life in a very loyal manner, in fact, it is widely used for applications in industry. In addition, the chemostat model is very interesting as a mathematical object as well and constitutes an active branch of biomathematics, under framework called *the theory of the chemostat*.

The classical chemostat device consists of three tanks, called the *feed bottle*, the *culture vessel* and the *collection vessel*, respectively, that are interconnected by means of pumps. The substrate is stored in the first tank and pumped to the culture vessel, where the interactions with the species take place. Moreover, another flow is pumped from the culture vessel to the collection vessel such that the volume in the culture vessel remains constant. Figure 1 provides a representation of the classical chemostat device. [5]
<figure><img src="image_2.png" /><figcaption>Fig. 1: Diagram of the chemostat device</figcaption></figure>

This dynamics is described by the following differential system (see [25])

$$ \frac{ds}{dt} = \frac{Q}{V}(s^{in} - s) - \mu(s)x, \qquad (1.1) $$

$$ \frac{dx}{dt} = -\frac{Q}{V}x + \mu(s)x, \qquad (1.2) $$

where s and x denote the concentration of substrate (or nutrient) and species (or microorganisms), respectively, $s^{in}$ is the input concentration of substrate, Q represents the input flow rate, V describes the volume of liquid media in the culture vessel and $\mu$ is a function describing the kinetics of the substrate is being transformed by the microbial species. In this paper we consider the classical Monod consumption function

$$ \mu(s) = \frac{\bar{\mu}_0 s}{k + s}, \qquad (1.3) $$

where $\bar{\mu}_0$ is the maximal growth rate coefficient of the consumer species and k represents the half-saturation constant. We recall that every parameter is positive.

We shall later on denote $D = Q/V$ the dilution rate.

System (1.1)-(1.2) has been widely studied in the last decades from both the biological and the mathematical points of view. Among the different assumptions made when considering system (1.1)-(1.2), one is to assume that the liquid media inside the culture vessel is perfectly mixed. This is possible by using a perfect agitator when having culture vessels with small volumes but it becomes difficult to achieve if the volume of the culture vessel is large, as it often happens in industrial devices or natural reservoirs such as lakes.

This encourages scientists to think about modeling non ideally mixed chemostats. To this end, some of the most common ways consists on considering systems of partial differential equations (see [29, 18, 35]) or a finite number of interconnected chemostats with different flow conditions described by systems of ordinary differential equations, as in the general gradostat model (see [38, 41, 24]).
Concerning the existing literature, most of them consider spatial heterogeneity only in the axial dimension of the chemostat as in tubular devices, where the culture vessel is very high and thin and the input flow is considerable large, (see [23, 17, 19, 45, 16]) and simple gradostats (see [32, 33, 43, 39]). Nevertheless, configurations in parallel have been much less treated, apart from simple considerations in chemical reaction engineering (see [31, 20]). However, it has been reported in the literature that tanks often present dead zones that can be simply represented by a part of the volume connected to the rest of the volume by diffusion (see [14] and reference herein).

Hence, in this paper our aim is to consider two classical chemostats where the culture vessels are connected in parallel by Fickian diffusion, as illustrated in Figure 2. This amounts to have a bi-directional flow rate between the tanks, preserving constant volumes in the tanks.

<figure><img src="image_3.png" /><figcaption>Fig. 2: Diagram of two chemostat devices connected by Fickian diffusion</figcaption></figure>

Therefore, the resulting model is given by the following differential system

$$ \frac{ds_1}{dt} = \frac{Q_1}{V_1}(s_1^{in} - s_1) - \mu(s_1)x_1 + \frac{d}{V_1}(s_2 - s_1), \quad (1.4) $$

$$ \frac{dx_1}{dt} = -\frac{Q_1}{V_1}x_1 + \mu(s_1)x_1 + \frac{d}{V_1}(x_2 - x_1), \quad (1.5) $$

$$ \frac{ds_2}{dt} = \frac{Q_2}{V_2}(s_2^{in} - s_2) - \mu(s_2)x_2 + \frac{d}{V_2}(s_1 - s_2), \quad (1.6) $$

$$ \frac{dx_2}{dt} = -\frac{Q_2}{V_2}x_2 + \mu(s_2)x_2 + \frac{d}{V_2}(x_1 - x_2), \quad (1.7) $$

where $s_i$ and $x_i$ denote the concentration of substrate and species in the culture vessel $i = 1, 2$, respectively, $s_i^{in}$ is the input concentration of substrate in the culture vessel $i = 1, 2$, $Q_i$ represents the input flow rate in the culture vessel $i = 1, 2$, $V_i$
describes the volume of the culture vessel $i = 1, 2$ and $d$ is the diffusion parameter. We recall that every parameter is positive.

Even though system (1.4)-(1.7) could provide a good approximation of real devices, it is a deterministic model and we know that phenomena in real life are subject to random fluctuations. Then, in order to obtain a more realistic model, we go further in this paper and consider random perturbations of the diffusion parameter.

There are many different ways to introduce randomness or stochasticity in deterministic systems (see, for instance, [6]). However, in this work we use a bounded function of the well-known Ornstein-Uhlenbeck process to ensure that the random perturbations remain bounded in a certain interval, as it happens in real devices. Then, we obtain the following system of random differential equations

$$ \frac{ds_1}{dt} = \frac{Q_1}{V_1}(s_1^{in} - s_1) - \mu(s_1)x_1 + \frac{d + \Phi(z^*(\theta_t \omega))}{V_1}(s_2 - s_1), \quad (1.8) $$

$$ \frac{dx_1}{dt} = -\frac{Q_1}{V_1}x_1 + \mu(s_1)x_1 + \frac{d + \Phi(z^*(\theta_t \omega))}{V_1}(x_2 - x_1), \quad (1.9) $$

$$ \frac{ds_2}{dt} = \frac{Q_2}{V_2}(s_2^{in} - s_2) - \mu(s_2)x_2 + \frac{d + \Phi(z^*(\theta_t \omega))}{V_2}(s_1 - s_2), \quad (1.10) $$

$$ \frac{dx_2}{dt} = -\frac{Q_2}{V_2}x_2 + \mu(s_2)x_2 + \frac{d + \Phi(z^*(\theta_t \omega))}{V_2}(x_1 - x_2), \quad (1.11) $$

where $\Phi(z^*(\theta_t \omega))$ denotes the noise, which will be presented in more detail in Section 2.

One may wonder about the choice of such a noise. Indeed, it has been shown to be a relevant way of modeling real random fluctuations in chemostat models (see, for instance [6, 44, 3, 11, 12, 9, 10], where the authors perturb different parameters in several chemostat models and other models in population dynamics).

We remark that it would be also possible to consider, for instance, the standard Wiener process. However, its realizations are unbounded, in fact it can take arbitrary large values, both positive or negative, and this leads into several drawbacks from the biological point of view. We refer to interested readers to [5, 6, 4], where the authors explain in more detail the inconvenience found when perturbing different chemostat models by the standard Brownian motion.

The rest of the paper is organized as follows. In Section 2 we present some preliminaries about the Ornstein-Uhlenbeck process and how to model bounded random fluctuations. Then, in Section 3 we prove the existence and uniqueness of non-negative global solution of the random system (1.8)-(1.11). After that, in Section 4 we prove the existence of deterministic compact absorbing and attracting sets for the solutions of the random system (1.8)-(1.11). Then, we study the internal structure of the attracting set in Section 5 and provide conditions under which the
extinction of the species cannot be avoided and conditions to ensure the persistence of the species, the main goal pursued by practitioners. Next, in Section 6 we illustrate the theoretical results with some numerical simulations. Finally, in Section 7 we present conclusions to sum up our contributions in this work.

## 2. Preliminaries

In this section we present some preliminaries about the Ornstein-Uhlenbeck process, including several properties which will be useful in the rest of the paper and we also explain the way in which we model the random bounded fluctuations.

### 2.1. The Ornstein-Uhlenbeck process

In this section we present the Ornstein-Uhlenbeck (OU) process and some of its properties. We refer readers to [3, 6, 7] for details.

Let $W$ be a two sided standard Wiener process. Kolmogorov's theorem ensures that $W$ has a continuous version, denoted usually by $\omega$, whose canonical interpretation is as follows. Let $\Omega$ be defined by

$$ \Omega := \{\omega \in \mathcal{C}(\mathbb{R}, \mathbb{R}) : \omega(0) = 0\}, $$

the space of continuous functions that are zero at zero, $\mathcal{F}$ be the Borel $\sigma$-algebra of measurable subsets of $\Omega$ generated by the compact open topology (see [1]) and $\mathbb{P}$ the Wiener measure on $\mathcal{F}$. Then, we have the probability space $(\Omega, \mathcal{F}, \mathbb{P})$.

In addition, to track the effect of the noise in time, i.e., each state $\omega$ at $t=0$ with its state after some time $t$, we define the mapping $\theta_t: \Omega \to \Omega$, which is called the Wiener shift flow, given by

$$ \theta_t \omega(\cdot) = \omega(\cdot + t) - \omega(t), \quad t \in \mathbb{R}. $$

The OU process is a stationary mean-reverting Gaussian stochastic process defined as the random variable

$$ (t, \omega) \mapsto z(t, \omega) := z^*(\theta_t \omega) = -\beta \nu \int_{-\infty}^{0} e^{\beta s} \theta_t \omega(s) ds, \quad (2.1) $$

for all $t \in \mathbb{R}$ and $\omega \in \Omega$, where $\beta > 0$ is the mean reversion constant representing the strength with which the process is attracted by the mean and $\nu > 0$ is the volatility constant describing the variation or the size of the noise.

We note that the OU process (2.1) can be obtained as the stationary solution of the Langevin equation

$$ dz + \beta z dt = \nu d\omega. \quad (2.2) $$
Usually, the OU process (2.1) can be used to 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 ($\beta = 0, \nu = 1$) and no noise at all ($\beta = 1, \nu = 0$).

We recall in the next proposition some of its properties.

**Proposition 2.1 (See [1, 8]).** There exists a $\theta_t$-invariant set $\tilde{\Omega} \in \mathcal{F}$ of $\Omega$ of full measure such that for $\omega \in \tilde{\Omega}$ and $\beta, \nu > 0$, we have

(i) the random variable $|z^*(\omega)|$ is tempered with respect to $\{\theta_t\}_{t \in \mathbb{R}}$, i.e., for a.e. $\omega \in \tilde{\Omega}$,

$$ \lim_{t \to \infty} e^{-\eta t} \sup_{t \in \mathbb{R}} |z^*(\theta_{-t}\omega)| = 0, \quad \text{for all } \eta > 0. $$

(ii) this mapping is a stationary solution of (2.2) with continuous trajectories

$$ (t, \omega) \to z^*(\theta_t \omega) = -\beta \nu \int_{-\infty}^{0} e^{\beta s} (\theta_t \omega)(s) ds. $$

(iii) for any $\omega \in \tilde{\Omega}$ one has

$$ \lim_{t \to \pm \infty} \frac{|z^*(\theta_t \omega)|}{t} = 0; \quad \lim_{t \to \pm \infty} \frac{1}{t} \int_0^t z^*(\theta_s \omega) ds = 0; $$

$$ \lim_{t \to \pm \infty} \frac{1}{t} \int_0^t |z^*(\theta_s \omega)| ds = \mathbb{E}[|z^*|] < \infty. $$

From now on we consider $\beta$ and $\nu$ fixed.

## 2.2. Modeling random bounded fluctuations

In this section we present a way to model bounded random fluctuations which has been proved to be very useful when modeling real random disturbances in biological models (see [3, 9, 10, 12]) fitting the real devices.

Let us consider an interval $[\underline{d}, \bar{d}]$, where $0 < \underline{d} < d < \bar{d} < \infty$.

Define now the function $\Phi: \mathbb{R} \to [-\alpha, \alpha]$ given by

$$ \Phi(z) = \frac{2\alpha}{\pi} \arctan(z), \qquad (2.3) $$

where $\alpha = \bar{d} - d = d - \underline{d} > 0$.

Then, the random diffusion parameter satisfies

$$ \underline{d} < d + \Phi(z^*(\theta_t \omega)) < \bar{d} \quad \text{for all } t \in \mathbb{R}, \qquad (2.4) $$
i.e., $d + \Phi(z^*(\theta_t\omega))$ remains positive and bounded in the interval $[\underline{d}, \bar{d}]$ for every time.

In addition, we have the following ergodic property which will be very useful.

**Proposition 2.2 (See [12]).** Let be $\Phi(z) = \frac{2\alpha}{\pi} \arctan(z)$ defined as in (2.3). Then

$$ \lim_{t \to +\infty} \frac{1}{t} \int_{0}^{t} \Phi(z^*(\theta_s \omega)) ds = 0, \quad \text{a.s. in } \Omega. \qquad (2.5) $$

**Remark 2.1.** We would like to notice that (2.5) remains true as long as $\Phi$ is an odd measurable function such that

$$ \lim_{z \to +\infty} \Phi(z) = \alpha > 0. $$

**Remark 2.2.** It worth mentioning that the random systems obtained when using the previous way to model bounded random perturbations, generate random dynamical systems and then we could use the very well-known theory of random dynamical systems and pullback attractors (see [7]) to study the corresponding random models. However, in this paper we investigate the random chemostat (1.8)-(1.11) for every fixed realization of the noise $\omega \in \Omega$ since it allows us to prove every result in forward sense, which is more natural than the pullback one in this case.

### 3. Existence and uniqueness of non-negative global solution

In this section we prove that the random chemostat model (1.8)-(1.11) has a unique global solution which remains non-negative for every non-negative initial condition.

From now on, we denote $\mathcal{X} := \{(s_1, x_1, s_2, x_2) \in \mathbb{R}^4 : s_1 \ge 0, x_1 \ge 0, s_1 \ge 0, s_2 \ge 0\}$ the positive cone in $\mathbb{R}^4$.

For the sake of clarity, we present in Table 1 some new notations which will allow us to simplify the expressions in what follows. We posit $D_i = Q_i/V_i$ and $d_i = d/V_i$ for $i \in \{1, 2\}$.

$$ D^{\max} = \max\{D_1, D_2\} \qquad d^{\max} = \max\{d_1, d_2\} \qquad V_{\text{inv}}^{\max} = \max\left\{\frac{1}{V_1}, \frac{1}{V_2}\right\} $$

$$ D^{\min} = \min\{D_1, D_2\} \qquad d^{\min} = \min\{d_1, d_2\} \qquad V_{\text{inv}}^{\min} = \min\left\{\frac{1}{V_1}, \frac{1}{V_2}\right\} $$

<caption>Table 1: Definition of some new parameters</caption>
**Theorem 3.1.** For any initial value $u_0 := (s_{10}, x_{10}, s_{20}, x_{20}) \in \mathcal{X}$ and $\omega \in \Omega$, the random chemostat model (1.8)-(1.11) possesses a unique global solution

$$u(\cdot; 0, \omega, u_0) := (s_1(\cdot; 0, \omega, u_0), x_1(\cdot; 0, \omega, u_0), s_2(\cdot; 0, \omega, u_0), x_2(\cdot; 0, \omega, u_0)),$$

where $u(\cdot; 0, \omega, u_0) \in C^1([0, +\infty), \mathcal{X})$, with $u(0; 0, \omega, u_0) = u_0$, $s_{10} = s_1(0; 0, \omega, u_0)$, $x_{10} = x_1(0; 0, \omega, u_0)$, $s_{20} = s_2(0; 0, \omega, u_0)$ and $x_{20} = x_2(0; 0, \omega, u_0)$.

**Proof.** To prove this theorem, we first rewrite system (1.8)-(1.11) as

$$\frac{du}{dt} = L(\theta_t \omega) u + F(u),$$

where $u := (s_1, x_1, s_2, x_2)^t$, $L(\theta_t \omega)$ is given by

$$\left( \begin{array}{cccc} -D_1 - d_1(\theta_t \omega) & 0 & d_1(\theta_t \omega) & 0 \\ 0 & -D_1 - d_1(\theta_t \omega) & 0 & d_1(\theta_t \omega) \\ d_2(\theta_t \omega) & 0 & -D_2 - d_2(\theta_t \omega) & 0 \\ 0 & d_2(\theta_t \omega) & 0 & -D_2 - d_2(\theta_t \omega) \end{array} \right),$$

where $d_i(\theta_t \omega) := (d + \Phi(z^*(\theta_t \omega))) / V_i$, $i \in \{1, 2\}$, and

$$F(u) = \begin{pmatrix} D_1 s_1^{in} - \mu(s_1) x_1 \\ \mu(s_1) x_1 \\ D_2 s_2^{in} - \mu(s_2) x_2 \\ \mu(s_2) x_2 \end{pmatrix}.$$

Since $z^*(\theta_t \omega)$ is continuous with respect to $t$, thanks to Proposition 2.1 (ii), and $\Phi$ is continuous, then $L(\theta_t \omega)$ is continuous with respect to $t$ and generates an evolution system on $\mathbb{R}^4$. Moreover, $F(u)$ is continuously differentiable with respect to $s_1, x_1, s_2$ and $x_2$, then it is locally Lipschitz with respect to $s_1, x_1, s_2$ and $x_2$. Hence, system (1.8)-(1.11) possesses a unique local solution.

Now, we prove that every solution of system (1.8)-(1.11) remains in $\mathcal{X}$ for every initial value in $\mathcal{X}$. To this end, we consider $s_1 \ge 0, x_1 \ge 0, s_2 \ge 0$ and $x_2 \ge 0$.
Then, we have

$$
\begin{align*}
\left. \frac{ds_1}{dt} \right|_{s_1=0} &= D_1 s_1^{in} + \frac{d + \Phi(z^*(\theta_t \omega))}{V_1} s_2 > 0, \\
\left. \frac{dx_1}{dt} \right|_{x_1=0} &= \frac{d + \Phi(z^*(\theta_t \omega))}{V_1} x_2 \ge 0, \\
\left. \frac{ds_2}{dt} \right|_{s_2=0} &= D_2 s_2^{in} + \frac{d + \Phi(z^*(\theta_t \omega))}{V_2} s_1 > 0, \\
\left. \frac{dx_2}{dt} \right|_{x_2=0} &= \frac{d + \Phi(z^*(\theta_t \omega))}{V_2} x_1 \ge 0,
\end{align*}
$$

whence we have that every solution starting in $\mathcal{X}$ remains there for every $t \ge 0$.

To finish the proof, we consider the variable $y := s_1 + x_1 + s_2 + x_2$ which satisfies
the following random differential equation

$$
\begin{equation}
\begin{split}
\frac{dy}{dt} ={}& D_1 s_1^{in} + D_2 s_2^{in} - D_1 s_1 - D_1 x_1 - D_2 s_2 - D_2 x_2 \\
& + \frac{d + \Phi(z^*(\theta_t \omega))}{V_1} (s_2 + x_2) + \frac{d + \Phi(z^*(\theta_t \omega))}{V_2} (s_1 + x_1) \\
& - \frac{d + \Phi(z^*(\theta_t \omega))}{V_1} (s_1 + x_1) - \frac{d + \Phi(z^*(\theta_t \omega))}{V_2} (s_2 + x_2).
\end{split}
\tag{3.1}
\end{equation}
$$

From (3.1), it is not difficult to obtain

$$
\frac{dy}{dt} \leq D_1 s_1^{in} + D_2 s_2^{in} + ((V_{\text{inv}}^{\text{max}} - V_{\text{inv}}^{\text{min}})(d + \Phi(z^*(\theta_t \omega))) - D^{\text{min}})y, \quad (3.2)
$$

where $V_{\text{inv}}^{\text{max}}$, $V_{\text{inv}}^{\text{min}}$ and $D^{\text{min}}$ are defined as in Table 1.

We note that, in order to obtain (3.2), we used (a) $D_i \ge D^{\min}$, where $i = 1, 2$, (b) $\frac{1}{V_i} \le V_{\text{inv}}^{\max}$ for $i = 1, 2$, (c) $\frac{1}{V_i} \ge V_{\text{inv}}^{\min}$, with $i = 1, 2$, and (d) $d + \Phi(z^*(\theta_t \omega)) > 0$ thanks to (2.4).

As $d + \Phi(z^*(\theta_t \omega))$ is bounded from above by $\bar{d}$, we can directly use the Gronwall
Lemma to conclude that

$$
\begin{equation}
\begin{split}
y(t; 0, \omega, y_0) \le {}& y_0 e^{-(D^{\min}+d(V_{\text{inv}}^{\min}-V_{\text{inv}}^{\max})t-(V_{\text{inv}}^{\min}-V_{\text{inv}}^{\max})\int_0^t \Phi(z^*(\theta_\tau \omega))d\tau)} \\
& + (D_1 s_1^{\text{in}} + D_2 s_2^{\text{in}}) \int_0^t e^{-(D^{\min}+d(V_{\text{inv}}^{\min}-V_{\text{inv}}^{\max})(t-s)-(V_{\text{inv}}^{\min}-V_{\text{inv}}^{\max})\int_s^t \Phi(z^*(\theta_\tau \omega))d\tau)} ds
\end{split}
\tag{3.3}
\end{equation}
$$

for every $\omega \in \Omega$, $y_0 := s_{10} + x_{10} + s_{20} + x_{20} > 0$ and $t \ge 0$.

From (3.3), we notice that $y \ge 0$ does not blow up at finite time, then $s_1, x_1, s_2$ and $x_2$ do not blow up at finite time either. In addition, every state variable remain non-negative for every time. Then, the unique local solution is a global one. $\square$
### 4. Existence of absorbing and attracting sets

In this section we to study the existence of compact absorbing and attracting sets for the solutions of the random system (1.8)-(1.11).

**Theorem 4.1.** Assume that

$$D^{\min} + \bar{d}(V_{inv}^{\min} - V_{inv}^{\max}) > 0 \qquad (4.1)$$

holds true. Then, for every $\varepsilon > 0$, the deterministic set

$$\mathcal{B}_{\varepsilon} := \{(s_1, x_1, s_2, x_2) \in \mathcal{X} : y_l^* - \varepsilon \le s_1 + x_1 + s_2 + x_2 \le y_u^* + \varepsilon\} \qquad (4.2)$$

is a compact absorbing set for the solutions of the random chemostat model (1.8)-(1.11), where $y_l^*$ and $y_u^*$ are given by

$$y_l^* := \frac{D_1 s_1^{in} + D_2 s_2^{in}}{D^{\max} + \bar{d}(V_{inv}^{\max} - V_{inv}^{\min})} \quad \text{and} \quad y_u^* := \frac{D_1 s_1^{in} + D_2 s_2^{in}}{D^{\min} + \bar{d}(V_{inv}^{\min} - V_{inv}^{\max})}. \qquad (4.3)$$

**Proof.** We split this proof into two parts. The first one focuses on obtaining the upper bound $y_u^*$. To this end, we notice that $V_{inv}^{\min} - V_{inv}^{\max} < 0$ whence $-(V_{inv}^{\min} - V_{inv}^{\max}) > 0$. In addition, thanks to (2.4), we know that $d + \Phi(z^*(\theta_t \omega)) \le \bar{d}$ for every $t \in \mathbb{R}$. By taking into account these considerations, from (3.3) we deduce that

$$y(t; 0, \omega, y_0) \le y_0 e^{-(D^{\min} + \bar{d}(V_{inv}^{\min} - V_{inv}^{\max}))t} + \frac{D_1 s_1^{in} + D_2 s_2^{in}}{D^{\min} + \bar{d}(V_{inv}^{\min} - V_{inv}^{\max})} \left[1 - e^{-(D^{\min} + \bar{d}(V_{inv}^{\min} - V_{inv}^{\max}))t}\right] \qquad (4.4)$$

for every $t \ge 0$, $\omega \in \Omega$ and $y_0 > 0$.

Hence, thanks to (4.4), as long as (4.1) is satisfied, for any $\varepsilon > 0$, $\omega \in \Omega$ and every initial value $u_0 \in F$, where $F$ denotes a bounded set, there exists some time $T_F(\omega, \varepsilon) > 0$ such that

$$y(t; 0, \omega, y_0) \le \frac{D_1 s_1^{in} + D_2 s_2^{in}}{D^{\min} + \bar{d}(V_{inv}^{\min} - V_{inv}^{\max})} + \varepsilon$$

for every $t \ge T_F(\varepsilon, \omega)$.

The second part of the proof focuses on obtaining the lower bound $y_l^*$. To this end, from (3.1), similarly to the first part, we obtain

$$\frac{dy}{dt} \ge D_1 s_1^{in} + D_2 s_2^{in} + ((V_{inv}^{\min} - V_{inv}^{\max})(d + \Phi(z^*(\theta_t \omega))) - D^{\max})y, \qquad (4.5)$$

where $V_{inv}^{\min}$, $V_{inv}^{\max}$ and $D^{\max}$ are given as in Table 1.

We notice that, in order to have (4.5), we used (a) $D_i \le D^{\max}$, then $-D_i \ge -D^{\max}$, where $i = 1, 2$, (b) $\frac{1}{V_i} \ge V_{inv}^{\min}$ for $i = 1, 2$, (c) $\frac{1}{V_i} \le V_{inv}^{\max}$, whence $-\frac{1}{V_i} \ge -V_{inv}^{\max}$ with $i = 1, 2$, and (d) $d + \Phi(z^*(\theta_t \omega)) > 0$ thanks to (2.4).
As $d + \Phi(z^*(\theta_t\omega))$ is bounded from below by $\underline{d}$, we can directly use the Gronwall Lemma to conclude that

$$y(t; 0, \omega, y_0) \geq y_0 e^{-D^{\max}t - (V_{\text{inv}}^{\max} - V_{\text{inv}}^{\min}) \int_0^t (d+\Phi(z^*(\theta_\tau\omega))) d\tau} \\ + (D_1 s_1^{\text{in}} + D_2 s_2^{\text{in}}) \int_0^t e^{-D^{\max}(t-s) - (V_{\text{inv}}^{\max} - V_{\text{inv}}^{\min}) \int_s^t (d+\Phi(z^*(\theta_\tau\omega))) d\tau} ds \quad (4.6)$$

for every $t \geq 0$, $\omega \in \Omega$ and $y_0 > 0$.

Recall now that $V_{\text{inv}}^{\max} - V_{\text{inv}}^{\min} > 0$ then $-(V_{\text{inv}}^{\max} - V_{\text{inv}}^{\min}) < 0$. Moreover, thanks to (2.4), $d + \Phi(z^*(\theta_t\omega)) \leq \bar{d}$ for every $t \in \mathbb{R}$. Then, from (4.6), it is not difficult to check that

$$y(t; 0, \omega, y_0) \geq y_0 e^{-(D^{\max}+\bar{d}(V_{\text{inv}}^{\max}-V_{\text{inv}}^{\min}))t} + \frac{D_1 s_1^{\text{in}} + D_2 s_2^{\text{in}}}{D^{\max} + \bar{d}(V_{\text{inv}}^{\max}-V_{\text{inv}}^{\min})} \left[1 - e^{-(D^{\max}+\bar{d}(V_{\text{inv}}^{\max}-V_{\text{inv}}^{\min}))t}\right] \quad (4.7)$$

for every $t \geq 0$, $\omega \in \Omega$ and $y_0 > 0$.

Then, thanks to (4.7), for every given $\varepsilon > 0$, $\omega \in \Omega$ and $y_0 > 0$, there exists some time $T_F(\varepsilon, \omega) > 0$ such that

$$y(t; 0, \omega, y_0) \geq \frac{D_1 s_1^{\text{in}} + D_2 s_2^{\text{in}}}{D^{\max} + \bar{d}(V_{\text{inv}}^{\max} - V_{\text{inv}}^{\min})} - \varepsilon$$

for every $t \geq T_F(\varepsilon, \omega)$.

Hence, the deterministic set $B_\varepsilon$, given as in (4.2), is a compact absorbing set for the solutions of system (1.8)-(1.11) as long as (4.1) holds true. $\square$

**Remark 4.1.** We remark that the absorbing set $B_\varepsilon$ in (4.2) is deterministic, i.e., it does not depend on the realization of the noise, which is a nice property.

**Remark 4.2.** Let us highlight that the absorbing set $B_\varepsilon$ in (4.2) is well defined, i.e.,

$$y_l^* = \frac{D_1 s_1^{\text{in}} + D_2 s_2^{\text{in}}}{D^{\max} + d^{\max} - d^{\min}} \leq \frac{D_1 s_1^{\text{in}} + D_2 s_2^{\text{in}}}{D^{\min} - d^{\max} + d^{\min}} = y_u^*. \quad (4.8)$$

Indeed, notice that $D^{\max} + 2\bar{d}V_{\text{inv}}^{\max} \geq D^{\min} + 2\bar{d}V_{\text{inv}}^{\min}$, i.e., $D^{\max} + \bar{d}V_{\text{inv}}^{\max} + \bar{d}V_{\text{inv}}^{\min} \geq D^{\min} + \bar{d}V_{\text{inv}}^{\min} + \bar{d}V_{\text{inv}}^{\max}$, whence $D^{\max} + \bar{d}V_{\text{inv}}^{\max} - \bar{d}V_{\text{inv}}^{\min} \geq D^{\min} + \bar{d}V_{\text{inv}}^{\min} - \bar{d}V_{\text{inv}}^{\max}$ and then $D^{\max} + \bar{d}(V_{\text{inv}}^{\max} - V_{\text{inv}}^{\min}) \geq D^{\min} + \bar{d}(V_{\text{inv}}^{\min} - V_{\text{inv}}^{\max})$, whence (4.8) follows easily.

Thanks to Theorem 4.1, we have the following corollary ensuring the existence of a deterministic compact attracting set for the solutions of the random chemostat
model (1.8)-(1.11).

**Corollary 4.1.** As long as (4.1) holds true, where $y_l^*$ and $y_u^*$ are defined as in (4.3), the compact deterministic set

$$ \mathcal{A} := \{(s_1, x_1, s_2, x_2) \in \mathcal{X} : y_l^* \le s_1 + x_1 + s_2 + x_2 \le y_u^*\} \quad (4.9) $$

is an attracting set for the solutions of system (1.8)-(1.11) in forward sense, i.e.,

$$ \lim_{t \to +\infty} \sup_{u_0 \in F} \inf_{b_0 \in \mathcal{A}} |u(t; 0, \omega, u_0) - b_0|_{\mathcal{X}} = 0. $$

The proof of this corollary is straightforward and follows from Theorem 4.1.

**5. Internal structure of the attracting set: extinction and persistence of the species**

In this section we study the internal structure of the attracting set (4.9) in order to obtain more detailed information concerning the random dynamics of system (1.8)-(1.11). In this way, we are able to provide conditions under the extinction of the species cannot be avoided and conditions which ensure the persistence of the species.

We present first the case of the extinction of the species.

**Theorem 5.1.** Provided

$$ \bar{\mu}_0 + d(V_{inv}^{\max} - V_{inv}^{\min}) < D^{\min}, \quad (5.1) $$

the following limit holds true for every $\omega \in \Omega$ and every initial value $u_0 \in F$

$$ \lim_{t \to +\infty} (x_1(t; 0, \omega, u_0) + x_2(t; 0, \omega, u_0)) = 0, $$

which means the extinction of the species in both culture vessels.

**Proof.** Define a new state variable $x := x_1 + x_2$ describing the total concentration of species, which satisfies the following differential inequalities

$$ \begin{aligned} \frac{dx}{dt} &= -D_1 x_1 + \mu(s_1)x_1 + \frac{d + \Phi(z^*(\theta_t \omega))}{V_1} x_2 - \frac{d + \Phi(z^*(\theta_t \omega))}{V_1} x_1 \\ &\quad - D_2 x_2 + \mu(s_2)x_2 + \frac{d + \Phi(z^*(\theta_t \omega))}{V_2} x_1 - \frac{d + \Phi(z^*(\theta_t \omega))}{V_2} x_2 \\ &\le -D^{\min}(x_1 + x_2) + \bar{\mu}_0(x_1 + x_2) + V_{inv}^{\max}(d + \Phi(z^*(\theta_t \omega)))(x_1 + x_2) \\ &\quad - V_{inv}^{\min}(d + \Phi(z^*(\theta_t \omega)))(x_1 + x_2) \\ &= (\bar{\mu}_0 + d(V_{inv}^{\max} - V_{inv}^{\min}) - D^{\min} + (V_{inv}^{\max} - V_{inv}^{\min})\Phi(z^*(\theta_t \omega)))x, \quad (5.2) \end{aligned} $$

where $V_{inv}^{\min}$, $V_{inv}^{\max}$ and $D^{\min}$ are defined in Table 1.
Recall that we used (a) $D_i \ge D^{\min}$ then $-D_i \le -D^{\min}$ for $i = 1, 2$, (b) $\frac{1}{V_i} \le V_{\text{inv}}^{\max}$ for $i = 1, 2$, (c) $\frac{1}{V_i} \ge V_{\text{inv}}^{\min}$ then $-\frac{1}{V_i} \le -V_{\text{inv}}^{\min}$ for $i = 1, 2$ and (d) the ergodic property in Theorem 2.5.

Now, we multiply by

$$e^{-(\bar{\mu}_0+d(V_{\text{inv}}^{\max}-V_{\text{inv}}^{\min})-D^{\min})t-(V_{\text{inv}}^{\max}-V_{\text{inv}}^{\min})\int_0^t \Phi(z^*(\theta_\tau\omega))d\tau}$$

and then we integrate over $[0, t]$. Then, we obtain

$$x(t; 0, \omega, x_0) \le x_0 e^{(\bar{\mu}_0+d(V_{\text{inv}}^{\max}-V_{\text{inv}}^{\min})-D^{\min})t+(V_{\text{inv}}^{\max}-V_{\text{inv}}^{\min})\int_0^t \Phi(z^*(\theta_\tau\omega))d\tau} \quad (5.3)$$

for every $\omega \in \Omega$, $x_0 := x_{10} + x_{20} \ge 0$ and $t \ge 0$.

From (5.3), as long as (5.1) holds true, we have that

$$\lim_{t \to +\infty} (x_1(t; 0, \omega, u_0) + x_2(t; 0, \omega, u_0)) = 0. \quad \square$$

**Proposition 5.1.** Condition (5.1) can be replaced by

$$\mu(y_u^*) + d(V_{\text{inv}}^{\max} - V_{\text{inv}}^{\min}) < D^{\min}, \quad (5.4)$$

which is more refined than (5.1), where $y_u^*$ is given as in (4.3) and $V_{\text{inv}}^{\max}$, $V_{\text{inv}}^{\min}$ and $D^{\min}$ are defined in Table 1.

**Proof.** Recall that, for every time large enough, from Theorem (4.1) we have

$$s_i \le s_1 + x_1 + s_2 + x_2 \le y_u^* + \varepsilon$$

for $i \in \{1, 2\}$. Then, as $\mu$ is an increasing function, we obtain

$$\mu(s_i) \le \mu(y_u^* + \varepsilon)$$

for $i \in \{1, 2\}$. Hence, it is possible to redo the proof of Theorem 5.1. To be more precise, (5.2) becomes

$$\frac{dx}{dt} \le (\mu(y_u^* + \varepsilon) + d(V_{\text{inv}}^{\max} - V_{\text{inv}}^{\min}) - D^{\min} + (V_{\text{inv}}^{\max} - V_{\text{inv}}^{\min})\Phi(z^*(\theta_t\omega))) x,$$

whose integration gives

$$x(t; 0, \omega, x_0) \le x_0 e^{(\mu(y_u^*+\varepsilon)+d(V_{\text{inv}}^{\max}-V_{\text{inv}}^{\min})-D^{\min})t+(V_{\text{inv}}^{\max}-V_{\text{inv}}^{\min})\int_0^t \Phi(z^*(\theta_\tau\omega))d\tau} \quad (5.5)$$

for every $\omega \in \Omega$, $x_0 > 0$ and $t$ large enough.

From (5.5), as long as (5.4) holds true, it is possible to consider $\varepsilon > 0$ small enough such that

$$\mu(y_u^* + \varepsilon) + d(V_{\text{inv}}^{\max} - V_{\text{inv}}^{\min}) < D^{\min}$$

and then we also obtain

$$\lim_{t \to +\infty} (x_1(t; 0, \omega, u_0) + x_2(t; 0, \omega, u_0)) = 0. \quad \square$$
Now, we focus on providing conditions under which we can ensure the persistence of the species in both culture vessels in the random chemostat model (1.8)-(1.11).

**Theorem 5.2.** As long as

$$ \bar{\mu}_0 \frac{D_i s_i^{in}}{D_i + \bar{d}_i} - (D_i + \bar{d}_i)(k + y_u^*) > 0 \quad (5.6) $$

holds true for some $i \in \{1, 2\}$, where $\bar{d}_i := \bar{d}/V_i$, we have that

$$ \liminf_{t \to +\infty} x_i(t; 0, \omega, u_0) \ge \bar{x}_i > 0 \quad (5.7) $$

for any $\omega \in \Omega$ and every $u_0 \in \mathcal{X}$, where

$$ \bar{x}_i := \frac{D_i s_i^{in}}{D_i + \bar{d}_i} - \frac{1}{\bar{\mu}_0} (D_i + \bar{d}_i)(k + y_u^*), \quad (5.8) $$

which means the persistence of the species in the corresponding tank $i \in \{1, 2\}$.

**Proof.** We first study the evolution of the sum of the concentrations of substrate and the microbial biomass in each culture vessel.

Let us start with the first one. To this end, we consider (1.8) and (1.9), whence we have the following expression

$$ \begin{aligned} \frac{d(s_1 + x_1)}{dt} &= D_1(s_1^{in} - s_1 - x_1) + \left(d_1 + \frac{1}{V_1}\Phi(z^*(\theta_t\omega))\right)(s_2 + x_2) \\ &\quad - \left(d_1 + \frac{1}{V_1}\Phi(z^*(\theta_t\omega))\right)(s_1 + x_1) \\ &\ge D_1 s_1^{in} - \left(D_1 + d_1 + \frac{1}{V_1}\Phi(z^*(\theta_t\omega))\right)(s_1 + x_1), \end{aligned} \quad (5.9) $$

where we used that $(d_1 + \frac{1}{V_1}\Phi(z^*(\theta_t\omega)))(s_2 + x_2) \ge 0$.

By integrating now (5.9), we obtain

$$ \begin{aligned} s_1(t; 0, \omega, u_0) + x_1(t; 0, \omega, u_0) &\ge (s_{10} + x_{10})e^{-(D_1+d_1)t-\frac{1}{V_1}\int_0^t \Phi(z^*(\theta_\tau\omega))d\tau} \\ &\quad + D_1 s_1^{in} \int_0^t e^{-(D_1+d_1)(t-r)-\frac{1}{V_1}\int_r^t \Phi(z^*(\theta_\tau\omega))d\tau} dr \\ &\ge (s_{10} + x_{10})e^{-(D_1+d_1)t-\frac{1}{V_1}\int_0^t \Phi(z^*(\theta_\tau\omega))d\tau} \\ &\quad + D_1 s_1^{in} \int_0^t e^{-(D_1+\bar{d}_1)r} dr \\ &= (s_{10} + x_{10})e^{-(D_1+d_1)t-\frac{1}{V_1}\int_0^t \Phi(z^*(\theta_\tau\omega))d\tau} \\ &\quad + \frac{D_1 s_1^{in}}{D_1 + \bar{d}_1} [1 - e^{-(D_1+\bar{d}_1)t}] \end{aligned} \quad (5.10) $$
for every $t \ge 0$, $\omega \in \Omega$ and $u_0 \in F$.

Then, from (5.10), for any $\varepsilon > 0$, $\omega \in \Omega$ and every initial value $u_0 \in F$, there exists some time $T_F(\omega, \varepsilon) > 0$ such that

$$s_1(t; 0, \omega, u_0) + x_1(t; 0, \omega, u_0) \ge \frac{D_1 s_1^{in}}{D_1 + \bar{d}_1} - \varepsilon \quad (5.11)$$

for every $t \ge T_F(\omega, \varepsilon)$.

Repeating the same steps with (1.10) and (1.11) we conclude that, for every $\varepsilon > 0$, $\omega \in \Omega$ and every initial value $u_0 \in F$, there exists some time $T_F(\omega, \varepsilon) > 0$ such that

$$s_2(t; 0, \omega, u_0) + x_2(t; 0, \omega, u_0) \ge \frac{D_2 s_2^{in}}{D_2 + \bar{d}_2} - \varepsilon \quad (5.12)$$

for every $t \ge T_F(\omega, \varepsilon)$.

Once reached this point, we can investigate the equation describing the dynamics of the species in the first culture vessel $x_1$. To this end, from (1.9), for every time large enough we have

$$\frac{dx_1}{dt} = -D_1x_1 + \mu(s_1)x_1 - \frac{d + \Phi(z^*(\theta_t\omega))}{V_1}x_1 \quad (5.13)$$

$$\ge -D_1x_1 + \mu(s_1)x_1 - \bar{d}_1x_1$$

$$= -(D_1 + \bar{d}_1)x_1 + \frac{\bar{\mu}_0 s_1}{k + s_1}x_1$$

$$\ge \left[ -(D_1 + \bar{d}_1) + \frac{\bar{\mu}_0 \left( \frac{D_1 s_1^{in}}{D_1 + \bar{d}_1} - \varepsilon - x_1 \right)}{k + y_u^* + \varepsilon} \right] x_1, \quad (5.14)$$

where we used (5.11).

It is easy to check that the sign of the last term in (5.14) (depending on $x_1$) coincides with the sign of the parabolic function given by

$$f(x) = \alpha_2 x^2 + \alpha_1(\varepsilon)x,$$

where $\alpha_2$ and $\alpha_1(\varepsilon)$ are given by

$$\alpha_2 := -\bar{\mu}_0,$$

$$\alpha_1(\varepsilon) := \bar{\mu}_0 \frac{D_1 s_1^{in}}{D_1 + \bar{d}_1} - (D_1 + \bar{d}_1)(k + y_u^*) - \varepsilon(D_1 + \bar{d}_1 + \bar{\mu}_0). \quad (5.15)$$

Define

$$x_1^*(\varepsilon) := \frac{D_1 s_1^{in}}{D_1 + \bar{d}_1} - \frac{1}{\bar{\mu}_0}(D_1 + \bar{d}_1)(k + y_u^*) - \frac{1}{\bar{\mu}_0}(\bar{\mu}_0 + D_1 + \bar{d}_1)\varepsilon. \quad (5.16)$$
On the one hand, $\alpha_2 < 0$. On the other hand, thanks to (5.6), one can find $\varepsilon > 0$ small enough such that $\alpha_1(\varepsilon) > 0$. Then, $f(x) > 0$ for $x \in (0, x_1^*(\varepsilon))$, $f(x) < 0$ for $x > x_1^*(\varepsilon)$ and $f(0) = f(x_1^*(\varepsilon)) = 0$, where $\varepsilon > 0$ small enough.

As a result, we deduce that

$$ \left. \frac{dx_1}{dt} \right|_{x_1=\tilde{x}_1} > 0 $$

for every $\tilde{x}_1 \in (0, x_1^*(\varepsilon))$, where $\varepsilon > 0$ small enough.

Thus, by taking limit when $\varepsilon$ goes to zero, we have that

$$ \liminf_{t \to +\infty} x_1(t; \omega, u_0) \ge \bar{x}_1 > 0 \qquad (5.17) $$

for every $\omega \in \Omega$ and any $u_0 \in \mathcal{X}$, where

$$ \bar{x}_1 := \frac{D_1 s_1^{in}}{D_1 + \bar{d}_1} - \frac{1}{\bar{\mu}_0} (D_1 + \bar{d}_1)(k + y_u^*). \qquad (5.18) $$

From (5.17), $x_2$ cannot goes to zero and we have an analogous result when condition (5.6) is fulfilled for $i = 2$ instead of $i = 1$. More precisely, we obtain that

$$ \liminf_{t \to +\infty} x_2(t; \omega, u_0) \ge \bar{x}_2 > 0 $$

for any $\omega \in \Omega$ and every $u_0 \in \mathcal{X}$, where

$$ \bar{x}_2 := \frac{D_2 s_2^{in}}{D_2 + \bar{d}_2} - \frac{1}{\bar{\mu}_0} (D_2 + \bar{d}_2)(k + y_u^*). \qquad (5.19) $$

**Remark 5.1.** It is worth mentioning that, if a condition that ensures persistence in one tank is fulfilled, one has immediately persistence in the other tank since, as $\bar{d} > 0$, the washout in the other tank is repulsive.

Thanks to the previous theorem, we provide practical conditions under which the persistence of species in both culture vessels can be ensured, which is often a goal pursued by practitioners.

**Remark 5.2.** It is not difficult, but laborious, to prove that the lower positive bounds $\bar{x}_1$ and $\bar{x}_2$ found in Theorem 5.2 are consistent with the previous theoretical results, that is,

$$ \bar{x}_1 + \bar{x}_2 \le y_u^*. \qquad (5.20) $$

To prove that, it is enough to write the definition of $\bar{x}_1$ in (5.11), $\bar{x}_2$ in (5.12) and $y_u^*$ in (5.20) and expand every term. After some lengthy calculations, the reader can check that this inequality is fulfilled.
## 6. Numerical simulations

In this section we present some numerical simulations to illustrate the theoretical results proved in the previous sections.

In Fig. 3 we depict a realization (in blue color) of the perturbed diffusion parameter $d$, where the red dashed line represents its deterministic value. As we explained in the preliminary section, it remains positive and bounded in a certain interval for every time. In this case we consider $d = 0.5$, $\alpha = 0.4$ and we take $\beta = 1$ and $\nu = 1$.

<figure><img src="image_4.png" /><figcaption>Fig. 3: Realization of the perturbed diffusion parameter. Values of the parameters: $d = 0.5$, $\alpha = 0.4$, $\beta = 1$ and $\nu = 1$.</figcaption></figure>

Now, we present in Fig. 4 the case in which the species in both culture vessels become extinct. To this end, we consider $s_1^{in} = 4$, $s_2^{in} = 1$, $k = 4$, $\bar{\mu}_0 = 0.1$, $Q_1 = 1$, $V_1 = 1.1$, $Q_2 = 1.4$, $V_2 = 1.4$, $\alpha = 0.4$, $d = 1$, we set $\beta = 1$ and $\nu = 1$ and we take $x_{10} = 5$, $s_{10} = 4$, $x_{20} = 2$ and $s_{20} = 5$ as initial values.
<figure><img src="image_5.png" /><figcaption>Fig. 4: **Extinction of the species in both culture vessels.** Values of the parameters: $s_1^{in} = 4$, $s_2^{in} = 1$, $k = 4$, $\bar{\mu}_0 = 0.1$, $Q_1 = 1$, $V_1 = 1.1$, $Q_2 = 1.4$, $V_2 = 1.4$, $\alpha = 0.4$, $d = 1$, $\beta = 1$, $\nu = 1$, $x_{10} = 5$, $s_{10} = 4$, $x_{20} = 2$ and $s_{20} = 5$</figcaption></figure>

It is easy to check that condition (4.1) to ensure the existence of the deterministic compact absorbing set $B_\varepsilon$, given by (4.2), holds true. In addition, since (5.1) is also verified, then we observe the extinction of the species in both culture vessels, as proved in Theorem 5.1.

In Fig. 5 we plot the case in which we have persistence of the species in both culture vessels. We present two different panels. The one in the top illustrates the dynamics of the species in the first culture vessel and the one in the bottom focuses on the dynamics of the species in the second culture vessel. In each one we simulate different realizations. In addition, in both panels we include another little one with a zoom of the dynamics around the final time.

In this case we consider $s_1^{in} = 0.5$, $s_2^{in} = 0.4$, $k = 0.04$, $\bar{\mu}_0 = 10$, $Q_1 = 0.3$, $V_1 = 1.7$, $Q_2 = 0.5$, $V_2 = 1.9$, $\alpha = 0.4$, $d = 0.5$, we take $\beta = 1$, $\nu = 1$ and we choose $x_{10} = 4$, $s_{10} = 1$, $x_{20} = 5$ and $s_{20} = 3$ as initial values.
<figure><img src="image_6.png" /><figcaption>Fig. 5: Persistence of the species in both culture vessels. Values of the parameters: $s_1^{in} = 0.5$, $s_2^{in} = 0.4$, $k = 0.04$, $\bar{\mu}_0 = 10$, $Q_1 = 0.3$, $V_1 = 1.7$, $Q_2 = 0.5$, $V_2 = 1.9$, $\alpha = 0.4$, $d = 0.5$, $\beta = 1$, $\nu = 1$, $x_{10} = 4$, $s_{10} = 1$, $x_{20} = 5$ and $s_{20} = 3$</figcaption></figure>

Notice that condition (4.1) is satisfied in this case as well, whence we have existence of a deterministic compact absorbing set. In addition, as (5.6) in Theorem 5.2 holds true for $i = 1, 2$, then we can guarantee the persistence of the species in both culture vessels, as observed in the figure.

Finally, in Fig. 6 we depict again some numerical simulations in the case of persistence of the species in both culture vessels. Nevertheless, in this case we set $s_1^{in} = 3.5$, $s_2^{in} = 4.4$, $k = 1$, $\bar{\mu}_0 = 5$, $Q_1 = 0.7$, $V_1 = 4.7$, $Q_2 = 0.4$, $V_2 = 4.9$, $\alpha = 0.5$, $d = 0.8$, we take $\beta = 1$, $\nu = 1$ and we choose $x_{10} = 2$, $s_{10} = 8$, $x_{20} = 1$ and $s_{20} = 10$ as initial values.
<figure><img src="image_7.png" /><figcaption>Fig. 6: **Persistence of the species in both culture vessels.** Values of the parameters: $s_1^{in} = 3.5$, $s_2^{in} = 4.4$, $k = 1$, $\bar{\mu}_0 = 5$, $Q_1 = 0.7$, $V_1 = 4.7$, $Q_2 = 0.4$, $V_2 = 4.9$, $\alpha = 0.5$, $d = 0.8$, $\beta = 1$, $\nu = 1$, $x_{10} = 2$, $s_{10} = 8$, $x_{20} = 1$ and $s_{20} = 10$</figcaption></figure>

In this case, for the sake of clarity, we present in Fig. 7 a zoom around the final time of the previous figure.

<figure><img src="image_8.png" /><figcaption>Fig. 7: **Persistence of the species in both culture vessels.** Values of the parameters: $s_1^{in} = 03.5$, $s_2^{in} = 4.4$, $k = 1$, $\bar{\mu}_0 = 5$, $Q_1 = 0.7$, $V_1 = 4.7$, $Q_2 = 0.4$, $V_2 = 4.9$, $\alpha = 0.5$, $d = 0.8$, $\beta = 1$, $\nu = 1$, $x_{10} = 2$, $s_{10} = 8$, $x_{20} = 1$ and $s_{20} = 10$</figcaption></figure>
We notice that conditions (4.1) holds true, whence we have existence of the deterministic compact absorbing set. Moreover, since (5.6) is also satisfied for $i = 1, 2$, from Theorem 5.2, we have that the species in both culture vessels persist.

It worth mentioning an interesting fact concerning the numerical simulations in Figs. 5, 6 and 7. Notice that, even though the system in study is random (see, for instance, Fig. 3 to observe how the random diffusion parameter behaves), the dynamics of the concentration of the species in both culture vessel seems to be almost deterministic even though the noise amplitude is not so small. Let us underline that, differently to disturbances on the input term $Q$ as considered in [9, 10, 11, 12], the noise acts on an *internal* term $d$ of the dynamics.

## 7. Conclusions

We recall that we consider two chemostat models connected by Fickian diffusion in this paper, where the diffusion parameter is perturbed by a bounded noise which essentially is a function of the well-known Ornstein-Uhlenbeck process.

As a first step, we proved the existence and uniqueness of non-negative global solution of the resulting random system in Theorem 3.1 and we also proved the existence of compact absorbing and attracting sets in Theorem 4.1 which, in addition, were deterministic. We notice that this is a nice property since such a sets do not depend on the realizations of the noise.

After that, we are interested in studying the internal structure of the deterministic compact attracting set in order to obtain more detailed information about the long-time dynamics of the random system in study.

To this end, Theorem 5.1 proves that the extinction of the species in both culture vessels cannot be avoided as long as condition 5.1, which is given by

$$ \bar{\mu}_0 + d(V_{\text{inv}}^{\text{max}} - V_{\text{inv}}^{\text{min}}) < D^{\text{min}}, $$

holds true.

Notice that condition 5.1 is quite logical since it means that the minimum of the dilution rates $D^{\text{min}}$ is greater than the growth rate coefficient of the consumer species $\bar{\mu}_0$ and the diffusion parameter $d$.

Apart from that, we are also interested in investigating conditions to ensure the persistence of the species in both culture vessels, which is the main goal pursued by practitioners. To this end, we present Theorem 5.2 where we ensure the persistence of the species $i = 1, 2$ as long as condition (5.6), which is given by

$$ \bar{\mu}_0 \frac{D_i s_i^{in}}{D_i + \bar{d}_i} - (D_i + \bar{d}_i)(k + y_u^*) > 0, $$

holds true.
Notice that condition (5.6) can be verified when increasing the growth rate coeffi-
cient of the consumer species $\bar{\mu}_0$ and decreasing both the dilution rates $D_i, i = 1, 2$,
the half-saturation constant $k$ and the parameters $\bar{d}_i, i = 1, 2$. This is not surprising
since it means that the consumption of species should be large, in comparison with
the dilution rate and the diffusion parameter.

Finally, we provide in Theorem 5.2 positive lower bounds for the concentration of
the species in both culture vessels.

In this work we have considered noise on the diffusion term $d$ that can be consid-
ered as *internal* to the dynamics of the overall chemostat seen a *spatialized* system
composed of the two tanks. Future work might consider the effects of having noise
simultaneously on the input and the diffusion.

Acknowledgments

This work has been partially supported by Project PGC2018-096540-B-I00 from
FEDER and Ministerio Ministerio de Ciencia, Innovación 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).

References

1. L. Arnold, *Random Dynamical Systems*, Springer Berlin Heidelberg, 1998.

2. J. Barlow, F. de Noyelles, B. Peterson, J. Peterson and W. Schaffner, Continuous flow nutrient bioassays with natural phytoplankton populations, G. Glass (Editor): Bioassay Techniques and Environmental Chemistry, John Wiley & Sons Ltd., 1973.

3. T. Caraballo, R. Colucci, J. López-de-la-Cruz and A. Rapaport, A way to model stochastic perturbations in population dynamics models with bounded realizations, *Communications in Nonlinear Science and Numerical Simulation*, **77** (2019), 239–257.

4. T. Caraballo, M. J. Garrido-Atienza and J. López-de-la-Cruz, *Some Aspects Concerning the Dynamics of Stochastic Chemostats*, vol. 69, chapter 11, 227–246, Springer International Publishing, Cham, 2016.

5. T. Caraballo, M. J. Garrido-Atienza and J. López-de-la-Cruz, Dynamics of some stochastic chemostat models with multiplicative noise, *Communications on Pure and Applied Analysis*, **16** (2017), 1893–1914.

6. T. Caraballo, M. J. Garrido-Atienza, J. López-de-la-Cruz and A. Rapaport, Modeling and analysis of random and stochastic input flows in the chemostat model, *Discrete & Continuous Dynamical Systems - Series B*, **24** (2018), 3591–3614.

7. T. Caraballo and X. Han, *Applied Nonautonomous and Random Dynamical Systems, Applied Dynamical Systems*, Springer International Publishing, 2016.

8. T. Caraballo, P. E. Kloeden and B. Schmalfuss, Exponentially stable stationary solutions for stochastic evolution equations and their perturbation, *Applied Mathematics and Optimization*, **50** (2004), 183–207.

9. T. Caraballo and J. López-de-la-Cruz, Bounded random fluctuations on the input flow in chemostat models with wall growth and non-monotonic kinetics, Mathematics, **6** (2021), 4025–4052.
24 Tomás Caraballo, Javier López-de-la-Cruz and Alain Rapaport

10. T. Caraballo and J. López-de-la-Cruz, Survey on chemostat models with bounded random input flow, *Mathematical Modelling and Control*, **1** (2021), 52-78.

11. T. Caraballo, J. López-de-la-Cruz and A. Rapaport, Modeling bounded random fluctuations in biological systems: application to the chemostat model with two species, *IFAC-PapersOnLine*, **52** (2019), 187-192.

12. T. Caraballo, J. López-de-la-Cruz and A. Rapaport, Study of the chemostat model with non-monotonic growth under random disturbances on the removal rate, *Mathematical Biosciences and Engineering*, **17** (2021), 7480-7501.

13. I. F. Creed, D. M. McKnight, B. A. Pellerin, M. B. Green, B. A. Bergamaschi, G. R. Aiken, D. A. Burns, S. E. G. Findlay, J. B. Shanley, R. G. Striegl, B. T. Aulenbach, D. W. Clow, H. Laudon, B. L. McGlynn, K. J. McGuire, R. A. Smith and S. M. Stackpoole, The river as a chemostat: fresh perspectives on dissolved organic matter flowing down the river continuum, *Canadian Journal of Fisheries and Aquatic Sciences*, **72** (2015), 1272-1285.

14. M. Crespo and A. Rapaport, Analysis and optimization of the chemostat model with a lateral diffusive compartment, *Journal of Optimization Theory and Applications*, **185** (2020), 597-621.

15. G. D'Ans, P. Kokotovic and D. Gottlieb, A nonlinear regulator problem for a model of biological waste treatment, *IEEE Transactions on Automatic Control*, **16** (1971), 341-347.

16. S. Diehl, J. Zambrano and B. Carlsson, Steady-state analyses of activated sludge processes with plug-flow reactor, *Journal of Environmental Chemical Engineering*, **5** (2017), 795-809.

17. P. Doran, Design of mixing systems for plan cell suspensions in stireed reactors, *Biotechnology Progress*, **15** (1999), 319-335.

18. A. Dramé, A semilinear parabolic boundary-value problem in bioreactors theory, *Electronic Journal of Differential Equations*, **129** (2004), 1-13.

19. A. Dramé, C. Lobry, J. Harmand, A. Rapaport and F. Mazenc, Multiple stable equilibriums profiles in tubular bioreactors, *Mathematical and Computer Modelling*, **48** (2008), 1840-1853.

20. S. Fogler, *Elements of Chemical Reaction Engineering, 4<sup>rd</sup> edition*, Prentice Hall, New York, 2008.

21. R. Freter, *Human Intestinal Microflora in Health and Disease*, chapter Mechanisms that control the microflora in the large intestine, 33-54, Academic Press, New York, 1983.

22. R. Freter, An understanding of colonization of the large intestine requires mathematical analysis, *Microecology and Therapy*, **16** (1986), 147-155.

23. R. B. Grieves, W. O. Pipes, W. F. Milbury and R. K. Wood, Piston-flow reactor model for continuous industrial fermentations, *Journal of Applied Chemistry*, **14** (1964), 478-486.

24. I. Haidar, A. Rapaport and F. Gérard, Effects of spatial structure and diffusion on the performances of the chemostat, *Mathematical Biosciences and Engineering*, **8** (2011), 953-971.

25. J. Harmand, C. Lobry, A. Rapaport and T. Sari, *The Chemostat: Mathematical Theory of Micro-organisms Cultures*, Wiley, Chemical Engineering Series, John Wiley & Sons, Inc., 2017.

26. G. E. Hutchinson, *A treatise on limnology. Vol. II: Introduction to lake biology and the limnoplankton*, John Willery & Sons, 1967.

27. H. W. Jannasch, Steady state and the chemostat in ecology, *Limnology and Oceanography*, **19** (1974), 716-720.
28. J. Kalff and R. Knoechel, Phytoplankton and their dynamics in oligotrophic and eutrophic lakes, *Annual Review of Ecology and Systematics*, **9** (1978), 475-495.

29. C. M. Kung and B. Baltzis, The growth of pure and simple microbial competitors in a moving and distributed medium, *Mathematical Biosciences*, **111** (1992), 295-313.

30. J. W. M. La Rivière, Microbial ecology of liquid waste treatment, in *Advances in Microbial Ecology*, vol. 1, Springer US, 1977, 215-259.

31. O. Levenspiel, *Chemical reaction engineering, 3rd edition*, Wiley, New York, 1999.

32. R. Lovitt and J. Wimpenny, The gradostat: a tool for investigating microbial growth and interactions insolute gradients, *The Journal of General Microbiology*, **6**.

33. R. Lovitt and J. Wimpenny, The gradostat: a bidirectional compound chemostat and its applications in microbial research, *The Journal of General Microbiology*, **127**.

34. J. Monod, La technique de culture continue: Théorie et applications, *Annales de l'Institute Pasteur*, **79** (1950), 390-410.

35. H. Nie and J. Wu, Positive solutions of a competition model for two resources in the unstirred chemostat, *Journal of Mathematical Analysis and Applications*, **355** (2009), 231-242.

36. A. Novick and L. Szilard, Experiments with the chemostat on spontaneous mutations of bacteria, *Proceedings of the National Academy of Sciences*, **36** (1950), 708-719.

37. E. Rurangwa and M. C. J. Verdegem, Microorganisms in recirculating aquaculture systems and their management, *Reviews in Aquaculture*, **7** (2015), 117-130.

38. H. Smith, B. Tang and P. Waltman, Competition in a *n*-vessel gradostat, *Journal of Applied Mathematics*, **91** (1991), 1451-1471.

39. H. Smith and P. Waltman, The gradostat: A model of competition along a nutrient gradient, *Microbial Ecology*, **22** (1991), 207-226.

40. G. Stephanopoulos, R. Aris and A. Fredrickson, A stochastic analysis of the growth of competing microbial populations in a continuous biochemical reactor, *Mathematical Biosciences*, **45** (1979), 99-135.

41. G. Stephanopoulos and A. Fredrickson, Effect of inhomogeneities on the coexistence of competing microbial populations, *Biotechnology and Bioengineering*, **21** (1979), 1491-1498.

42. F. Stewart and B. Levin, The population biology of bacterial plasmids: A priori conditions for the existence of conjugationally transmitted factors, *Genetics*, **87** (1977), 209-28.

43. B. Tang, Mathematical investigations of growth of microorganisms in the gradostat, *Journal of Mathematical Biology*, **23**.

44. J. López-de-la-Cruz, Random and stochastic disturbances on the input flow in chemostat models with wall growth, *Stochastic Analysis and Applications*, **37** (2019), 668-698.

45. J. Zambrano, B. Carlsson and S. Diehl, Optimal steady state design of zone volumes of bioreactors with monod growth kinetics, *Biotechnology and Bioengineering*, **100** (2015), 59-66.