R Fekih-Salem

A Rapaport
email: alain.rapaport@montpellier.inra.fr

T Sari

Emergence of coexistence and limit cycles in the chemostat model with flocculation for a general class of functional responses

Keywords: Chemostat, Flocculation, Bi-stability, Coexistence, Hopf bifurcation, Limit cycle

We consider a model of two microbial species in a chemostat competing for a single-resource, involving the flocculation of the most competitive species which is present in two forms: isolated and attached. We first show that the model with one species and a non-monotonic growth rate of isolated bacteria may exhibit bi-stability and allows the appearance of unstable limit cycles through a sub-critical Hopf bifurcations due to the joined effect of inhibition and flocculation. We then show that the model with two species presents an even richer set of possible behaviors: coexistence, bi-stability and occurrence of stable limit cycles through a super-critical Hopf bifurcations. All these features cannot occur in the classical chemostat model, where generically at most one competitor can survive on a single resource.

Introduction

In the mathematical model of competition of n species for a single growth-limiting nutrient in a chemostat, a classical result, well-known as the Competitive Exclusion Principle (CEP), asserts that generically at most one species can survive to the competition [START_REF] Grover | Resource competition[END_REF][START_REF] Hansen | Single-nutrient microbial competition: qualitative agreement between experimental and theoretically forecast outcomes[END_REF][START_REF] Hardin | The competitive exclusion principale[END_REF][START_REF] Sari | A Lyapunov function for the chemostat with variable yields[END_REF][START_REF] Sari | Global dynamics of the chemostat with different removal rates and variable yields[END_REF][START_REF] Smith | The Theory of the Chemostat: Dynamics of Microbial Competition[END_REF]. The dynamical equations of the model are

     Ṡ = D(S in -S) - n i=1 1 y i f i (S)x i , ẋi = [f i (S) -D]x i , i = 1, . . . , n (1)
where S(t) denotes the concentration of the substrate at time t, x i (t) denotes the concentration of the species i at time t and n represents the number of species. The operating parameters S in and D denote, respectively, the concentration of the substrate in the feed device and the dilution rate of the chemostat. For i = 1, . . . , n, the function f i (•) represents the per-capita growth rate of the species i (or its functional response) and y i is the yield yield constant which can be chosen equal to one, without loss of generality. A crucial assumption in this classical chemostat model [START_REF] Ad | Kinetic models of suspension flocculation by polymers[END_REF] is that the specific growth rates only depend upon the substrate concentration and that are independent of the concentration of microbial species.

The model [START_REF] Ad | Kinetic models of suspension flocculation by polymers[END_REF] has been extensively studied in the literature, see for example Smith and Waltman [START_REF] Smith | The Theory of the Chemostat: Dynamics of Microbial Competition[END_REF]. In [START_REF] Butler | A mathematical model of the chemostat with a general class of functions describing nutrient uptake[END_REF], Butler and Wolkowicz have studied the model [START_REF] Ad | Kinetic models of suspension flocculation by polymers[END_REF] for a general class of growth rates including monotonic and non-monotonic growth functions such as the Monod and Haldane laws. This last law takes into account the growth-limiting for low concentrations of substrate and the growth-inhibiting for high concentrations. For distinct break-even concentrations, they have demonstrated that the competitive exclusion holds where at most one competitor avoids extinction. In some cases, the species that wins the competition depends on the initial condition. Extension with inhibition by reaction product has been also considered in [START_REF] Patnaik | Dynamic sensitivity of a chemostat for a microbial reaction with substrate and product inhibition[END_REF].

Essajee and Tanner [START_REF] Essajee | The effect of extracellular variables on the stability of the continuous baker's yeast-ethanol fermentation process[END_REF] have introduced in their study variable yields that depend linearly on the substrate concentration and they have showed the occurrence of limit cycles. See also [START_REF] Crooke | Hopf bifurcations for a variable yield continuous fermentation model[END_REF][START_REF] Crooke | The effect of the specific growth rate and yield expressions on the existence of oscillatory behaviour of a continuous fermentation model[END_REF] for a numerical solutions and other growth models. Recently, Sari and Mazenc, in [START_REF] Sari | Global dynamics of the chemostat with different removal rates and variable yields[END_REF], were able to construct a Lyapunov function to study the global dynamics of the model [START_REF] Ad | Kinetic models of suspension flocculation by polymers[END_REF] with a general class of growth rates, distinct removal rates for each species and variable yields, depending on the concentration of substrate. They showed that at most one competitor can survive, the one with the lowest break-even concentration, that is the species that consumes less substrate at steady state. Sari, in [START_REF] Sari | A Lyapunov function for the chemostat with variable yields[END_REF][START_REF] Sari | Competitive exclusion for chemostat equations with variable yields[END_REF], proposed a new Lyapunov function for the study of global asymptotic behavior of model [START_REF] Ad | Kinetic models of suspension flocculation by polymers[END_REF], which is an extension of the Lyapunov functions used by Hsu [START_REF] Hsu | Limiting behavior for competing species[END_REF] and by Wolkowicz and Lu [START_REF] Wolkowicz | Global dynamics of a mathematical model of competition in the chemostat: general response functions and differential death rates[END_REF].

However, the CEP contradicts the biodiversity that is observed, for instance in aquatic ecosystems where phytoplankton species competing for same resources can coexist (see [START_REF] Hutchinson | The paradox of the plankton[END_REF][START_REF] Scheffer | Why plankton communities have no equilibrium: solutions to the paradox[END_REF]). Such a biodiversity is also observed in laboratory, with mixed cultures including at least two competitors for a single resource (see [START_REF] Heßeler | Coexistence in the chemostat as a result of metabolic by-products[END_REF][START_REF] Schmidt | Characterization of a three bacteria mixed culture in a chemostat: Evaluation and application of a quantitative Terminal-Restriction Fragment Length Polymorphism (T-RFLP) analysis for absolute and species specific cell enumeration[END_REF]).

To construct mathematical models that are more consistent with real-world observations, several improvements of the idealized model of competition have been proposed. Typically, adding terms of inter-specific competition between populations of microorganisms and/or intra-specific competition between individuals of the same species in the classical chemostat model leads to dynamics where species can coexist at the equilibrium, [START_REF] Leenheer | Crowding effects promote coexistence in the chemostat[END_REF][START_REF] Fekih-Salem | Sur un modèle de compétition et de coexistence dans le chemostat[END_REF][START_REF] Wolkowicz | Direct interference on competition in a chemostat[END_REF]. Many papers in the literature have proposed extensions of the classical chemostat model that present periodic solutions or limit cycles due to impulsive effect [START_REF] Li | Periodic solution of a chemostat model with variable yield and impulsive state feedback control[END_REF][START_REF] Zhou | Analysis of competitive chemostat models with the Beddington-DeAngelis functional response and impulsive effect[END_REF]. The predator-prey models show the existence of periodic solutions due to Hopf bifurcations [START_REF] Smith | The Theory of the Chemostat: Dynamics of Microbial Competition[END_REF]. Research on such models dates back to Drake et al. [START_REF] Drake | The food chain[END_REF]. Recently, a predator-prey model with three different simultaneous time delays and diffusion shows the existence of a periodic solution and Hopf bifurcation [START_REF] Zhang | Hopf bifurcation analysis in a Monod-Haldane predator-prey model with delays and diffusion[END_REF]. Fowler [START_REF] Fowler | Starvation kinetics of oscillating microbial populations[END_REF] studied a competition model of several species on a single resource which exhibits oscillations when the competition is not entirely antagonistic but is partly syntrophic. In a starvation situation, these oscillations are extreme and the deterministic model becomes inappropriate and must be replaced by the stochastic model that permits the extinction of species in finite time.

Flocculation is a physical and chemical process in which the isolated or planktonic bacteria naturally aggregate, reversibly, to one another to form macroscopic flocs. This mechanism of attachement could be to a wall like biofilms [START_REF] Costerton | Overview of microbial biofilms[END_REF][START_REF]Mathematical modeling of biofilms[END_REF] or simply a formation of flocs or aggregates [START_REF] Thomas | Flocculation modelling: a review[END_REF]. Jones et al. [START_REF] Jones | The Freter model: A simple model of biofilm formation[END_REF] studied the Freter model of biofilm formation (that represents the functioning of intestine) where the parameter values used for the simulations have been chosen from the experimental data of Freter et al. [START_REF] Freter | Survival and implantation of Escherichia coli in the intestinal tract[END_REF].

