# Multiplicity of neutrally stable periodic orbits with coexistence in the chemostat subject to periodic removal rate

Thomas Guilmeau, Alain Rapaport

► To cite this version:

Thomas Guilmeau, Alain Rapaport. Multiplicity of neutrally stable periodic orbits with coexistence in the chemostat subject to periodic removal rate. SIAM Journal on Applied Mathematics, 2024, 84 (1), pp.39-59. 10.1137/23M1552450. hal-04248438

HAL Id: hal-04248438

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

Submitted on 18 Oct 2023

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
MULTIPLICITY OF NEUTRALLY STABLE PERIODIC ORBITS
WITH COEXISTENCE IN THE CHEMOSTAT SUBJECT TO
PERIODIC REMOVAL RATE

THOMAS GUILMEAU* AND ALAIN RAPAPORT†

**Abstract.** We identify a taxonomic property on the growth functions in the multi-species chemostat model which ensures the coexistence of a subset of species under periodic removal rate. We show that proportions of some powers of the species densities are periodic functions, leading to an infinity of distinct neutrally stable periodic orbits depending on the initial condition. This condition on the species for neutral stability possesses the feature to be independent of the shape of the periodic signal for a given mean value. We give also conditions allowing the coexistence of two distinct subsets of species. Although these conditions are non-generic, we show in simulations that when these conditions are only approximately satisfied, then the behavior of the solutions is close from the non-generic case over a long time interval, justifying the interest of our study.

**Key words.** Chemostat model, Periodic removal rate, Coexistence, Poincaré map, Multiplicity of periodic orbits, Neutral stability.

**MSC codes.** 34C25, 37C25, 92-10, 92D25, 92D40

**1. Introduction.** The mathematical model of the chemostat can represent a vast array of natural phenomena where different living species compete for a common limiting resource. In particular, it is widely used to model waste water bio-processes or to capture competition for a single resource in ecological modeling [15, 5]. Recall that the chemostat is originally an experimental device in which bacterial species grow in a perfectly stirred vessel of constant volume, continuously removed and fed with fresh substrate. Of course, the mathematical model of the chemostat reflects also the behavior of this device.

When several species (of densities $x_i$, $i = 1, 2, ...$) all consume a single substrate that is fed in the ecosystem at a constant rate, then the mathematical theory of the chemostat claims that the *Competitive Exclusion Principle* holds, meaning that all the species will asymptotically disappear from the system, except the fittest one [15, 5]. However, in more complex environments, this property does not always hold. This is for instance the case of periodic operation of bioreactors [13] or temporal fluctuations in ecology [6, 4, 10]. If the environment favors for some time one species and then another one in a balanced way, then one may expect the two species to coexist. A large part of the literature is dedicated to the study of the asymptotic behavior of these dynamics, with two or more species, under periodic removal rate [16, 1, 7, 9] or periodic nutrient input [6, 14, 4, 17] or both [18, 11].

Indeed, when a vessel with two species of concentrations $x_i(t)$ at time $t$ ($i = 1, 2$) is fed with a periodic removal rate, the system can exhibit the following different behaviors :

* The two species are washed out, that is $\lim_{t \to +\infty} x_i(t) = 0$, for $i = 1, 2$.

* One species only is washed out while the other survives.

* The two species coexist, that is $x_i(t) > \epsilon > 0$, for $t \ge 0$ and $i = 1, 2$, for some $\epsilon > 0$.

Those behaviors can be predicted using the Floquet theory. Namely, the Flo-

<footnote>* CVN, Univ. Paris-Saclay, Inria, CentraleSupélec, Gif-sur-Yvette, France (thomas.guilmeau@inria.fr).</footnote>

<footnote>†MISTEA, Univ. Montpellier, INRAE, Institut Agro, Montpellier, France (alain.rapaport@inrae.fr).</footnote>
quet exponents allow to determine if the periodic solutions with only one species are unstable. If both of these single-species periodic solutions are hyperbolic unstable, then the theory of competitive planar systems [2] applied to the chemostat model shows that any positive solution converges asymptotically to a periodic solution with species coexistence (a complete description of this theory is given for instance in [15]). However, the possibility of having several attracting periodic orbits remains an open problem. Although there is no theoretical obstruction for this, no such example has been yet exhibited in the literature. Let us underline that for a given set of species, the condition for the single-species periodic solutions to be unstable depends on the periodic removal rate function, and not only on its mean value. Differently to the stationary environment for which the single break-even concentrations determine which species can survive, the shapes of the growth functions come into play in periodic environments.

In the present work, we investigate how having species whose growth functions share similar shapes allow their coexistence and the multiplicity of neutrally stable periodic orbits (with two or more species). We do not assume the single-species periodic solutions to be hyperbolic, and generalize in a functional way the neutral stability condition in constant environment (which is given by the equality of break-even concentrations). The construction that we propose here is non-generic, relying on a relatively strong condition on the growth functions that is deemed as unfeasible in real life. However, one may face practical situations close from this non-generic case, which may provide a better understanding of the periodic chemostat over a long duration.

Our construction lies on the concept of "taxon", which in this paper, denotes a group of species whose growth functions share the same shape, in a sense which is made precise later on. We show that the fittest species from each taxon behave like one species to some extent, and when they can coexist, then there exists an infinite number of periodic orbits. Moreover, we also study the case when two taxa are present in the chemostat. In this situation, we show that the fittest species from each taxon can coexist all together, under an additional condition of instability of single-species periodic solutions chosen in each taxon.

The paper is organized as follows. In the next Section 2, we recall the equations of the model with the usual assumptions and give some useful definitions and preliminary results. In Section 3, we present our new conditions with a concept of "taxon" and show some relevant properties. Sections 4 and 5 give our main results, first for a single taxon and then for competition between taxa. Finally, we presents and discusses several numerical simulations to illustrate our results Section 6, before drawing conclusions in Section 7.

**2. The setting.** We consider the multi-species chemostat model written as follows

$$ (2.1) \qquad \begin{aligned} \dot{s} &= u(t)(s_{in} - s) - \sum_{i=1}^{n} \frac{\mu_i(s)}{Y_i} x_i, \\ \dot{x}_i &= (\mu_i(s) - u(t))x_i, \quad 1 \le i \le n, \end{aligned} $$

(with $n \ge 2$) where $s$ is the substrate concentration, $s_{in}$ the input substrate concentration and $x_i$ ($1 \le i \le n$) are the respective concentrations of the $n$ populations. We recall the well-known fact that the yield coefficients $Y_i$, $1 \le i \le n$ can be taken equal to one without loss of generality, by a change of variables ($x_i/Y_i$ replaced by $x_i$). As usual, the growth functions $\mu_i(\cdot)$ satisfy the following properties.
Assumption 2.1. For any $i = 1, \dots, n$, $\mu_i$ is a $C^1$ increasing function from $\mathbb{R}_+$ to $\mathbb{R}_+$ with $\mu_i(0) = 0$.

Under this assumption, we define classically the break-even concentration for each species $i = 1, \dots, n$ as the function

$$ \lambda_i(v) := \sup \{s \in \mathbb{R}_+; \mu_i(s) < v\}, \quad v > 0. $$

Note that when $\lambda_i(v)$ is finite, one has necessarily $\mu_i(\lambda_i(v)) = v$.

Assumption 2.2. The removal rate $u(\cdot)$ is a measurable function from $[0, +\infty)$ to $[u_-, u_+]$ with $0 < u_- \le u_+$, which is $T$-periodic (with $0 < T < +\infty$). We posit

$$ \bar{u} = \frac{1}{T} \int_{0}^{T} u(t) \, dt. $$

Consider now the variable $b := s + \sum_{i=1}^{n} x_i$, whose dynamics is given by

$$ (2.2) \qquad \dot{b} = u(t)(s_{in} - b). $$

Under Assumption 2.2, one has $u(t) \ge u_- > 0$ at any $t \ge 0$ and one can deduce that the solution of (2.2) converges exponentially to $s_{in}$, independently of the initial condition of (2.1). With Assumption 2.1, the solutions of (2.1) are uniquely defined and bounded for any non-negative initial condition. Therefore the asymptotic behavior of the $n+1$ dimensional system (2.1) is determined by the $n$ dimensional dynamics

$$ (2.3) \qquad \dot{x}_i = \left[ \mu_i \left( s_{in} - \sum_{j=1}^{n} x_j \right) - u(t) \right] x_i, \quad 1 \le i \le n, $$

which leaves the set

$$ \Delta_0 := \left\{ x \in \mathbb{R}_+^n ; \text{ such that } \sum_{i=1}^{n} x_i \le s_{in} \right\} $$

forwardly invariant. In the remaining, we shall consider the asymptotic dynamic (2.3) on the set $\Delta_0$ only (which is biology relevant in the chemostat framework).

