# Stability of the chemostat system with a mutation factor

T. Bayen , H. Cazenave-Lacroutz , J. Coville Avignon Université, Laboratoire
de Mathématiques d’Avignon (EA 2151), F-84018 Avignon, France.
<EMAIL_ADDRESS>Université, Laboratoire de Mathématiques
d’Avignon (EA 2151), F-84018 Avignon, France. henri.cazenave-lacroutz@univ-
avignon.frUR 546 Biostatistique et Processus Spatiaux, INRAE, Domaine St Paul
Site Agroparc, F-84000 Avignon, France<EMAIL_ADDRESS>
###### Abstract

In this paper, we consider a resource-consumer model taking into account a
mutation effect between species (with constant mutation rate). The
corresponding mutation operator is a discretization of the Laplacian in such a
way that the resulting dynamical system can be viewed as a regular
perturbation of the classical chemostat system. We prove the existence of a
unique locally stable steady-state for every value of the mutation rate and
every value of the dilution rate not exceeding a critical value. In addition,
we give an expansion of the steady-state in terms of the mutation rate and we
prove a uniform persistence property of the dynamics related to each species.
Finally, we show that this equilibrium is globally asymptotically stable for
every value of the mutation rate provided that the dilution rate is with small
enough values.

Keywords : chemostat system, population dynamics, dynamical system, regular
perturbation, global stability.

##  1 Introduction

The chemostat system was introduced in the fifties to model the behavior of
bacteria competing for a same substrate (see [33, 34, 35]). It has now become
a reference model for the modeling of ecosystems (lakes, rivers,
microalgae,…), see, e.g., [23], and it is widely used in biotechnology, for
instance, for the control of the production of microalgae of interest or in
waste water treatment (see, e.g., [3, 18, 7, 4] and references herein). The
chemostat system with $n$ species competing for one same resource writes

${}\left|\begin{array}[]{cl}\dot{x}_{i}&=(\mu_{i}(s)-D)x_{i},\quad\quad 1\leq
i\leq n,\vspace{0.1cm}\\\
\dot{s}&=\displaystyle-\sum_{j=1}^{n}\frac{\mu_{j}(s)}{Y_{i}}x_{j}+D(s_{in}-s),\end{array}\right.$
(1.1)

where $x_{i}$ is the concentration of species $i$ (the consumers) and $s$
denotes the substrate concentration (the resource). The numbers $Y_{i}$ are
the yield coefficients, the parameter $s_{in}$ is the input substrate
coefficient, the functions $\mu_{i}$ are the kinetics, and $D$ is the dilution
rate. Properties of (1.1) has been studied a lot [18, 41, 19, 22, 25, 27, 28,
29, 39, 46], and one essential feature is the famous competitive exclusion
principle (CEP) which asserts that, asymptotically, only one species survives
[27, 46, 41, 25]. Many extensions of the CEP have been studied in presence of
delay, external inhibitors, or variable yields (see, e.g., [22, 28, 39] among
others). It is also worth mentioning that the CEP predicts exclusion of the
less competitive species and not coexistence in contrast with observations in
several ecosystems. That is why, extensions of the chemostat system were also
developed (such as in [19]) to cope with this reality. In this paper, we
consider another extension of the chemostat system related to the possibility
for a species to produce mutants or to appear through mutation (see, e.g.,
[35, 13, 14]). It turns out that mutation will modify the behavior of the
system leading to coexistence. There exist various approaches to model this
phenomenon: each species may convert into other species with a mutation rate
depending on various parameters such as the kinetics (see, e.g., [29] or
[21]). Throughout this paper, we shall assume that the dispersion is such that
each species $i$ converts into neighbor species $i+1$ and $i-1$ with a
constant mutation rate. This amounts to add a linear term $Tx$ in the sub-
system satisfied by the concentration vector $x$ in (1.1), where $T$ is the
mutation matrix. Our objective in this paper is to provide a thorough study of
asymptotic stability properties of the resulting system. Surprisingly, to our
best knowledge, few papers addressed this question apart [1, 13] and [5, 6]
which study a minimal time control problem to select optimally species of
interest (see also [31]).