In this paper, we consider a flocculation mechanism and show how it can lead also to oscillations and non intuitive phenomena of the dynamics. This mechanism is different than the ones previously considered in the literature for explaining the oscillations that are observed experimentally. Indeed, understanding and exploiting the flocculation process appears to be a major challenge to tackle contemporary issues in the fields of wastewater treatments and development of renewable energy, and to improve next future bioprocesses. In [START_REF] Haegeman | Modeling bacteria flocculation as density-dependent growth[END_REF], the effect of flocculation on the growth dynamics was analyzed with an arbitrary number of bacteria in flocs. Haegeman and Rapaport [START_REF] Haegeman | How flocculation can explain coexistence in the chemostat[END_REF] proposed a competition model of two microbial species on a single nutrient with monotonic increasing uptake functions, where attached bacteria or flocs of bacteria do not grow and are subject to the same dilution rate than isolated biomass. Assuming that the most competitive species inhibits its growth by the formation of flocs, they could explain the coexistence between two species. An extension of this model was studied in [START_REF] Fekih-Salem | La floculation et la coexistence dans le chemostat[END_REF] without neglecting the substrate consumption of attached bacteria, but assuming that they consume less substrate than the isolated bacteria, since the bacteria at the surface of flocs have easier access to the substrate than the bacteria inside the flocs. More recently, Fekih-Salem et al. [START_REF] Fekih-Salem | Extensions of the chemostat model with flocculation[END_REF] proposed a model of flocculation of n species that generalizes several models [START_REF] Jones | The Freter model: A simple model of biofilm formation[END_REF][START_REF] Pilyugin | The simple chemostat with wall growth[END_REF][START_REF] Tang | Population dynamics and competition in chemostat models with adaptive nutrient uptake[END_REF] that have been considered in the literature. Assuming that the flocculation and deflocculation dynamics are fast compared to the growth dynamics, Haegeman and Rapaport [START_REF] Haegeman | How flocculation can explain coexistence in the chemostat[END_REF] could build a density-dependent model with the same dilution rate that is studied in [START_REF] Lobry | A new hypothesis to explain the coexistence of n species in the presence of a single resource[END_REF][START_REF] Lobry | Sur un modèle densité-dépendant de compétition pour une ressource[END_REF]. More precisely, the specific growth rate of each species i does not depend only upon the substrate concentration but depends also on the concentration of the same species. In [START_REF] Lobry | Sur un modèle densité-dépendant de compétition pour une ressource[END_REF], the authors determines a sufficient conditions for coexistence of several species in competition for a single resource by introducing the concept of steady-state characteristic. When the dilution rates are identical, Lobry et al. [START_REF] Lobry | A new hypothesis to explain the coexistence of n species in the presence of a single resource[END_REF][START_REF] Lobry | Sur un modèle densité-dépendant de compétition pour une ressource[END_REF] were able to show the existence and the global stability of the coexistence equilibrium. From experimental data, Harmand and Godon [START_REF] Harmand | Density-dependent kinetics models for a simple description of complex phenomena in macroscopic mass-balance modeling of bioreactors[END_REF] have shown that the bioprocesses with attached biomass is better described using ratio-dependent kinetics. Moreover, the study of a flocculation model [START_REF] Fekih-Salem | Extensions of the chemostat model with flocculation[END_REF] with different dilution rates leads also to density-dependent dilution rates for the overall biomass [START_REF] Fekih-Salem | Extensions of the chemostat model with flocculation[END_REF], which is a new feature. In the present work, we revisit the flocculation model proposed in [START_REF] Haegeman | How flocculation can explain coexistence in the chemostat[END_REF], but considering that the attached bacteria consumes also the substrate. Moreover, we consider a general class of growth rates to study the effects on coexistence of two competitors, to be compared with the results obtained by Butler and Wolkowicz [START_REF] Butler | A mathematical model of the chemostat with a general class of functions describing nutrient uptake[END_REF], where at most one competitor can survive on a single resource in absence of flocculation.

In order to generalize the flocculation modelling in the literature, we consider the following flocculation model where the three first equations have been introduced in [START_REF] Fekih-Salem | Extensions of the chemostat model with flocculation[END_REF]:

       Ṡ = D(S in -S) -f (S)u -g(S)v -f 2 (S)x 2 u = (f (S) -D 0 )u -α(•)u + β(•)v v = (g(S) -D 1 )v + α(•)u -β(•)v ẋ2 = (f 2 (S) -D 2 )x 2 (2)
where u(t) and v(t) denote, respectively, the concentrations of isolated and attached bacteria of the first species at time t; f (•) and g(•) represent, respectively, the per-capita growth rates of the isolated and attached bacteria; D 0 , D 1 and D 2 represent, respectively, the removal rates of isolated and attached bacteria of the first species, and of the second species; α(•)u and β(•)v, denote, respectively, the flocculation and deflocculation rates.

Exactly as in [START_REF] Haegeman | How flocculation can explain coexistence in the chemostat[END_REF], we assume that the second species which is the less competitive does not require to inhibit its growth by a flocculation mechanism in order to coexist with the most competitive species. In fact, the flocs consume less substrate than isolated bacteria since they have less access to the substrate, given that this access to the substrate is proportional to the outside surface of the floc. Hence, a flocculation mechanism can be interpreted as an inhibition of the growth of species. In the literature [START_REF] Ad | Kinetic models of suspension flocculation by polymers[END_REF][START_REF] Yang | A new method for calculation of flocculation kinetics combining Smoluchowski model with fractal theory[END_REF][START_REF] Yu | The role of mixing conditions on floc growth, breakage and re-growth[END_REF], the attachment/detachment rates have been reported to be quite variable depending on the mixing conditions. In a former work [START_REF] Fekih-Salem | Extensions of the chemostat model with flocculation[END_REF], we have shown that when attachment and detachment rates are fast one can build a reduced model without distinction between isolated and attached bacteria, but the resulting growth rate is density-dependent as well as the dilution rate, due to the fact that attached and isolated bacteria have different removal rates. In the present paper, we do not assume that attachment and detachment rates are fast. The attachment/detchament processes and their consequences on growth and mortality are quite complex and not yet thoroughly understood. Nevertheless, many studies report that planktonic cells have often a better kinetics and also a higher mortality than attached ones (see [START_REF] Harms | Influence of substrate diffusion on degradation of dibenzofuran and 3chlorodibenzofuran by attached and suspended bacteria[END_REF][START_REF] Heffernan | Comparison of planktonic and biofilm cultures of Pseudomonas fluorescens DSM 8341 cells grown on fluoroacetate[END_REF][START_REF] Wentland | Spatial variations in growth rate within Klebsiella pneumoniae colonies and biofilm[END_REF]).

Table 1 summarizes the modelling assumptions and describes the flocculation and deflocculation rates used in the literature. Note that W = v/v max where v max denotes the maximum areal biomass density of adherent bacteria and G(•) is a decreasing function. The terms a and b are positive constants.

Modelling assumptions Flocculation and deflocculation rates References

Di = D, i = 0, 1, x2 = 0 α(•) = a(1 -W ), β(•) = b + g(S)(1 -G(W )) Jones et al. [29] Di = D, i = 0, 1, x2 = 0 α(•) = α(S), β(•) = β(S) Tang et al. [44] D0 = D, D1 = 0, x2 = 0 α(•) = a, β(•) = b Pilyugin and Waltman [37] Di = D, i = 0, 1, 2, g(S) = 0 α(•) = au, β(•) = b Haegeman and Rapaport [17] Di = D, i = 0, 1, 2 α(•) = au, β(•) = b Fekih-Salem et al. [12] Di = D, i = 0, 1, x2 = 0 α(•) = α(S, v, u), β(•) = β(S, u, v)
Fekih-Salem et al. [START_REF] Fekih-Salem | Extensions of the chemostat model with flocculation[END_REF] Table 1: Modelling assumptions and the description of flocculation and deflocculation rates. All growth rates are monotonic increasing.

However, in the literature it has not yet been studied the effect of a substrate inhibition on the growth of the planktonic bacteria, which is the matter of the present work.

The paper is organized as follows. We first introduce in Section 2 the flocculation model (3) along with assumptions and preliminary results. In Section 3 the model is studied with only one species, non-monotonic growth rate of isolated bacteria and monotonic growth rate of aggregated bacteria. We show that this model can present unstable limit cycles through a sub-critical Hopf bifurcation. An analysis of the two species model is derived in Section 4 with a monotonic growth rate of the second species. The numerical simulations show the occurrence of stable limit cycles through a super-critical Hopf bifurcations. In Section 5, we study the particular case of the flocculation model where the non-monotonic growth rate of attached bacteria is directly related to the isolated ones. Finally, we draw conclusions in the last Section 6. The mathematical proofs are given in Appendix A and the parameter values used in simulations are provided in Appendix B.

Mathematical model

The model of two competitors on one resource in a chemostat, taking into account the flocculation of the most competitive species, can be written as follows

         Ṡ = D(S in -S) -f (S)u -g(S)v -f 2 (S)x 2 u = f (S)u -au 2 + bv -Du v = g(S)v + au 2 -bv -Dv ẋ2 = f 2 (S)x 2 -Dx 2 (3)
We assume that two isolated bacteria can stick together to form a new floc with rate au, and that a floc can split and liberate isolated bacteria with rate b. For the simplicity of the mathematical analysis, we have assumed the same dilution rate D for all reactants, and we have considered the simple expressions α(S, u, v) = au and β(S, u, v) = b as in references [START_REF] Fekih-Salem | Modèles mathématiques pour la compétition et la coexistence des espèces microbiennes dans un chémostat[END_REF][START_REF] Fekih-Salem | La floculation et la coexistence dans le chemostat[END_REF][START_REF] Haegeman | How flocculation can explain coexistence in the chemostat[END_REF]. In [START_REF] Fekih-Salem | Modèles mathématiques pour la compétition et la coexistence des espèces microbiennes dans un chémostat[END_REF], it is shown that the two species flocculation model (3) with a non-monotonic growth rate only for isolated bacteria exhibits the emergence of stable limit cycle. In this work, we demonstrate that this model with only one species can have also a limit cycle but which is unstable.