Note that solutions of (2.3) with $x_i(0) = 0$ for some $i$ are such that $x_i(t) = 0$ for any $t$. Therefore, by uniqueness of solutions of (2.3), we deduce that a solution $x(\cdot)$ of (2.3) with a positive initial condition has to stay positive for any time.

We give below a Lemma that will be useful in the rest of the paper.

LEMMA 2.3. There exists a number $\underline{s} \in (0, s_{in})$ such that the subset

$$ \Delta := \left\{ x \in \Delta_0 ; s_{in} - \sum_{i=1}^{n} x_i > \underline{s} \right\} $$

is forward invariant and attractive by the dynamics (2.3).
*Proof.* Posit $s = s_{in} - \sum_{i=1}^{n} x_i$. Since $x$ is a solution of (2.3) and the set $\Delta_0$ is forward invariant, it comes

$$ 
\begin{aligned}
\dot{s} &= -\sum_{i=1}^{n} \mu_i(s)x_i(t) + u(t)(s_{in} - s) \\
&\geq -\left(\max_{1 \leq i \leq n} \mu_i(s)\right) \sum_{i=1}^{n} x_i(t) + u(t)(s_{in} - s) \\
&\geq -\left(\max_{1 \leq i \leq n} \mu_i(s)\right) s_{in} + u(t)(s_{in} - s) \\
&\geq g(s) := -\left(\max_{1 \leq i \leq n} \mu_i(s)\right) s_{in} + u_{-}(s_{in} - s).
\end{aligned}
$$

The function $g$ is continuous and decreasing with $g(0) = u_{-}s_{in} > 0$ and $g(s_{in}) < 0$. By the intermediate value Theorem, there exists a number $s_m \in (0, s_{in})$ such that $g(s_m) = 0$ with $g(s) > 0$ for $s < s_m$. Therefore, for any $s \in (0, s_m)$, the domain $\{s > s\}$ is forwardly invariant and attractive by the dynamics of $s$, which amounts to claim that the subset $\Delta$ is forward invariant and attractive by the dynamics (2.3). □

**3. A taxonomic assumption and its consequences.** We shall consider subsets $I \subset \{1, \dots, n\}$ of at least two species, whose growth functions $\mu_i$ ($i \in I$) share a common property, defining what we propose to call a *taxon* in the present context.

**DEFINITION 3.1.** A subset of populations $I \subset \{1, \dots, n\}$ belong to a same taxon if there exists a $C^1$ increasing function $\varphi$ with $\varphi(s) < u_{-}$ and $\varphi(s_{in}) > u_{+}$, and numbers $\alpha_i > 0$, $\beta_i \in \mathbb{R}$ for $i \in I$ such that

$$ (3.1) \qquad \mu_i(s) = \alpha_i \varphi(s) + \beta_i, \quad s \in [s, s_{in}], \quad i \in I. $$

We shall say that such a function $\varphi$ is a *generating growth function of the taxon*. We also define the subset $I^* \subset I$, which is such that

$$ I^* := \{i \in I; \lambda_i(\bar{u}) = \min_{j \in I} \lambda_j(\bar{u})\}. $$

In this definition, the choice of the generating function $\varphi$ is not unique but it can typically represent a *canonical growth function* that verifies Assumption 2.1 so that the growth functions $\mu_i$ among a taxon differ only by an affine transformation away for 0 i.e. on the interval $[s, s_{in}]$. Note that condition (3.1) cannot be imposed for any $s > 0$ if $\beta_i \neq 0$, because growth functions have to be equal to 0 at the origin (Assumption 2.1). Several examples of growth functions that satisfy Assumption 2.1 and condition (3.1) will be given in Section 6. Let us point out that this condition does not prevent the graphs of the functions $\mu_i$ to cross on the domain $[s, s_{in}]$.

The second part of Definition 3.1 concerns the subset $I^* \subset I$, which is made of all the species that share the same minimal break-even concentration. In some sense, the set $I^*$ contains the fittest species, as these are the surviving species in the autonomous chemostat model [5]. Generically, the set $I^*$ is reduced to a singleton, but we shall study in this work the non-generic situation where more than one species belong to $I^*$. Let us recall from [15] that sharing the same break-even concentration does not ensure in general the existence of periodic orbits with coexistence: the instability of periodic solutions with single species is required, which relies on an interplay between the shapes of the growth functions and the time-varying removal rate.
For convenience, we shall define some auxiliary numbers. The function $\varphi$ being increasing with $\varphi(\bar{s}) < u_-$ and $\varphi(s_{in}) > u_+$, and as $\bar{u} \in [u_-, u_+]$, there exists an unique $\bar{s} \in (\bar{s}, s_{in})$ such that

$$ (3.2) \qquad \varphi(\bar{s}) = \bar{u}, $$

and for any $i \in I$ we define the numbers

$$ (3.3) \qquad \gamma_i := \alpha_i^{-1}(\mu_i(\bar{s}) - \bar{u}) = \frac{\beta_i}{\alpha_i} + \bar{u}(1 - \alpha_i^{-1}). $$

*Remark 3.2.* We have for any $i \in I$ that $\mu_i(\lambda_i(\bar{u})) = \bar{u}$, which is equivalent to write

$$ \varphi(\lambda_i(\bar{u})) = \alpha_i^{-1}(\bar{u} - \beta_i) = \bar{u} - \gamma_i, = \varphi(\bar{s}) - \gamma_i. $$

Since the generating growth function $\varphi$ is increasing, a higher value $\gamma_i$ corresponds to a lower break-even concentration $\lambda_i$. In particular, we have that $i \in I^*$ exactly when $\gamma_i = \max_{j \in I} \gamma_j$.

The key point in our study will be to consider particular ratios of powers of species concentrations, for species that belong to $I$, defined as follows

$$ \rho_{ij} := \frac{\alpha \sqrt{x_i}}{\alpha_j \sqrt{x_j}}, \quad i, j \in I. $$

The dynamics of these variables present some particular properties, as shown in the next Proposition.

**PROPOSITION 3.3.** Let $I$ be a subset of populations that belong to a same taxon. For any solution positive solution in $\Delta$, the dynamics of the ratios $\rho_{ij}$ are as follows

$$ (3.4) \qquad \dot{\rho}_{ij} = (\gamma_i - \gamma_j + (\bar{u} - u(t))(\alpha_i^{-1} - \alpha_j^{-1}))\rho_{ij}, \quad i, j \in I. $$

*Proof.* Let us first differentiate the equality $\alpha \sqrt{x_j} \rho_{ij} = \alpha \sqrt{x_i}$ with respect to $t$:

$$ \frac{1}{\alpha_j} x_j^{\frac{1}{\alpha_j} - 1} \dot{x}_j \rho_{ij} + x_j^{\frac{1}{\alpha_j}} \dot{\rho}_{ij} = \frac{1}{\alpha_i} x_i^{\frac{1}{\alpha_i} - 1} \dot{x}_i, \quad i, j \in I $$

and replace $\dot{x}_k$ for $k = i, j$ by the expression $(\alpha_k \varphi(s(t)) + \beta_k - u(t))x_k$. One obtains

$$ x_j^{\frac{1}{\alpha_j}} \left( \varphi(s(t)) + \frac{\beta_j}{\alpha_j} - \alpha_j^{-1} u(t) \right) \rho_{ij} + x_j^{\frac{1}{\alpha_j}} \dot{\rho}_{ij} = x_i^{\frac{1}{\alpha_i}} \left( \varphi(s(t)) + \frac{\beta_i}{\alpha_i} - \alpha_i^{-1} u(t) \right) $$

for $i, j \in I$. Multiplying by $x_j^{\alpha_j}$, one can write

$$ \dot{\rho}_{ij} = \left( \frac{\beta_i}{\alpha_i} - \frac{\beta_j}{\alpha_j} - u(t)(\alpha_i^{-1} - \alpha_j^{1}) \right) \rho_{ij}, \quad i, j \in I. $$

Finally, from the definition (3.3) of numbers $\gamma_k$, $k \in i, j$, one can also write

$$ \dot{\rho}_{ij} = (\gamma_i - \gamma_j + (\bar{u} - u(t))(\alpha_i^{-1} - \alpha_j^{-1})) \rho_{ij}, \quad i, j \in I. \quad \square $$
The dynamics (3.4) of the ratios $\rho_{ij}$ within a taxon presents thus the remarkable feature that the time evolution of each ratio depends only on its initial value and the function $u$, i.e. their dynamics are decoupled. As a consequence, one obtains the following properties of the solutions of (2.3).

**PROPOSITION 3.4.** Let $I$ be a subset of populations that belongs to a same taxon. For any positive solution in $\Delta$, one has

1. For $i, j$ in $I^*$, $\rho_{ij}$ are periodic functions.

2. For $i \in I \setminus I^*$, $x_i$ converges asymptotically to 0.