Let us give a quick overview of [13] that introduced the chemostat system with
a mutation. The main result is a global stability property of the coexistence
steady-state provided that the kinetics are sufficiently close to a nominal
one $\mu_{0}$ as well as yield coefficients which also should be close to a
nominal value $Y_{0}$. This means that the quantities
$\|\mu_{i}-\mu_{0}\|_{L^{\infty}}$ and $|Y_{i}-Y_{0}|$ should be small enough
for every $1\leq i\leq n$ to ensure the global stability property. This result
(in the spirit of [1]) is interesting in itself but it does not predict the
behavior of the system whenever kinetics are not necessarily close to a common
one. In this paper, we consider the more general situation where kinetics are
of Monod type, but not necessary close to a nominal one. Our aim is to address
stability properties of the corresponding system. Based on experimental
studies (see, e.g., [35], we shall assume that the yield coefficients are
equal to one. As in [13], we shall see that mutation implies coexistence in
contrast with the CEP for the chemostat model.

The paper is structured as follows. In Section 2, we introduce the chemostat
system with mutation and we recall the CEP. In Section 3, we show in
Proposition 3.1 that there is exactly one locally asymptotically stable (LAS)
steady-state provided that the dilution rate does not exceed a certain value
(for which extinction of species would occur). This result extends the
analysis of [13] and relies on eigenvalue properties of a rank one
perturbation of a symmetric non-positive matrix (see [9]). In Section 4, we
compute an expansion of the steady-state in terms of the mutation factor. We
obtain that way an interesting result asserting that, at steady-state, few
species dominate, namely the one that wins the selection in absence of
mutation, and its neighbors (see, Proposition 4.1). We also study the converse
case, i.e., when the mutation factor becomes large (w.r.t. the kinetics of the
system). In Section 5, we show a uniform persistence property (see [40]). This
property asserts that, asymptotically, each species is present in the system
(and not only the total biomass [13]). This uniform persistence property
highlights the difference of the chemostat system with mutation w.r.t. the
classical chemostat system (leading to exclusion of every less competitive
species). In Section 6, we give our main result (Theorem 6.2) about global
stability of the steady-state for small enough dilution rates. To show that
the equilibrium is GAS, we proceed in three steps. First, we study stability
properties of the system without dilution rate (with mutation). Next, we show
that the GAS property is valid on an invariant attractive manifold associated
with the system for small enough values of the dilution rate. This requires to
prove a robust persistence property (in line with [13]) and to use
perturbation results of [42] (see also [38, 42, 44]). We conclude by using the
theory of asymptotically autonomous systems (see, e.g., [43]).

##  2 Recap on the chemostat model and preliminary properties

Throughout this paper, we consider a chemostat system with $n\geq 1$ species
including a mutation effect between species. We suppose that each species $i$
is able to convert into species $i+1$ and $i-1$ with a constant mutation rate.
This yields the following dynamical system

${}\left|\begin{array}[]{cl}\dot{x}_{i}&=(\mu_{i}(s)-u)x_{i}+\varepsilon(Tx)_{i},\quad\quad
1\leq i\leq n,\vspace{0.1cm}\\\
\dot{s}&=\displaystyle-\sum_{j=1}^{n}\mu_{j}(s)x_{j}+u(1-s),\end{array}\right.$
(2.1)

where:

  * $\bullet$

For every $1\leq i\leq n$, $x_{i}$ denotes the concentration of species $i$
and $s$ the substrate concentration.

  * $\bullet$

For every $1\leq i\leq n$, the kinetics $\mu_{i}$ of species $i$ is supposed
to be of Monod type, i.e., $\mu_{i}(s)=\frac{m_{i}s}{a_{i}+s}$ ($a_{i}$,
$m_{i}$ are positive numbers such that
$i\not=j\;\Rightarrow(m_{i},a_{i})\not=(m_{j},a_{j})$).

  * $\bullet$

The input substrate concentration has been renormalized to $1$ and the
dilution rate is $u\in\mathbb{R}_{+}$.

  * $\bullet$

The mutation parameter is $\varepsilon\geq 0$, $x:=(x_{1},...,x_{n})^{\top}$
denotes the column vector of the species concentrations (the symbol ⊤ is the
transpose operator), and the mutation matrix111As usual, matrices are named
using capital letters and coefficients are represented by lower case letters.
$T\in\mathbb{R}^{n\times n}$ is:

${}T:=\left[\begin{array}[]{ccccc}-1&1&0&\cdots&0\\\ 1&-2&1&\cdots&0\\\
\vdots&\ddots&\ddots&\ddots&\vdots\\\ 0&\cdots&1&-2&1\\\
0&\cdots&0&1&-1\end{array}\right].$ (2.2)

Note that the symmetric matrix $T$ corresponds to the discretization of the
one-dimensional Laplace equation (Poisson problem) with Neumann boundary
conditions. We recall that it is quasi-positive (i.e., for $i\not=j$,
$t_{i,j}\geq 0$) and irreducible (since $(T+rI_{n})^{k}$ is with positive
entries for $r=3$ and $k$ large enough). From Perron-Frobenius’s Theorem (see,
e.g., [8]), the largest eigenvalue of $T$ (called the Perron root) is simple
and the Perron vector (i.e., the corresponding unitary eigenvector) is
positive. Finally, the sum of the coefficients of $T$ on a row is always
zero222This property is also essential for proving the invariance of the set
$\Delta$ (Lemma 2.1)., so, $0$ is necessarily the Perron root of $T$ and
$a/\sqrt{n}$ is the Perron vector where $a:=(1,...,1)$. Note also that $T$ is
non-positive.

######  Remark 2.1.

More complex mutation terms between species can be also considered in the
chemostat model as for instance in [1, 29] where the mutation factor involves
the kinetics of the species. Mutation could also involve a pool of species
close to some index $i$ (not only the two closest indexes of $i$), but, in
this paper, we restrict our attention to a mutation term $\varepsilon Tx$
where $T$ is given by (2.2) (see [13]) and $\varepsilon\geq 0$ is eventually a
small parameter. System (2.1) can be viewed as an approximation of a
population dynamics model involving a phenotypic trait, see, e.g., [16, 17,
32, 36] (among others).

When $\varepsilon=0$, we retrieve the classical chemostat model with $n\geq 1$
species described by the system

${}\left|\begin{array}[]{cl}\dot{x}_{i}&=(\mu_{i}(s)-u)x_{i},\quad\quad 1\leq
i\leq n,\vspace{0.1cm}\\\
\dot{s}&=\displaystyle-\sum_{j=1}^{n}\mu_{j}(s)x_{j}+u(1-s),\end{array}\right.$
(2.3)

in such a way that (2.1) can be viewed as a regular perturbation333For the
concept of regular perturbation of a dynamical system, we refer to [2, 42, 30]
(see also references herein). of (2.3) for small values of $\varepsilon$. When
dealing with the chemostat system, it is usual to introduce the so-called
break-even concentrations $\lambda_{i}(u)\in[0,+\infty)$ that play a key role
in the chemostat system:

${}\lambda_{i}(u):=\left\\{\begin{array}[]{lll}\mu_{i}^{-1}(u)&\mathrm{if}\;\;\mu_{i}(1)>u,&\vspace{0.1cm}\\\
+\infty&\mathrm{otherwise},&\end{array}\right.\quad 1\leq i\leq n.$ (2.4)

In order to study asymptotic stability properties of (2.1), it will be helpful
to recall the global stability properties of (2.3). Doing so, set,

$E_{i}=(0,...,0,1-\lambda_{i}(u),0,...,0,\lambda_{i}(u))\in\mathbb{R}^{n+1},$

for $1\leq i\leq n$ and observe that $E_{i}$ is a steady-state of (2.3)
provided that $\lambda_{i}(u)<+\infty$. In addition, the point

$E_{wo}:=(0,...,0,1)\in\mathbb{R}^{n+1},$

is also an equilibrium of (2.3) (called washout steady-state). Thus, (2.3) has
at most $n+1$ steady-states. The well-known competitive exclusion principle
(CEP) can be now stated as follows.

######  Theorem 2.1.

$\mathrm{(i)}$. Let $u>0$. If there is a unique $1\leq i_{0}\leq n$ such that
$\lambda_{i_{0}}(u)=\min_{1\leq i\leq n}\lambda_{i}(u)<+\infty$, then, for
every initial condition $(x^{0},s^{0})\in[0,+\infty)^{n}\times[0,1]$ such that
$x^{0}_{i_{0}}>0$, the unique solution of (2.3) starting at $(x^{0},s^{0})$ at
time $0$ converges to $E_{i_{0}}$.

$\mathrm{(ii)}$. Let $u>0$. If $\min_{1\leq i\leq n}\lambda_{i}(u)=+\infty$,
then, for every initial condition
$(x^{0},s^{0})\in[0,+\infty)^{n}\times[0,1]$, the unique solution of (2.3)
starting at $(x^{0},s^{0})$ at time $0$ converges to $E_{wo}$.

######  Remark 2.2.

In case (i) of the previous theorem, if the minimum is non-unique, then,
coexistence may occur [41, 25], but, we do not develop this point here because
it is non-generic. Throughout the paper, we shall assume (if necessary) that
$u$ is such that the minimum is unique.

The competitive exclusion principle asserts a global stability property of one
species for (2.1) initially present in the vessel, i.e., only one species
survives generically (namely the one with the least break-even concentration).

######  Remark 2.3.

There are various proofs of this result (see, e.g., [25, 41, 37] among
others). When kinetics are of Monod type, a direct way is to use a Lyapunov
function. Doing so, write $E_{i_{0}}=(x_{1}^{*},...,x_{n}^{*},s^{*})$ and
(2.3) as

${}\left|\begin{array}[]{cl}\dot{\tilde{x}}_{i}&=(\mu_{i}(s^{*})-u)\tilde{x}_{i}+(\mu_{i}(s)-\mu_{i}(s^{*}))x_{i},\quad\quad
1\leq i\leq n,\vspace{0.1cm}\\\
\dot{\tilde{s}}&=-u\tilde{s}-\sum_{j=1}^{n}\mu_{j}(s)\tilde{x}_{j}-\sum_{j=1}^{n}(\mu_{j}(s)-\mu_{j}(s^{*}))x_{j}^{*},\end{array}\right.$
(2.5)

where $\tilde{x}_{i}:=x_{i}-x_{i}^{*}$ and $\tilde{s}:=s-s^{*}$. Next, it can
be verified that the function

${}V(\tilde{x},\tilde{s}):=\phi_{s^{*}}(\tilde{s})+\sum_{j=1}^{n}\frac{a_{j}+s^{*}}{a_{j}}\phi_{x_{j}^{*}}(\tilde{x}_{j})+\frac{1}{2}\Big{[}\tilde{s}+\sum_{j=1}^{n}\tilde{x}_{j}\Big{]}^{2},$
(2.6)

is a strict Lyapunov function for (2.5) where
$\phi_{a}(\sigma):=\sigma-a\ln(1+\sigma/a)$, $a>0$, see, e.g., [22, 27] and
references herein. However, even if (2.1) is a regular perturbation of (2.3),
it is an open question how to construct a Lyapunov function for (2.1) based on
(2.6) (see [22, 27]) or on relative entropy identities (see [11, 12]).
Besides, global stability property may fail to hold under small perturbations
of a dynamical system444As an example, consider the system
$\dot{x}=-x/(1+x^{2})+\varepsilon x$ for which $0$ is GAS for $\varepsilon=0$
and LAS for every $\varepsilon\in[0,1)$. But $0$ is never GAS for every
$\varepsilon>0$. We thank F. Mazenc for indicating to us such an example..

Going back to (2.1), observe that solutions to (2.1) are defined globally over
$\mathbb{R}_{+}$ and that the dynamics of $x$ can be rewritten

${}\dot{x}=B(s,u,\varepsilon)x,$ (2.7)

where

${}B(s,u,\varepsilon):=M(s)-uI_{n}+\varepsilon T\in\mathbb{R}^{n\times n},$
(2.8)

$I_{n}\in\mathbb{R}^{n\times n}$ denotes the identity matrix, and $M(s)$
stands for the diagonal matrix
$M(s):=\mathrm{diag}(\mu_{1}(s),...,\mu_{n}(s))$. Note that the matrix
$B(s,u,\varepsilon)$ is quasi-positive for every $s\in\mathbb{R}$, so
$\mathbb{R}_{+}^{n}$ is forward invariant by (2.7) (see, e.g., [13]). In
contrast with (2.3), it is enough to suppose that only one species is present
at time $0$ to ensure that for every time $t>0$, one has $x_{i}(t)>0$ for
every $1\leq i\leq n$, as we now show.

######  Property 2.1.

Let $(\varepsilon,u)\in\mathbb{R}_{+}^{*}\times\mathbb{R}_{+}$ and let
$x(\cdot)$ be a solution to (2.1). If there is $1\leq i\leq n$ such that
$x_{i}(0)>0$, then, for every time $t>0$, one has $x_{j}(t)>0$ for every
$1\leq j\leq n$.

###### Proof.

Recall that $\mathbb{R}_{+}^{n}$ is forward invariant by (2.7). We claim that
for every time $t\geq 0$, one has $x_{i}(t)>0$. Indeed, let
$t_{0}:=\inf\\{t>0\;;\;x_{i}(t)=0\\}$ and suppose that $t_{0}<+\infty$. Since
$x_{i}(t)>0$ for $t\in[0,t_{0})$, one has

$x_{i}(t_{0})=0\quad\mathrm{and}\quad\dot{x_{i}}(t_{0})=\varepsilon(x_{i+1}(t_{0})+x_{i-1}(t_{0}))\leq
0,$

which implies $x_{i+1}(t_{0})=x_{i-1}(t_{0})=0$ Observe now that
$\dot{x}_{i+1}(t_{0})=\varepsilon x_{i+2}(t_{0})$ and since $x_{i+1}(\cdot)$
vanishes at $t=t_{0}$, we deduce that

$\dot{x}_{i+1}(t_{0})=\varepsilon x_{i+2}(t_{0})\leq 0.$

Thus, one must have $x_{i+2}(t_{0})=0$. In the same way, we get that
$x_{i-2}(t_{0})=0$. By induction over $j$, we deduce that for every $1\leq
j\leq n$, one has $x_{j}(t_{0})=0$. By Cauchy-Lipschitz’s Theorem, one must
have $x\equiv 0$ over $\mathbb{R}_{+}$ which is a contradiction since
$x_{i}(0)>0$. This proves our claim.

Let us now show that $x_{i+1}$ never vanishes over $(0,+\infty)$. If there is
$t_{1}>0$ such that $x_{i+1}(t_{1})=0$, then, we would have
$\dot{x}_{i+1}(t_{1})\leq 0$, thus

$\dot{x}_{i+1}(t_{1})=\varepsilon(x_{i}(t_{1})+x_{i+2}(t_{1}))>0,$

since $x_{i}$ is positive over $\mathbb{R}_{+}$. This is a contradiction,
therefore, one must have $x_{i+1}(t)>0$ for every time $t>0$. We can repeat
this argument step by step for every species, which proves the desired
property. ∎

Note also that if $s(0)\in[0,1]$, then one has $s(t)\in[0,1]$ for every $t\geq
0$. Hence, we shall consider initial conditions in the set

$\mathcal{D}:=(\mathbb{R}_{+}^{n}\backslash\\{0\\})\times[0,1],$

when dealing with (2.1). The next property is related to the quantity

$b:=s+\sum_{j=1}^{n}x_{j},$

and it is crucial in the rest of the paper.

######  Lemma 2.1.

For every $(\varepsilon,u)\in\mathbb{R}_{+}\times\mathbb{R}_{+}^{*}$, the set

$\Delta:=\Big{\\{}(x,s)\in\mathcal{D}\;;\;\sum_{j=1}^{n}x_{j}+s=1\Big{\\}},$
(2.9)

is an invariant and attractive manifold for (2.1).

###### Proof.

From (2.1), $b$ satisfies $\dot{b}=u(1-b)$, hence $b(t)=1+(b(0)-1)e^{-tu}$ for
$t\geq 0$, whence the result. ∎

This lemma makes possible (if necessary) to reduce the stability properties of
(2.1) to the system

$\dot{x}=B\Big{(}1-\sum_{j=1}^{n}x_{j},u,\varepsilon\Big{)}x,$

obtained from (2.1) by considering conditions in $\Delta$. Note that if
$(x^{0},s^{0})\in\Delta$, then $\sum_{j=1}^{n}x_{j}^{0}\leq 1$, that is why,
it is also useful to introduce the set

${}\mathcal{D}^{\prime}:=\Big{\\{}x\in[0,+\infty)^{n}\;;\;\sum_{j=1}^{n}x_{j}\leq
1\Big{\\}},$ (2.10)

when dealing with initial conditions in $\Delta$. The next property is well-
known for (2.3) (see, e.g., [25, 41]) and it remains unchanged for (2.1).

######  Property 2.2.

For every $u>0$, there is $c_{u}>0$ such that for every $\varepsilon\geq 0$
and for every initial condition in $\mathcal{D}$, the unique corresponding
solution to (2.1) satisfies:

${}\liminf_{t\rightarrow+\infty}s(t)\geq c_{u}.$ (2.11)

###### Proof.

Let $\delta_{u}:=\sup\\{s\in[0,1]\;;\;\max_{1\leq j\leq
n}\mu_{j}(s)\leq\frac{u}{8}\\}$, Since $\mu_{i}(0)=0$ for every $1\leq i\leq
n$, $\delta_{u}$ is well-defined. and we can set
$c_{u}:=\min(\delta_{u},1/2)$. From Lemma 2.1, there is $t_{0}\geq 0$ such
that $\sum_{j=1}^{n}x_{j}(t)\leq 2$ for every $t\geq t_{0}$. If now there is
$t\geq t_{0}$ such that $s(t)\leq c_{u}$, one has

$\dot{s}(t)=-\sum_{j=1}^{n}\mu_{j}(s(t))x_{j}(t)+u(1-s(t))\geq-\frac{u}{8}\sum_{j=1}^{n}x_{j}(t)+\frac{u}{2}\geq-\frac{u}{4}+\frac{u}{2}=\frac{u}{4}.$

From the preceding inequality, $[0,c_{u}]$ is a repelling set for the dynamics
of $s$ which then implies (2.11). ∎

##  3 Local asymptotic stability

###  3.1 Existence of a locally stable equilibrium

Throughout the paper, given a symmetric matrix $A\in\mathbb{R}^{n\times n}$,
we denote by $\lambda(A)$ its largest eigenvalue.

######  Lemma 3.1.

Let $(\varepsilon,u)\in\mathbb{R}_{+}\times\mathbb{R}_{+}^{*}$. Then, one has
$\lambda(B(0,u,\varepsilon))<0$ and the mapping
$s\mapsto\lambda(B(s,u,\varepsilon))$ is increasing over $[0,1]$.

###### Proof.

Observe that $B(0,u,\varepsilon)=-uI_{n}+\varepsilon T$ which implies
$\lambda(B(0,u,\varepsilon))=-u<0$. Now, thanks to the Perron-Frobenius
Theorem [8], for every $s\in[0,1]$, $\lambda(B(s,u,\varepsilon))$ exists and
is of multiplicity one. It is also the unique eigenvalue associated with a
positive eigenvector. Recall now that given two quasi-positive irreducible and
symmetric matrices $C,D\in\mathbb{R}^{n\times n}$ such that $c_{i,j}\leq
d_{i,j}$ for every $1\leq i,j\leq n$ (with a strict inequality for at least
one coefficient), one has $\lambda(C)<\lambda(D)$, see [8]. Since for every
$1\leq i\leq n$, $\mu_{i}$ is increasing, so is
$s\mapsto\lambda(B(s,u,\varepsilon))$. ∎

Next, we study the existence of a locally stable equilibrium point for (2.1).
Doing so, we shall use a result of [9] about rank one perturbations of a
singular $M$-matrix $A=\rho(H)I_{n}-H\in\mathbb{R}^{n\times n}$ where
$\rho(H)$ denotes the spectral radius of a given matrix
$H\in\mathbb{R}^{n\times n}$. Let us recall the concept of $M$-matrix.

######  Definition 3.1.

Given $A=(a_{i,j})\in\mathbb{R}^{n\times n}$, we say that $A$ is an $M$-matrix
if there exists a matrix $H\in\mathbb{R}^{n\times n}$ and $s>\rho(H)$ such
that

$A=sI_{n}-H\quad\mathrm{and}\quad s>\rho(H).$

Notice that if $A$ is an $M$-matrix, then its eigenvalues are with nonnegative
real parts and $a_{i,j}\leq 0$ for $i\not=j$. Theorem 2.7 of [9] provides
sufficient conditions for a matrix $A+vw^{\top}$ (where
$v,w\in\mathbb{R}^{n}$) to be positive stable if $0$ is a geometrically simple
eigenvalue of $A$. We refer to [9] for the precise statement of those
conditions. The next Proposition (point (ii) only) extends the analysis of
[13] showing that, depending on the values of $(\varepsilon,u)$, (2.1) has a
unique locally stable equilibrium. The notation $\|\cdot\|$ stands for the
euclidean norm in $\mathbb{R}^{n}$.

######  Proposition 3.1.

(i) If $(\varepsilon,u)\in\mathbb{R}_{+}^{*}\times\mathbb{R}_{+}^{*}$ is such
that $\lambda(B(1,u,\varepsilon))\leq 0$, then, the washout steady-state
$E_{wo}$ is the only equilibrium of (2.1) and it is stable. If
$\lambda(B(1,u,\varepsilon))<0$, it is globally asymptotically stable.
(ii) If $(\varepsilon,u)\in\mathbb{R}_{+}^{*}\times\mathbb{R}_{+}^{*}$ is such
that $\lambda(B(1,u,\varepsilon))>0$, then, (2.1) admits a unique locally
stable equilibrium
$E_{\varepsilon,u}:=(x^{\varepsilon,u},s^{\varepsilon,u})\in(0,+\infty)^{n}\times(0,1)$
called coexistence steady-state and $E_{wo}$ is unstable.

###### Proof.

For sake of completeness, we give the proof of (i) which can also be found in
[13]. If $(x,s)$ is a steady-state of (2.1), then

${}\left\\{\begin{array}[]{ll}B(s,u,\varepsilon)x&=0,\vspace{0.1cm}\\\
\sum_{j=1}^{n}\mu_{j}(s)x_{j}&=u(1-s).\end{array}\right.$ (3.1)

The equation $B(s,u,\varepsilon)x=0$ with $x\not=0$ implies that 0 is an
eigenvalue, hence $\lambda(B(s,u,\varepsilon))\geq 0$. It is possible only if
$s=1$. Indeed, otherwise, since
$s<1\;\Rightarrow\;\lambda(B(s,u,\varepsilon))<\lambda(B(1,u,\varepsilon))$,
we would have $\lambda(B(s,u,\varepsilon))<0$ and a contradiction. It follows
from (3.1) that any equilibrium verifies $x=0$, so, the only equilibrium point
is the washout. The Jacobian of (2.1) at $E_{wo}$ is the block matrix

$\left[\begin{array}[]{cc}B(1,u,\varepsilon)&0\\\
-\mu_{1}(1)\cdots-\mu_{n}(1)&-u\end{array}\right]\in\mathbb{R}^{(n+1)\times(n+1)}.$

If $\lambda(B(1,u,\varepsilon))\leq 0$, then, $\dot{x}=B(s,u,\varepsilon)x\leq
B(1,u,\varepsilon)x$, thus there is $C_{1}>0$ such that $\|x(t)\|\leq
C_{1}\|x(0)\|$ for every $t\geq 0$. Using that $\Delta$ is attractive for
(2.1) (recall Lemma 2.1), we deduce that $E_{wo}$ is stable. If, in addition,
$\lambda(B(1,u,\varepsilon))<0$, $B(1,u,\varepsilon)$ is a Hurwitz matrix,
and, thanks to the inequality $\dot{x}\leq B(1,u,\varepsilon)x$, we deduce
that $x(t)\rightarrow 0$ as $t\rightarrow+\infty$ which proves the desired
property using Lemma 2.1.

In case (ii), we find two equilibria depending if $s=1$ or $s<1$. If $s=1$,
then $x=0$ and the corresponding steady-state is the washout that is unstable
since $\lambda(B(1,u,\varepsilon))>0$. The other possible steady-states
satisfy (3.1) with $s<1$, so,
$x\in\mathrm{Ker}(B(s,u,\varepsilon))\backslash\\{0\\}$. But, the largest
eigenvalue of $B(s,u,\varepsilon)$ is the only one with a positive eigenvector
(thanks to the Perron-Frobenius Theorem). So, we necessarily have
$\lambda(B(s,u,\varepsilon))=0$ which has a unique solution
$s^{\varepsilon,u}$ (because of the monotonicity of
$\lambda(B(\cdot,u,\varepsilon))$ and the fact that
$\lambda(B(0,u,\varepsilon))\lambda(B(1,u,\varepsilon))<0$). Hence, zero is
the Perron root of $B(s^{\varepsilon,u},u,\varepsilon)$ (it is a simple
eigenvalue) and we denote by $a^{\varepsilon,u}\in\mathbb{R}^{n}$ its Perron
vector. We deduce that $x$ necessarily satisfies

$x=\nu
a^{\varepsilon,u}\quad\mathrm{and}\quad\sum_{j=1}^{n}x_{j}+s^{\varepsilon,u}=1,$

where $\nu\in\mathbb{R}_{+}^{*}$. These two equalities define a unique point
$x^{\varepsilon,u}\in\mathbb{R}^{n}$ such that $x^{\varepsilon,u}_{i}>0$ for
every $1\leq i\leq n$. Let us now turn to the local asymptotic stability
property. Observe that (2.1) is equivalent to

$\left|\begin{array}[]{cl}\dot{x}&=B(b-\sum_{j=1}^{n}x_{j},u,\varepsilon)x,\vspace{0.1cm}\\\
\dot{b}&=u(1-b),\end{array}\right.$

(recall that $b=s+\sum_{j=1}^{n}x_{j}$). The Jacobian matrix of the preceding
system at $(x^{\varepsilon,u},1)$ is

$J_{\varepsilon,u}:=\left[\begin{array}[]{cc}A_{\varepsilon,u}&d_{\varepsilon,u}\\\
0&-u\end{array}\right]\in\mathbb{R}^{(n+1)\times(n+1)},$

where
$d_{\varepsilon,u}:=M^{\prime}(s^{\varepsilon,u})x^{\varepsilon,u}\in\mathbb{R}^{n}$
and
$A_{\varepsilon,u}:=B(s^{\varepsilon,u},u,\varepsilon)-d_{\varepsilon,u}a^{\top}\in\mathbb{R}^{n\times
n}$ is a rank-one perturbation of $B(s^{\varepsilon,u},u,\varepsilon)$. For
proving our claim, it is then enough to show that $A_{\varepsilon,u}$ is a
Hurwitz matrix. Observe that one has
$-B(s^{\varepsilon,u},u,\varepsilon)=\rho(H)I_{n}-H$ where
$H:=\rho(B(s^{\varepsilon,u},u,\varepsilon))I_{n}+B(s^{\varepsilon,u},u,\varepsilon)$.
Hence, $-B(s^{\varepsilon,u},u,\varepsilon)$ can be written as a singular
$M$-matrix which is thus non-negative. We can now apply Theorem 2.7 (v) of [9]
with the matrix $H$ (for this, note that $d_{\varepsilon,u}$ and $a$ have
positive coefficients) and deduce that $-A_{\varepsilon,u}$ is strictly
positive stable (which means that all eigenvalues of $-A_{\varepsilon,u}$ are
with positive real parts). We can thus conclude that $J_{\varepsilon,u}$ is a
Hurwitz matrix which ends the proof. ∎

In the rest of the paper, we keep the notation $a^{\varepsilon,u}$ for the
Perron vector associated with the $0$ eigenvalue of the matrix
$B(s^{\varepsilon,u},u,\varepsilon)$.

###  3.2 Occurrence of the washout and coexistence steady-states

In this part, we make more explicit the condition about
$\lambda(B(1,u,\varepsilon))$ which separates washout and coexistence
equilibria in Proposition 3.1. It is convenient to introduce the functions

$\bar{\mu}(s):=\frac{1}{n}\sum_{j=1}^{n}\mu_{j}(s)\;;\;\hat{\mu}(s):=\max(\mu_{1}(s),...,\mu_{n}(s)),$

which are increasing over $\mathbb{R}_{+}$. Also, we set $m:=\hat{\mu}(1)$.

######  Proposition 3.2.

For every $(\varepsilon,u)\in\mathbb{R}_{+}^{*}\times\mathbb{R}_{+}^{*}$ such
that $\lambda(B(1,u,\varepsilon))>0$, one has:

${}\bar{\mu}(s^{\varepsilon,u})\leq u\leq\hat{\mu}(s^{\varepsilon,u})\leq
u+2\varepsilon.$ (3.2)

In addition, for every
$(\varepsilon,u)\in\mathbb{R}_{+}^{*}\times\mathbb{R}_{+}^{*}$, the quantity
$\lambda(B(1,u,\varepsilon))$ satisfies the following inequalities:

${}\max\Big{(}\hat{\mu}(1)-u-2\varepsilon,\bar{\mu}(1)-u\Big{)}\leq\lambda(B(1,u,\varepsilon))\leq\hat{\mu}(1)-u.$
(3.3)

###### Proof.

First, observe that $B(s,u,\varepsilon)\leq\hat{B}(s)$ (componentwise) where
$\hat{B}(s):=\max_{1\leq j\leq n}(\mu_{j}(s)-u)I_{n}+\varepsilon
T\in\mathbb{R}^{n\times n}$, and that both matrices are quasi-positive and
irreducible. We can thus deduce that
$\lambda(B(s,u,\varepsilon))\leq\lambda(\hat{B}(s))=\max_{1\leq j\leq
n}\mu_{j}(s)-u$. This gives the second inequality in (3.3).

Now, for $s\in[0,1]$, we set $\zeta:=\lambda(B(s,u,\varepsilon))$ and let $x$
be an eigenvector of $B(s,u,\varepsilon)$ with eigenvalue $\zeta$. Recall from
Perron-Frobenius’s Theorem that $x>0$. The equality $B(s,u,\varepsilon)x=\zeta
x$ rewrites

${}\left\\{\begin{array}[]{ll}(\mu_{1}(s)-u-\varepsilon)x_{1}+\varepsilon
x_{2}&=\zeta x_{1},\\\
(\mu_{i}(s)-u-2\varepsilon)x_{i}+\varepsilon(x_{i+1}+x_{i-1})&=\zeta
x_{i},\quad 2\leq i\leq n-1,\\\ (\mu_{n}(s)-u-\varepsilon)x_{n}+\varepsilon
x_{n-1}&=\zeta x_{n},\end{array}\right.$ (3.4)

Summing those $n$ equalities with $s=s^{\varepsilon,u}$ and
$x=x^{\varepsilon,u}$ gives

$\sum_{j=1}^{n}\mu_{j}(s^{\varepsilon,u})x_{j}^{\varepsilon,u}=u\sum_{j=1}^{n}x_{j}^{\varepsilon,u},$

thus, one obtains
$u=\sum_{j=1}^{n}p_{j}^{\varepsilon,u}\mu_{j}(s^{\varepsilon,u})$ (where for
$1\leq i\leq n$,
$p_{i}^{\varepsilon,u}:=\frac{x_{i}^{\varepsilon,u}}{\sum_{j=1}^{n}x_{j}^{\varepsilon,u}}$).
Because $p_{i}^{\varepsilon,u}\leq 1$, we deduce that
$u\leq\hat{\mu}(s^{\varepsilon,u})$, which gives the second inequality in
(3.2). From (3.4) with $s\in[0,1]$, we also get

${}\zeta=\mu_{i}(s)-u-a_{i}\varepsilon+\varepsilon\frac{x_{i+1}+x_{i-1}}{x_{i}},$
(3.5)

with $a_{i}=2$ for $2\leq i\leq n-1$, $a_{1}=a_{n}=1$, and the convention that
$x_{0}=x_{n+1}=0$. Since $x_{i}>0$ for $1\leq i\leq n$, this equality entails

$\zeta\geq\max_{1\leq j\leq n}\mu_{j}(s)-u-2\varepsilon.$

From the preceding inequality, we can deduce the last inequality in (3.2)
(taking $s=s^{\varepsilon,u})$ and also the inequality
$\lambda(B(1,u,\varepsilon))\geq\hat{\mu}(1)-u-2\varepsilon$ in (3.3) (taking
$s=1$).

To conclude, we need to prove the two inequalities in (3.2)-(3.3) involving
the mean values of the kinetics. Doing so, we sum equalities (3.5) which gives

$\zeta=\frac{1}{n}\sum_{j=1}^{n}\mu_{j}(s)-u-\frac{2(n-1)}{n}\varepsilon+\frac{\varepsilon}{n}\left[\sum_{j=1}^{n-1}\frac{x_{j+1}}{x_{j}}+\sum_{j=2}^{n}\frac{x_{j-1}}{x_{j}}\right].$

Applying the arithmetic-geometric mean inequality yields

$\frac{1}{n-1}\sum_{j=1}^{n-1}\frac{x_{j+1}}{x_{j}}\geq\left(\prod_{j=1}^{n-1}\frac{x_{j+1}}{x_{j}}\right)^{\frac{1}{n-1}}=\left(\frac{x_{n}}{x_{1}}\right)^{\frac{1}{n-1}}\;;\;\frac{1}{n-1}\sum_{j=1}^{n-1}\frac{x_{j-1}}{x_{j}}\geq\left(\frac{x_{1}}{x_{n}}\right)^{\frac{1}{n-1}},$

which implies

$\zeta\geq\frac{1}{n}\sum_{j=1}^{n}\mu_{j}(s)-u-\frac{2(n-1)}{n}\varepsilon+\frac{(n-1)\varepsilon}{n}\left[\left(\frac{x_{n}}{x_{1}}\right)^{\frac{1}{n-1}}+\left(\frac{x_{1}}{x_{n}}\right)^{\frac{1}{n-1}}\right].$

Using that $y^{\frac{1}{n-1}}+y^{-\frac{1}{n-1}}\geq 2$ for every $y>0$, we
obtain the inequality $\zeta\geq\bar{\mu}(s)-u$. Specializing this inequality
with $s=s^{\varepsilon,u}$ and $s=1$ gives us the left inequalities in (3.2)
and (3.3) which concludes the proof. ∎

Thanks to this proposition, we can make the following observations:

  * $\bullet$

From (3.2)-(3.3), the steady-state $E_{\varepsilon,u}$ occurs whenever
$\bar{\mu}(1)>u$ and this condition does not depend on the parameter
$\varepsilon$.

  * $\bullet$

Observe also that, thanks to those inequalities, we recover the fact that if
$u\geq m$, then, only washout occurs.

For every $\varepsilon\geq 0$, we can also uniquely define a critical value
for the dilution rate

$u_{c}(\varepsilon):=\lambda(M(1)+\varepsilon T),$

which is such that only washout occurs if the dilution rate is such that
$u\geq u_{c}(\varepsilon)$ (according to Proposition 3.1 (i)). From (3.3) and
the previous remarks, the value $u_{c}(\varepsilon)$ satisfies:

${}\forall\varepsilon\geq 0,\;u_{c}(\varepsilon)\in[\hat{u}(\varepsilon),m],$
(3.6)

where $\hat{u}(\varepsilon):=\max(m-2\varepsilon,\bar{\mu}(1))$.
Interestingly, the presence of mutation in the system implies occurrence of
the washout for values of the dilution rate in the interval
$[m-2\varepsilon,m)$ for which the species with the least break-even
concentration would survive (without mutation). In addition, we can observe
that the larger the mutation rate is, the lower the dilution need to be to
avoid washout.

We now recall a result related to the differentiability of $\lambda(\cdot)$
that will be applied several times in this paper. Given a symmetric quasi-
positive matrix $A$, the largest eigenvalue of $A$ is simple and thus
$\lambda(\cdot)$ is analytic as a function of its $n^{2}$ coefficients in some
neighborhood of $A$ in the symmetric matrices space (see, e.g., [47, 15]). In
addition, the first derivative of $\lambda(\cdot)$ (i.e., the matrix whose
$(i,j)$ entry is $\frac{\partial\lambda}{\partial a_{i,j}}(A))$ is given by

${}D_{1}\lambda(A)=ww^{\top},$ (3.7)

where $w$ denotes the Perron vector associated with $\lambda(A)$ (see [15,
24]).

######  Proposition 3.3.

The function $\varepsilon\mapsto u_{c}(\varepsilon)$ is non-increasing over
$\mathbb{R}_{+}$. In addition, one has $u_{c}(0)=m$, and
$u_{c}(\varepsilon)\rightarrow\bar{\mu}(1)$ as
$\varepsilon\rightarrow+\infty$.

###### Proof.

Applying the previous property with the symmetric quasi-positive matrix
$B(1,0,\varepsilon)=M(1)+\varepsilon T$ gives

$u_{c}^{\prime}(\varepsilon)=\sum_{1\leq i,j\leq
n}\frac{\partial\lambda(B(1,0,\varepsilon))}{\partial b_{i,j}}\frac{\partial
b_{i,j}(1,0,\varepsilon)}{\partial\varepsilon}=(v^{\varepsilon})^{\top}Tv^{\varepsilon}\leq
0,$

where $v^{\varepsilon}$ is the Perron vector associated with
$B(1,0,\varepsilon)$ and $b_{i,j}(1,0,\varepsilon)$ denote the $n^{2}$ entries
of $B(1,0,\varepsilon)$. This shows that $u_{c}(\cdot)$ is non-increasing over
$\mathbb{R}_{+}$. Now, from the CEP, we have immediately $u_{c}(0)=m$.
Finally, using (3.7), we have the expansion

${}u_{c}(\varepsilon)=\varepsilon\lambda\left(T+\frac{1}{\varepsilon}M(1)\right)=\varepsilon\left[\lambda(T)+\frac{1}{\varepsilon}\frac{a^{\top}M(1)a}{a^{\top}a}+o\left(\frac{1}{\varepsilon}\right)\right]=\bar{\mu}(1)+o(1),$
(3.8)

as $\varepsilon\rightarrow+\infty$, which concludes the proof. ∎

###  3.3 Global stability in the two species case

In this part, we prove that $E_{\varepsilon,u}$ is GAS for (2.1) when $n=2$
and we also give explicit expressions for the steady-state and the critical
value of the dilution rate. We start by addressing the global stability
property.

######  Proposition 3.4.

For $n=2$ and
$(\varepsilon,u)\in\mathbb{R}_{+}^{*}\times(0,u_{c}(\varepsilon))$,
$E_{\varepsilon,u}$ is globally asymptotically stable in $\mathcal{D}$.

###### Proof.

For those initial conditions in the set $\Delta$, (2.1) is equivalent to

${}\dot{x}=h(x),$ (3.9)

where $h:\mathcal{D}^{\prime}\rightarrow\mathbb{R}^{2}$ is given by (recall
(2.10))

$h(x):=\left(\begin{array}[]{c}(\mu_{1}(1-x_{1}-x_{2})-u)x_{1}+\varepsilon(x_{2}-x_{1})\\\
(\mu_{2}(1-x_{1}-x_{2})-u)x_{2}+\varepsilon(x_{1}-x_{2})\end{array}\right).$

If we set $\varphi(x_{1},x_{2}):=\frac{1}{x_{1}x_{2}}$ for $x_{1},x_{2}>0$, a
direct computation shows that the quantity

$\frac{\partial(\varphi h_{1})}{\partial x_{1}}+\frac{\partial(\varphi
h_{2})}{\partial
x_{1}}=-\varepsilon\frac{x_{1}^{2}+x_{2}^{2}}{x_{1}^{2}x_{2}^{2}}-\frac{\mu^{\prime}_{1}(1-x_{1}-x_{2})}{x_{2}}-\frac{\mu^{\prime}_{2}(1-x_{1}-x_{2})}{x_{2}},$

is negative in the interior of $\mathcal{D}^{\prime}$. It follows from the
Bendixson-Dulac Theorem that no periodic orbits occurs in
$\mathcal{D}^{\prime}$ for (3.9). Now, Proposition 3.1 implies that only two
equilibria occur, namely the washout $(0,0)$ which is unstable and the point
$x^{\varepsilon,u}$ in the interior of $\mathcal{D}^{\prime}$ which is locally
asymptotically stable. Since there are no periodic orbits, we deduce that
$x^{\varepsilon,u}$ is globally asymptotically stable for (2.1) restricted to
$\Delta$. Now, coming back to (2.1) for $n=2$, the sub-system satisfied by $x$
reads

${}\dot{x}=\tilde{h}(t,x),$ (3.10)

where $\tilde{h}:\mathbb{R}\times\mathbb{R}^{2}\rightarrow\mathbb{R}^{2}$ is
defined as

$\tilde{h}(t,x):=\left(\begin{array}[]{c}(\mu_{1}(b(t)-x_{1}-x_{2})-u)x_{1}+\varepsilon(x_{2}-x_{1})\\\
(\mu_{2}(b(t)-x_{1}-x_{2})-u)x_{2}+\varepsilon(x_{1}-x_{2})\end{array}\right).$

Clearly, $\tilde{h}(t,x)\rightarrow h(t,x)$ uniformly locally w.r.t. $x$, thus
(3.10) is a non-autonomous perturbation of (3.9). Since every solution to
(3.9) converges to $x^{\varepsilon,u}$, we deduce from [43] that every
solution to (3.10) also converges to this point. To conclude, let us given a
solution $(x(\cdot),s(\cdot))$ of (2.1). We have proved that $x(t)\rightarrow
x^{\varepsilon,u}$ when $t$ goes to infinity. Since
$s(t)+x_{1}(t)+x_{2}(t)\rightarrow 1$ as $t\rightarrow+\infty$, we deduce that
$s(t)\rightarrow s^{\varepsilon,u}$ as $t$ goes to infinity which ends the
proof. ∎

We now turn to explicit expressions involving the steady-state. The sub-system
$(x_{1},x_{2})$ of (2.1) rewrites

$\dot{x}=B(s,u,\varepsilon)x\quad\mathrm{with}\quad
B(s,u,\varepsilon)=\left[\begin{array}[]{cc}\mu_{1}(s)-(u+\varepsilon)&\varepsilon\\\
\varepsilon&\mu_{2}(s)-(u+\varepsilon)\end{array}\right].$

The largest eigenvalue of $B(s,u,\varepsilon)$ can thus be explicitly
computed:

${}\lambda(B(s,u,\varepsilon))=-\varepsilon-u+\frac{\mu_{1}(s)+\mu_{2}(s)}{2}+\frac{1}{2}\sqrt{(\mu_{1}(s)-\mu_{2}(s))^{2}+4\varepsilon^{2}}.$
(3.11)

Hence, the coexistence steady-state $(x^{\varepsilon,u},s^{\varepsilon,u})$
exists provided that $\lambda(B(1,u,\varepsilon))>0$ (see Proposition 3.1)
which amounts to saying that the dilution rate fulfills the inequality

$u<u_{c}(1)=-\varepsilon+\frac{\mu_{1}(1)+\mu_{2}(1)}{2}+\frac{1}{2}\sqrt{(\mu_{1}(1)-\mu_{2}(1))^{2}+4\varepsilon^{2}}.$

For a given $u$ satisfying the previous inequality, we can compute
$s^{\varepsilon,u}$ numerically solving $\lambda(B(s,u,\varepsilon))=0$ w.r.t.
$s$, thanks to (3.11) (see Fig. 1). It follows that the dilution rate $u$ is
related to $s^{\varepsilon,u}$ via the equality

$u=-\varepsilon+\frac{\mu_{1}(s^{\varepsilon,u})+\mu_{2}(s^{\varepsilon,u})}{2}+\frac{1}{2}\sqrt{(\mu_{1}(s^{\varepsilon,u})-\mu_{2}(s^{\varepsilon,u}))^{2}+4\varepsilon^{2}}.$

Using that $x^{\varepsilon,u}_{1}+x^{\varepsilon,u}_{2}+s^{\varepsilon,u}=1$,
one also obtains

$x^{\varepsilon,u}_{1}=\frac{(1-s^{\varepsilon,u})(u-\mu_{2}(s^{\varepsilon,u}))}{\mu_{1}(s^{\varepsilon,u})-\mu_{2}(s^{\varepsilon,u})}\;;\;x^{\varepsilon,u}_{2}=\frac{(1-s^{\varepsilon,u})(u-\mu_{1}(s^{\varepsilon,u}))}{\mu_{2}(s^{\varepsilon,u})-\mu_{1}(s^{\varepsilon,u})}.$

######  Remark 3.1.

The previous expressions of $x^{\varepsilon,u}_{i}$ are valid if the kinetics
do not intersect. If there is a (unique) $\bar{s}\in(0,1)$ such that
$\mu_{1}(\bar{s})=\mu_{2}(\bar{s})$, then, these expressions are valid only if
$u\not=\mu_{1}(\bar{s})$. If $u=\bar{u}:=\mu_{1}(\bar{s})$, then, one has
$s^{\varepsilon,\bar{u}}=\bar{s}$ and
$x^{\varepsilon,\bar{u}}=\big{(}\frac{1-\bar{s}}{2},\frac{1-\bar{s}}{2}\big{)}$.

Fig. 1 depicts $\varepsilon\mapsto x^{\varepsilon,u}$, $\varepsilon\mapsto
s^{\varepsilon,u}$, and $\varepsilon\mapsto u_{c}(\varepsilon)$ for a fixed
$u>0$ such that species $2$ survives when $\varepsilon=0$ (see the plof of
$\mu_{1}$ and $\mu_{2}$ below). We verify numerically that
$E_{\varepsilon,u}\rightarrow E_{2}$ as $\varepsilon\downarrow 0$ (see Section
4) and that $u_{c}(0)=\max(\mu_{1}(1),\mu_{2}(1))$ and
$u_{c}(+\infty)=\frac{\mu_{1}(1)+\mu_{2}(1)}{2}$ (see (3.8)).

Figure 1: Plot of $\mu_{1}(s):=\frac{s}{s+1}$ (in red) and
$\mu_{2}(s):=\frac{0.7s}{0.5+s}$ (in blue) on Fig. left. For $u=0.2$ and
$\varepsilon=0$, species $2$ survives and $\lambda_{2}(0.2)=0.2$ (recall
(2.4)). Next, from left to right, plot of $s^{\varepsilon,u}$,
$x^{\varepsilon,u}_{1}$, and $x^{\varepsilon,u}_{2}$ as a function of
$\varepsilon\in[0,2]$. Fig. right depics $u_{c}(\varepsilon)$ illustrating
(3.8).

###  3.4 Illustration of the global stability property for $n\geq 3$

In view of the local stability property of $E_{\varepsilon,u}$ and the global
stability of this equilibrium for $n=2$, one can wonder if this property
remains valid for $n\geq 3$, $\varepsilon>0$, and
$u\in(0,u_{c}(\varepsilon))$. Although we know the behavior of (2.1) for
$\varepsilon=0$, it turns out that this question is delicate even if
$\varepsilon$ is arbitrarily small (see also Remark 2.3). In Section 6, we
address this question when $\varepsilon>0$ is fixed and $u$ is with small
enough values.

We present below numerical simulations of solutions to (2.1) for $n=20$,
$\varepsilon=1$, and $u=5$, see Fig. 2. The kinetics associated with the
species are arbitrary functions of Monod type. Our observations are as
follows:

  * $\bullet$

First, we observe convergence of the system to the coexistence equilibrium for
a large set of initial conditions.

  * $\bullet$

Interestingly, we also see that even though the system converges to the
coexistence equilibrium, very few species have a significant concentration
asymptotically. We shall give an explanation of this phenomenon in Section 4
for small values of the parameter $\varepsilon>0$.

Figure 2: Fig. left : plot of $x_{i}(\cdot)$. Fig. middle: plot of the
$L^{\infty}$-norm error $\|x(\cdot)-x^{\varepsilon,u}\|_{L^{\infty}}$ between
the solution and the steady-state. Fig. right: plot of the value of species
concentrations when $t$ goes to infinity. Data are such that $n=20$,
$\varepsilon=1$, $u=5$ and $\mu_{i}(s)=\frac{20s}{a_{i}+s}$ with
$a_{i}=1+\frac{1}{2}(10-i)^{2}$, $1\leq i\leq 20$.

##  4 Behavior of the coexistence steady state

In this part, we study the behavior of the coexistence equilibrium w.r.t. the
parameter $\varepsilon$. Based on the implicit function theorem, we give an
expansion of $E_{\varepsilon,u}$ up to the first order as
$\varepsilon\downarrow 0$ (for a fixed dilution rate $u$) and we also study
its limit as $\varepsilon\rightarrow+\infty$. Recall that $m=\hat{\mu}(1)$.

######  Proposition 4.1.

Suppose that $u\in(0,m)$ and that there is a unique $1\leq i_{0}\leq n$ such
that $\lambda_{i_{0}}(u)=\min_{1\leq i\leq n}\lambda_{i}(u)<+\infty$. Then,
there exist $\xi\in\mathbb{R}^{n}$ and $\sigma>0$ such that when
$\varepsilon\downarrow 0$, the following expansion is fulfilled:

${}E_{\varepsilon,u}=E_{i_{0}}+\varepsilon(\xi,\sigma)+o(\varepsilon).$ (4.1)

In addition, the vector $\xi$ and $\sigma$ are given by (with the convention
that $\xi_{0}=\xi_{n+1}=0$):

${}\left\\{\begin{array}[]{cll}\xi_{j}&=0&\mathrm{if}\;j\in\\{1,...,n\\}\backslash\\{i_{0}-1,i_{0},i_{0}+1\\},\vspace{0.1cm}\\\
\xi_{j}&=\frac{1-\lambda_{i_{0}}(u)}{u-\mu_{j}(\lambda_{i_{0}}(u))}&\mathrm{if}\;j\in\\{i_{0}-1,i_{0}+1\\}\backslash\\{0,n+1\\},\vspace{0.1cm}\\\
\xi_{i_{0}}&=-\xi_{i_{0}-1}-\xi_{i_{0}+1}-\sigma,\vspace{0.1cm}&\\\
\sigma&=-\frac{t_{i_{0},i_{0}}}{\mu^{\prime}_{i_{0}}(\lambda_{i_{0}}(u))}.\end{array}\right.$
(4.2)

###### Proof.

Since $0<u<m$, one has $m-2\varepsilon>u$ for $\varepsilon$ small enough,
thus, (3.6) implies that $u<\hat{u}(\varepsilon)\leq u_{c}(\varepsilon)$ so
that the steady-state $E_{\varepsilon,u}$ exists for every $\varepsilon>0$
small enough.

Now, for convenience, we write $E_{i_{0}}$ as $E_{i_{0}}=(x^{*},s^{*})$. We
start by proving that the mapping $\varepsilon\mapsto s^{\varepsilon,u}$ is of
class $C^{1}$ in some right neighborhood of $\varepsilon=0$. Doing so, let
$\theta>0$ and let us define the open set
$D_{\theta}:=(-\theta,1)\times(0,1)$. Consider the $C^{1}$ mapping
$\tilde{B}:D_{\theta}\rightarrow\mathbb{R}^{n\times n}$ given by
$\tilde{B}(\varepsilon,s):=B(s,u,\varepsilon)$ for $(\varepsilon,s)\in
D_{\theta}$ (here $u=\mu_{i_{0}}(s^{*})>0$ is fixed). Note that for every
$(\varepsilon,s)\in D_{\theta}$, the matrix $\tilde{B}(\varepsilon,s)$ is
symmetric quasi-positive and that for $\varepsilon=0$, zero is the largest and
simple eigenvalue of $\tilde{B}(0,s^{*})=M(s^{*})-uI_{n}$ (observe that
$\tilde{B}(0,s^{*})$ is diagonal with $n$ distinct eigenvalues). It follows
that $\lambda(\cdot)$ is analytic as a function of its $n^{2}$ coefficients in
some neighborhood of $\tilde{B}(0,s^{*})$ in the space of symmetric matrices.
Since $\tilde{B}(\cdot,\cdot)$ is of class $C^{1}$ w.r.t. $(\varepsilon,s)$,
there are $\theta^{\prime}>0$ and $\nu>0$ small enough such that the
composition

$(\varepsilon,s)\mapsto\varphi(\varepsilon,s):=\lambda(\tilde{B}(\varepsilon,s))$

is of class $C^{1}$ over
$(-\theta^{\prime},\theta^{\prime})\times(s^{*}-\nu,s^{*}+\nu)$. For
$\varepsilon=0$, the unitary eigenvector of $\tilde{B}(0,s^{*})$ for the zero
eigenvalue is the $i_{0}$-th vector $w=e_{i_{0}}$ of the canonical basis of
$\mathbb{R}^{n}$. It follows from (3.7) that $\frac{\partial\varphi}{\partial
s}(0,s^{*})=\mu^{\prime}_{i_{0}}(s^{*})>0$. So, we can apply the implicit
function theorem locally around $(0,s^{*})$. Hence, the mapping
$\varepsilon\mapsto s^{\varepsilon,u}$ is of class $C^{1}$ over
$(-\varepsilon_{0},\varepsilon_{0})$, and in particular in some right
neighborhood of $\varepsilon=0$. In addition, one has:

$\forall\varepsilon\in(-\varepsilon_{0},\varepsilon_{0}),\quad\varphi(\varepsilon,s^{\varepsilon,u})=0.$

By differentiating the preceding equality w.r.t. $\varepsilon$ and letting
$\varepsilon\downarrow 0$, we find

${}\sum_{1\leq i,j\leq
n}\frac{\partial\lambda(\tilde{B}(0,s^{*}))}{\partial\tilde{b}_{i,j}}\frac{\partial\tilde{b}_{i,j}(0,s^{*})}{\partial\varepsilon}+\mu^{\prime}_{i_{0}}(s^{*})\frac{ds^{\varepsilon,u}}{d\varepsilon}_{|_{\varepsilon=0}}=0,$
(4.3)

where $\tilde{b}_{i,j}(0,s^{*})$ denote the $n^{2}$ entries of
$\tilde{B}(0,s^{*})$. Combining (3.7) and (4.3), we obtain

$t_{i_{0},i_{0}}+\mu^{\prime}_{i_{0}}(s^{*})\frac{ds^{\varepsilon,u}}{d\varepsilon}_{|_{\varepsilon=0}}=0,$

which implies
$\frac{ds^{\varepsilon,u}}{d\varepsilon}_{|_{\varepsilon=0}}=\sigma$ and the
desired expansion of $s^{\varepsilon,u}$ up to the first order as in
(4.1)-(4.2).

Let us now turn to the expansion of $x^{\varepsilon,u}$ w.r.t. $\varepsilon$.
Doing so, let us consider the $C^{1}$ mapping
$\psi:(-\theta,\theta)\times\mathbb{R}^{n}\rightarrow\mathbb{R}^{n+1}$ defined
as

$\psi(\varepsilon,x):=\Big{(}\tilde{B}(\varepsilon,s^{\varepsilon,u})x,\sum_{j=1}^{n}x_{j}-1\Big{)},$

whose differential w.r.t. $x$ at $(0,x^{*})$ satisfies

$D_{x}\psi(0,x^{*})h=\Big{(}\tilde{B}(0,s^{*})h,\sum_{j=1}^{n}h_{j}\Big{)}.$

We can check that the kernel of $D_{x}\psi(0,x^{*})$ is reduced to $\\{0\\}$,
so, $D_{x}\psi(0,x^{*})$ is invertible. Hence, by the implicit function
theorem, we can conclude that $\varepsilon\mapsto x^{\varepsilon,u}$ is of
class $C^{1}$ in some neighborhood of $\varepsilon=0$. To obtain the desired
expansion of $x^{\varepsilon,u}$, let us write
$x^{\varepsilon,u}=x^{*}+\varepsilon d+o(\varepsilon)$. Since
$\sum_{j=1}^{n}x^{\varepsilon,u}_{j}+s^{\varepsilon,u}=1$, one has

${}\sigma+\sum_{j=1}^{n}d_{j}=0.$ (4.4)

Expanding $B(s^{\varepsilon,u},u,\varepsilon)$ w.r.t. $\varepsilon$ up to the
first order, we get:

$\displaystyle 0=B(s^{\varepsilon,u},u,\varepsilon)x^{\varepsilon,u}$
$\displaystyle=B(s^{*}+\sigma\varepsilon+o(\varepsilon),u,\varepsilon)x^{\varepsilon,u}$
$\displaystyle=(M(s^{*}+\sigma\varepsilon+o(\varepsilon))-uI_{n}+\varepsilon
T)x^{\varepsilon,u}$ $\displaystyle=(M(s^{*})-uI_{n}+\varepsilon T+\sigma
M^{\prime}(s^{*})\varepsilon+o(\varepsilon))(x^{*}+\varepsilon
d+o(\varepsilon))$ $\displaystyle=\big{(}Tx^{*}+\sigma
M^{\prime}(s^{*})x^{*}+(M(s^{*})-uI_{n})d\big{)}\varepsilon+o(\varepsilon),$

using the relation $(M(s^{*})-uI_{n})x^{*}=0$ in the last equality. Hence, we
deduce that

$Tx^{*}+\sigma M^{\prime}(s^{*})x^{*}+(M(s^{*})-uI_{n})d=0,$

which gives

$(Tx^{*})_{j}+\sigma\mu^{\prime}_{j}(s^{*})x^{*}+(\mu_{j}(s^{*})-u)d_{j},\quad
1\leq j\leq n.$

In the case where $1<i_{0}<n$, we obtain (4.2) for
$j\in\\{1,...,n\\}\backslash\\{i_{0}\\}$ from the preceding equation. For
$j=i_{0}$, (4.4) gives (4.2). A similar computation gives (4.2) whenever
$i_{0}=1$ or $i_{0}=n$, which concludes the proof. ∎

From Proposition 4.1, species of index $i_{0}$ is the only one with a positive
value at the zero order. Observe that it satisfies the inequality
$\xi_{i_{0}}<\lambda_{i_{0}}(u)$. In addition, only neighbors of $i_{0}$
(i.e., species with index $i_{0}-1$ or $i_{0}+1$) are significant up to the
first order. Species with index $j\notin\\{i_{0}-1,i_{0},i_{0}+1\\}$ are
(asymptotically) not significant w.r.t. species with index $i_{0}-1$ and
$i_{0}+1$. We now turn to the case where $\varepsilon$ tends to $+\infty$.

######  Proposition 4.2.

For every $u\in(0,\bar{\mu}(1))$, the point $E_{\varepsilon,u}$ has a limit
when $\varepsilon\rightarrow+\infty$ and

${}\lim_{\varepsilon\rightarrow+\infty}E_{\varepsilon,u}=\Big{(}\frac{1-\bar{\mu}^{-1}(u)}{n}a,\bar{\mu}^{-1}(u)\Big{)}.$
(4.5)

###### Proof.

Using the implicit function theorem locally around each $\varepsilon>0$, we
deduce that the derivative of $\varepsilon\mapsto s^{\varepsilon,u}$ w.r.t.
$\varepsilon$ exists and is non-negative (see the proof of Proposition 3.3).
Hence, $s^{\varepsilon,u}$ is non-increasing, and thus it admits a limit
$s^{u}$ as $\varepsilon$ goes to infinity because $s^{\varepsilon,u}\in[0,1]$,
for every $\varepsilon>0$. By definition of $s^{\varepsilon,u}$ we have:

$\displaystyle 0$
$\displaystyle=\lambda\left(B(s^{\varepsilon,u},u,\varepsilon)\right)=\varepsilon\lambda\left(T+\frac{M(s^{\varepsilon,u})}{\varepsilon}\right)-u$
$\displaystyle=\varepsilon\left[\lambda\left(T+\frac{M(s^{\varepsilon,u})}{\varepsilon}\right)-\lambda\left(T+\frac{M(s^{u})}{\varepsilon}\right)\right]+\varepsilon\lambda\left(T+\frac{M(s^{u})}{\varepsilon}\right)-u$
$\displaystyle=\varepsilon\left[\lambda\left(T+\frac{M(s^{\varepsilon,u})}{\varepsilon}\right)-\lambda\left(T+\frac{M(s^{u})}{\varepsilon}\right)\right]+\frac{1}{n}\sum\limits_{j=1}^{n}\mu_{j}(s^{u})-u+o(1),$

using the same expansion as in (3.8). As $\lambda(\cdot)$ is of class $C^{1}$
in some neighborhood of $T$, it is in particular locally Lipschitz, so, the
first term goes to 0 as $\varepsilon$ goes to infinity (using that $M(\cdot)$
is also of class $C^{1}$ and that $s^{\varepsilon,u}\rightarrow s^{u}$ as
$\varepsilon\rightarrow+\infty$). Hence, one must have $\bar{\mu}(s^{u})=u$,
that is $s^{u}=\bar{\mu}^{-1}(u)$. Let us now turn to the limit of
$x^{\varepsilon,u}$ as $\varepsilon\rightarrow+\infty$. From the proof of
Proposition 3.1, the vector $x^{\varepsilon,u}$ satisfies the system

${}M(s^{\varepsilon,u})x^{\varepsilon,u}-ux^{\varepsilon,u}+\varepsilon
Tx^{\varepsilon,u}=0,$ (4.6)

and it is proportional to $a^{\varepsilon,u}$:

${}x^{\varepsilon,u}=\frac{1-s^{\varepsilon,u}}{\sum_{j=1}^{n}a^{\varepsilon,u}_{j}}a^{\varepsilon,u}.$
(4.7)

Since for every $\varepsilon>0$, $\|a^{\varepsilon,u}\|=1$, there is
$\tilde{a}^{u}\in\mathbb{R}^{n}$ with $\|\tilde{a}^{u}\|=1$ such that, up to a
sub-sequence, one has $a^{\varepsilon,u}\rightarrow\tilde{a}^{u}$. By passing
to the limit in (4.6), we find that $T\tilde{a}^{u}=0$, thus
$\tilde{a}^{u}=\frac{a}{\sqrt{n}}$. Now, $\tilde{a}^{u}$ is also the limit of
every converging sub-sequence of $(a^{\varepsilon,u})_{\varepsilon}$, hence
$\tilde{a}^{u}$ is the limit of $(a^{\varepsilon,u})_{\varepsilon}$. Letting
$\varepsilon\rightarrow+\infty$ in (4.7) then gives (4.5), which ends the
proof. ∎

Even if the case $\varepsilon\rightarrow+\infty$ may have no meaning from an
application point of view, this result shows that species are asymptotically
uniformly distributed.

##  5 Persistence of all the species

In this section, we give an extension of [13, Theorem 3] showing that each
species (individually) is persistent. We refer to [40, 44, 45] for the
mathematical theory of persistence. The persistence result in [13] is related
to the total biomass (i.e., the sum of the concentrations of the species). In
our setting, it can be stated as follows.

######  Theorem 5.1 ([13]).

There is $c>0$ such that for every
$(\varepsilon,u)\in\mathbb{R}_{+}^{*}\times(0,u_{c}(\varepsilon))$ and every
initial condition in the set $\mathcal{D}$, the unique solution of (2.1)
associated with this initial condition satisfies

${}\liminf_{t\rightarrow+\infty}\;\sum_{j=1}^{n}x_{j}(t)\geq\beta_{\varepsilon,u}:=c\frac{\min_{1\leq
j\leq n}v_{j}^{\varepsilon,u}}{\max_{1\leq j\leq n}v_{j}^{\varepsilon,u}},$
(5.1)

where $v^{\varepsilon,u}$ is the Perron vector associated with the matrix
$B(1,u,\varepsilon)$.

######  Remark 5.1.

It follows that for every
$(\varepsilon,u)\in\mathbb{R}_{+}^{*}\times(0,u_{c}(\varepsilon))$, one has
(recall (2.11)):

${}c_{u}\leq\liminf_{t\rightarrow+\infty}s(t)\leq\limsup_{t\rightarrow+\infty}s(t)\leq
1-\beta_{\varepsilon,u}<1,$ (5.2)

for every solution to (2.1) starting in $\mathcal{D}$. Section 6.1
investigates the particular case where $\varepsilon>0$ and $u=0$.

Before proving that each species is uniformly persistent, let us recall some
definitions of [10, 20] about the notion of persistence. Hereafter, the
interior, resp. the boundary of a set $A\subset\mathbb{R}^{n}$ is denoted by
$\mathrm{Int}(A)$, resp. $\partial A$, $B(x,r)$ denotes the open ball of
center $x\in\mathbb{R}^{n}$ and radius $r>0$. Finally, for every $r>0$, we
define $S(A,r):=\\{x\in\mathbb{R}^{n}\;;\;d(x,A)\leq r\\}$ where $d$ is a
distance over $\mathbb{R}^{n}$ and $d(x,A):=\inf_{a\in A}d(x,a)$. Consider now
a differential equation $\dot{x}=f(x)$ where
$f:\mathbb{R}^{n}\rightarrow\mathbb{R}^{n}$ is smooth and such that every
solution to this equation is global. Let us denote by $\mathcal{F}$ the
associated flow.

######  Definition 5.1.

Given two non-empty subsets $\mathcal{Y},\mathcal{Z}\subset\mathbb{R}^{n}$,
the sets $W^{\pm}(\mathcal{Y})$ stand respectively for

$W^{+}(\mathcal{Y}):=\\{x\in\mathcal{Z}\;;\;\omega(x)\subset\mathcal{Y}\\}\;;\;W^{-}(\mathcal{Y}):=\\{x\in\mathcal{Z}\;;\;\alpha(x)\subset\mathcal{Y}\\},$

where $\omega(x)$ and $\alpha(x)$ denote respectively the $\omega$-limit and
$\alpha$-limit sets of some point $x\in\mathbb{R}^{n}$ for the flow
$\mathcal{F}$.

######  Definition 5.2.

Let $E$ be a non-empty closed subset of $\mathbb{R}^{n}$ that is forward
invariant by $\mathcal{F}$. We say that $\mathcal{F}$ is uniformly persistent
related to $E$ if there is $\kappa>0$ such that for every initial condition in
$\mathrm{Int}(E)$, the corresponding solution $x(\cdot)$ satisfies

${}\liminf_{t\rightarrow+\infty}d(x(t),\partial E)>\kappa,$ (5.3)

where $d$ is a distance over $\mathbb{R}^{n}$

In the next theorem, we show that each species is uniformly persistent. The
proof is based on (5.1).

######  Theorem 5.2.

For every $(\varepsilon,u)\in\mathbb{R}_{+}^{*}\times(0,u_{c}(\varepsilon))$,
there exists $\gamma_{\varepsilon,u}>0$ such that for every initial condition
in the set $\mathcal{D}$, the unique solution of (2.1) associated with this
initial condition satisfies

${}\liminf_{t\rightarrow+\infty}x_{i}(t)\geq\gamma_{\varepsilon,u},$ (5.4)

for every $1\leq i\leq n$.

###### Proof.

Fix $\varepsilon>0$, $u\in(0,u_{c}(\varepsilon))$, and consider the sets
$E_{\eta}$ given by

$E_{\eta}:=\Big{\\{}(x,s)\in\mathbb{R}_{+}^{n}\times[0,1]\;;\;s+\sum_{j=1}^{n}x_{j}\leq
1+\eta\Big{\\}},$

where $\eta>0$. Obviously, $E_{\eta}$ is a closed subset of
$\mathbb{R}_{+}^{n}\times[0,1]$ that is positively invariant by (2.1). In
addition, it is easily seen that its boundary satisfies

$\partial
E_{\eta}=\Big{\\{}(x,s)\in\mathbb{R}_{+}^{n}\times[0,1]\;;\;\exists\,i\in\\{1,...,n\\},\;x_{i}=0\;\;\mathrm{or}\;\;s=0\;\;\mathrm{or}\;\;\sum_{j=1}^{n}x_{j}+s=1+\eta\Big{\\}}.$

Let $(x^{0},s^{0})\in\mathcal{D}$ be an initial condition and let us denote by
$(x(\cdot),s(\cdot))$ the corresponding solution of (2.1). From Property 2.1,
one has $x_{i}(t)>0$ for every $t>0$ and $1\leq i\leq n$. In addition,
$b(t)\rightarrow 1$ as $t\rightarrow+\infty$ and $s(\cdot)$ cannot approach
$0$ because of (5.2). We deduce in particular that $S(\partial
E_{\eta},1)\cap\mathrm{Int}(E_{\eta})$ is point dissipative (see [10, 20]),
which means the following:

$\forall(x^{0},s^{0})\in S(\partial
E_{\eta},1)\cap\mathrm{Int}(E_{\eta}),\;\forall
t>0,\;(x(t),s(t))\in\mathrm{Int}(E_{\eta}).$

We now show that the maximal invariant subset $N$ of $\partial E_{\eta}$ by
(2.1) is acyclic555This property amounts to verify that $N^{c}\cap
W^{-}(N)\cap W^{+}(N)=\varnothing$ where $N^{c}$ is the complement of $N$ in
$E_{\eta}$, see [10, 20] or [26, 45] for a more detailed definition. and
isolated. First, observe that $N=\\{0_{\mathbb{R}^{n}}\\}\times[0,1+\eta]$
using Property 2.1. Considering now the distance $d$ over
$\mathbb{R}_{+}\times[0,1]$ defined as

$d((x,s),(x^{\prime},s^{\prime})):=\sum_{j=1}^{n}|x_{j}-x^{\prime}_{j}|+|s-s^{\prime}|,$

for $(x,s),(x^{\prime},s^{\prime})\in\mathbb{R}_{+}\times[0,1]$, one has using
(5.1)

$\liminf_{t\rightarrow+\infty}d((x(t),s(t)),N)=\liminf_{t\rightarrow+\infty}\sum_{j=1}^{n}x_{j}(t)\geq\beta_{\varepsilon,u}>0,$

for every solution of (2.1) starting in $\mathcal{D}$. If now $N$ and
$\partial E_{\eta}$ stand respectively for $\mathcal{Y}$ and $\mathcal{Z}$ in
Definition 5.1, the previous inequality implies that $W^{+}(N)=N$. Hence $N$
is necessarily acyclic.

Finally, the set $N$ is isolated because for every initial condition
$(x^{0},s^{0})\in\mathcal{D}\backslash N$, (5.1) implies the existence of
$t_{0}\geq 0$ such that

$\forall t\geq
t_{0},\;d((x(t),s(t)),N)=\sum_{j=1}^{n}x_{j}(t)\geq\frac{\beta_{\varepsilon,u}}{2}>0.$

We are now in a position to use [20, Theorem 4.3] which asserts that the flow
defined by (2.1) is uniformly persistent related to the set $E_{\eta}$
provided that there is $\delta>0$ such that

${}W^{+}(N)\cap S(\partial
E_{\eta},\delta)\cap\mathrm{Int}(E_{\eta})=\varnothing.$ (5.5)

But, (5.5) is clearly verified with $\delta:=1$ because $W^{+}(N)\cap
S(\partial E_{\eta},\delta)\subset N\subset\partial E_{\eta}$, so we have
proved that for every $\eta>0$, the flow defined by (2.1) is uniformly
persistent related to the set $E_{\eta}$. To conclude the proof, fix $\eta>0$
and apply (5.3) with $E_{\eta}$ in place of $E$. Note that Property 2.1 and
Lemma 2.1 imply that every solution is necessarily with values in
$\mathrm{Int}(E_{\eta})$ over $\mathbb{R}_{+}^{*}$. Hence, we deduce that
there exists $\kappa>0$ such that

$\liminf_{t\rightarrow+\infty}d((x(t),s(t)),\partial E_{\eta})\geq\kappa,$

for every solution starting in $\mathcal{D}$. In view of the definition of
$\partial E_{\eta}$, we can write $\partial
E_{\eta}=\bigcup_{i=1}^{n}F_{i}\cup\tilde{F}$ where
$F_{i}:=\\{(x,s)\in\mathbb{R}_{+}^{n}\times[0,1]\;;\;x_{i}=0\\}$ and
$\tilde{F}$ is the complement of $\bigcup_{i=1}^{n}F_{i}$ in $E_{\eta}$. It
follows that for every $1\leq i\leq n$ and every initial condition in
$\mathrm{Int}(E_{\eta})$, one has

${}\kappa\leq\liminf_{t\rightarrow+\infty}d((x(t),s(t)),\partial
E_{\eta})\leq\liminf_{t\rightarrow+\infty}d((x(t),s(t)),F_{i})=\liminf_{t\rightarrow+\infty}x_{i}(t).$
(5.6)

Finally, for every initial condition $(x^{0},s^{0})\in\mathcal{D}$, there is a
time $t^{\prime}_{0}\geq 0$ such that for every time $t\geq t^{\prime}_{0}$,
the associated solution to (2.1) satisfies
$(x(t),s(t))\in\mathrm{Int}(E_{\eta})$. Combining this property with (5.6)
then yields the desired property (5.4) with $\gamma_{\varepsilon,u}:=\kappa$.
∎

##  6 Global stability property of (2.1)

###  6.1 Asymptotic behavior of (2.1) with $u=0$

We start by studying (2.1) in batch mode, i.e., we take $u=0$. This will be
useful to prove Theorem 6.2. The dynamics of $x$ then becomes

${}\dot{x}=M(s(t))x+\varepsilon Tx.$ (6.1)

If $\varepsilon=0$, the solution $(x(\cdot),s(\cdot))$ of (2.3) converges to
some point $(x^{\infty},0)\in\Delta$ such that
$\sum_{j=1}^{n}x^{\infty}_{j}=b(0)$. So, we suppose in what follows that
$\varepsilon>0$.

######  Proposition 6.1.

If $u=0$ and $\varepsilon>0$, every solution $(x(\cdot),s(\cdot))$ of (2.1)
starting in $\mathcal{D}$ satisfies

$\lim_{t\rightarrow+\infty}(x(t),s(t))=\Big{(}\frac{b(0)}{n}a,0\Big{)}.$

###### Proof.

Observe that the mapping $t\mapsto s(t)$ decreases over $\mathbb{R}_{+}$, and
that $s\geq 0$. Thus, $s(\cdot)$ necessarily converges to some value
$\bar{s}$. By Barbalat’s Lemma, $\lim_{t\rightarrow+\infty}\dot{s}(t)$ exists
and is zero. Suppose now by contradiction that $\bar{s}>0$. It follows that
$\sum_{j=1}^{n}x_{j}(t)\rightarrow b(0)-\bar{s}$ and that $b(0)-\bar{s}>0$
because $x(0)\in\mathcal{D}$. Hence, there is $t_{0}\geq 0$ such that for
every $t\geq t_{0}$, one has

$\dot{s}(t)=-\sum_{j=1}^{n}\mu_{j}(s(t))x_{j}(t)\leq-\upsilon\sum_{j=1}^{n}x_{j}(t)\leq-\upsilon\frac{b(0)-\bar{s}}{2},$

where $\upsilon:=\frac{1}{2}\min_{1\leq j\leq n}(\mu_{j}(\bar{s}))$. We have
thus obtained a contradiction with the fact that $\dot{s}(t)\rightarrow 0$ as
$t\rightarrow+\infty$. Let us now come back to (6.1) which is a non-autonomous
perturbation of the linear system

${}\dot{x}=\varepsilon Tx.$ (6.2)

In order to apply the theory of asymptotically autonomous system [43], we need
to rewrite (6.2) in the orthogonal of $\mathbb{R}a$ in such a way that the
corresponding autonomous dynamics possesses a unique globally asymptotically
stable equilibrium (this is not the case with (6.2) since zero is an
eigenvalue of $T$). Doing so, we know that there exists an invertible matrix
$P\in\mathbb{R}^{n\times n}$ such that $P^{-1}TP=D$ where
$D:=\mathrm{diag}(0,\alpha_{2},...,\alpha_{n})$ with $\alpha_{i}<0$ for $2\leq
i\leq n$. In addition, without any loss of generality, we may assume that the
first column of $P$ is exactly equal to the vector $a$ (that is collinear to
the Perron vector of $T$), and we also set
$\tilde{D}:=\mathrm{diag}(\alpha_{2},...,\alpha_{n})\in\mathbb{R}^{(n-1)\times(n-1)}$.
Multiplying (6.1) on the left by $P^{-1}$ then gives666Given
$v\in\mathbb{R}^{n}$, the notation $v_{-1}$ indicates that $v=(v_{1},v_{-1})$
with $v_{-1}\in\mathbb{R}^{n-1}$ ($v_{-1}$ is the vector obtained from $v$ by
removing the first component).

${}\left|\begin{array}[]{cl}\dot{y}&=(P^{-1}M(s(t))x)_{1},\\\
\dot{z}&=(P^{-1}M(s(t))x)_{-1}+\varepsilon\tilde{D}z,\end{array}\right.$ (6.3)

where $x=P\left(\begin{array}[]{c}y\\\ z\end{array}\right)$. Next, the ODE
satisfied by $z$ can be rewritten

$\dot{z}=F(t,z)+\varepsilon\tilde{D}z,$

where $F:\mathbb{R}_{+}\times\mathbb{R}^{n-1}\rightarrow\mathbb{R}^{n-1}$ is
defined by

$F(t,z):=\left(P^{-1}M(s(t))P\left(\begin{array}[]{c}y(t)\\\
z\end{array}\right)\right)_{-1}.$

Since $s(t)\rightarrow 0$ when $t\rightarrow+\infty$ and $y(\cdot)$ is
bounded, the preceding system is a non-autonomous perturbation of the linear
system

${}\dot{z}=\varepsilon\tilde{D}z.$ (6.4)

Now, one has $F(t,z)\rightarrow 0$ when $t\rightarrow+\infty$ uniformly
locally w.r.t. $z$ and observe that every solution to (6.4) converges to zero.
We deduce from [43] that every solution $(y(\cdot),z(\cdot))$ to (6.3) is such
that $z(t)\rightarrow 0$ when $t\rightarrow+\infty$. Coming back to the
original variable $x$, the solution $x(\cdot)$ can be written

$x(t)=y(t)a+o(1).$

To conclude, observe that $t\mapsto b(t)$ is constant. Thus, for every $t\geq
0$,

$b(0)=b(t)=\sum_{j=1}^{n}x_{j}(t)+s(t)=ny(t)+o(1).$

Hence, $y(t)\rightarrow b(0)/n$ as $t\rightarrow+\infty$, which ends the
proof. ∎

######  Remark 6.1.

This proposition shows that if $\varepsilon>0$ and $u=0$, then, every species
concentration converges to the same value $\frac{b(0)}{n}$ as
$t\rightarrow+\infty$. In that case, any solution to (2.1) converges to the
point $\Big{(}\frac{b(0)}{n}a,0\Big{)}$ that depends on the initial condition.

###  6.2 Global stability for $\varepsilon>0$ and $u$ small enough

Let us first recall Corollary 2.3 of [42] which is a fundamental result about
global stability of a perturbed steady-state. Let $k\geq 1$, and $G$, $U$ two
closed subsets of $\mathbb{R}^{n}$ and $\mathbb{R}^{k}$ respectively. Consider
a continuous function $g:G\times U\rightarrow\mathbb{R}^{n}$, $(x,u)\mapsto
g(x,u)$ where $u\in U$ is a parameter. Suppose that $D_{x}g(x,u)$ exists and
is continuous over $G\times U$ and that solutions $x(\cdot,x_{0},u)$ to the
Cauchy Problem

${}\left|\begin{array}[]{cl}\dot{x}&=g(x,u),\vspace{0.1cm}\\\
x(0)&=x^{0},\end{array}\right.$ (6.5)

are unique and remain in $G$ for every time $t\geq 0$ and every $(x^{0},u)\in
G\times U$.

######  Theorem 6.1 ([42]).

Let $(x^{\star},u^{\star})\in G\times U$ be such that
$x^{\star}\in\mathrm{Int}(G)$ and $g(x^{\star},u^{\star})=0$. Suppose that the
matrix $D_{x}g(x^{\star},u^{\star})$ is Hurwitz and that $x^{\star}$ is
globally attracting for solutions to (6.5) with $u=u^{\star}$. If there is a
non-empty compact set $K\subset G$ such that for each $(x^{0},u)\in G\times
U$, $x(t,x^{0},u)\in K$ for $t$ large enough, then, there are $r>0$ and a
unique point $x^{\star}(u)\in G$ for every $u\in B(u^{\star},r)$ such that
$g(x^{\star}(u),u)=0$ and:

${}\forall u\in B(u^{\star},r),\;\forall x^{0}\in
G,\;x(t,x^{0},u)\xrightarrow[t\rightarrow+\infty]{}x^{\star}(u).$ (6.6)

######  Remark 6.2.

This result also applies if $x^{\star}$ is on the boundary of $G$ provided
that the dynamics $g$ can be extended to a $C^{1}$ mapping in some convex
neighborhood of $x^{\star}$ (see [42, Corollary 2.3]).

The next lemma is based on Proposition 5.1 (see [13]) and it will be useful to
prove Theorem 6.2.

######  Lemma 6.1.

For every
$(\varepsilon,u_{1})\in\mathbb{R}_{+}^{*}\times(0,u_{c}(\varepsilon))$, there
is $\beta_{\varepsilon}>0$ such that for every $u\in[0,u_{1}]$ and every
initial condition in $\mathcal{D}$, the unique solution of (2.1) associated
with this initial condition satisfies

$\liminf_{t\rightarrow+\infty}\sum_{j=1}^{n}x_{j}(t)\geq\beta_{\varepsilon}.$

###### Proof.

Let $(\varepsilon,u_{1})\in\mathbb{R}_{+}^{*}\times(0,u_{c}(\varepsilon))$.
Given $u\in[0,u_{1}]$, the Perron vector $v^{\varepsilon,u}$ associated with
the greatest eigenvalue of the matrix $B(1,u,\varepsilon)$ is the unique
solution to the system

$\left\\{\begin{array}[]{rl}B(1,u,\varepsilon)w-\lambda(B(1,u,\varepsilon))w&=0,\\\
\|w\|-1&=0.\end{array}\right.$

Now, consider the $C^{1}$-mapping

$\tilde{\psi}:(u,w)\in[0,u_{1}]\times\mathbb{R}^{n}\mapsto\tilde{\psi}(u,w):=\Big{(}B(1,u,\varepsilon)w-\lambda(B(1,u,\varepsilon))w,\frac{\|w\|^{2}-1}{2}\Big{)}\in\mathbb{R}^{n+1}.$

Its partial derivative w.r.t. $w$ at the point
$(u,v^{\varepsilon,u})\in[0,u_{1}]\times\mathbb{R}^{n}$ is given by:

$\frac{\partial\tilde{\psi}}{\partial
w}(u,v^{\varepsilon,u})=\Big{(}B(1,u,\varepsilon)-\lambda(B(1,u,\varepsilon))I_{n},v^{\varepsilon,u}\Big{)}.$

Hence, if $w$ is in the kernel of $\frac{\partial\tilde{\psi}}{\partial
w}(u,v^{\varepsilon,u})$, it satisfies
$B(1,u,\varepsilon)w=\lambda(B(1,u,\varepsilon))w$ and $v^{\varepsilon,u}\cdot
w=0$ (here, $\cdot$ is the scalar product in $\mathbb{R}^{n}$). The first
equality implies that there is $\nu\in\mathbb{R}$ such that $w=\nu
v^{\varepsilon,u}$. Using the second equality, we find that $\nu=0$, thus
$w=0$ and $\frac{\partial\tilde{\psi}}{\partial w}(u,v^{\varepsilon,u})$ is
invertible. Thanks to the implicit function theorem, we obtain that way that
$u\mapsto v^{\varepsilon,u}$ is locally continuous around every
$u\in[0,u_{1}]$, thus it is continuous over $[0,u_{1}]$. Now, Proposition 5.1
of [13] implies that

$\liminf_{t\rightarrow+\infty}\sum_{j=1}^{n}x_{j}(t)\geq c\min_{1\leq j\leq
n}v_{j}^{u,\varepsilon}.$

Since the mapping $u\mapsto v^{\varepsilon,u}$ is continuous over $[0,u_{1}]$,
so is $u\mapsto\min_{1\leq j\leq n}v_{j}^{u,\varepsilon}$, hence,

$\liminf_{t\rightarrow+\infty}\sum_{j=1}^{n}x_{j}(t)\geq\beta_{\varepsilon}:=c\min_{u\in[0,u_{1}]}\min_{1\leq
j\leq n}v^{u,\varepsilon}_{j}.$

Because $u\mapsto\min_{1\leq j\leq n}v_{j}^{u,\varepsilon}$ is positive and
continuous over $[0,u_{1}]$, we get that $\beta_{\varepsilon}>0$ which ends
the proof. ∎

We now give our main result about the global stability of the steady-state
$E_{\varepsilon,u}$ when $\varepsilon>0$ is fixed and $u$ is with small enough
values.

######  Theorem 6.2.

For every $\varepsilon>0$, there is
$u_{s}(\varepsilon)\in(0,u_{c}(\varepsilon)]$ such that for every
$u\in(0,u_{s}(\varepsilon))$, the steady-state $E_{\varepsilon,u}$ is globally
asymptotically stable.

###### Proof.

First, Proposition 3.2 implies that for every $u>0$ such that
$u<u_{c}(\varepsilon)$, then one has $\lambda(B(1,u,\varepsilon))>0$.
Therefore, for every $u\in(0,u_{c}(\varepsilon))$, the point
$E_{\varepsilon,u}$ is the unique locally asymptotically stable point of (2.1)
in $\mathcal{D}$. We start by proving the result for those initial conditions
that are in the set $\Delta$. Fix $u_{1}\in(0,u_{c}(\varepsilon))$. The
dynamics of $x$ can be then written $\dot{x}=g(x,u)$ where $g:G\times
U\rightarrow\mathbb{R}^{n}$ is defined as

$g(x,u):=B\Big{(}1-\sum_{j=1}^{n}x_{j},u,\varepsilon\Big{)}x,$

with $G:=\mathcal{D}^{\prime}$ (recall (2.10)) and $U:=[0,u_{1}]$. We set
$x^{\star}:=\frac{a}{n}$, $u^{\star}=0$. We are then in a position to verify
the hypotheses of Theorem 6.1:

  * $\bullet$

At $(x^{\star},u^{\star})$, one has $g(x^{\star},u^{\star})=0$ and
$x^{\star}\in\partial G$ since $\sum_{j=1}^{n}x^{\star}_{j}=1$ ;

  * $\bullet$

The Jacobian matrix of $g$ w.r.t. $x$ at $(x^{\star},u^{\star})$ writes

$D_{x}g(x^{\star},u^{\star})=B\Big{(}1-\sum_{j=1}^{n}x_{j}^{\star},u^{\star},\varepsilon\Big{)}-d_{0}\,a^{\top}=\varepsilon
T-d_{0}\,a^{\top},$

where $d_{0}:=M^{\prime}(0)x^{\star}$ (recall the proof of Proposition 3.1
(ii)). It is a rank one perturbation of $\varepsilon T$. By using a similar
argumentation as in the proof of Proposition 3.1 (ii), we can check that it is
Hurwitz ;

  * $\bullet$

By extending $\mu_{i}$ as a $C^{1}$ function over $\mathbb{R}$, the dynamics
$g$ can be extended to a $C^{1}$ function in $(G\cap
B(x^{\star},\delta))\times U$ for every $\delta>0$ ;

  * •

The set $K:=\\{x\in G\;;\;\beta_{\varepsilon}\leq\sum_{j=1}^{n}x_{j}\leq 1\\}$
is compact, and according to Lemma 6.1, for every $u\in[0,u_{1}]$ and for $t$
large enough, one has $x(t)\in K$ for every solution to $\dot{x}=g(x,u)$.

We can then apply Theorem 6.1 which implies the existence of
$u_{0}\in(0,u_{1}]$ such that for every $u\in[0,u_{0}]$, the point
$x^{\varepsilon,u}$ is GAS for the dynamics restricted to $\Delta$.

We now consider initial conditions in $\mathcal{D}$ and let $u\in(0,u_{0}]$ be
fixed. The first $n$ equations in (2.1) write

$\dot{x}=B\Big{(}b(t)-\sum_{j=1}^{n}x_{j},u,\varepsilon\Big{)}x.$

This system is a non-autonomous perturbation of the autonomous system
$\dot{x}=g(x,u)$ since $b(t)\rightarrow 1$ when $t\rightarrow+\infty$. Using a
similar argumentation as in the proof of Proposition 6.1 (see [43]), we deduce
that for every initial condition $x^{0}\in[0,+\infty)^{n}\backslash\\{0\\}$,
one has $x(t)\rightarrow x^{\varepsilon,u}$ as $t\rightarrow+\infty$. Now,
given some initial condition $(x^{0},s^{0})\in\mathcal{D}$ for system (2.1),
one has $b(t)\rightarrow 1$ as $t\rightarrow+\infty$. So, one has
$s(t)\rightarrow s^{\varepsilon,u}$ as $t\rightarrow+\infty$ which shows that
for every $u\in[0,u_{0}]$, then $E_{\varepsilon,u}$ is GAS.

We now argue that the same reasoning can be employed starting from the point
$x^{\varepsilon,u_{0}}$ (which is GAS) in place of $x^{\star}$. We obtain that
way the existence of $u^{\prime}_{0}>u_{0}$ such that $E_{\varepsilon,u}$ is
GAS for every $u\in[u_{0},u^{\prime}_{0}]$. Repeating this argumentation, one
can define

$u_{s}(\varepsilon):=\sup\\{u_{0}\in(0,u_{c}(\varepsilon))\;;\;\forall
u\in[0,u_{0}],\;E_{\varepsilon,u}\;\mathrm{is}\;\mathrm{GAS}\\}.$

This concludes the proof. ∎

Showing that $E_{\varepsilon,u}$ is GAS for every
$(\varepsilon,u)\in\mathbb{R}_{+}^{*}\times(0,u_{c}(\varepsilon))$ seems a
difficult question that could deserve further investigations based on results
of Section 5 (Theorem 5.2). We can make the following observations:

  * $\bullet$

If $u_{s}(\varepsilon)=u_{c}(\varepsilon)$, then we have the desired property.
But, at this step, we only know that $u_{s}(\varepsilon)\leq
u_{c}(\varepsilon)$. If $u_{s}(\varepsilon)<u_{c}(\varepsilon)$, note that
$E_{\varepsilon,u}$ remains LAS for every
$u\in[u_{s}(\varepsilon),u_{c}(\varepsilon))$, i.e., no bifurcation occurs at
$u=u_{s}(\varepsilon)$. So, one can wonder if in this setting, such a loss of
global stability is possible or not.

  * $\bullet$

Another approach consists in showing that $E_{\varepsilon,u}$ is GAS for every
$u>0$ provided that $\varepsilon>0$ is small enough using a similar result as
in Lemma 6.1, and proceeding as in the proof of Theorem 6.2. One should prove
that for every $u\in(0,\bar{\mu}(1))$, there is a constant
$\gamma^{\prime}_{u}>0$ (that does not depend on $\varepsilon$) such that

$\liminf_{t\rightarrow+\infty}x_{i_{0}}(t)\geq\gamma^{\prime}_{u},$

for every $\varepsilon$ small enough and every solution of (2.1), where
species $i_{0}$ wins the competition in absence of mutation.

In the next table, we summarize asymptotic properties about (2.1) that have
been established in this paper (including also the case without mutation,
under the hypotheses of Theorem 2.1777As in Theorem 2.1, we do not mention
here the (non-generic) cases where the dilution rate $u$ would be such that
$u=\mu_{i}(s)=\mu_{j}(s)$ for some indexes $i\not=j$ and $s\in(0,1)$. ).

$u$ $\varepsilon$ | $\varepsilon=0$ | $\varepsilon>0$
---|---|---
$u=0$ | Convergence into $\hat{\Delta}$ | Convergence to $\Big{(}\frac{b(0)}{n}a,0\Big{)}$
$0<u<u_{s}(\varepsilon)$ |  $E_{i_{0}}$ GAS in $\mathcal{D}_{i_{0}}$ |  $E_{\varepsilon,u}$ GAS in $\mathcal{D}$
$u_{s}(\varepsilon)\leq u<u_{c}(\varepsilon)$ |  $E_{i_{0}}$ GAS in $\mathcal{D}_{i_{0}}$ |  $E_{\varepsilon,u}$ LAS
$u=u_{c}(\varepsilon)$ |  $E_{wo}$ S in $\mathbb{R}_{+}^{n}\times[0,1]$ |  $E_{wo}$ S in $\mathbb{R}_{+}^{n}\times[0,1]$
$u>u_{c}(\varepsilon)$ |  $E_{wo}$ GAS in $\mathbb{R}_{+}^{n}\times[0,1]$ |  $E_{wo}$ GAS in $\mathbb{R}_{+}^{n}\times[0,1]$
Table 1: Summary of asymptotic properties of (2.1). Here,
$\hat{\Delta}:=\\{(x,s)\in\Delta\;;\;s=0\\}$,
$\mathcal{D}_{i_{0}}:=\mathcal{E}_{i_{0}}\times[0,1]$,
$\mathcal{E}_{i_{0}}:=\\{x\in\mathbb{R}_{+}^{n}\;;\;x_{i_{0}}\not=0\\}$ and
$\mathcal{D}=(\mathbb{R}_{+}^{n}\backslash\\{0\\})\times[0,1]$. The
abbreviations $S$, LAS, and GAS stand respectively for stable, locally
asymptotically stable, globally asymptotically stable.

##  7 Conclusion and perspectives

In this paper, we could extend some results of [13] showing that the
coexistence steady-state of (2.1) is always LAS and in particular GAS provided
that the dilution rate is small enough (assuming only that kinetics are of
Monod type). Let us emphasize that in contrast with the chemostat system,
mutation implies coexistence, i.e., each species is present asymptotically.
Future works could investigate global stability via a Lyapunov approach at
least for $\varepsilon>0$ small enough taking into account the knowledge of a
Lyapunov function for $\varepsilon=0$. Asymptotic stability properties could
be also addressed with more complicated mutation terms such as in [1, 29]. As
well, most properties proved in this paper are still valid if the kinetics are
only increasing, hence, one can wonder if such properties remain valid with
more sophisticated growth functions such as Haldane’s kinetics. Finally, it
could be also interesting to study continuous models describing the growth of
a population structured by a phenotypical trait living in a limited substrate
environment (see [36]).

## Acknowledgment

This research benefited from the support of Avignon Université (AAP
Agro&Sciences) and from the support of the FMJH Program PGMO and from the
support to this program from EDF-THALES-ORANGE. The authors would also like to
thank Francis Mairet, Pedro Gajardo, and Frédéric Mazenc for helpful
discussions about Lyapunov functions. The authors are grateful to P. De
Leenheer and A. Rapaport for fruitful exchanges.

## References

  * [1] S.S. Arkin, Microbial evolution in the chemostat, PhD Thesis, 2010, http://hdl.handle.net/10044/1/11305
  * [2] Z.S. Athanassov, Perturbation Theorems for Nonlinear Systems of Ordinary Differential Equations, J. Math. Anal. Appl., vol. 86, pp. 194–207, 1982.
  * [3] G. Bastin, D. Dochain, On-line estimation and adaptive control of bioreactors, Elsevier, New York, 1990.
  * [4] T. Bayen, F. Mairet, P. martinon, M. Sebbah, Analysis of a periodic optimal control problem connected to microalgae anaerobic digestion, Optimal Control Appl. Methods, vol. 36, 6, pp. 750–773, 2015.
  * [5] T. Bayen, F. Mairet, Optimization of the separation of two species in a chemostat, Automatica J. IFAC, Vol. 50, 4, pp. 1243–1248, 2014.
  * [6] T. Bayen, F. Mairet, Optimization of strain selection in evolution experiments in chemostat, Internat. J. Control, vol. 90, 12 , pp. 2748–2759, 2017.
  * [7] T. Bayen, J. Harmand, M. Sebbah, Time-optimal control of concentration changes in the chemostat with one single species, Appl. Math. Model., vol. 50, pp. 257–278, 2017.
  * [8] A. Berman, R. Plemmons, Nonnegative Matrices in the Mathematical Sciences, SIAM, Philadelphia, PA, 1994.
  * [9] J. Bierkens, A. Ran, A singular M-matrix perturbed by a nonnegative rank one matrix has positive principal minors; is it D-stable?, Linear Algebra Appl, vol. 457, pp. 191–208, 2014.
  * [10] G. Buttler, P. Waltman, Persistence in Dynamical Systems, J. Differential Equations, 63, pp. 255–263, 1986.
  * [11] J. Coville, Convergence to equilibrium for positive solutions of some mutation-selection model, preprint (2013). Available at arXiv:1308.6471.
  * [12] J. Coville, F. Fabre, Convergence to the equilibrium in a Lotka-Volterra ODE competition system with mutations, preprint (2013). Available at arXiv:1301.6237.
  * [13] P. De Leenheer, J. Dockery, T. Gedeon, S. Pilyugin, The chemostat with lateral gene transfer, J. Biol. Dyn., vol. 4, 6, pp. 607–620, 2010.
  * [14] P. De Leenheer, S.S. Pilyugin, Multistrain virus dynamics with mutations: A global analysis, Math. Med. Biol., vol. 25, 4, pp. 285–322, 2008.
  * [15] E. Deutsch, M. Neumann, Derivatives of the Perron Root at an Essentially Nonnegative Matrix and the Group Inverse of an AA-Matrix, J. Math. Anal. Appl., 102, pp. 1–29, 1984.
  * [16] O. Diekmann, A beginners guide to adaptive dynamics, Banach Center Publ., vol. 63, pp. 47–86, 2004.
  * [17] O. Diekmann, P.-E. Jabin, S. Mischler, B. Perthame, The dynamics of adaptation: an illuminating example and a Hamilton–Jacobi approach, Theoretical population biology, vol. 67, pp. 257–271, 2005.
  * [18] D. Dochain, P. Vanrolleghem, Dynamical modelling and estimation in wastewater treatment processes, IWA Publishing, vol. 4, London, 2001.
  * [19] R. Fekih-Salem, J. Harmand, C. Lobry, A. Rapaport, T. Sari., Extensions of the chemostat model with floculation. J. Math. Anal. Appl. 397, vol. 1, pp. 292–305, 2013.
  * [20] H. L Freedman, S. Ruan, M. Tang, Uniform Persistence and Flows Near a Closed Positively Invariant Set, J. Dynam. Differential Equations, vol. 6, 4, 1994.
  * [21] C. Fritsch, F. Campillo, O. Ovaskainen, A numerical approach to determine mutant invasion fitness and evolutionary singular strategies, Theoretical Population Biology, vol. 115, pp. 89–99, 2017.
  * [22] P. Gajardo, F. Mazenc, H. Ramirez, Competitive exclusion principle in a model of chemostat with delays, Dyn. Contin. Discrete Impuls. Syst. Ser. A Math. Anal., vol. 16, pp. 253–272, 2009.
  * [23] A. Gaudy, E. Gaudy, Microbiology of waste waters, Annu. Rev. Microbiol., vol. 20, pp. 319–36, 1966.
  * [24] P.T. Harker, Derivatives of the Perron Root of a Positive Reciprocal Matrix: With Application to the Analytic Hierarchy Process, Appl. Math. Comput., 22, pp. 217–232, 1987.
  * [25] J. Harmand, C. Lobry, A. Rapaport, T. Sari, The Chemostat: Mathematical Theory of Microorganism Cultures, Wiley-ISTE, 2017.
  * [26] M. W. Hirsch, H. L. Smith, X.-Q. Zhao, Chain transitivity, attractivity and strong repellors for semidynamical systems, J. Dynam. Differential Equations, vol. 13, 1, pp. 107–131, 2000.
  * [27] S.-B. Hsu, Limiting behavior for competing species, SIAM J. Appl. Math., vol. 34, pp.760–763, 1978.
  * [28] S.B.-Hsu, P. Waltman, A survey of mathematical models of competition with an inhibitor, Math. Biosci., vol. 187, pp.53–91, 2004.
  * [29] C. Lobry La compétition dans le chémostat, Travaux En Cours 81 : Des Nombres et des Mondes, pp. 119–187, édition Herman, Paris, 2013.
  * [30] P. Magal, Perturbation of a Globally Stable Steady State and Uniform Persistence, J. Dynam. Differential Equations, vol. 21, pp., 1–20, 2009.
  * [31] P. Masci, O. Bernard, F. Grognard, Continuous selection of the fastest growing species in the chemostat, IFAC Proceedings Volumes, vol. 41, 2, pp. 9707–9712, 2008.
  * [32] S. Mirrahimi, B. Perthame, J.Y. Wakano, Evolution of species trait through resource competition, J. Math. Biol., vol. 64, 7, pp. 1189–1223, 2012.
  * [33] J. Monod, Recherches sur la Croissance des Cultures Bactériennes, Hermann, Paris 1942.
  * [34] J. Monod, La technique de culture continue théorie et applications, Ann. Inst. Pasteur, 79, pp. 390–410, 1950.
  * [35] A. Novick, L. Szilard, Experiments with the chemostat on spontaneous mutations of bacteria, PNAS 36: pp.708–719, 1950.
  * [36] B. Perthame Transport equations in biology, Birkhäuser Verlag, Berlin, 2007.
  * [37] A. Rapaport, M. Veruete, A new proof of the competitive exclusion principle in the chemostat, Discrete Contin. Dyn. Syst. Ser. B, vol. 24, pp. 3755–3764, 2019.
  * [38] P.L. Salceanu, Robust uniform persistence in discrete and continuous dynamical systems using Lyapunov exponents, Math. Biosci. Eng., vol. 8, 3, pp. 807–825, 2011.
  * [39] T. Sari, A Lyapunov function for the chemostat with variable yields, C. R. Math. Acad. Sci. Paris, vol. 348, 13–14, pp. 747–751, 2010.
  * [40] H.L. Smith, H.R. Thieme, Dynamical systems and population persistence, Providence, R.I: American Mathematical Society, 2011.
  * [41] H.L. Smith, P. Waltman, The theory of the chemostat, Dynamics of microbial competition, Cambridge University Press, 1995.
  * [42] H.L. Smith, P. Waltman, Perturbation of a globally stable steady state, Proc. Amer. Math. Soc., vol. 127, 2, pp. 447–453, 1999.
  * [43] H. Thieme, Convergence results and a Poincaré Bendixson trichotomy for asymptotically autonomous differential equations, J. Math. Biol, vol. 30, pp. 755–763, 1992.
  * [44] H. Thieme, Persistence under relaxed point-dissipativity (with application to an endemic model), SIAM J. Math. Anal., vol. 24, 2, pp.407–435, 1993.
  * [45] H. Thieme, Uniform weak implies uniform strong persistence for non-autonomous semiflows, Proc. Amer. Math. Soc., vol. 127, 8, pp. 2395–2403, 1999.
  * [46] G.S.K. Wolkowicz, Z. Lu, Global dynamics of a mathematical model of competition in the chemostat: general response functions and differential death rates, SIAM J. Appl. Math., vol. 52, pp. 222–233, 1992.
  * [47] Y. Xu, Y. Lai, Derivatives of functions of eigenvalues and eigenvectors for symmetric matrices, J. Math. Anal. Appl., 444, pp. 251–274, 2016.