Monod and Haldane growth functions are quite popular in microbial models. The Haldane expression can been seen as a generalization of the Monod one, allowing an inhibition for large values of substrate. When the planktonic bacteria exhibits such inhibition, the isolated bacteria, that are more "protected", can exhibit or not an inhibition, depending on the strains and the environment. For sake of simplicity of the mathematical analysis, we have considered more general class of growth functions than the precise expressions of Monod and Haldane ones. As a first step in the analysis of (3), we consider a non-monotonic growth rate only for isolated bacteria. Then, we consider a non-monotonic growth rates for isolated and attached bacteria with a certain correlation. More precisely, the growth rates satisfy the following assumptions: H1: The function f : R + → R + is continuously differentiable, f (0) = 0 and there exist two positive real numbers λ 0 and µ 0 , such that λ 0 < µ 0 and

f (S) > D if S ∈]λ 0 , µ 0 [ f (S) < D if S / ∈ [λ 0 , µ 0 ].
We shall allow λ 0 and/or µ 0 to be equal to +∞, and our results can then be applied for any monotonic growth rate for the isolated bacteria. More precisely, there are two very common types of rate expression for the function f (S): monotonic, with f > 0, and non-monotonic, with typically f > 0 on an interval (0, c) and f < 0 on (c, +∞) where c is a positive constant. The popular Monod and Haldane growth functions are particular instances of such functions.

H2: g(0) = 0 and g (S) > 0 for all S > 0.

H3: f 2 (0) = 0 and f 2 (S) > 0 for all S > 0.

When equations g(S) = D and f 2 (S) = D have solutions, they are unique and we define the usual break-even concentrations

λ 1 = g -1 (D) and λ 2 = f -1 2 (D). Otherwise, we put λ k = +∞ (k = 1, 2).
The following mathematically result is a consequence of the fact that we assume that the microbial species in the model (3) do not die of natural causes and hence the bacteria and the substrate have the same dilution rate D. This assumption allows to reduce the dimension of the model (3) by one and simplifies analysis. In future perspectives, it would be interesting to study the robustness of our results to small perturbations in the bacteria decay rate away from zero. One has the following property. Proposition 2.1. For any non-negative initial condition, the forward solution of (3) remains non-negative and positively bounded. The set

Ω = (S, u, v, x 2 ) ∈ R 4 + : Z = S + u + v + x 2 = S in
is positively invariant and is a global attractor for the dynamics (3).

In the following, we shall use for convenience the abbreviation LES for Locally Exponentially Stable equilibria.

3 Study of the one species model

In order to better understand the qualitative behaviors of trajectories of the flocculation model ( 3), we propose in this section to study first the one species case, that is the 3d model

     Ṡ = D(S in -S) -f (S)u -g(S)v u = f (S)u -au 2 + bv -Du v = g(S)v + au 2 -bv -Dv. (4)
This study has led us to define several functions, given below, that are also useful for the analysis of the two species model studied in Section 4. Note that Prop. 2.1 will allow the analysis of the system of three equations ( 4) to be reduced to the analysis of a planar system.

Existence of equilibria

The equilibria of ( 4) are given by the solutions of the following system

     D (S in -S) = f (S)u + g(S)v 0 = f (S)u -au 2 + bv -Du 0 = g(S)v + au 2 -bv -Dv. (5)
At equilibrium, if u = 0 then one has necessarily v = 0 and vice versa. Thus, one cannot observe a steady state with extinction of isolated or attached bacteria only. Denote

ϕ(S) = f (S) -D and ψ(S) = g(S) -D. (6)
The sum of the second and the third equation of ( 5) gives the following equation

ϕ(S)u + ψ(S)v = 0. (7)
Under the assumptions H1-H3, three possible cases may arise:

λ 0 < λ 1 < µ 0 , λ 0 < µ 0 < λ 1 or λ 1 < λ 0 < µ 0 .
The latter case appears to be unrealistic from the biological point of view because the attached bacteria are expected to have a less easy access to substrate than the isolated bacteria (this access being proportional to the external surface of flocs); see e.g. these papers [START_REF] Harms | Influence of substrate diffusion on degradation of dibenzofuran and 3chlorodibenzofuran by attached and suspended bacteria[END_REF][START_REF] Heffernan | Comparison of planktonic and biofilm cultures of Pseudomonas fluorescens DSM 8341 cells grown on fluoroacetate[END_REF][START_REF] Wentland | Spatial variations in growth rate within Klebsiella pneumoniae colonies and biofilm[END_REF]. So, we shall consider only the two first possible cases in the following. Equation ( 7) has a positive solution (u, v) if and only if ϕ(S) and ψ(S) at steady state have opposite signs (see Fig. 1). Let I and J be the sets given in Table 2. If ϕ(S) < 0 and ψ(S) > 0, then S has to belong to the set I, and when ϕ(S) > 0 and ψ(S) < 0, S has to belong to the set J. Those conditions can be summarized by the single condition S ∈ I ∪ J, with

I =]λ 0 , min(µ 0 , λ 1 )[ and J =] max(µ 0 , λ 1 ), +∞[. (8)
Then we have to seek solutions (S, u, v) of ( 5) with S ∈ I ∪ J. As one has ψ(S) = 0, equation ( 7) can be rewritten as

v = - ϕ(S) ψ(S) u. ( 9
) Case I J λ 1 < µ 0 ]λ 0 , λ 1 [ ]µ 0 , +∞[ µ 0 < λ 1 ]λ 0 , µ 0 [ ]λ 1 , +∞[
Table 2: Intervals of existence of positive equilibria according to the case.

Replacing v by its expression [START_REF] Fekih-Salem | Modèles mathématiques pour la compétition et la coexistence des espèces microbiennes dans un chémostat[END_REF] in the second equation of ( 5), we obtain (for positive u, v) By replacing u by ( 10) in ( 9), we obtain

u = U (S) with U (S) := ϕ(S) aψ(S) [ψ(S) -b] . (10)
v = V (S) with V (S) := - ϕ 2 (S) aψ 2 (S) [ψ(S) -b] . (11)
If the equation ψ(S) = b has solution, it is unique and we set

λ b = ψ -1 (b).
Otherwise, we let λ b = +∞. Note that u and v defined by ( 10) and ( 11), respectively, are positive if and only if 5), we deduce that S in -S = u + v. Replacing u and v by ( 10) and ( 11), we obtain S in -S = H(S) where