3. If $\liminf_{t>0} x_i(t) > 0$ for some $i \in I^*$, then $\liminf_{t>0} x_j(t) > 0$ for any other $j \in I^*$.

*Proof.* From (3.4), one obtains the expression

$$ \frac{d}{dt} \log \rho_{ij} = (\gamma_i - \gamma_j) + (\bar{u} - u(t))(\alpha_i^{-1} - \alpha_j^{-1}), \quad t \ge 0 $$

that we integrate between $t$ and $t+T$:

$$ \log \rho_{ij}(t+T) = \log \rho_{ij}(t) + (\gamma_i - \gamma_j)T, \quad t \ge 0 $$

which gives equivalently

$$ \rho_{ij}(t+T) = \rho_{ij}(t)e^{(\gamma_i - \gamma_j)T}, \quad t \ge 0. $$

Then, for $i, j \in I^*$, $\gamma_i = \gamma_j$ and the function $\rho_{ij}$ is thus periodic. If $i \notin I^*$, for any $j \in I^*$, one has $\gamma_j > \gamma_i$ (see Remark 3.2) and thus $x_i(t) \to 0$ for $t \to +\infty$. The last point of the Proposition is a straightforward consequence of point 1. $\square$

This result states that the Competitive Exclusion Principle occurs within a taxon in the periodic chemostat. It also means that when one or several species of a same taxon persist in a periodic chemostat, it can be invaded by a new one belonging to the same taxon preserving the coexistence of all resident populations, under the condition that all species have the same minimal break even concentration (for the average removal rate). Diversity can be then (theoretically) augmented within a same taxon. Once one has shown that ratios are periodic functions, it is not surprising that one could obtain coexistence of species. However, one has to study the asymptotic behavior of the total biomass within a taxon to show the effective convergence of the solutions to a periodic orbit in $\mathbb{R}_+^n$, which is the matter of the next sections.

The (non-generic) property of having identical break-even concentrations that implies coexistence of species is already known in the classical chemostat model with constant removal rate (see for instance [5]). However, as we have recalled earlier, this property does not guarantee the coexistence under periodic removal rate. Instead, integral conditions which depends on the periodic function $u$ have to be fulfilled [15]. Here, the remarkable feature within a taxon is that under the simple condition of equal break-even concentrations, coexistence can be guaranteed whatever is the periodic function $u$ (provided that its average value is equal to the fixed value $\bar{u}$). This property presents thus a robustness with respect to the removal rate $u(\cdot)$ fluctuating about its mean value.

Now and for the rest of the paper, we shall assume that each species population can persist alone, which is ensured by the following hypothesis.

**Assumption 3.5.** One has $\mu_i(s_{in}) > \bar{u}$ for any $i = 1, \dots, n$.
**4. Multiplicity of periodic solutions within a single taxon.** In this section, we consider that all the species belong to a same taxon. Let us underline that we deal with dynamics in $\mathbb{R}_+^n$ (with $n$ possibly larger than 2), for which the mathematical tools for studying asymptotic periodic solutions that are specific to planar dynamics [2] do not apply here. In the proof of the coming results, we shall exploit the results of Proposition 3.4 to rewrite the dynamics of each species as a non-autonomous dynamical system and use the theory of asymptotic periodic semi-flows [19].

**PROPOSITION 4.1.** Assume that the whole set $I = \{1, \dots, n\}$ of populations belong to a same taxon. Then for any positive initial condition in $\Delta$, the solution of (2.3) converges asymptotically to a periodic solution composed of all species in $I^*$, the other species being washout. Moreover, if $I^*$ is not reduced to a singleton, the system (2.3) admits an infinite number of periodic solutions in $\Delta$ with distinct orbits, which depend on the initial proportions of species in $I^*$ only.

*Proof.* Let $x(\cdot)$ be a positive solution of (2.3) in $\Delta$ and $\rho_{ij}(\cdot)$, $i, j \in I$, be the corresponding solutions of (3.4). Take $i^*$ in $I^*$ such that

$$ \alpha_{i^*} = \min_{i \in I^*} \alpha_i. $$

We show that $x_{i^*}$ converges asymptotically to a positive periodic solution.

The variable $x_{i^*}$ can be written as the solution of the non-autonomous scalar dynamics

$$ \dot{x}_{i^*} = (f(t, x_{i^*}) - u(t))x_{i^*} $$

where

$$ f(t, y) := \mu_{i^*} \left( s_{in} - \sum_{i \in I^*} \rho_{ii^*}(t)^{\alpha_i} y^{\frac{\alpha_i}{\alpha_{i^*}}} - \sum_{i \notin I^*} x_i(t) \right). $$

Let $\tilde{\mu}_{i^*}$ be a $C^1$ increasing extension of the function $\mu_{i^*}$ for negative arguments, and consider the dynamics

$$ (4.1) \qquad \dot{y} = (\tilde{f}(t, y) - u(t))y $$

on $\mathbb{R}$, where $\tilde{f}$ consists in replacing $\mu_{i^*}$ by $\tilde{\mu}_{i^*}$ in the expression of the function $f$. Note that the choice of $x_{i^*}$ implies that one has $\alpha_i/\alpha_{i^*} \ge 1$ for any $i \in I^*$. Therefore, $\tilde{f}$ is Lipschitz with respect to $y$ and the solutions of (4.1) are well defined. At $y = s_{in}$, the argument of $\tilde{\mu}_{i^*}$ is negative, but as $\mu_i(0) = 0$ and $\tilde{\mu}_{i^*}$ is increasing, one has necessarily $\tilde{\mu}_{i^*} < 0$ for negative arguments, and then $\tilde{f}(t, s_{in}) \le 0$ for any $t \ge 0$. The set $[0, s_{in}]$ is thus forwardly invariant. Clearly, $x_{i^*}$ is the solution of (4.1) for the initial value $y(0) = x_{i^*}(0)$ (which belongs to $[0, s_{in}]$).

According to Proposition 3.4, the functions $\rho_{ii^*}$ for $i \in I^*$ are $T$-periodic and variables $x_i$ for $i \notin I^*$ converge asymptotically to 0. We then consider the limiting function

$$ \tilde{f}_l(t, y) := \tilde{\mu}_{i^*} \left( s_{in} - \sum_{i \in I^*} \rho_{ii^*}(t)^{\alpha_i} y^{\frac{\alpha_i}{\alpha_{i^*}}} \right) $$
which is T-periodic and verifies

$$ \lim_{t \to +\infty} |\tilde{f}(t, y)y - \tilde{f}_l(t, y)y| = 0 \quad \text{uniformly for } y \in [0, s_{in}]. $$

Consequently, by Proposition 3.2 in [19], the non-autonomous semi-flow of (4.1) in $[0, s_{in}]$ is asymptotically periodic with limit periodic semi-flow of

$$ (4.2) \qquad \dot{y} = (\tilde{f}_l(t, y) - u(t))y $$

(for which $[0, s_{in}]$ is also forwardly invariant).

We follow now the approach exposed in [15] for one dimensional periodic dynamics, but adapted here to our context. Let us consider the Poincaré map $P$ associated to the periodic dynamics (4.2)

$$ P: y_0 \in [0, s_{in}] \mapsto y(T, y_0) \in [0, s_{in}] $$

where $y(\cdot, y_0)$ denotes the solution of (4.2) with $y(0) = y_0$. One has clearly $P(0) = 0$, and from the Theorem of continuous dependency of the solution of ordinary differential equation with respect to the initial condition, $P$ is continuously differentiable with $P'(y_0) = z(T)$, where $z(\cdot)$ is solution of

$$ \dot{z} = (\partial_y \tilde{f}_l(t, y(t, y_0))y(t, y_0) + \tilde{f}_l(t, y(t, y_0)) - u(t))z, \quad z(0) = 1 $$

that is

$$ z(T) = \exp \left( \int_{0}^{T} \partial_{y} \tilde{f}_{l}(t, y(t, y_{0})) y(t, y_{0}) + \tilde{f}_{l}(t, y(t, y_{0})) - u(t) dt \right) > 0. $$

The map $P$ is thus increasing and one has $P'(0) = \exp(T(\mu_{i^*}(s_{in}) - \bar{u})) > 1$ (by Assumption 3.5). So 0 is a repulsive fixed point of the map $P$, and for any $y_0 > 0$, the sequence $\{P^k y_0\}_{k \in \mathbb{N}}$ is strictly monotonic and bounded, thus converging to a positive fixed point $y^*$ of $P$. Moreover, as $y(\cdot, y^*)$ is periodic, one has

$$ \int_{0}^{T} \tilde{f}_{l}(t, y(t, y^{*})) - u(t) dt = 0 $$

and thus

$$ P'(\bar{y}) = \exp \left( \int_{0}^{T} \partial_{y} \tilde{f}_{l}(t, y^{*}) y(t, y^{*}) dt \right). $$

