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

2007 International Conference in Honor of Claude Lobry

Practical coexistence in the chemostat
with arbitrarily close growth functions

Alain Rapaporta,\* — Denis Dochainb — Jérôme Harmandc

a UMR Analyse des systèmes et Biométrie et
Equipe-Projet INRA-INRIA "Modélisation et Ressources en Eau"
2 pl. Viala,
34060 Montpellier, FRANCE

b CESAME, Université Catholique de Louvain,
4-6 avenue G. Lemaître,
1348 Louvain-la-Neuve, BELGIUM
denis.dochain@uclouvain.be

c Laboratoire de Biotechnologie de l'Environnement et
Equipe-Projet INRA-INRIA "Modélisation et Ressources en Eau"
route des étangs,
11100 Narbonne, FRANCE
harmand@supagro.inra.fr

\* Corresponding author
rapaport@supagro.inra.fr

**RÉSUMÉ.** Nous montrons que la coexistence entre différentes espèces en compétition sur une même ressource peut durer sensiblement, lorsque leurs courbes de croissance sont arbitrairement proches. Le comportement transitoire est analysé en termes de dynamiques lente-rapide. Nous prouvons que des espèces non dominantes peuvent d'abord croître avant de décroître, en fonction de leurs proportions initiales.

**ABSTRACT.** We show that the coexistence of different species in competition for a common resource may be substantially long when their growth functions are arbitrarily closed. The transient behavior is analyzed in terms of slow-fast dynamics. We prove that non-dominant species can first increase before decreasing, depending on their initial proportions.

**MOTS-CLÉS** : chémostat, compétition, persistance, dyamiques lente-rapide

**KEYWORDS** : chemostat, competition, persistence, slow-fast dynamics
# 1. Introduction

The Competitive Exclusion Principle is well known in microbial ecology. It expresses the fact that when several microbial species grow on the same substrate in a chemostat, generically at most one species will eventually survive [5, 1, 12, 2]. This is an asymptotic property and the transient behavior of such dynamics has not yet been investigated, to our knowledge.

In the present paper, we study transient properties of multiple species growing on the same substrate, depending on the initial species distribution. We exhibit a concept of *practical coexistence*, in the sense that even if the best species will finally be the winner, the transient before the other species disappear may eventually be substantially long. It appears that the different species may coexist for a long period of time before the competitive exclusion practically applies. More precisely, some of species may be first increasing (before finally decreasing) depending on the initial conditions. Our analysis is based on a slow-fast characterization of the system dynamics, and provides an evaluation of the bound from below of the time at which each species stops increasing and therefore starts decreasing.

In Section 2, we give some preliminaries about the dynamics and recall the Competitive Exclusion Principle. Section 3 presents a slow-fast description of the dynamics. Then, the reduced order model is analyzed in Section 4. The results are then illustrated with numerical simulations in Section 5.

# 2. Preliminaries

Let us consider the chemostat model with one limited resource and several species

$$ egin{cases} \dot{x}_i = \mu_i(s)x_i - Dx_i, & (i = 1 \cdots n) \\ \dot{s} = -\sum_{i=1}^{n} \frac{\mu_i(s)}{y_i}x_i + D(S_{in} - s) \end{cases} \quad (1) $$

The growth functions $\mu_i(\cdot)$ are assumed to be $C^1$ non-negative functions such that $\mu_i(0) = 0$.

Without any loss of generality, we shall assume in the following that all yield factors $y_i$ have been taken equal to one (one can easily check that this amounts to replace $x_i$ by $x_i/y_i$ or to change the unit measuring each stock $x_i$). Let us first recall the following Lemma.

**Lemma 2.1** *The domain*

$$ \mathcal{D} = \left\{ (x, s) \in \mathbb{R}_{+}^{n+1} \mid \sum_{i=1}^{n} x_i + s \leq S_{in} \right\} $$

is invariant and attractive by the dynamics (1) in the positive cone $\mathbb{R}_{+}^{n+1}$.

*Proof.* When $x_i = 0$, one has $\dot{x}_i = 0$. Consequently, the trajectories cannot cross the axes $x_i = 0$.
When $s = 0$, one has $\dot{s} = DS_{in} > 0$. The trajectories cannot approach the axis $s = 0$.

From these two facts, one concludes that the cone $\mathbb{R}_+^{n+1}$ is an invariant domain. Consider now the variable

$$z = \sum_{i=1}^{n} x_i + s$$