S ∈ I ∪ J b with J b = J ∩ [0, λ b [. We remark that the interval J b is empty if b < ψ(µ 0 ) in the case λ 1 < µ 0 (see Fig. 1(a)) and is empty if b = 0 in the case µ 0 < λ 1 (see Fig. 1(b)). From (
H(S) := ϕ(S) aψ 2 (S) [ψ(S) -b] [ψ(S) -ϕ(S)] . (12)
We can then state the following result.

Proposition 3.1. The system (4) has the following equilibria:

1. the washout equilibrium E 0 = (S in , 0, 0), that always exists.

2. a positive equilibrium E 1 = S, ū, v with S solution of the equation

H(S) = S in -S, ū = U ( S), v = V ( S), that exists if and only if S ∈ I ∪ J b .
A straightforward calculation gives the following expression of the derivative of H(•).

H = f (ψ -b)(ψ -2ϕ) aψ 2 + g ϕ -ϕψ + 2ϕ(ψ -b) + bψ aψ 3 (13)
whose sign can be positive or negative at S ∈ I ∪ J b (see Figs. 2 and3). Moreover, the function H(•) is defined and positive on this interval. It vanishes at λ 0 , µ 0 , λ b and tends to infinity as S tends to λ 1 .

We have obtained the following results.

Proposition 3.2.

• If S in λ 0 , then there is no positive equilibrium.

• If λ 0 < S in < µ 0 or S in > λ b , then there exists at least one positive equilibrium. Generically, there is an odd number of positive equilibria.

• If µ 0 < S in < λ b , then the system has generically an even number of positive equilibria. In this case, there exists at least two positive equilibria if λ 1 < µ 0 , otherwise, the system can have no positive equilibrium. In Fig. 2(c), the line of equation x 1 = S in -S seems to be quasi-horizontal due to the scale differences on the axes. In the same manner, Fig. 3 illustrates the case µ 0 < λ 1 where the function H(•) is decreasing (a) or non-monotonic (b) on J b . In all figures, we have chosen the red color for LES equilibria, blue color for unstable equilibria and green color when an equilibrium can change its stability.

(a) x1 H Sin Sin Sin E 1 1 E 2 1 E 3 1 E0 λ0 I λ1 µ0 J b λ b Sin S (b) x1 Sin Sin E 1 1 E 2 1 E 3 1 E 4 1 λ0 λ1 µ0 λ b S (c) x1 Sin H E 1 1 E 2 1 E 3 1 λ0 λ1 S
(a) x1 H Sin Sin Sin E 1 1 E 2 1 E 3 1 E0 λ0 I µ0 λ1 J b λ b Sin S (b) x1 Sin Sin E 1 1 E 2 1 E 3 1 E 4 1 E0 λ0 µ0 λ1 λ b S

Stability of equilibria

In the next Proposition, we give a condition for which the washout is the unique globally asymptotically stable equilibrium of (4). When this condition is not fulfilled, dynamics (4) admits multi-equilibria and we focus then on the study of their local asymptotic stability. Indeed, we are interested by conditions for which the washout is the single stable equilibrium.

Proposition 3.3. If S in < min(λ 0 , λ 1 ), then the washout equilibrium E 0 is globally asymptotically stable for (4) with respect to initial conditions in R 3 + .

To study the local stability of each equilibrium point of dynamics (4), we consider the density of total mass in the chemostat z = S + u + v and the vector y = (S, u) . It is easy to see that system (4) possesses a cascade structure (see e.g. [START_REF] Isidori | Nonlinear Control Systems II[END_REF] and Appendix of [START_REF] Smith | The Theory of the Chemostat: Dynamics of Microbial Competition[END_REF]) in the (z, y) coordinates:

ż = D(S in -z) ẏ = φ 2 (z, y),
where

φ 2 (z, y) = D(S in -S) -f (S)u -g(S) (z -S -u) (f (S) -au -D) u + b (z -S -u) .
Using results from [START_REF] Mischaikow | Asymptotically autonomous semiflows: chain recurrence and Lyapunov functions[END_REF], the three order system (4) can be reduced (for the local stability) to the two-dimensional system which is simply the projection on the plane (S, u)

Ṡ = D(S in -S) -f (S)u -g(S)(S in -S -u) u = (f (S) -au -D)u + b(S in -S -u) (14)
with z = S in the equilibrium of the first dynamics. Recall that the function ϕ(•) was defined in [START_REF] Leenheer | Crowding effects promote coexistence in the chemostat[END_REF]. Then we have proved the following results.

Proposition 3.4. The washout equilibrium E 0 is LES if and only if ϕ(S in ) < 0 and S in < λ b .

Note that when Prop. 3.3 is applicable so the weaker Prop. 3.4 also.

Proposition 3.5. A positive equilibrium E 1 is LES if and only if

H ( S) > -1 if S ∈ I H ( S) < -1 and tr A 1 < 0 if S ∈ J b . ( 15
)
where A 1 is the Jacobian matrix of ( 14) at the equilibrium E 1 = ( S, ū) corresponding to the equilibrium E 1 of (4).

We summarize these results in Table 3 Equilibria Existence condition Stability condition

E 0 always exists ϕ(Sin) < 0 and Sin < λ b E 1 H(S) = Sin -S has solution S ∈ I ∪ J b condition (15)
Table 3: Existence and local stability of equilibria of system (4).

One-parameter bifurcations of equilibria

Our aim in this section is to study the behavior of system ( 14) when the parameter S in is varying and all other parameters are fixed. Denote Sin the critical value of S in for which the curve of the function H(•) is tangent to the line of equation x 1 = S in -S. For S in ∈]λ b , Sin [, there exist a positive equilibrium denoted E 1 1 that is LES, the washout equilibrium E 0 and the positive equilibrium E 2 1 that are unstable (saddles) while the positive equilibrium E 3 1 can change its stability (see Fig. 2(a)). We define

D 1 (S in ) := det A 3 1 and T 1 (S in ) := tr A 3 1
where A 3 1 is the Jacobian matrix of ( 14) at E 3 1 = ( S, ū). Indeed, this latter equilibrium satisfies the condition H ( S) < -1 whereas T 1 (S in ) can change its sign as S in increases. Fig. 4(a) illustrates this change of sign of T 1 (S in ) where D 1 (S in ) is positive or, equivalently, H ( S) < -1. Furthermore, the Jacobian matrix A 3 1 has one pair of complex-conjugate eigenvalues λj (S in ) = α(S in ) ± iβ(S in ), j = 1, 2 which becomes purely imaginary for a particular value S in = S c in such that α(S c in ) = 0, with β(S c in ) = 0. We assume the following property (that has been checked numerically)

dα dS in (S c in ) > 0 .
Therefore, E 3 1 changes its stability through a Hopf bifurcation. To illustrate the change of asymptotic behavior of E 3 1 , we represent the variations of the eigenvalues as the parameter S in increases from S in = 5 to 14.5 (see Fig. 4(b)) where the pair of complex-conjugate eigenvalues crosses the imaginary axis at S c in ≈ 9.9117 from negative half plane to positive half plane. Fig. 4(c) illustrates the one-parameter bifurcation diagram for system [START_REF] Freter | Survival and implantation of Escherichia coli in the intestinal tract[END_REF] showing the effect on the S variable of all non-negative equilibria as S in varies. The unstable cycle appears about E 3 1 and defines its basin of attraction. More precisely, for S in < λ 0 , E 0 is stable. E 1 1 bifurcates from E 0 into the positive quadrant when S in = λ 0 . For S in ∈]λ 0 , µ 0 [, the stability is transferred to E 1 1 , while E 0 becomes a saddle, via a transcritical bifurcation. Similarly, E 2 1 bifurcates from E 0 into the positive quadrant when S in = µ 0 . For S in ∈]µ 0 , λ b [, the stability is transferred to E 0 , while E 2 1 becomes a saddle. Then, (a) u in [, the unstable manifold of E 0 has E 3 1 as omega limit set, while the stable manifold of E 2 1 intersects the axis u = 0. The stable manifold of E 2 1 (red curve in Fig. 5(a)) divides the phase plane in the two basins of attraction of E 1

E 1 1 E 2 1 E 3 1 E0 S (b) u E 1 1 E 2 1 E 3 1 E0 S (c) u E 1 1 E 2 1 E 3 1 E0 S (d) u E 1 1 E 2 1 E 3 1 E0 S (e) u E 1 1 E 2 1 E 3 1 E0 S (f ) u E 1 1 E 2 1 E 3 1 E0 S
1 and E 3 1 .

ii. For S in = S 1 in (heteroclinic bifurcation): the unstable manifold of E 0 and the stable manifold of E 2 1 merge. Once again, the two basins of attraction of E 1

1 and E 3 1 are separated by the stable manifold of E 2 1 (the green and red curves in Fig. 5(b)).

iii. For S in ∈]S 1 in , S 2 in [: the unstable manifold of E 0 has E 1 1 as omega limit set, while the stable manifold of E 2 1 intersects the line S + u = S in . This stable manifold of E 2 1 (red curve in Fig. 5(c)) divides the phase plane into the two basins of attraction of E 1 1 and E 3 1 .

iv. For S in = S 2 in (homoclinic bifurcation): the two stable and unstable manifolds of E 2 1 merge together (the green curve in Fig. 5(d)). This homocline orbit delimits the domain of attraction of E 3 1 .

v. For S in ∈]S 2 in , S c in [ (Andronov-Leontovich bifurcation, see [START_REF] Kuznetsov | Elements of Applied Bifurcation Theory[END_REF], Theorem 6.1 and Fig. 6.8): the unstable cycle bifurcates from the homoclinic cycle (the cyan curve in Fig. 5(e)).

vi. For S in = S c in (sub-critical Hopf bifurcation, see [START_REF] Kuznetsov | Elements of Applied Bifurcation Theory[END_REF]): the unstable limit cycle disappears. For S in ∈]S c in , Sin [, E 3 1 is unstable (see Fig. 5(f)). For S in = Sin (saddle-node bifurcation): E 2 1 and E 3 1 collide, forming a nonhyperbolic saddle-node point, and disappear.

We conclude that the one species model (4) with a non-monotonic functional response presents a richness of behaviors with possibly multiple positive equilibria, bi-stability and existence of unstable limit cycle resulting from a sub-critical Hopf bifurcation, all of them being not possible without flocculation. Such richness of behavior is impossible for the flocculation model (4) based upon monotone kinetics (see [START_REF] Fekih-Salem | La floculation et la coexistence dans le chemostat[END_REF] where there exists a unique coexistence equilibrium which is locally exponentially stable once there exists). In the next section, we study the effects on the asymptotic behavior of system (3) when another species that does not aggregate and with monotonic growth is present.

Study of the two species model

In this section, we extend the previous analysis to the full model (3) with two species.

Existence and local stability of equilibria

Under the assumptions H1-H3, we have proved the following result: Proposition 4.1. The system (3) admits the following equilibria:

1. The washout equilibrium F 0 = (S in , 0, 0, 0), that always exists.

2. The equilibrium F 2 = (λ 2 , 0, 0, S in -λ 2 ) of extinction of species 1, that exists if and only if λ 2 < S in .

3. One or more equilibria F 1 = ( S, ū, v, 0) of extinction of species 2, where S is solution of the equation

H(S) = S in -S, ū = U ( S), v = V ( S), that exist if and only if S ∈ I ∪ J b .
4. The positive equilibrium

F * = (λ 2 , u * , v * , x * 2 ) with u * = U (λ 2 ), v * = V (λ 2 ), x * 2 = S in -λ 2 -H(λ 2 ), that exists if and only if λ 2 ∈ I ∪ J b and S in -λ 2 > H(λ 2 ),
where the function U (•), V (•) and H(•) are defined in [START_REF] Fekih-Salem | Extensions of the chemostat model with flocculation[END_REF], [START_REF] Fekih-Salem | Sur un modèle de compétition et de coexistence dans le chemostat[END_REF] and [START_REF] Fekih-Salem | La floculation et la coexistence dans le chemostat[END_REF].

In a similar fashion as the proof of Prop. 3.3, the following result can be obtained: Proposition 4.2. If S in < min(λ 0 , λ 1 , λ 2 ), then the washout equilibrium F 0 is globally asymptotically stable for (3) with respect to initial conditions in R 4 + . As for the study of the one species model, a change of coordinates that reveals a cascade structure of the dynamical system is convenient for the stability analysis. We consider the total mass density Z = S + u + v + x 2 and the vector Y = (S, u, x 2 ) . The system (3) is then equivalent to the system

Ż = D(S in -Z) Ẏ = φ 3 (Z, Y ),
where

φ 3 (Z, Y ) =    D(S in -S) -f (S)u -g(S) (Z -S -u -x 2 ) -f 2 (S)x 2 [f (S) -au -D] u + b (Z -S -u -x 2 ) [f 2 (S) -D] x 2    .
Thus, the fourth-order system (3) can be reduced to the three order system which is simply the projection on the three-dimensional space (S, u, x 2 )

     Ṡ = D(S in -S) -f (S)u -g(S) (S in -S -u -x 2 ) -f 2 (S)x 2 u = [f (S) -au -D] u + b (S in -S -u -x 2 ) ẋ2 = [f 2 (S) -D] x 2 (16)
with Z = S in . Our results about the local stability conditions of the equilibria F 0 , F 1 and F 2 are summarized in the following proposition:

Proposition 4.3.
1. F 0 is LES if and only if ϕ(S in ) < 0 and S in < min(λ b , λ 2 ).

2. F 2 is LES if and only if ϕ(λ 2 ) < 0 and λ 2 < λ b .

3. F 1 is LES if and only if S < λ 2 and the condition ( 15) is satisfied.

Thus we find the stability condition of the one species model (4) and when S < λ 2 it recalls the CEP where the species that consumes less substrate at steady state wins the competition in the classical chemostat model. The next proposition gives one of the key results of the paper -sufficient conditions for the stability or instability of the coexistence steady-state F * = (λ 2 , u * , v * , x * 2 ), depending on the break-even concentration λ 2 .

Proposition 4.4.

1. When λ 2 ∈ I, F * is LES if H (λ 2 ) > -1.
2. When λ 2 ∈ J b , F * is always unstable.

Note that if λ 2 ∈ I and H (λ 2 ) < -1, then F * can be stable or unstable as it will be illustrated in Section 4.2 .

We summarize the results of this section in Table 4 Equilibria Table 4: Existence and local stability of equilibria of system (3).

The result of Prop. 4.4 is illustrated on Fig. 6, with the existence of a unique positive equilibrium F * when λ 1 < µ 0 .

On this figure, one can see that there are three equilibria F 1 , that we denote F k 1 with k = 1, 2, 3.

-when λ 2 ∈ I (case a), F * is LES and all the equilibria F 1 1 , F 2 1 , F 3 1 , F 2 and F 0 are unstable. Numerical simulations can show the global convergence towards the coexistence equilibrium F * from any positive initial condition.

-when λ 2 ∈ J b (case b), F * is unstable as F 2 1 , F 3 1 and F 0 while F 1 1 and F 2 are LES. Numerical simulations can show a bi-stability with convergence either to F 1 1 (and consequently the exclusion of the second species) or to F 2 (and consequently the exclusion of the first species).

Occurrence of limit cycles

In the following, we focus on the change of stability of the positive equilibrium F * through a Hopf bifurcation. More precisely, we analyze the bifurcations according to the parameters D and S in , when the conditions µ 0 < λ 1 , λ 2 ∈ I and H (λ 2 ) < -1 are satisfied (since the conditions of the Routh-Hurwitz criterion (A.10) are not fulfilled).

To illustrate the Hopf bifurcation, we consider the same growth rates f (•) and g(•) for the first species than in Section 3, and for the second species we consider a growth rate f 2 (•) of Monod-type

f 2 (S) = m 2 S K 2 + S ,
where m 2 is the maximum growth rate and K 2 the Michaelis-Menten constant. All the values of the parameters required for the simulation are given in Table B.2. Recall from the previous section that Sin denotes the critical value of S in for which the graph of the function H(•) is tangent to the graph of the line S → S in -S. We distinguish two main cases depending on the position of S in relatively to this critical value.

Pictures when S

in ∈]µ 0 , Sin [.
Depending on the value λ 2 , the break-even concentration of the second species, the following change of stability occur -For λ 2 ∈]λ 0 , S1 [, the equilibria F 0 , F 2 , F 1 1 and F 2 1 are unstable and F * is LES (see Fig. 7(a)).

-For λ 2 = S1 , F * coalesces with F 1 1 .

-For λ 2 ∈] S1 , S2 [, F * disappears in a saddle-node bifurcation and transfers stability to F 1 1 .

-For The Jacobian matrix of reduced system (16) at the equilibrium F * = (λ 2 , u * , x *

λ 2 = S2 , F 2 1 coalesces with F * . -For λ 2 ∈] S2 , λ c1 2 [, F * is unstable. -For λ 2 ∈]λ c1 2 , µ 0 [, F * is LES. (a) Sin x1 1 1 H F 1 F * F2 F * F2 F0 λ0 I S1 S2 λ c 1
2 ) has one negative eigenvalue and one pair of complex-conjugate eigenvalues λj (D) = α(D) ± iβ(D), j = 1, 2 which becomes purely imaginary for a particular value D = D c1 such that α(D c1 ) = 0 with β(D c1 ) = 0. One can check (numerically) that the following inequality is fulfilled

dα dD (D c1 ) > 0 .
We have first considered numerically the two following cases.

-For D = 0.8993 (then one has λ 2 2.042), F * is a saddle-focus on the blue curve since the Jacobian matrix has one negative eigenvalue and one pair of complex-conjugate eigenvalues with positive real part λ3 -4.913 and λ1,2 0.0035 ± 0.0536 i .

-For D = 0.9004 (then one has λ 2 = 2.047), F * changes its stability and becomes a stable focus. The eigenvalues are given by λ3 -4.853 and λ1,2 -0.966e -4 ± 0.0567 i .

To better illustrate the change of stability of the positive equilibrium F * , we have drawn the variations of the eigenvalues as D increases (see Fig. 7(b)) where the pair of complex-conjugate eigenvalues crosses the imaginary axis at D c2 0.90037 from negative half plane to positive half plane. Fig. 7(c) gives the corresponding real part of the complex-conjugate eigenvalues.

Fig. 8 shows the ω-limit set projected in coordinates S, x 1 and x 2 depending on the value of the parameter D, which shows the existence of a limit cycle for a certain range of the values of D. Finally, one has the following possible pictures.

1. For D < D c1 0.899275, F * is a saddle-focus and all trajectories converge to the stable node F 1 1 where there is competitive exclusion of the second species: Fig. 9(a) shows the convergence to the equilibrium F 1 1 from any positive initial condition.

2. For D c1 < D < D c2 , the system exhibits a bi-stability with two basins of attraction, one to the limit cycle and the other one to the stable node F 1 1 (see Fig. 9(b)). When D decreases, the radius of limit cycle increases until the critical value D c1 when the limit cycle crosses the saddle point F 2 1 with a Andronov-Leontovich bifurcation [START_REF] Kuznetsov | Elements of Applied Bifurcation Theory[END_REF].

3. For D > D c2 , F * changes its stability and becomes a stable-focus with a super-critical Hopf bifurcation: the trajectories converge to the stable focus F * or to the stable node F 1 1 : Fig. 9(c) shows the bi-stability, with either convergence to F * (2.055, 0.686, 0.432, 0.326) or F 1 1 (1.235, 1.144, 1.121, 0). Fig. 10(a) illustrates the time course of system (3) in the case of exclusion of the second species and the convergence to the equilibrium F 1 1 . Figs. 10(b-c) illustrate a positive, periodic, solution representing coexistence of the two species and show that a limit cycle disappears through a homoclinic bifurcation where the period of the solution changes as a function of D. In fact, the period asymptotes to infinity at a finite value of the bifurcation parameter D = D c1 . Finally, Fig. 10(c) illustrates the convergence the positive equilibrium F * which becomes a stable focus.

Pictures when S in > Sin .

In order to detect the two super-critical Hopf bifurcations, the parameter m 2 is considered as variable and all other parameters are fixed. Indeed, depending on the value of λ 2 , one has the following change of stability.

-For λ 2 ∈]λ 0 , λ c2 2 [, the equilibria F 0 , F 2 are unstable whereas the equilibrium F * is LES (see Fig. 11(a)). (a) -For

x2 x1 S F2 • F0 • F * • F 2 1 • F 1 1 • (b) x2 x1 S F2 • F0 • F * • F 2 1 • F 1 1 • (c) x2 x1 S F2 • F0 • F * • F 2 1 • F 1 1 • Figure
λ 2 ∈]λ c2 2 , λ c1 2 [, F * is unstable. -For λ 2 ∈]λ c1 2 , µ 0 [, F * is LES.
The Jacobian matrix of reduced system (16) at the equilibrium F * = (λ 2 , u * , x * 2 ) has one negative eigenvalue and one pair of complex-conjugate eigenvalues that crosses the imaginary axis at m c1 2 1.864 from negative half plane to positive half plane as m 2 increases. Then it returns to the negative half plane by crossing the imaginary axis at m c2 2 2.085 (see Fig. 11(b)). Fig. 11(c) illustrates the corresponding real part of the complex-conjugate eigenvalues, as a function of m 2 . Fig. 12 illustrates the ω-limit set of the trajectories depending on the parameter m 2 , projected on the axes S, x 1 and x 2 . Finally, the following possible pictures occur.

(a) Sin x1 H F * F2 F * F2 F * F2 λ0 I λ c 2 2 λ c 1 2 µ0 λ1 F0 S (b) β(m2) λ1 λ2 α(m2) ¢ ¢ f f x (c) α(m2) m2 m c 1 2 m c 2 2
1. For m 2 < m c1
2 , all solutions of (3) converge to the stable focus F * . Fig. 13(a) shows the global convergence to F * (2, 0.757, 0.524, 0.37) in three dimensions for m 2 = 1.8, where F 0 and F 2 are unstable.

For m 2 ∈]m c1

2 , m c2 2 [, there is a super-critical Hopf bifurcation with appearance of limit cycles where F * becomes a saddle-focus. Fig. 13(b) shows the convergence to a limit cycle from any positive initial condition when m 2 = 1.95, for which F * (1.714, 0.970, 0.844, 0.121) is a saddle-focus and F 0 and F 2 are unstable. Fig. 13(c) shows the oscillatory coexistence with constant amplitude and frequency over the time. As m 2 increases, the size of the limit cycle get greater, then shrinks and finally disappears in a second super-critical Hopf bifurcation.

For m 2 > m c2

2 , the equilibrium F * returns to be a stable focus and the numerical simulations can show the global convergence towards the coexistence equilibrium of two species from any positive initial condition. Hence the importance of flocculation dynamics to avoid the washout of species and prevent their extinction. Fig. 14 shows the appearance and disappearance of limit cycles in the three-dimensional space (S, x 1 , x 2 ) through two super-critical Hopf bifurcations for different values of m 2 . One can observe the following behaviors.

-For m 2 = 1.83, all trajectories converge to F * .

-For m 2 = 1.875, all trajectories converge to a limit cycle whose size becomes greater until the parameter m 2 reaches the value m 2 = 1.95.

-For m 2 = 2.05, the limit cycle shrinks.

-For m 2 = 2.1, the second Hopf bifurcation occurs when F * changes its stability and becomes a stable focus. Then, all trajectories converge to F * from any positive initial condition.

Correlated growth rates of isolated and attached bacteria

In this section, we investigate the particular case of isolated and attached bacteria growth rates that are correlated in a specific way: we consider that the growth rate of attached bacteria is the same than the isolated one but for a modified substrate concentration that can be interpreted as an "apparent" substrate concentration for the attached biomass:

g(S) = f (pS)
where p is a parameter that belongs to ]0, 1]. The modified concentration pS represents the proportion of the substrate that is accessible or "apparent" for the attached biomass. The multiplicative parameter p measures the effect of being attached depending for a given bacteria strain. Such a modeling consideration is motivated 1. by the fact that the nutrient access is less easy for the attached bacteria compared to isolated ones, 2. by a choice of distinguishing isolated and attached bacteria growth rates with a single parameter p, 3. by a way to obtain the model (3) from a continuous transformation of the classical chemostat model (1) that do not distinguish isolated and attached bacteria (when p is equal to one).

Our objective is to focus on the bifurcations related to the parameter p, under Assumptions H1 and H3. One has λ 1 = λ 0 p and µ 1 = µ 0 p that verify g(λ 1 ) = g(µ 1 ) = D. Hence, the inequalities λ 0 < λ 1 , µ 0 < µ 1 and λ 1 < µ 1 are fulfilled. Consequently, there are two possible cases Note that the interval J b is empty if b = 0 in both cases. Moreover, the function H(•) is defined and positive on the interval I ∪ J b . It vanishes at λ 0 , µ 0 , λ b and µ b , and tends to infinity as S tends to λ 1 or µ 1 . Fig. 15(a) illustrates the interval of existence of equilibria F 1 and F * for the case λ 1 < µ 0 . Fig. 15(b) shows the existence of a unique positive equilibrium F * which is LES on the red curve, while F 0 , F 2 , F 1 1 , F 2 1 , F 3 1 and F 4 1 are unstable. On the blue curve, F * exists and is unstable like F 0 , F 2 1 , F 3 1 and F 4 1 while F 2 and F 1 1 are LES. Fig. 16(a) illustrates the interval of existence of equilibria F 1 and F * in the case µ 0 < λ 1 . Fig. 16(b) shows the existence of a unique positive equilibrium F * which is LES on the red curve, while F 0 , F 2 , F 1 1 , F 2 1 , F 3 1 and F 4 1 are unstable. On the black curve, F * exists and can change its stability with the occurrence of a limit cycle where F 0 , F 2 , F 2 1 , F 3 1 and F 4 1 are unstable while F 1 1 is LES. On the blue curve, F * exists and is unstable as F 0 , F 2 1 , F 3 1 and F 4 1 while F 2 and F 1 1 are LES.

λ 0 < λ 1 < µ 0 < µ 1 or λ 0 < µ 0 < λ 1 < µ 1 .
F * F2 F 1 1 p c 1 p c 2 p c 3 p (b) x1 F 1 1 F * F 2 1 F0 F2 p c 1 p c 2 p c 3 p (c) x2 F2 F * F0 F 1 1 F 2 1 F 1 1 p c
Condition F 0 F 2 F * F 1 1 F 2 1 p ∈ [0, p c 1 [ U U S p ∈]p c 1 , p c 2 [ U U S U p ∈]p c 2 , p c 3 [ U U S U U p ∈]p c 3 , 1] U U S U
Table 5: Existence and local stability of equilibria of (3) corresponding to steady-state diagram in Fig. 17.

Fig. 17 illustrates the steady-state diagram for system (3) in the case λ 2 ∈]λ 0 , µ 0 [ that shows the effect on the variables of all non-negative equilibria as p varies. Table 5 summarizes existence and local stability of equilibria of (3) according to the parameter p where the letter S (resp. U) means stable (resp. unstable). Absence of letter means that the corresponding equilibrium does not exist.

In this case, the positive equilibrium F * disappears with a saddle-node bifurcation with F 1 1 which becomes LES while F 2 1 , F 2 and F 0 are unstable for p > p c 3 . Indeed, the curves f (•) and g(•) coalesce as the parameter p tends to zero 1 and the graph of the curve H on I and J b converges to the graphs of S = λ 0 and S = µ 0 , respectively. The case λ 2 > µ 0 can be treated similarly where the system exhibits a bi-stability of F 1 1 and F 2 for p large enough. Thus, we find the same result than for the classical chemostat model with a non-monotonic growth rate [START_REF] Butler | A mathematical model of the chemostat with a general class of functions describing nutrient uptake[END_REF], where the system may exhibit bi-stability and (generically) at most one species wins the competition according to the initial condition.

We conclude that making the parameter p varying from one to zero, the coexistence equilibrium can appear to be LES while all the other equilibria are unstable, which show that even under the stronger assumption of isolated and attached bacteria growth rates being correlated, the flocculation phenomenon can inhibit the growth of the most competitive species and allow then another species to coexist.

Discussion and conclusion

In this work, we have investigated mathematically and through numerical simulations a model of competition of two species when one of them flocculates, considering general classes of growth functions.

In a first step, we have studied this model when the species that do not flocculate is absent. We have shown the multiplicity of positive equilibria with the possibility of bi-stability of two positive equilibria. That is, the behavior of system depends on the initial condition and the coexistence equilibria can promote either isolated bacteria and/or attached bacteria. The bistability is already known for the chemostat model without flocculation when the growth function is non-monotonic. It can be explained by the fact that for large values of substrate the bacteria get "lazy" and are then washed-out, while it is not the case for small values of substrate concentration. Here, the consideration of attached bacteria plays the role of an "ecological niche" that protect bacteria from washed-out. Therefore, the washed-out equilibrium is replaced by a positive equilibrium with a relatively large value of the concentration, while the equilibrium with small concentration of substrate is preserved. This first result is new and non intuitive: it cannot occur in the classical chemostat model (1) even with a non-monotonic growth rate (in this case, the system may exhibit a bi-stability but with at most one positive equilibrium, the other equilibrium being the washout). Here, the proposed flocculation phenomenon allows to avoid the attraction basin of the washout equilibrium and consequently prevent the species extinction. Flocculation has thus an effect of "protection" of the species. Moreover and even more surprisingly, flocculation and substrate inhibition together allow the occurrence of an unstable limit cycle through a sub-critical Hopf bifurcations. As the input concentration of nutrient S in is increasing, it is expected that the total biomass for large times increases, and there the tradeoff between attached and isolated bacteria could be less favorable to the isolated one, explaining the instability of the positive equilibrium for large values of substrate concentration, for which the attached bacteria are playing a more important role to avoid the washed-out.

Our mathematical analysis of the complete two species model has revealed even richer possible behaviors. We have first shown the existence of a unique positive equilibrium that may be locally exponentially stable while all other equilibria are unstable. Then, the study of bifurcations according to the concentration of substrate in the feed bottle and the dilution rate shows the appearance of stable limit cycles due to Hopf bifurcations and the disappearance through homoclinic bifurcations. Numerical simulations show that the system may exhibit bi-stability with convergence either to a limit cycle or to the exclusion of the second species. In some cases, a coexistence exists about either a (locally exponentially) stable positive equilibrium or a limit cycle, depending on the initial condition. We deduce that the flocculation mechanism could be responsible of a coexistence, and moreover that a growth inhibition of the most competitive species could lead to the occurrence of limit cycles with coexistence of species.

It is possible to get periodic solutions without flocculation if we have a variable yield coefficient. However, the justification of variable yield is questionable since the experimental evidences seem to be rare. We show here that a completely different mechanism could explain periodic solutions, based on the flocculation phenomenon that is experimentally observed with microscopes.

For the one species model, the inhibition and flocculation are antagonist phenomena: inhibition tends to wash out the bacteria while flocculation protects them. These mechanisms can create oscillations (as trajectories of the system are necessarily bounded) or compensate at bi-stability. For the two species model, it is already well known that when one species has a density-dependant growth rate, coexistence is possible with another species, even when this last one has a growth rate that is not density-dependant (see Lobry et al. [START_REF] Lobry | Persistence in ecological models of competition for a single resource[END_REF]). Here, the species in its both form behave similarly to a virtual species with a density-dependant growth rate (see our paper [START_REF] Fekih-Salem | Extensions of the chemostat model with flocculation[END_REF] when attachement/detachment is fast). Moreover, if one considers the competition between both species without flocculation, it is known that there is exclusion (see Butler and Wolkowicz [START_REF] Butler | A mathematical model of the chemostat with a general class of functions describing nutrient uptake[END_REF]) but the winner could depend on the initial condition (due to the inhibition of the first species). Then, the aggregated form of the first species is playing a role of a mediator between the two possible extinctions leading to possible oscillations and/or bi-stability. In any cases, the distinction between bi-stability and limit cycle can be obtained only thru a mathematical analysis that depends on precise values of the parameters. Thus, it suffices to prove that E 0 = (S in , 0) is a LES equilibrium of the reduced system [START_REF] Freter | Survival and implantation of Escherichia coli in the intestinal tract[END_REF]. From (A.2), the Jacobian matrix of ( 14) at E 0 is given by

A 0 = ψ(S in ) ψ(S in ) -ϕ(S in ) -b ϕ(S in ) -b . Hence, tr A 0 = [ψ(S in ) -b] + ϕ(S in ), det A 0 = [ψ(S in ) -b]ϕ(S in ).
Thus, the two eigenvalues of A 0 are ψ(S in ) -b and ϕ(S in ). They are negative, that is, E 0 is a stable node if and only if S in < λ b and ϕ(S in ) < 0. This completes the proof.

Proof of Prop. 3.5. At the equilibrium E 1 = ( S, ū, v), it is easy to see from Prop. 3.1 that z = S + ū + v = S in .

Then it suffices to prove that E 1 = ( S, ū) is a LES equilibrium of the reduced system [START_REF] Freter | Survival and implantation of Escherichia coli in the intestinal tract[END_REF]. From expression (A. According to [START_REF] Fekih-Salem | Modèles mathématiques pour la compétition et la coexistence des espèces microbiennes dans un chémostat[END_REF], we can replace ϕ/ψ by -v/u in the function H (•) given by ( 13) and so we obtain Consequently, if H ( S) > -1, then E 1 is LES. Finally, if H ( S) < -1, then D 1 ( S) < 0 and we deduce that E 1 is unstable.

In case S ∈ J b , we have ψ( S) > 0. From expression (A.6), if H ( S) > -1, it follows that D 1 ( S) < 0 and so E 1 is unstable. This completes the proof. For convenience, we denote ϕ 2 (S) = f 2 (S) -D for the second species. Let J be the Jacobian matrix of ( 16) at (S, u, x 2 ), that we write with the following notation 1. At the washout equilibrium F 0 = (S in , 0, 0, 0), we have Z = S in . Thus, it suffices to prove that F 0 = (S in , 0, 0) is a LES equilibrium of the reduced system [START_REF] Grover | Resource competition[END_REF]. From (A.9), the Jacobian matrix of ( 16) at F 0 is given by

Proof of

J 0 =    A 0 g(S in ) -f 2 (S in ) -b 0 0 ϕ 2 (S in )   
where A 0 is the 2 × 2 Jacobian matrix ( 14) at E 0 . The eigenvalues of J 0 are the eigenvalues of A 0 and ϕ 2 (S in ). From Prop. 3.4, one can conclude that the equilibrium F 0 is LES if and only if ϕ(S in ) < 0 and S in < min(λ b , λ 2 ).

2. At the equilibrium F 2 = (λ 2 , 0, 0, x2 ) where x2 = S in -λ 2 , we have Z = S in . It suffices to prove that F 2 = (λ 2 , 0, x2 ) is a LES equilibrium of the reduced system [START_REF] Grover | Resource competition[END_REF]. From (A.9), the Jacobian matrix of ( 16) at F 2 , is given by

J 2 =    ψ(λ 2 ) -f 2 (λ 2 )x 2 ψ(λ 2 ) -ϕ(λ 2 ) ψ(λ 2 ) -b ϕ(λ 2 ) -b -b f 2 (λ 2 )x 2 0 0    .
The characteristic polynomial is given by P (λ) = det(J 2 -λI), where I is the 3 × 3 identity matrix. Denote C i and L i the columns and lines of the matrix J 2 -λI. The replacements of C 1 by C 1 -C 3 and L 3 by L 3 + L 1 preserve the determinant and one has

P (λ) = -f 2 (λ 2 )x 2 -λ ψ(λ 2 ) -ϕ(λ 2 ) ψ(λ 2 ) 0 ϕ(λ 2 ) -b -λ -b 0 ψ(λ 2 ) -ϕ(λ 2 ) ψ(λ 2 ) -λ      

Figure 1 :

 1 Figure 1: Growth function f of Haldane type and g of Monod type in two cases: λ 1 < µ 0 and µ 0 < λ 1 . Intervals of existence of positive equilibria.

Figure 2 :

 2 Figure 2: The case λ 1 < µ 0 : Multiplicity of positive equilibria: H is increasing (a) and non-monotonic (b) on I. (c) Numerical example of existence of three positive equilibria. Parameter values are given in Table B.1.

Figure 3 :

 3 Figure 3: The case µ 0 < λ 1 : Multiplicity of positive equilibria: H is decreasing (a) and non-monotonic (b) on J b . Parameter values are given in Table B.1.

SinFigure 4 :

 4 Figure 4: Sub-critical Hopf bifurcation: (a) change of sign of T 1 (S in ) where D 1 (S in ) > 0 for λ b = 4 < S in < Sin ≈ 14.9124, (b) variation of a pair of complex-conjugate eigenvalues as S in increases. (c) Steady-state diagram : In blue the unstable equilibria, in red the stable equilibria and in cyan the unstable cycle. Parameter values are given in Table B.1.

Figure 5 :

 5 Figure 5: Case (a) : case S in = 7 < S 1 in ≈ 8.8763. Case (b): case S in = S 1 in (Heteroclinic bifurcation). Case (c): case S 1 in < S in = 8.9 < S 2 in ≈ 9.1162. Case (d) : case S in = S 2 in (Homocline bifurcation). Case (e): case S 2 in < S in = 9.4 < S c in (Unstable limit cycle). Case (f): case S c in < S in = 11 < Sin . On this picture, in red: stable manifold, in blue: unstable manifold and in green: (hetero-)cycle. Parameter values are given in Table B.1.

F 2 F 1 H

 21 ) < 0 and Sin < min(λ b , λ2) Sin > λ2 ϕ(λ2) < 0 and λ2 < λ b (S) = Sin -S has solution S ∈ I ∪ J b S < λ2 and condition (15) F * λ2 ∈ I ∪ J b and Sin -λ2 > H(λ2) λ2 ∈ I and Routh-Hurwitz condition (A.10)

Figure 6 :

 6 Figure 6: When λ 1 < µ 0 : Case a: λ 2 ∈ I; Case b: λ 2 ∈ J b . Parameter values are given in Table B.2.

Figure 7 :

 7 Figure 7: (a) Change of stability: F * is LES on the red curve and unstable on the blue curve. Variation of a pair of complex-conjugate eigenvalues (b) and the corresponding real part (c) as D increases. Parameter values of Table B.2 were used.

Figure 8 :

 8 Figure 8: Projections of the ω-limit set in variables S, x 1 and x 2 , as a function of D, reveal the occurrence of limit cycles through a super-critical Hopf bifurcation and disappearance through a homoclinic bifurcation. Parameter values of Table B.2 were used.

9 :Figure 10 :

 910 Figure 9: Super-critical Hopf bifurcation and homoclinic bifurcation: bi-stability, coexistence and limit cycle: (a) D = 0.8992 < D c1 , (b) D c1 < D = 0.9 < D c2 , (c) D = 0.9020 > D c2 . Parameter values of Table B.2 were used.

Figure 11 :

 11 Figure 11: (a) Change of asymptotic behavior for S in > Sin and λ 2 ∈ Two super-critical Hopf bifurcations: Variation of a pair of complex-conjugate eigenvalues (b) and the corresponding real part (c) as m 2 increases. Parameter values are given in Table B.2.

Figure 12 :

 12 Figure 12: Projections of the ω-limit set in coordinates S, x 1 and x 2 depending on m 2 : Occurrence and disappearance of limit cycles through two super-critical Hopf bifurcations. Parameter values are given in Table B.2.

Figure 13 :

 13 Figure 13: Coexistence: global convergence to the positive equilibrium F * or to the limit cycle. Parameter values are given in Table B.2.

m2 = 1 Figure 14 :

 114 Figure 14: Occurrence and disappearance of limit cycles in the three-dimensional space (S, x 1 , x 2 ) for different values of m 2 . Parameter values are given in Table B.2.

  Generically, equation ψ(S) = b has two positive solutions λ b and µ b (with λ b < µ b ) or no solution and we then put λ b = +∞ and µ b = +∞. Consider the interval I defined in (8) and define J b = J ∩ (]0, λ b [∪]µ b , +∞[) where J =] max(µ 0 , λ 1 ), µ 1 [.

Figure 15 :Figure 16 :

 1516 Figure 15: (a) Interval of existence of equilibria F 1 and F * in the case λ 1 < µ 0 . (b) Existence of a unique positive equilibrium F * . (a)

3 pFigure 17 :

 317 Figure 17: Steady-state diagram of S, x 1 and x 2 as the parameter p varies. Parameter values are given in Table B.2. F 1 1 coalesces with F 0 when p = p c 1 ≈ 0.1251, F 2 1 coalesces with F 0 when p = p c 2 ≈ 0.25 and F 1 1 coalesces with F * when p = p c 3 ≈ 0.4128.

2 ) evaluated at E 1 ,

 21 D 1 ( S) := det A 1 gives D 1 ( S) = a 11 m 22 + a 21 a 12 = -m 22 ψ( S) + f ( S)ū + g ( S)v m 22 + f ( S)ūa 12 -ba 12 . Since m 22 ψ( S) + ba 12 = ϕ( S) ψ( S) -b = aūψ( S), (A.3) it follows that D 1 ( S) = f ( S)ū(a 12 + m 22 ) + g ( S)vm 22 -aūψ( S). (A.4) Replacing u by its expression (10), m 22 can be rewritten as m 22 = -b ϕ( S) ψ( S) + aū + b (A.5) which is positive. Therefore, m 22 = -ϕψ + 2ϕ(ψ -b) + bψ ψ and a 12 + m 22 = -(ψ -b)(ψ -2ϕ) ψ .

H = -f a 12

 12 + m 22 aψ -g v m 22 auψ .Multiplying (H + 1) by auψ finally yields the following expressionD 1 ( S) = -aūψ( S) 1 + H ( S) . (A.6)From the expression (A.5) of m 22 and the equations (9-10), a straightforward calculation of T 1 ( S) := tr A 1 givesT 1 ( S) = -f ( S)ū + g ( S)v + aū + a ū2 v + b v ū .In the case S ∈ I, we have a 12 > 0 and ψ( S) < 0. From expression (A.6), D 1 ( S) > 0 if and only if H ( S) > -1. From (A.4), it follows that D 1 ( S) > 0 ⇐⇒ -aūψ( S) + g ( S)vm 22 > -f ( S)ū(a 12 + m 22 ). (A.7) Dividing this last inequality by a 12 + m 22 , we obtain -f ( S)ū < -aū a 12 + m 22 ψ( S) + g ( S)v m 22 a 12 + m 22 < -ψ( S) + g ( S)v since a 12 + m 22 = aū + P with P and -ψ( S) are positive. Thus a 11 > 0 and so one has T 1 ( S) = -(a 11 + m 22 ) < 0.

Prop. 4 . 1 .(A. 8 )

 418 The equilibria of (3) are solutions of the following system of equations in-S) = f (S)u + g(S)v -f 2 (S)x 2 0 = f (S)u -au 2 + bv -Du 0 = g(S)v + au 2 -bv -Dv 0 = [f 2 (S) -D] x 2 .The fourth equation gives x 2 = 0 or S = λ 2 . If x 2 = 0 then, according to Prop. 3.1, we deduce the existence of equilibria F 0 and F 1 . If x 2 = 0 and u = 0, then one has S = λ 2 , v = 0 and the first equation leads to the equality x 2 = S in -λ 2 . Thus, the equilibrium F 2 = (λ 2 , 0, 0, S in -λ 2 ) exists if and only if λ 2 < S in . If x 2 = 0 and u = 0, then S = λ 2 , v = 0 and from the proof of Prop. 3.1, it follows that u = U (λ 2 ) and v = V (λ 2 ) which are positive if and only if λ 2 ∈ I ∪ J b . From the first equation of (A.8), we obtain x 2 = S in -λ 2 -H(λ 2 ) which is positive if and only if S in -λ 2 > H(λ 2 ).