As the functions $\tilde{\mu}_i$ are assumed to be increasing, one has $\tilde{\mu}_{i^*} > 0$ which implies $\partial_y \tilde{f}_l < 0$, and thus $P'(\bar{y}) < 1$. Therefore the map $y \mapsto P(y) - y$ is decreasing at each root, which implies that it cannot have more than one root. We conclude that the (positive) fixed point $y^*$ is unique.

As $P$ admits a finite number of fixed points (indeed only one), one can apply the results about asymptotically autonomous discrete dynamical systems (Theorem 2.4 in [19]), from which one gets

$$ \lim_{k \to +\infty} x_{i^*}(t + kT) = y(t, y^*), \quad t \in [0, T]. $$
and we conclude that the solution $x(\cdot)$ converges asymptotically to a periodic solution of (2.3) in $\Delta$, for which all species in $I^*$ are present, the other species being excluded:

$$ \lim_{k \to +\infty} |x(t + kT) - x^p(t)| = 0, \quad t \in [0, T] $$

where $x^p(\cdot)$ is the periodic solution given by

$$ x_i^p(t) = \begin{cases} \rho_{ii^*}(t)^{\alpha_i} y(t, y^*)^{\frac{\alpha_i}{\alpha_{i^*}}} > 0, & i \in I^*, \\ 0, & i \notin I^*, \end{cases} $$

for $t \in [0, T]$.

Consider now another positive initial condition in $\Delta$ but with the same initial ratios $\rho_{ii^*}(0)$ for $i \in I^*$. According to (3.4), the functions $\rho_{ii^*}(i \in I^*)$ are identical and consequently the limiting periodic dynamics (4.2) is also identical. As this later one admits an unique periodic solution, we conclude that the solution of (2.3) converges asymptotically to the same periodic solution $x^p(\cdot)$.

We now show how to construct an infinity of distinct periodic solutions, when $I^*$ is not reduced to a singleton. Consider a sequence $\{\rho_0^k\}_{k \in \mathbb{N}}$ of positive vectors in $\mathbb{R}^n$ such that

$$ (4.3) \qquad \max_{i \in I^* \setminus \{i^*\}} (\rho_0^{k+1})_i > \max_{t \in [0, T]} \max_{i \in I^* \setminus \{i^*\}} \rho_{ii^*}^k(t), \quad k \in \mathbb{N} $$

where $\rho_{ii^*}^k(\cdot)$ are the periodic solutions of (3.4) with $\rho_{ii^*}^k(0) = (\rho_0^k)_i$ for $i \in I^* \setminus \{i^*\}$. Condition (4.3) imposes that for each $i \in I^* \setminus \{i^*\}$, the orbits $\gamma^+(\rho_{ii^*}^k) = \{\rho_{ii^*}^k(t), t \in [0, T]\}$ of (3.4) are all disjoint for $k \in \mathbb{N}$. Moreover, for each $k$, there exists an unique periodic solution $x^{pk}(\cdot)$ of (2.3) in $\Delta$ for which all species in $I^*$ are present with ratios given by the functions $\rho_{ii^*}^k$. This implies that the periodic orbits $\gamma^+(x^{pk}) = \{x^{pk}(t), t \in [0, T]\}$, $k \in \mathbb{N}$, of (2.3) are all disjoint. Indeed, if $\gamma^+(x^{pk}) = \gamma^+(x^{pl})$ for some $k \neq l$, there should exist $\tau \ge 0$ such that $x^{pk}(t) = x^{pl}(t + \tau)$ for any $t \ge 0$. In particular, one should have $x_i^{pk}(t) = x_i^{pl}(t + \tau)$ for $i \neq i^*$ in $I^*$ and any $t \ge 0$, but as the orbits $\gamma^+(\rho_{ii^*}^k)$, $\gamma^+(\rho_{ii^*}^l)$ are disjoint, one should have $\rho_{ii^*}^k(t + \tau) \neq \rho_{ii^*}^l(t)$ for some $t$, that is $x_i^{pk}(t) \neq x_i^{pl}(t + \tau)$ and thus a contradiction with the fact that the orbits are non distinct. $\square$

The results of Proposition 4.1 are twofold. First, they complete those of Proposition 3.4, since we now have the persistence of the species of $I^*$, while the less fit species of $I \setminus I^*$ are washed-out. This stronger form of the competitive exclusion principle comes from Assumption 3.5 and the fact that there is no species outside the taxon $I$. This result allows coexistence in the periodic setting under the non-generic assumption that the surviving species belong to $I^*$. It generalizes known similar results in the non-autonomous setting [5]. Then, the second aspect of this result concerns the number of distinct periodic coexistence solutions of the system (2.3), which are a continuum of neutrally stable periodic solutions. Proposition 4.1 shows that if the growth functions are close enough to each other, in the sense that the corresponding species belong to $I^*$, then there are infinitely many periodic coexistence solutions. We may expect that eventual sufficient conditions for uniqueness forbid growth functions to be too close in a sense close to ours.
**5. Multiplicity of periodic solutions with more than one taxon.** In this section, we show that it is possible to have coexistence of two taxa in competition, each of them preserving the proportions of species having the same break-even concentrations, leading to an infinite number of periodic orbits. Roughly speaking, the idea of the proof of Proposition 5.2 below is to consider a non-autonomous planar system that is asymptotically periodic, and to revisit the results about periodic planar dynamics in this framework, with the help of the theory of asymptotic periodic semi-flows. For technicalities, we need in this section the following additional hypothesis.

*Assumption 5.1.* The functions $\mu_i$, $i \in \{1, \dots, n\}$, are analytic at any $s > \underline{s}$, and $u$ is an analytic function of $t$.

For convenience, we shall denote for any integrable scalar function $\zeta(\cdot)$ the average quantity by

$$ \langle \zeta \rangle_T := \frac{1}{T} \int_0^T \zeta(t) dt. $$

**PROPOSITION 5.2.** Assume that one has $\{1, \dots, n\} = I_a \cup I_b$, where species in $I_a$, resp. $I_b$, belong to a same taxon. For any fixed positive initial condition in $\Delta$, let $(s_a^p(\cdot), x_a^p(\cdot))$, resp. $(s_b^p(\cdot), x_b^p(\cdot))$ be the asymptotic periodic solution of (2.1) when only species in $I_a$, resp. $I_b$, are initially present. If the conditions

$$ (5.1) \qquad \lambda_{ba} := \langle \mu_i(s_b^p) \rangle_T - \bar{u} > 0, \quad i \in I_a^*, $$

$$ (5.2) \qquad \lambda_{ab} := \langle \mu_i(s_a^p) \rangle_T - \bar{u} > 0, \quad i \in I_b^*, $$

are fulfilled, then the solution of (2.3) converges asymptotically to a periodic solution for which all species in $I_a^*$ and $I_b^*$ are present, the other species being washed out. Moreover, when there exists at least one positive initial condition satisfying the above conditions, and at least one of the subsets $I_a^*$, $I_b^*$ is not reduced to a singleton, there exists an infinity of distinct periodic orbits of (2.3) for which all the species in $I_a^* \cup I_b^*$ are present.

*Proof.* The taxa are characterized by generating functions $\varphi_a$, $\varphi_b$ and numbers $\alpha_i > 0$, $\beta_i \ge 0$ such that

$$ \mu_i(s) = \alpha_i \varphi_a(s) + \beta_i, \quad i \in I_a, \quad \mu_i(s) = \alpha_i \varphi_b(s) + \beta_i, \quad i \in I_b. $$

Take $i_a^* \in I_a^*$ and $i_b^* \in I_b^*$. Let $x(\cdot)$ be a positive solution of (2.3) in $\Delta$ and define the functions, for $i \in I_a^* \cup I_b^*$

$$ (5.3) \qquad r_i(t) = \begin{cases} \rho_{ii_a^*}(t), & i \in I_a^*, \\ \rho_{ii_b^*}(t), & i \in I_b^*, \end{cases} \quad t \ge 0. $$

Then, variables $x_{i_a^*}$, $x_{i_b^*}$ are solutions of the non-autonomous planar dynamics

$$ \dot{y}_a = (f_a(t, y_a, y_b) - u(t))y_a, $$

$$ \dot{y}_b = (f_b(t, y_a, y_b) - u(t))y_b, $$

with

$$ f_a(t, y_a, y_b) := \mu_{i_a^*}(s(t, y_a, y_b)), \quad f_b(t, y_a, y_b) := \mu_{i_b^*}(s(t, y_a, y_b)), $$
where

$$s(t, y_a, y_b) = s_{in} - \sum_{i \in I_a^*} r_i(t)^{\alpha_i} y_a^{\frac{\alpha_i}{\alpha_i^*}} - \sum_{i \in I_b^*} r_i(t)^{\alpha_i} y_b^{\frac{\alpha_i}{\alpha_i^*}} - \sum_{i \notin I_a^* \cup I_b^*} x_i(t).$$