which is solution of the ordinary differential equation $\dot{z} = D(S_{in} - z)$. One immediately concludes that the domain $\mathcal{D} = \mathbb{R}_+^{n+1} \cap \{z \le S_{in}\}$ is invariant and attractive. $\square$

Let us now introduce the following assumption.

**Assumption A0.** Functions $\mu_i(\cdot)$ are increasing for any $i = 1 \cdots n$.

Under Assumption A0, it is usual to define the break-even concentrations :

$$\lambda_i(D) = \begin{cases} s_i & \text{such that } \mu_i(s_i) = D, \\ +\infty & \text{if } \mu_i(s) < D \text{ for any } s \ge 0, \end{cases} \quad (2)$$

for each $i = 1 \cdots n$. Let us recall the Competitive Exclusion Principle (CEP), (see for instance Theorem 3.2 in [11]), for which the following assumption is required.

**Assumption A1.** There exists an unique $i^* \in \{1 \cdots n\}$ such that

$$\lambda_{i^*}(D) = \min_{i=1 \cdots n} \lambda_i(D) < +\infty.$$

**Proposition 2.2 (CEP)** Under Assumptions A0 and A1, any trajectory of (1) with initial condition in $\mathcal{D}$ such that $x_{i^*}(0) > 0$ fulfills the following properties :

- the substrate concentration $s(\cdot)$ converges asymptotically toward the steady state value :

$$s^* = \min(\lambda_{i^*}(D), S_{in}),$$

- the species concentration $x_{i^*}(\cdot)$ converges asymptotically toward $S_{in} - s^*$,

- any species concentration $x_i(\cdot)$ with $i \ne i^*$ converges asymptotically toward zero.

**Corollary 2.3** When $s^* < S_{in}$, the convergence given by Proposition 2.2 is exponential.

*Proof.* One can easily check the $n + 1$ eigenvalues of the Jacobian matrix at the non-null equilibrium are $-D < 0$, $-\mu_{i^*}'(s^*)(S_{in} - s^*) < 0$ and $\mu_i(s^*) - D < 0$ for any $i \ne i^*$. $\square$

The CEP provides information about the asymptotic behavior of solutions of (1). In the present work, we focus on transient stages of the trajectories of system (1), when functions $\mu_i(\cdot)$ are “close” to each other.

In the following, we shall assume that A0 and A1 are fulfilled with $s^* < S_{in}$.
### 3. A slow-fast characterization

We consider the average growth function, defined as follows.

$$ \bar{\mu}(s) = \frac{1}{n} \sum_{i=1}^{n} \mu_i(s), $$

and define the number

$$ \epsilon = \max_{i \le n} \max_{s \in [0, S_{in}]} |\mu_i(s) - \bar{\mu}(s)|. \quad (3) $$

Notice that $\epsilon$ is positive under Assumption A1 (functions $\mu_i(\cdot)$ cannot coincide on the whole interval $[0, S_{in}]$). Then, consider the $C^1$ functions

$$ \nu_i(s) = \frac{\mu_i(s) - \bar{\mu}(s)}{\epsilon}, \quad (i = 1 \cdots n). \quad (4) $$

Then the growth functions $\mu_i(\cdot)$ can be expressed as follows.

$$ \mu_i(s) = \bar{\mu}(s) + \epsilon\nu_i(s), \quad (i = 1 \cdots n). $$

Let us now consider the total biomass $b$ and the proportions $p_i$, defined as follows.

$$ b = \sum_{i=1}^{n} x_i, \quad p_i = \frac{x_i}{b}. $$

Then the dynamics of the variables $b$, $s$ and $p_i$ are given by the following equations :