11 =

 11 -ψ(S) + f (S)u + g (S)(S in -S -u -x 2 ) + f 2 (S)x 2 , a 12 = f (S) -g(S), a 13 = f 2 (S) -g(S),a 21 = f (S)u -b, m 22 = -ϕ(S) + 2au + b, m 31 = f 2 (S)x 2 and a 33 = ϕ 2 (S). Proof of Prop. 4.3.

   . Table B.1: Parameter values and corresponding values λ 0 , µ 0 , λ 1 and λ b . Table B.2: Parameter values for the two species model.

	Parameter	m11	K11	Ki	m12	K12	D	a	b		Sin		λ0	µ0	λ1	λ b
		(h -1 )	(g/l)	(g/l)	(h -1 )	(g/l)	(h -1 )	(l/h/g)	(l/h/g)	(g/l)
	Fig. 2(a)	2.7	1	1.2	1.5	1	0.9	0.1	0.27		3		0.710	1.690	1.5	3.545
	Fig. 2(c)	3	0.1	0.01	1	2	0.1	2	2		15		0.003	0.286	0.222	+∞
	Fig. 3(a)	2.7	0.5	1.2	1.3	1	0.9	0.1	0.15	Variable	0.283	2.117	2.25	4.2
	Figs. 4 and 5	3	1	1.2	1.5	1	0.9	0.005	0.3	Variable	0.528	2.272	1.5	4
	Parameter		m11	K11	Ki	m12	K12	D		a		b		Sin	m2	K2
			(h -1 )	(g/l)	(g/l)	(h -1 )	(g/l)	(h -1 )	(l/h/g)	(l/h/g)	(g/l)	(h -1 )	(g/l)
	Fig. 6(a) Fig. 6(b)		2.7	1	1.2	1.5	1	0.9		0.1		0.27	3	3 1.8	2
	Figs. 7, 8, 9 and 10 Figs. 11, 12, 13 and 14	3	0.5	1.2	1.25	1.5	Variable 0.9	2		2		3.5 3.65	1.78 Variable	2
	Fig. 17		3.5	0.5	1			1.3		2		0.1		4	3.5	1.5