Let $\tilde{\mu}_{i_a^*}$, $\tilde{\mu}_{i_b^*}$ be $C^1$ increasing extensions of the functions $\mu_{i_a^*}$, $\mu_{i_b^*}$ for negative arguments, and consider the dynamics

$$\begin{aligned} \dot{y}_a &= (\tilde{f}_a(t, y_a, y_b) - u(t))y_a, \\ \dot{y}_b &= (\tilde{f}_b(t, y_a, y_b) - u(t))y_b, \end{aligned} \tag{5.4}$$

in the plane, where the functions $\tilde{f}_a$, $\tilde{f}_b$ are defined with the expressions of $f_a$, $f_b$ replacing the functions $\mu_{i_a^*}$, $\mu_{i_b^*}$ by their extensions $\tilde{\mu}_{i_a^*}$, $\tilde{\mu}_{i_b^*}$. As one has $\alpha_i/\alpha_{i_a^*} \ge 1$ for $i \in I_a^*$, and $\alpha_i/\alpha_{i_b^*} \ge 1$ for $i \in I_b^*$, this dynamics is Lipschitz in $(y_a, y_b)$. Moreover the set

$$S := [0, s_{in}] \times [0, s_{in}]$$

is forward invariant as one has $f_a(t, s_{in}, y_b) < 0$, $f_b(t, y_a, s_{in}) < 0$ for any $(y_a, y_b) \in S$ and $t \ge 0$. Solutions of (5.4) are thus well defined in $S$ and unique. The pair $(x_{i_a^*}(\cdot), x_{i_b^*}(\cdot))$ is such a solution for $(y_a(0), y_b(0)) = (x_{i_a^*}(0), x_{i_b^*}(0))$.

With Proposition 3.4, we know that the functions $r_i$ with $i \in I_a^* \cup I_b^*$ are T-periodic, and variables $x_i$ with $i \notin I_a^* \cup I_b^*$ converge asymptotically to 0. We thus consider the limiting dynamics

$$\begin{aligned} \dot{y}_a &= (\tilde{f}_a^l(t, y_a, y_b) - u(t))y_a, \\ \dot{y}_b &= (\tilde{f}_b^l(t, y_a, y_b) - u(t))y_b, \end{aligned} \tag{5.5}$$

where

$$\tilde{f}_a^l(t, y_a, y_b) := \tilde{\mu}_{i_a^*}(s^l(t, y_a, y_b)), \quad \tilde{f}_b^l(t, y_a, y_b) := \tilde{\mu}_{i_b^*}(s^l(t, y_a, y_b))$$

with

$$s^l(t, y_a, y_b) := s_{in} - \sum_{i \in I_a^*} r_i(t)^{\alpha_i} y_a^{\frac{\alpha_i}{\alpha_i^*}} - \sum_{i \in I_b^*} r_i(t)^{\alpha_i} y_b^{\frac{\alpha_i}{\alpha_i^*}},$$

which are time periodic functions. One has also

$$\lim_{t \to +\infty} |(\tilde{f}_a(t, y_a, y_b) - \tilde{f}_a^l(t, y_a, y_b))y_a| = \lim_{t \to +\infty} |(\tilde{f}_b(t, y_a, y_b) - \tilde{f}_b^l(t, y_a, y_b))y_b| = 0$$

uniformly for $(y_a, y_b) \in S$. Therefore, the non-autonomous semi-flow of (5.4) in $S$ is asymptotically periodic with limit periodic semi-flow of (5.5) (see Proposition 3.2 in [19]). The system (5.5) is competitive and we can apply the results of the literature about periodic competitive planar systems, which states that any bounded solution converges to a periodic solution $(y_a^p(\cdot), y_b^p(\cdot))$ (see for instance Theorem 4.2 in [15]). Let $P$ be the Poincaré map associated to this dynamics

$$P: Y_0 \in S \mapsto Y(T, Y_0) \in S$$

where $Y(\cdot, Y_0)$ denotes the solution $(y_a(\cdot), y_b(\cdot))$ of (5.5) with $(y_a(0), y_b(0)) = Y_0$. On the axis $y_a = 0$ or $y_b = 0$, the dynamics is with a single taxon. One can then reproduce the arguments of the proof of Proposition 4.1 to show that there are unique fixed points $Y_a^* = (y_a^*, 0)$, $Y_b^* = (0, y_b^*)$ of $P$ in $S$ with $y_a^* > 0$, $y_b^* > 0$. Moreover,
one has $y_a^* = (x_a^p)_{i_a^*}(0)$, $y_b^* = (x_b^p)_{i_b^*}(0)$, where $x_a^p(\cdot)$, $x_b^p(\cdot)$ are the asymptotic periodic solutions of (2.3) for the initial conditions $\xi, \zeta$

$$ \xi_i = \begin{cases} x_i(0), & i \in I_a, \\ 0, & i \in I_b, \end{cases} \qquad \zeta_i = \begin{cases} 0, & i \in I_a, \\ x_i(0), & i \in I_b \end{cases} $$

(remind from Proposition 3.4 that functions $\rho_{ii^*}$ ($i \in I_a$) or $\rho_{ii^*}$ ($i \in I_b$) remain the same).

The linearized dynamics $\dot{Y} = M(t)Y$ of (5.5) is given by the matrix

$$ M(t) = \begin{bmatrix} a(t) & b(t) \\ c(t) & d(t) \end{bmatrix} $$

with

$$ a(t) = \tilde{f}_a^l(t, y_a(t), y_b(t)) - u(t) + \partial_{y_a} \tilde{f}_a^l(y_t, y_a(t), y_b(t))y_a(t), $$

$$ b(t) = \partial_{y_b} \tilde{f}_a^l(t, y_a(t), y_b(t))y_a(t), $$

$$ c(t) = \partial_{y_a} \tilde{f}_b^l(t, y_a(t), y_b(t))y_b(t), $$

$$ d(t) = \tilde{f}_b^l(t, y_a(t), y_b(t)) - u(t) + \partial_{y_b} \tilde{f}_b^l(t, y_a(t), y_b(t))y_b(t). $$

Along the periodic solution $(x_a^p(\cdot), 0)$, one has

$$ M(t) = \begin{bmatrix} \star & \star \\ 0 & f_b^l(t, ((x_a^p)_{i_a^*}(t), 0) - u(t)) \end{bmatrix} $$

for which the characteristics multiplier $\exp \int_0^T f_b^l(t, ((x_a^p)_{i_a^*}(t), 0) - u(t)) dt$ is equal to $\exp(T\lambda_{ba})$ and larger than 1 under condition (5.1). The fixed point $Y_a^*$ is thus hyperbolic repulsive. In a similar way, $Y_b^*$ is an hyperbolic repulsive fixed point under condition (5.2). This implies that $Y_a^*$ and $Y_b^*$ are isolated fixed points of P.

For $Y_0 = 0$, the solution of (5.5) is identically null and one has

$$ M(t) = \begin{bmatrix} \mu_{i_a^*}(s_{in}) - u(t) & 0 \\ 0 & \mu_{i_b^*}(s_{in}) - u(t) \end{bmatrix}. $$

The characteristics multipliers are thus $\exp \int_0^T \mu_{i_a^*}(s_{in}) - y(t) dt$, $\exp \int_0^T \mu_{i_b^*}(s_{in}) - y(t) dt$ which are larger than one under Assumption 3.5. The zero solution is thus repulsive.

Along any positive solution $y_a(\cdot), y_b(\cdot)$, note that one has $b(t) < 0$ and $c(t) < 0$ at any $t \in [0, T]$. Then, one has $\dot{Y}_1 > 0$ for $Y_1 = 0$ and $Y_2 < 0$, and $\dot{Y}_2 < 0$ for $Y_1 > 0$ and $Y_2 = 0$. Therefore, the second and fourth quadrant are invariant by the linear dynamics $\dot{Y} = M(t)Y$, which implies that the matrix $P'(Y_0)$ has strictly positive diagonal elements and strictly negative off-diagonal elements for a positive $Y_0 \in S$. Let $S' = \{Y \in S; s^l(0, Y) \in (\underline{s}, s_{in}]\}$. By Lemma 2.3, $S'$ is invariant by P and any fixed point of P belongs to $S'$, including $Y_a^*$ and $Y_b^*$. Following the arguments given in [2], the positive fixed points of P lie on a continuous curve $\Gamma$ in $S'$, which connects the fixed points $Y_a^*, Y_b^*$. Under Assumption 5.1, $x(\cdot)$ is analytic and the functions $r_i(\cdot)$ as well. Therefore, the map P is analytic on $S'$. Then, the curve $\Gamma$ is also analytic (see [2, 4]). If there were an infinite number of fixed points of P in $S$ then all the points of the curve $\Gamma$ will be fixed points by analyticity, which contradicts the fact that $Y_a^*$ and $Y_b^*$ are isolated fixed points.
Finally, as P has a finite number of fixed points on S, we can apply the results about asymptotically autonomous discrete dynamical systems (Theorem 2.4 in [19]), from which one gets