$$ \left\{ \begin{aligned} \dot{b} &= \bar{\mu}(s)b - Db + \epsilon \left( \sum_{i=1}^{n} \nu_i(s)p_i \right) b \\ \dot{s} &= -\bar{\mu}(s)b + D(S_{in} - s) - \epsilon \left( \sum_{i=1}^{n} \nu_i(s)p_i \right) b \\ \dot{p}_i &= \epsilon \left( \sum_{j=1}^{n} (\nu_i(s) - \nu_j(s))p_j \right) p_i \quad (i = 1 \cdots n) \end{aligned} \right. \quad (5) $$

Let us consider the change of time variable $\tau = \epsilon t$. System (5) can then be equivalently written as follows :

$$ \left\{ \begin{aligned} \epsilon \frac{db}{d\tau} &= \bar{\mu}(s)b - Db + \epsilon \left( \sum_{i=1}^{n} \nu_i(s)p_i \right) b \\ \epsilon \frac{ds}{d\tau} &= -\bar{\mu}(s)b + D(S_{in} - s) - \epsilon \left( \sum_{i=1}^{n} \nu_i(s)p_i \right) b \\ \frac{dp_i}{d\tau} &= \left( \sum_{j=1}^{n} (\nu_i(s) - \nu_j(s))p_j \right) p_i \quad (i = 1 \cdots n) \end{aligned} \right. \quad (6) $$

When $\epsilon$ is small i.e. the growth functions $\mu_i(\cdot)$ are all close to the average $\bar{\mu}(\cdot)$, system (6) is in the “slow-fast” form, for which the vector :

$$ \xi = \begin{pmatrix} b \\ s \end{pmatrix} $$
corresponds to the “fast” variables, and the “boundary-layer” dynamics is given by the system :

$$
\left\{
\begin{array}{l}
\dot{\tilde{b}} = \bar{\mu}(\tilde{s})\tilde{b} - D\tilde{b} \\
\dot{\tilde{s}} = -\bar{\mu}(\tilde{s})\tilde{b} + D(S_{in} - \tilde{s})
\end{array}
\right.
\qquad (7)
$$

Notice that system (7) is nothing else that the simple chemostat model with a single (fictitious) species of growth $\bar{\mu}(\cdot)$. Denote $\bar{\lambda}(\cdot)$ the break-even concentration associated to the function $\bar{\mu}(\cdot)$ and assume the following hypothesis.

**Assumption A2.** $\bar{\lambda}(D) < S_{in}$.

Then, under this last assumption, dynamics (7) admits the positive equilibrium

$$
\bar{E} = \left( \begin{array}{c} S_{in} - \bar{\lambda}(D) \\ \bar{\lambda}(D) \end{array} \right)
$$

which is globally exponentially stable on $\mathbb{R}_+ \setminus \{0\} \times \mathbb{R}_+$ (Corollary 2.3).

We show now that fixing an arbitrary small neighborhood $\mathcal{V}$ of $\bar{E}$ and an arbitrary small positive number $\tau$, there exists $\bar{\epsilon} > 0$ such that, for any $\epsilon < \bar{\epsilon}$ the state vector $\xi(\cdot)$ enters and remains in $\mathcal{V}$ within the time $\tau$.

**Proposition 3.1** Assume that A2 is fulfilled. For any initial condition in $\mathcal{D}$ with $b(0) > 0$, there exist positive numbers $k_1, k_2$ such that for any $\epsilon > 0$ sufficiently small, one has :

$$
\|\xi(\tau) - \bar{E}\| \le \epsilon k_1 + k_2 e^{-\beta \tau / \epsilon}, \quad \forall \tau \ge 0. \quad (8)
$$

*Proof.* Let us fix a positive initial condition in $\mathcal{D}$ and denote :

$$
v(t) = \sum_{i=1}^{n} \nu_i(s(t))p_i(t)
$$

along the solution of system (5). Then, the fast part of the dynamical system (5) can be
written as the following non-autonomous system in variables $(b, z = b + s)$.

$$
\left\{
\begin{array}{l}
\dot{b} = \bar{\mu}(z-b)b - Db + \epsilon v(t)b \\
\dot{z} = D(S_{in}-z)
\end{array}
\right.
\qquad (9)
$$

Recall from Lemma 2.1 that the solutions of (5) remain in the bounded domain $\mathcal{D}$, and
notice that $|v(\cdot)|$ is bounded by 1 whatever is $\epsilon$. Then, $z(\cdot)$ being solution of the differential
equation $\dot{z} = D(S_{in} - z)$, the function $b(\cdot)$ is bounded by functions $b^{-}(\cdot)$, $b^{+}(\cdot)$, defined
as solutions of

$$
\begin{cases}
\dot{b}^{+} = (\bar{\mu}(z(t) - b^{+}) - D + \epsilon)b^{+} \\
\dot{b}^{-} = (\bar{\mu}(z(t) - b^{-}) - D - \epsilon)b^{-}
\end{cases}
$$

with $b^{+}(0) = b^{-}(0) = b(0)$. One has clearly :

$$
\xi(t) \in [b^{-}(t), b^{+}(t)] \times [s^{-}(t), s^{+}(t)], \quad \forall t \ge 0 \quad (10)
$$
where $s^{-}(t) = z(t) - b^{+}(t)$ and $s^{+}(t) = z(t) - b^{-}(t)$. Similarly, the solution $\tilde{\xi}(\cdot)$ of the boundary-layer dynamics (7) with initial condition $\tilde{\xi}(0) = \xi(0)$ fulfills the frame property.

$$ \tilde{\xi}(t) \in [b^{-}(t), b^{+}(t)] \times [s^{-}(t), s^{+}(t)], \quad \forall t \ge 0 \quad (11) $$

One can also check that $\xi^{+} = (b^{+}, s^{-})$ and $\xi^{-} = (b^{-}, s^{+})$ are solutions of the dynamics :

$$ \begin{cases} \dot{b}^{+} = (\bar{\mu}(s^{-}) + \epsilon)b^{+} - Db^{+} \\ \dot{s}^{-} = -(\bar{\mu}(s^{-}) + \epsilon)b^{+} + D(S_{in} - s^{-}) \end{cases} \quad (12) $$

and

$$ \begin{cases} \dot{b}^{-} = (\bar{\mu}(s^{+}) - \epsilon)b^{-} - Db^{-} \\ \dot{s}^{+} = -(\bar{\mu}(s^{+}) - \epsilon)b^{-} + D(S_{in} - s^{+}) \end{cases} \quad (13) $$

respectively. Let

$$ \begin{aligned} E^{+} &= (S_{in} - \bar{\lambda}(D - \epsilon), \bar{\lambda}(D - \epsilon)), \\ E^{-} &= (S_{in} - \bar{\lambda}(D + \epsilon), \bar{\lambda}(D + \epsilon)) \end{aligned} $$

be the two equilibria of (12) and (13), respectively. Similarly to equations (7), $\xi^{+}(\cdot)$ and $\xi^{-}(\cdot)$ converge asymptotically towards $E^{+}$ and $E^{-}$, respectively. Furthermore, equilibria $E^{+}$ and $E^{-}$ are exponentially stable. Then there exist positive numbers $\gamma$ and $\beta$ such that the inequalities :

$$ \begin{aligned} & ||\xi^{+}(t) - E^{+}|| \leq \gamma ||\xi(0) - E^{+}|| e^{-\beta t}, \\ & ||\tilde{\xi}(t) - \bar{E}|| \leq \gamma ||\tilde{\xi}(0) - \bar{E}|| e^{-\beta t}, \\ & ||\xi^{-}(t) - E^{-}|| \leq \gamma ||\xi(0) - E^{-}|| e^{-\beta t} \end{aligned} \quad (14) $$

are fulfilled for any $t \ge 0$.

The function $\bar{\mu}(\cdot)$ being regular and strictly increasing, there exists a constant $k_1$ such that

$$ \bar{\lambda}(D) - \epsilon k_1/2 < \bar{\lambda}(D - \epsilon) < \bar{\lambda}(D) < \bar{\lambda}(D + \epsilon) < \bar{\lambda}(D) + \epsilon k_1/2 \quad (15) $$

for any $\epsilon > 0$ small enough.

Finally from properties (10), (11), (14) and (15), one obtains

$$ \begin{aligned} ||\xi(t) - \bar{E}|| & \le ||\xi(t) - \tilde{\xi}(t)|| + ||\bar{\xi}(t) - \bar{E}|| \\ & \le ||\xi^{+}(t) - \xi^{-}(t)|| + ||\bar{\xi}(t) - \bar{E}|| \\ & \le \epsilon k_1 + \gamma(||\xi(0) - E^{+}|| + ||\xi(0) - E^{-}|| + ||\xi(0) - \bar{E}||)e^{-\beta t} \\ & = \epsilon k_1 + k_2 e^{-\beta t} \end{aligned} $$

where $k_2 = \gamma(||\xi(0) - E^{+}|| + ||\xi(0) - E^{-}|| + ||\xi(0) - \bar{E}||)$. Finally, property (8) is fulfilled, with time $\tau = \epsilon t$. $\square$

Let $p = (p_i)_{i=1\dots n}$ be the vector of the “slow” variables, and consider the reduced dynamics :

$$ \frac{dp_i}{d\tau} = \left( \sum_{j=1}^{n} (\nu_i(\bar{s}) - \nu_j(\bar{s})) p_j \right) p_i, \quad (i = 1 \dots n), \quad (16) $$
where $\bar{s} = \bar{\lambda}(D)$. In the next section, we shall study the solutions of system (16) and compare them with the solutions of the original system (6).

## 4. The reduced dynamics

The reduced dynamical equations (16) of the “slow” part is given by the bilinear dynamical equation

$$ \frac{dp_i}{d\tau} = \sum_{j=1}^{n} A_{ij}p_j p_i \quad (17) $$

where $A = [A_{ij}]$ is a skew symmetric matrix with $A_{ij} = \nu_i(\bar{s}) - \nu_j(\bar{s})$ (for simplicity, we have dropped the $\bar{s}$ notation over the variables $p_i$). Let us consider the generic case.

Let us first consider the following assumption.

**Assumption A3.** For any $i \neq j$, one has $\nu_i(\bar{s}) \neq \nu_j(\bar{s})$.

Without any loss of generality, we can also assume that species are numbered such that

$$ \nu_n(\bar{s}) > \nu_{n-1}(\bar{s}) > \cdots > \nu_1(\bar{s}). $$

Let us then define numbers $B_i = A_{ni}$, where

$$ B_1 > B_2 > \cdots > B_{n-1} > B_n = 0. \quad (18) $$

Since $\sum_j p_j = 1$, one has

$$ \sum_j A_{ij}p_i = \sum_j (\nu_i - \nu_j)p_j = -\nu_n + \nu_i + \sum_j (\nu_n - \nu_j)p_j = -B_i + \sum_j B_j p_j, $$

and one can write equivalently the dynamical equations (17) as follows :

$$ \frac{dp_i}{d\tau} = \left( -B_i + \sum_{j=1}^{n} B_j p_j \right) p_i, \quad i = 1 \cdots n. \quad (19) $$

**REMARQUE.** — Under Assumptions A2 and A3, one has $\lambda_n(D) = s^*$ for $\epsilon$ small enough. In accordance with the CEP (Proposition 2.2), the n-th species asymptotically wins the competition because it has the unique smallest break-even concentration.

Under Assumption A3, system (19) admits exactly $n$ distinct equilibriums, which are exactly the vertexes of the simplex :

$$ S = \left\{ p \in \mathbb{R}_+^n \mid \sum_{i=1}^n p_i = 1 \right\}, $$

One can easily check that $S$ is invariant by the dynamics (19) and the eigenvalues of the Jacobian matrix at an equilibrium $\bar{p} \in S$ such that $\bar{p}_i = 1$ are $B_i$ and $B_j - B_i$ for $j \neq i$. Consequently one obtains immediately the following properties for the dynamics defined on $S$.
- when $\bar{p}_1 = 1$, $\bar{p}$ is a source,

- when $\bar{p}_n = 1$, $\bar{p}$ is a sink,

- when $\bar{p}_i = 1$ with $i \in 2 \cdots n-1$, $\bar{p}$ is a saddle point with a stable manifold of dimension $i-1$ contained in the face :

$$ \mathcal{F}_i = \left\{ p \in \mathcal{S} \mid \sum_{j=1}^{i} p_j = 1 \right\}. $$

**REMARQUE.** — Notice that solutions $q_i(\cdot)$ of dynamics $\dot{q}_i = -B_i q_i (i = 1 \cdots n)$ fulfill $\frac{d}{dt}\bar{q}_i = (-B_i + \sum_j B_j \bar{q}_j)\bar{q}_i$ with $\bar{q}_i = q_i / \sum_j q_j$. We deduce that the solutions of system (19) are given by the analytical formula :

$$ p_i(\tau) = \frac{p_i(0)e^{-B_i\tau}}{\sum_{j=1}^{n} p_j(0)e^{-B_j\tau}}, \quad i = 1 \cdots n. \qquad (20) $$

Let $p^*$ be the equilibrium $(0, \cdots, 0, 1)' \in \mathcal{S}$. Its stability property is given by the following Lemma.

**Lemma 4.1** For any initial condition $p(0) \in \mathcal{S}$ with $p_n(0) > 0$, the solution $p(\cdot)$ of the reduced dynamics (19) converges exponentially toward the equilibrium $p^*$.

*Proof.* From equations (20), one has one has $p_n(\tau) \to 1$ and $p_i(\tau) \to 0$ for any $i = 1 \cdots n-1$, when $\tau \to +\infty$. The linearized dynamics of (17) about $p^*$ is simply $\dot{p}_i = -B_i p_i (i = 1 \cdots n)$. Consequently, each component $p_i$ for $i < n$ converges exponentially towards 0 and $p_n = 1 - \sum_{i<n} p_i$ converges exponentially towards 1. $\square$

Let us now compare solutions $p(\cdot)$ of the reduced dynamics (17) with solutions $p_\epsilon(\cdot)$ of the original dynamics (6), when $\epsilon$ is small. When $x_n(0) > 0$, we already know from Corollary 2.3 and Lemma 4.1 that both $p(\cdot)$ and $p_\epsilon(\cdot)$ converge exponentially towards $p^*$. We give now a result that compares these solutions during their transient stage.

**Corollary 4.2** Assume that A2 and A3 are fulfilled. For any initial condition of (1) in $\mathcal{D}$ with $x_n(0) > 0$ and any $T > 0$, there exists $\bar{\epsilon} > 0$ such that

$$ \epsilon < \bar{\epsilon} \Rightarrow p(\tau) - p_{\epsilon}(\tau) = O(\epsilon), \text{ uniformly for } \tau > T. \qquad (21) $$

*Proof.* Recalling the facts :

- the equilibrium $\bar{E}$ of the boundary layer dynamics (7) is exponentially stable (Corollary 2.3),

- the equilibrium $p^*$ of the reduced dynamics (17) is exponentially stable (Lemma 4.1),
the Tikhonov's theorem (see for instance Theorem 9.4 in [9]) gives the conclusion (21) for any initial condition close to $(\bar{E}, p^*)$ in $\mathbb{R}_+^2 \times S$, and then extended to larger initial conditions by Proposition 3.1. □

Consider now the time function :

$$ \pi(\tau) = \sum_{j=1}^{n} B_j p_j(\tau). $$

The transient behavior of the solutions of (19) can be characterized by $\pi(0)$, as described by the following result.

**Proposition 4.3** Under Assumption A3, for any initial condition $p(0)$ in $S$, the solution $p(\cdot)$ of the reduced dynamics (19) fulfills the following properties.

- for indexes $i$ such that $\pi(0) \le B_i$, $p_i(\cdot)$ is decreasing,

- for indexes $i$ such that $\pi(0) > B_i$, $p_i(\cdot)$ is increasing up to $T_i$ such that $p_i(T_i) = B_i$ and is then decreasing. Furthermore, one has :

$$ T_i \ge \frac{1}{B_1} \log \frac{\pi(0)(B_1 - B_i)}{B_i(B_1 - \pi(0))}. \quad (22) $$

*Proof.* One has immediately, from equations (19) :

$$ \frac{d\pi}{d\tau} = -\sum_{j=1}^{n} B_j^2 p_j + \left(\sum_{j=1}^{n} B_j p_j\right)^2 = -\sum_{j=1}^{n} \phi(B_j)p_j + \phi\left(\sum_{j=1}^{n} B_j p_j\right) \quad (23) $$

where $\phi(\cdot)$ is the square function. $\phi(\cdot)$ being a convex function, one deduces that $d\pi/d\tau \le 0$. The function $\tau \mapsto \pi(\tau)$ is non-increasing. Note that one has also :

$$ \frac{dp_i}{d\tau} = (-B_i + \pi(\tau))p_i. $$

Consequently, the function $\tau \mapsto p_i(\tau)$ is always decreasing when $\pi(0) \le B_i$. Otherwise, $\tau \mapsto p_i(\tau)$ is increasing up to $T_i$ such that $\pi(T_i) = B_i$ and then decreasing.

From (23), one can derive the inequality :

$$ \frac{d\pi}{d\tau} = \pi^2 - B_1\pi + \sum_{j=1}^{n} B_j p_j(B_1 - B_j) \ge \pi^2 - B_1\pi, $$

and deduce an estimation from below of the function $\pi(\cdot)$ :

$$ \pi(\tau) \ge \pi^-(\tau), \quad \tau \ge 0, $$

where $\pi^-(\cdot)$ is solution of the differential equation :

$$ \frac{d\pi^-}{d\tau} = \pi^{-2} - B_1\pi^-, \quad \pi^-(0) = \pi(0). $$
It is straightforward to check that $\pi^-(\cdot)$ is given by the following expression :

$$ \pi^-(\tau) = \frac{\pi(0)B_1}{\pi(0) + (B_1 - \pi(0))e^{B_1\tau}} \quad (24) $$

Let us fix an initial condition $p_1(0), \dots, p_n(0)$ and consider $i_0$ the smallest index $i = 1, \dots, n$ such that $B_i < \pi(0)$. Note that inequality $B_1 > \pi(0)$ is fulfilled exactly when $i_0 \le n-1$. For $i_0 \le n-1$, the following bound from below of time $T_i$ is obtained from (24).

$$ T_i \ge \frac{1}{B_1} \log \frac{\pi(0)(B_1 - B_i)}{B_i(B_1 - \pi(0))} \quad \square $$

REMARQUE. — Note that $\tau \mapsto p_1(\tau)$ is always a non-increasing function, and the function $\tau \mapsto p_n(\tau)$ is always non-decreasing.

Let $i_0$ be the smallest index in $1, \dots, n$ such that $B_i < \pi(0)$. One has necessarily $T_{i_0} < T_{i_0+1} < \dots < T_{n-1}$. From expression (22), one can deduce the following qualitative properties :

i. when $\pi(0)$ is closed to $B_1$ (i.e. species 1 is majority at initial time), all the species concentrations, except for the species 1, are increasing for a long time;

ii. when $\pi(0)$ is closed to $B_i$ with $i > 1$, the concentrations of species $j$ for $j \le i$ are quickly decreasing.

## 5. Simulation results

Numerical simulations have been performed in order to illustrate the concepts developed here above. Figures 1 and 2 illustrate two typical behaviors. In both figures, in order to provide figures that are easily readable, we have limited the number of species to 5.

The operating conditions of the chemostat have been selected as follows :

$$ D = 0.1 \text{ h}^{-1}, \quad s_{in} = 5 \text{ g/l} \quad (25) $$

Moreover we have considered, for the sake of simplicity, yield coefficients $y_i$ equal to 1.

<figure><img src="image_2.png" /><figcaption>**Figure 1.** *Classical Competitive Exclusion Principle situation*</figcaption></figure>
We have considered Monod specific growth rates :

$$ \mu_i = \frac{\mu_{max,i}s}{K_{s,i} + s} \qquad (26) $$

where $\mu_{max,i}$ and $K_{s,i}$ are the maximum specific rate (h⁻¹) and the affinity constant (g/l) associated to each species $x_i$. For the sake of the numerical simulations and in order to fulfill the fast-slow dynamics requirements, we have selected one unique value for all the parameters $\mu_{max,i}$ :

$$ \mu_{max,i} = 0.5 \text{ h}^{-1}, \quad i = 1, \dots, 5 \qquad (27) $$

and a mean value of the affinity constants $K_{s,i}$, i.e.

$$ K_{s,mean} = 1 \text{ g/l} \qquad (28) $$

such the steady-state values of the nutrient concentration that would correspond if only one out of the five species is present deviates by no more than 2 % from the mean steady-state value $\bar{s}$ (= 0.25 g/l). This results in the following values for the $K_{s,i}$:

$$ K_{s,1} = 1.02 \text{ g/l}, \quad K_{s,2} = 1.01 \text{ g/l}, \quad K_{s,3} = 1 \text{ g/l} \qquad (29) $$

$$ K_{s,4} = 0.99 \text{ g/l}, \quad K_{s,5} = 0.98 \text{ g/l} \qquad (30) $$

For the same reasons of readability of the numerical results, we have standardized the initial conditions such that the sum of the initial values of the species $x_i$ is equal to the steady-state value of $s_{in} - s$ (= 4.75 g/l):

$$ \sum_{i=1}^{5} x_{i}(0) = s_{in} - \bar{s} = 4.75 \text{ g/l} \qquad (31) $$

In both figures, we have also drawn the time evolution of the following biodiversity index : $\sum_{i=1}^{5} p_{i}^{2}$. If its final value will eventually tend to 1 (since only one species survives), it will follow a trajectory that is representative of the distribution of the different species.

In Figure 1, we have considered the following initial conditions :

$$ x_{1} = 3.1434 \text{ g/l}, \quad x_{2} = 0.0105 \text{ g/l}, \quad x_{3} = 0.0118 \text{ g/l} \qquad (32) $$

$$ x_{4} = 0.0126 \text{ g/l}, \quad x_{5} = 1.5717 \text{ g/l} \qquad (33) $$

i.e. the initial condition of the species $x_5$ that will eventually win the competition is larger than those of the other species except for the extreme species $x_1$ (whose initial condition is twice that of $x_5$) that, from our analysis, will decrease whatever its initial value. The species $x_3$ and $x_4$ have initial conditions for which the condition of Proposition 4.1, $\Pi(0) - B_i$, is positive :

$$ \Pi(0) - B_3 = 5.0375e-04, \quad \Pi(0) - B_4 = 0.0013 \qquad (34) $$

yet very low : the predicted values of the time when they start to decrease $T_i$ are 205 and 548, respectively (and indeed very close to those observed in simulation, i.e. 206 and 553, respectively). We are therefore in a rather classical CEP configuration : $x_5$ increases and $x_1$ decreases, while $x_2$ decreases from the beginning, and $x_3$ and $x_4$ have a rather short
period of increase before decreasing. This latter variations are not clearly visible, yet they are illustrated via the biodiversity index (Figure 1, right) which initially decreases before increasing toward the value of 1.

In Figure 2, we have considered the following initial conditions :

$$x_1 = 0.1052 \text{ g/l}, x_2 = 2.5242 \text{ g/l}, s_{,3} = 2.1035 \text{ g/l} \qquad (35)$$

$$x_4 = 0.0084 \text{ g/l}, x_5 = 0.0088 \text{ g/l} \qquad (36)$$

<figure><img src="image_3.png" /><figcaption>**Figure 2.** *Practical coexistence*</figcaption></figure>

The situation is now different : the initial conditions of the three intermediate species are now larger than that of the species $x_5$. The condition of Proposition 4.1, $\Pi(0) - B_i$, are positive indeed again for $x_3$ and $x_4$ :

$$\Pi(0) - B_3 = 4.5219e-04, \Pi(0) - B_4 = 0.0013 \qquad (37)$$

indeed not very different from those of Figure 1. Yet the initial conditions are sufficiently large to exhibit practical coexistence via a long increase period for $x_3$ ($t = 2217$) and $x_4$ ($t = 3448$). this is also illustrated the biodiversity graph (Figure 2, right) which remains low until around $t = 3500$, i.e. when $x_4$ finally decreases. Note that in the present instance, the predicted values of the time when they start to decrease $T_i$ deviate more from the observed ones : they are equal to 182 and 526, respectively. Note that the time scale has changed between Figure 1 and Figure 2.

## 6. Conclusion

In this paper we have studied the behavior of multiple species competing for the same substrate in a chemostat when initial conditions can substantially modify the transient behavior of the system. We have shown in particular by considering slow-fast dynamics that intermediate species can survive for a substantial period of time before starting to decrease and leave the room for the species that has the best affinity with the nutrient. This formalizes the practical situation when coexistence of multiple species can last for long period before substantial decrease of the non-dominant species takes place. The results have been illustrated in numerical simulation.

**Acknowledgments.** This work has been achieved within the INRA-INRIA project 'ME-RE', during the visit of the second author.
# 7. Bibliographie

[1] R. ARIS, A.E. HUMPHREY, « Dynamics of a chemostat in which two organisms compete for a common substrate », *Biotechnology and Bioengineering*, vol. 19, pp. 1375-1386, 1977.

[2] G.J. BUTLER, G.S.K. WOLKOWICZ, « A mathematical model of the chemostat with a general class of functions describing nutrient uptake », *SIAM Journal on Applied Mathematics*, vol. 45, n° 1, pp. 138–151, 1985.

[3] G.J. BUTLER, S.B.. HSU, P. WALTMAN, « A mathematical model of the chemostat with periodic washout rate », *SIAM Journal on Applied Mathematics*, vol. 45, n° 3, pp. 435–449, 1985.

[4] C. CENENS, I.Y. SMETS, J.F. VAN IMPE, « Modeling the competition between floc-forming, filamentous bacteria in actiavted sludge waste water treatment systems - II. A prototype mathematical model based on kinetic selection, filamentous backbone theory », *Water Research*, vol. 34, n° 9, pp. 2535-2541, 2000.

[5] G. HARDIN, « The competition exclusion principle », *Science*, vol. 131, pp. 1292-1298, 1960.

[6] C. LOBRY, F. MAZENC, A. RAPAPORT, « Persistence in ecological models of competition for a single resource », *C.R. Acad. Sci. Paris, Ser I* vol. 340, pp. 199-240, 2004.

[7] C. LOBRY, J. HARMAND, « A new hypothesis to explain the coexistence of n species in the presence of a single resource », *C.R. Biologies*, vol. 329, pp. 40-46, 2006.

[8] C. LOBRY, F. MAZENC, A. RAPAPORT, « Sur un modèle densité-dépendant de compétition pour une resource », *C.R. Biologies*, vol. 329, pp. 63-70, 2006.

[9] H.K. KHALIL, « Nonlinear Systems », *Prentice Hall, Second Edition*, 1996.

[10] H.L. SMITH, « Competitive coexistence in an oscillating chemostat », *SIAM Journal on Applied Mathematics*, vol. 40, n° 3, pp. 498–522, 1981.

[11] H. SMITH, P. WALTMAN, « The Theory of the Chemostat : Dynamics of Microbial Competition », *Cambridge Studies in Mathematical Biology*, Cambridge University Press, 1985.

[12] G. STEPHANOPOULOS, R. ARIS, A.G. FREDERICKSON, « A stochastic analysis of the growth of competing microbial populations in a continuous biochemical reactor », *Mathematical Biosciences*, vol. 45, pp. 99-135, 1979.