Acknowledgments

The authors wish to thank the financial support of TREASURE euro-Mediterranean research network (https: //project.inria.fr/treasure/). This work was partly done in the PhD thesis of the first author within the INRA/INRIA team MODEMIC, with the financial support of the Averroes program, the PHC UTIQUE project No. 13G1120 and the COADVISE project.

Appendix A: Proofs

Proof of Prop. 2.1. One has

Hence S(t) 0, v(t) 0 and x 2 (t) 0 for all t 0. One has also u = 0 ⇒ u = bv 0 , and then u(t) 0 for all t 0.

Let Z = S + u + v + x 2 . One obtain from (3) Ż = D(S in -Z) and thus the explicit solution

One can write Z(t) max(S in , Z(0)) for all t 0.

Therefore, the solution of (3) is positively bounded and is defined for all t 0. From (A.1), it can be deduced that the set Ω is positively invariant and is a global attractor for [START_REF] Costerton | Overview of microbial biofilms[END_REF].

Proof of Prop. 3.3. If (u + v)(0) 0, then from the first equation of (4), we deduce that the solution S(t) enters and remains in [0, S in ] in a finite time T . For t > T , one has S(t) ∈ [0, S in ] and therefore max(f (S(t)), g(S(t))) α

where α = max(f (S in ), g(S in )). Since S in < min(λ 0 , λ 1 ), we deduce that α < D. Thus,

It follows that u + v converges to 0. From (A.1), we conclude that z = S + u + v converges to S in . This completes the proof.

Let A be the Jacobian matrix of ( 14) at (S, u), that we write with the following notation

with

In the following, we denote E or F the equilibriums of the reduced systems associated to equilibriums E or F of the full dynamics.

Proof of Prop. 3.4. At washout E 0 = (S in , 0, 0), we have z = S in , that is a globally exponentially equilibrium of the z sub-system.

Moreover, the replacements of C 1 by C 1 -C 2 and L 2 by L 1 + L 2 lead to

are the eigenvalues of J 2 . We conclude that F 2 is LES if and only if ϕ(λ 2 ) < 0 and λ 2 < λ b .

3. At the equilibrium F 1 = ( S, ū, v, 0), it is easy to see from Prop. 4.1 that Z = S + ū + v = S in . Then, it suffices to prove that F 1 = ( S, ū, 0) is a LES equilibrium of the reduced system [START_REF] Grover | Resource competition[END_REF]. From (A.9), the Jacobian matrix of ( 16) at F 1 , is given by

The eigenvalues of this matrix are the eigenvalues of the 2 × 2 sub-matrix A 1 and ϕ 2 ( S). From Prop. 3.5, it follows that F 1 is LES if and only if S < λ 2 and the condition (15) holds.

Proof of Prop. 4.4. We denote with a * the quantities evaluating at the positive equilibrium F * = (λ 2 , u * , v * , x * 2 ). It is easy to see from Prop. 4.1 that Z * = λ 2 +u * +v * +x * 2 = S in . Then, it suffices to prove that F * = (λ 2 , u * , x * 2 ) is a LES equilibrium of the reduced system [START_REF] Grover | Resource competition[END_REF]. Let J * denote the Jacobian matrix of reduced system ( 16) at the equilibrium F * . The characteristic polynomial is given by P According to Routh-Hurwitz criterion, F * is LES if and only if

(A.10)

From equation (A.3), we obtain

If λ 2 ∈ J b , then c 3 > 0 and so the Routh-Hurwitz criterion implies that F * is unstable. In the following, we consider the case λ 2 ∈ I and H (λ