$$ \lim_{k \to +\infty} (x_{i_a^*}(t + kT), x_{i_b^*}(t + kT)) = Y(T, Y^*), \quad t \in [0, T]. $$

where $Y^*$ is a fixed point of P in S. As the fixed points on the axes 0, $Y_a^*$ and $Y_b^*$ are all repulsive, we conclude that $Y^*$ is positive, and that $x(\cdot)$ converges asymptotically to the periodic solution $x^p(\cdot)$ given by

$$ x_i^p(t) = \begin{cases} r_i(t)^{\alpha_i} y_a(t, Y^*)^{\frac{\alpha_i}{\alpha_{i_a^*}}} > 0, & i \in I_a^*, \\ r_i(t)^{\alpha_i} y_b(t, Y^*)^{\frac{\alpha_i}{\alpha_{i_b^*}}} > 0, & i \in I_b^*, \\ 0, & i \notin I_a^* \cup I_b^*, \end{cases} $$

for $t \in [0, T]$.

When $I_a^*$ is not reduced to a singleton, take $i_a^\dagger \neq i_a^*$ in $I_a^*$ and consider perturbations $x^\epsilon(\cdot)$ of the solution $x(\cdot)$, as solutions of (2.3) for the initial condition

$$ (5.6) \qquad x_i^\epsilon(0) = \begin{cases} x_i(0), & i \neq i_a^\dagger, \\ (r_{i_a^\dagger}(0) + \epsilon)^{a_{i_a^\dagger}} x_{i_a^*}(0)^{\frac{\alpha_i}{\alpha_{i_a^*}}}, & i = i_a^\dagger, \end{cases} $$

with $\epsilon > 0$. By continuity of solutions of (2.1) with respect to the initial condition, there exists $\bar{\epsilon} > 0$ such that for any $\epsilon \in (0, \bar{\epsilon})$ $x^\epsilon(0)$ belongs to $\Delta$ and conditions (5.1), (5.2) are fulfilled for this new initial condition. As before, we deduce that $x^\epsilon(\cdot)$ converges asymptotically to a periodic solution $x^{\epsilon,p}(\cdot)$ for which all species in $I_a^*$ and $I_b^*$ are present.

Let $\rho_{ij}^\epsilon$ be the ratio functions for the initial condition $x_i^\epsilon(0)$. Note from (5.6) that one has $\rho_{i_a^\dagger i_a^*}^\epsilon(0) = \rho_{i_a^\dagger i_a^*}^\epsilon(0) + \epsilon$. One gets from (3.4)

$$ \begin{aligned} \rho_{i_a^\dagger i_a^*}^\epsilon(t) &= \rho_{i_a^\dagger i_a^*}^\epsilon(0) \exp \int_0^t (\bar{u} - u(\tau))(\alpha_{i_a^\dagger}^{-1} - \alpha_{i_a^*}^{-1}) d\tau \\ &= \rho_{i_a^\dagger i_a^*}^\epsilon(t) + \epsilon \exp \int_0^t (\bar{u} - u(\tau))(\alpha_{i_a^\dagger}^{-1} - \alpha_{i_a^*}^{-1}) d\tau, \quad t \ge 0. \end{aligned} $$

Therefore, the orbits $\gamma^+(\rho_{i_a^\dagger i_a^*}^\epsilon)$ for $\epsilon \in (0, \bar{\epsilon})$ are all distinct and we deduce, as in the proof of Proposition 4.1 that the orbits of the periodic solutions $\gamma^+(x^{\epsilon,p})$ are all distinct. $\square$

**Remark 5.3.** Conditions (5.1), (5.2) are independent of the choice of $i \in I_a^*$, $I_b^*$. Indeed, let $\bar{s}_a = \lambda_i(\bar{u})$ which is identical for any $i \in I_a^*$, and one has

$$ \begin{aligned} \langle \mu_i(s_b^p) \rangle_T - \bar{u} &= \alpha_i \langle \varphi_a(s_b^p) \rangle_T + \beta_i - \bar{u} \\ &= \alpha_i \langle \varphi_a(s_b^p) - \varphi_a(\bar{s}_a) \rangle_T + \alpha_i \varphi_a(\bar{s}_a) + \beta_i - \bar{u} \\ &= \alpha_i \langle \varphi_a(s_b^p) - \varphi_a(\bar{s}_a) \rangle_T \end{aligned} $$

(using the property $\bar{u} = \mu_i(\bar{s}_a) = \alpha_i \varphi_a(\bar{s}_a) + \beta_i$ for $i \in I_a^*$). The sign of $\langle \mu_i(s_b^p) \rangle_T - \bar{u}$ is thus independent of $i \in I_a^*$, the numbers $\alpha_i$ being positive. One obtains symmetrically the same property for the sign of $\langle \mu_i(s_a^p) \rangle_T - \bar{u}$ with $i \in I_b^*$.
Similarly to Proposition 4.1, Proposition 5.2 generalizes known results providing a continuum of neutrally stable periodic solutions of system (2.3). Indeed, when each taxon contains only one species, we recover exactly the results in [15]. However, thanks to our taxonomic assumptions, we are able to extend it from two species to two taxa, under very similar conditions. Then, we are also able to establish the existence of an infinite number of periodic coexistence solutions, generalizing the result of Proposition 4.1 from one to two taxa.

**6. Numerical illustrations.** In this section, we illustrate numerically our results with the class of Hill functions [8]

$$ \mu(s) = \frac{\mu_{max}s^p}{K^p + s^p} $$

that are parameterized by the three positive numbers $\mu_{max}$, K and p. This class of increasing growth functions is quite popular in microbiology, pharmacology or biochemistry for its flexibility and effectiveness in fitting experimental data, and also because it enlarges the well known class of Monod functions

$$ \mu(s) = \frac{\mu_{max}s}{K + s} $$

as the particular case of p equal to 1. The parameter p measures in some way a distance from the Monod model. Quite often, practitioners conduct experiments first with large values of s to estimate the maximal growth rate $\mu_{max} = \lim_{s \to +\infty} \mu(s)$ and then look for the value of s for which the growth rate is equal to $\mu_{max}/2$. One can straightforwardly check that this value is equal to the affinity constant K (sometimes also called the half-saturation constant), whatever is p. This is enough to identify the Monod growth function, but without additional data, different candidates of growth functions in the class of Hill functions could also suit, depending on the parameter p.

To obtain a variety of different situations, we have considered three generating functions with p equal to 2 or 4

$$ (6.1) \qquad \varphi_a(s) = \frac{m_a s^2}{K_a^2 + s^2}, \qquad \varphi_b(s) = \frac{m_b s^4}{K_b^4 + s^4}, \qquad \varphi_c(s) = \frac{m_c s^2}{K_c^2 + s^2}. $$

with parameters given in Table 6.1, so that one of them ($\varphi_b$) has a significantly different shape. Graphs of the functions $\varphi_a$, $\varphi_b$, $\varphi_c$ are depicted on Figure 6.1. For $p > 1$, the graphs of these functions are convex up to K and then concave for larger values. The larger is p, the more the convexity/concavity is pronounced, as one can see on Figure 6.1. Mixing species from taxa a and b or from taxa a and c allows then to generate contrasted situations. Imposing the same constant K is a way to consider species having similar affinity for the resource and to focus on the impact of the different shapes of the growth functions away from this point. We have generated nine growth functions within these three taxa with characteristic numbers given in Table 6.2 and $\underline{s} = 0.5$. Let us denote the sets of indices of species belonging to a same taxon $I_a = \{1, 2, 3\}$, $I_b = \{4, 5, 6\}$, $I_c = \{7, 8, 9\}$. The graphs of these functions are depicted in Figure 6.2, where we have considered for each of these nine growth functions a $C^1$ extension for $s \le \underline{s}$ as a polynomial increasing on $[0, \underline{s}]$ and null at 0. Clearly, Assumptions 2.1, 2.2, 3.5, 5.1 are satisfied. Within a taxon, each species i is characterized by the parameters $\alpha_i$ and $\beta_i$ following Equation (3.1). Since $\alpha_i$ multiplies the generating growth function, it affects more the behavior of the species for large values of the growth rates, hence for larger values of s (the larger is $\alpha_i$,
the higher is the maximal growth rate). Since $\beta_i$ is added to the generating growth function, its impact is greater for smaller values of the growth rates hence at small values of $s$ (the lower is $\beta_i$, the smaller is the growth function). These effects can be seen on Table 6.2 and Figure 6.2.

The operating conditions have been chosen as follows

$$s_{in} = 4, \quad \bar{u} = 0.8, \quad T = 10$$

where the periodic removal rate is

$$(6.2) \qquad u(t) = \bar{u} + 0.2 \sin \left( \frac{2\pi t}{T} \right)$$

For these values, we have checked in all our simulations that the solutions remain in the set $\Delta$ with $s = 0.5$. For the chosen value of $\bar{u}$, numbers $\bar{s}$ and $\gamma_i$ defined in (3.2) and (3.3) are given in Table 6.3 for each taxon. Note that we do not impose the value $\bar{s}$ to be the same for each taxon. Accordingly to Remark 3.2, the fittest species within each taxon are given by the subsets of indices $I_a^* = \{1, 2\}$, $I_b^* = \{4, 5\}$, $I_c^* = \{8, 9\}$.

<table><caption>Table 6.1: Parameters defining the generating functions $\varphi_a, \varphi_b, \varphi_c$</caption><thead><tr><th>$i$</th><th>$a$</th><th>$b$</th><th>$c$</th></tr></thead><tbody><tr><td>$m_i$</td><td>2</td><td>2.95</td><td>1.8</td></tr><tr><td>$K_i$</td><td>$\sqrt{3}$</td><td>$\sqrt{3}$</td><td>$\sqrt{3}$</td></tr></tbody></table>

<table><caption>Table 6.2: Characteristic numbers of the nine growth functions with respect to their taxon</caption><thead><tr><th>$i$</th><th>1</th><th>2</th><th>3</th><th>4</th><th>5</th><th>6</th><th>7</th><th>8</th><th>9</th></tr></thead><tbody><tr><th>$\varphi$</th><td colspan="3">$\varphi_a$</td><td colspan="3">$\varphi_b$</td><td colspan="3">$\varphi_c$</td></tr><tr><th>$\alpha_i$</th><td>0.9</td><td>1.15</td><td>0.85</td><td>0.8</td><td>1.05</td><td>0.6</td><td>0.9</td><td>1.1</td><td>0.7</td></tr><tr><th>$\beta_i$</th><td>0.161</td><td>-0.0165</td><td>0.1455</td><td>0.232</td><td>0.0545</td><td>0.326</td><td>0.161</td><td>0.019</td><td>0.261</td></tr></tbody></table>

<table><caption>Table 6.3: For each taxon, numbers $\bar{s}$ and $\gamma_i$ (with $\bar{u} = 0.8$)</caption><thead><tr><th>species</th><th>1</th><th>2</th><th>3</th><th>4</th><th>5</th><th>6</th><th>7</th><th>8</th><th>9</th></tr></thead><tbody><tr><th>$\varphi$</th><td colspan="3">$\varphi_a$</td><td colspan="3">$\varphi_b$</td><td colspan="3">$\varphi_c$</td></tr><tr><th>$\bar{s}$</th><td colspan="3">1.28</td><td colspan="3">1.30</td><td colspan="3">1.40</td></tr><tr><th>$\gamma_i$</th><td>0.09</td><td>0.09</td><td>0.03</td><td>0.09</td><td>0.09</td><td>0.01</td><td>0.09</td><td>0.09</td><td>0.03</td></tr></tbody></table>
<figure><img src="image_2.png" /><figcaption>Fig. 6.1: Graphs of the generating growth functions: $\varphi_a$ and $\varphi_c$ are close to each other with $\varphi_a$ slightly above $\varphi_c$, while $\varphi_b$ has a different shape.</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>Fig. 6.2: Graphs of the growth functions within each taxon. The orange and blue growth functions reach the value $\bar{u}$ for lower values of $s$ than the green one, implying that the species associated with the latter is expected to be washed-out by the two other ones.</figcaption></figure>

**6.1. Simulations with species of a single taxon.** To illustrate Proposition 4.1, we have run simulations with all the three species in each taxon. For various initial conditions, we found the same asymptotic orbits with the same initial ratios of species, as expected. We observed that species that are not the fittest within their taxon i.e. that do not belong to $I_a^*$, $I_b^*$ or $I_c^*$ are washed-out in presence of all the species in $I_a$, $I_b$ or $I_c$ respectively. On the opposite, species with the largest values of $\gamma_i$ coexist, that are 1 and 2 for taxon a, 4 and 5 for taxon b, and 7, 8 for taxon c (see Table 6.3). To illustrate trajectories over time, we have chosen for simplicity uniform initial distribution between species. Figure 6.3 shows how the green species 3, 6 or 9 are washed-out within their taxon, while trajectories of the orange and blue species which coexist are more or less close to each other depending on the taxon. This is related to the distance of their graphs (see Figure 6.2) or how close are the values of $\alpha_i$ (see Table 6.2). Additionally, Figure 6.4 represents the solutions in the phase portrait for different initial distributions. It shows the multiplicity of periodic orbits with coexistence of two species, which is the main result of the present work. For different initial distributions, we obtained qualitatively the same behaviors but with
orbits that are more or less elongated. This is explained by the fact that proportions are periodically conserved since initial time (but not the total biomass). One can also observe on Figure 6.4 that the size of the orbits depends on the taxon. This is related to the amplitude of the periodic solutions: the more similar the species are (i.e. the closer their graphs are ), the largest amplitudes are. This is why orbits with taxon b are smaller than for the two other taxa.

<figure><img src="image_4.png" /><figcaption>Fig. 6.3: Simulations with species of the same taxon, with $x_i(0) = 0.5$ for each $i$. One species is washed-out by the two others. The curves for taxa a and c look alike, related to the similarity between $\phi_a$ and $\phi_c$.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Fig. 6.4: Multiplicity of periodic orbits among a single taxon. The shape of the periodic solutions varies according to the taxon.</figcaption></figure>

**6.2. Simulations with species of two different taxa.** In this section, we have simulated species from two different pairs of taxa. As explained formerly, we have considered two contrasted situations: significantly different taxa (a and b) and taxa relatively close to each other (a and c). Initial conditions have been chosen with uniform distribution between the six species (three from each taxon), for simplicity. We have first computed numerically numbers $\lambda_{ba}$, $\lambda_{ab}$ defined in (5.1), (5.2):

$$ \lambda_{ba} = 0.012576 > 0, \quad \lambda_{ab} = 0.006272 > 0 $$

According to Proposition 5.2, inter-taxa coexistence is thus possible as these two numbers are positive for the periodic function (6.2). On Figure 6.5, one can see that species that were washed-out when together with species of the same taxon (species 3 and 6, see Section 6.1) are also washed-out here, as expected. We have now four
species from $I_a^* \cup I_b^*$ that coexist. Additional information about the transients are revealed. One can observe that species from taxon b (4 and 5) are initially raising faster than those of taxon a (1 and 2). This is due to the fact that the initial level s of the resource is relatively large, which favors species from the taxon b as their growth functions take larger values (see Figure 6.2). One can also observe that the time to reach a quasi-periodic regime is much slower than in the previous simulations with a single taxon. Note that there are twice as many species that coexist than before, and that the transient dynamics is governed by a scalar dynamics (4.1) for the single taxon case while it is now ruled by a two dimensional one (5.4) for which the transients can be slower.

<figure><img src="image_6.png" /><figcaption>Fig. 6.5: Simulations with $x_i(0) = 0.5$, for $i \in I_a \cup I_b$. Species 5 and 6 are washed-out while the others (which do not belong to the same taxon) are able to coexist together.</figcaption></figure>

We have then considered species of $I_a$ in presence with those of $I_c$ (Figure 6.6). In this case, we computed the numbers

$$ \lambda_{ca} = 0.071 > 0, \quad \lambda_{ac} = -0.0639 < 0. $$

According to Proposition 5.2, inter-taxa coexistence is now no longer guaranteed. However, as the number $\lambda_{ca}$ is positive, we deduce that taxon a can settle. However, although taxon c seems closer from taxon a than b if one looks at their growth curves only (Figure 6.2), any species from taxon c loses the competition with taxon a. here, one needs to explicitly compute the sign of numbers $\lambda_{ba}$ and $\lambda_{cc}$ to predict the issue of the competition. Finally, only species 1 and 2 are asymptotically present, with a transient speed similar to the simulations with the single taxon a (cf Figure 6.3a).

<figure><img src="image_7.png" /><figcaption>Fig. 6.6: Simulations with $x_i(0) = 0.5$, for $i \in I_a \cup I_c$. Species 1 and 2 from the taxon a coexist but wash-out all the other species.</figcaption></figure>
**6.3. Simulations under approximate taxon condition.** The conditions for species to belong to a same taxon, and that some of them have identical numbers $\gamma_i$, are not generic among all increasing growth functions. In the spirit of former works [12, 3] for constant removal rate, we investigate here numerically cases where the condition (3.1) is only approximately satisfied.

For this purpose, we considered an additional species labeled 2' whose growth function $\mu_{2'}$ is closed to $\mu_2$ but that does not belong to the taxon $I_a$. For the illustration, we have simply taken $\mu_{2'} = \mu_2 + \varepsilon\eta$, where $\eta$ is a smooth function null at 0 with $\eta(\lambda_2(\bar{u})) > 0$ that is not proportionate to $\varphi_a$, and $\varepsilon$ is a small number.

*Remark 6.1.* In the particular case where the perturbation is such that $\mu_{2'}$ still satisfies the taxonomic assumption with $\alpha_{2'} = \alpha_2$ and $\gamma_{2'} = \gamma_2 - \varepsilon$ (recall that $\gamma_2 = \gamma_1$), it is possible to study how the perturbation propagates over time. Indeed, we have from Proposition 3.3 that

$$ \dot{\rho}_{12'} = (\varepsilon + (\bar{u} - u(t))(\alpha_1^{-1} - \alpha_2^{-1})) \rho_{12'}. $$

If $x_2(0) = x_{2'}(0)$ and the other initial conditions are the same, we can integrate this dynamics over $k \in \mathbb{N}$ periods from $t = 0$ and obtain the following comparison result:

$$ \rho_{12'}(kT) = \rho_{12}(kT)e^{k\varepsilon T}. $$

For the numerical illustration, we have taken $\eta(s) = \frac{s}{1+s}$ and first run simulations with species 1, 2 and 2'. For $\varepsilon < 0$, the species 2' is asymptotically conducted to wash-out. However, as one can see on Figure 6.7, the transients can be very long when $\mu_{2'}$ is very close to $\mu_2$ so that the three species coexist in an almost periodic manner during a long time horizon. On the opposite, for $\varepsilon > 0$, the species 2' is the final winner of the competition (Figure 6.8). However, it can take a long time for the other species to decline, so that here also the three species coexist in an almost periodic manner during a long time period.

<figure><img src="image_8.png" /><figcaption>Fig. 6.7: Simulations with species 1, 2, 2' together when $\varepsilon < 0$. Species 2' is washed-out after a transient regime whose length depend on $\varepsilon$.</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>Fig. 6.8: Simulations with species 1, 2, 2' together when $\varepsilon > 0$. Species 1 and 2 are both washed-out after a transient regime whose length depend on $\varepsilon$.</figcaption></figure>

Then, we have considered species of taxa $I_a$ and $I_b$ together, as before (Figure 6.5), but where species 2 is replaced by species 2'. Figures 6.9, 6.10 show that depending on the sign of $\varepsilon$, species 2' belongs to the final composition of the ecosystem or not, and that the time necessary to distinguish this issue can be very long, even for values of $\varepsilon$ not extremely small. Note that coexistence of several species is maintained in both cases. A message here is that the analysis of non generic situations, as we do here, could be of some interest when facing cases likely to be closed from the non-generic case, and we believe that this could be even more likely when considering a large number of species.

<figure><img src="image_10.png" /><figcaption>Fig. 6.9: Simulations with species 1, 2', 3, 4, 5, 6 together when $\varepsilon = -10^{-3}$. Species 2' does not belong to $I_a^*$ and is washed-out after a long transient regime.</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>Fig. 6.10: Simulations with species 1, 2', 3, 4, 5, 6 together when $\varepsilon = 10^{-3}$. Now $I_a^*$ is reduced to species 2', which is the only species of taxon *a* to persist.</figcaption></figure>

**7. Conclusion.** In this work, we have proposed a taxonomic condition on growth functions, that allows the coexistence under periodic removal rate of two or more species with the same break-even concentrations within a single taxon. This generalizes the neutral condition of identical break-even concentrations in the chemostat with constant removal rate, which ensures coexistence at steady state.

Our condition is defined relatively to a generating growth function that characterizes a taxon, independently of the period and the shape of the periodic signal. This makes the coexistence result robust with this respect, as well as the extinction of species in the taxon with larger break-even concentrations. Relying on the theory of asymptotic periodic semi-flows, we have shown that when persistence of a taxon is guaranteed, then any solution converges to a periodic orbit, and moreover that there exists an infinite number of distinct periodic orbits depending on the initial proportions of the species (within the taxon). There exists then a continuum of neutrally stable periodic orbits. In addition, we have generalized the existing results about coexistence of two species in the periodic chemostat model to the case of two taxa in competition, leading to a double infinity of periodic orbits with species of both taxon.

The condition of belonging to a taxon is non-generic and gives non-hyperbolic periodic orbits, but we believe that it could guide future constructions of sufficient conditions for the uniqueness of hyperbolic periodic orbits in the chemostat model with periodic removal rate, an open problem for general classes of growth functions.

Finally, we have shown with the help of numerical simulations that when these non-generic conditions are just "almost" fulfilled, what is likely to happen among a huge number of species, then one can observe many almost periodic solutions with coexistence on a long time window, justifying a posteriori the present analysis of neutrally stable periodic solutions.

## REFERENCES

[1] G.J. BUTLER, S.B. HSU AND P.E. WALTMAN, *A mathematical model of the chemostat with periodic washout rate*, SIAM Journal of Applied Mathematics, 45, 435-449, 1985.

[2] P. DE MOTTONI AND A. SCHIAFFINO, *Competition systems with periodic coefficients: a geometric approach*, Journal of Mathematical Biology, 11, 319-335, 1981.

[3] M. EL HAJJI AND A. RAPAPORT, *Practical coexistence of two species in the chemostat - a slow-fast characterization*, Mathematical Biosciences, 218(1), 33-39, 2009.

[4] J.K. HALE AND A.S. SOMOLINOS, *Competition for fluctuating nutrient*, Journal of Mathematical Biology, 18, 255-280, 1983.

[5] J. HARMAND, C. LOBRY, A. RAPAPORT, T. SARI, *The chemostat, mathematical theory of continuous culture of micro-organisms*, ISTE Press, London, 2017.
[6] S.B. HSU, *A Competition Model for a Seasonally Fluctuating Nutrient* Journal of Mathematical Biology, 9, 115-132, 1980.

[7] P. LENAS AND S. PAVLOU, *Coexistence of three competing microbial populations in a chemostat with periodically varying dilution rate*, Mathematical Biosciences 129, 111-142, 1995.

[8] H. MOSER, *The dynamics of bacterial populations maintained in the chemostat*, Carnegie Institution of Washington Publication, 1958.

[9] S. NAKAOKA AND Y. TAKEUCHI, *How can three species coexist in a periodic chemostat? Mathematical and Numerical Study*, Proceedings of the 9th International Conference "Difference Equations and Discrete Dynamical Systems", L. Allen, B. Aulbach, S. Elaydi, R. Sacker (Editors), Los-Angeles (USA), 2-7 Aug. 2004, World Scientific, 121-133, 2005.

[10] S. NOWACK AND I. KLAPPER, *Exclusion in a Temporally Varying Chemostat System: Dependence on Trade-Offs*, SIAM Journal on Applied Mathematics 78(5), 2819-2839, 2018.

[11] Q.-L. PENG AND H.I.FREEDMAN, *Global Attractivity in a Periodic Chemostat with General Uptake Functions*, Journal of Mathematical Analysis and Applications, 249(2), 300-323, 2000.

[12] A. RAPAPORT, D. DOCHAIN AND J. HARMAND, *Long run coexistence in the chemostat with multiple species*, Journal of Theoretical Biology, 257(2), 252-259, 2009.

[13] P.L. SILVERSON, R.R. HUDGINS, *Periodic Operations of Bioreactors*, Elsevier, 2013.

[14] H. L. SMITH, *Competitive coexistence in an oscillating chemostat*, SIAM Journal on Applied Mathematics, 40, 498-522, 1981.

[15] H. SMITH, P. WALTMAN, *The theory of chemostat, dynamics of microbial competition*, Cambridge Studies in Mathematical Biology, Cambridge University Press, 1995.

[16] G. STEPHANOPOULOS, A. FREDERICKSON, R. ARIS, *The growth of competing microbial populations in a cstr with periodically varying inputs*, AIChE Journal, 25(5), 863-872, 1979.

[17] F. WANG AND G. PANG, *Competition in a chemostat with Beddington DeAngelis growth rates and periodic pulsed nutrient*, Journal of Mathematical Chemistry, 44(3), 691-710, 2008.

[18] G. WOLKOWICZ AND X.-Q. ZHAO, *N-species competition in a periodic chemostat*, Differential Integral Equations 11, 465-491, 1998.

[19] X.-Q. ZHAO, *Asymptotic behavior for asymptotic periodic semiflows with applications*, Communications on Applied Nonlinear Analysis, 3(4), 43-66, 1996.