R Fekih-Salem

A Rapaport
email: rapaport@supagro.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 growth yield which can be chosen equal to one, without loss of generality. Several extensions of this model have been studied in the literature. 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] Butler | A mathematical model of the chemostat with a general class of functions describing nutrient uptake[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].

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] Butler | A mathematical model of the chemostat with a general class of functions describing nutrient uptake[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 breakeven 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] Butler | A mathematical model of the chemostat with a general class of functions describing nutrient uptake[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] Fekih-Salem | Sur un modèle de compétition et de coexistence dans le chemostat[END_REF][START_REF] Leenheer | Crowding effects promote coexistence in the chemostat[END_REF][START_REF] Wolkowicz | Direct interference on competition in a chemostat[END_REF]. Recently, 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]. Study of predator-prey in the chemostat model with time delays and diffusion shows the existence of periodic solutions due to Hopf bifurcation [START_REF] Zhang | Hopf bifurcation analysis in a Monod-Haldane predator-prey model with delays and diffusion[END_REF].

In this paper, we consider a mechanism due to flocculation of bacteria 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 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] Freter | An understanding of colonization resistance of the mammalian large intestine requires mathematical analysis[END_REF][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 [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]. 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 densitydependent 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.

The paper is organized as follows. We first introduce in Section 2 the flocculation model ( 2) 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 (2)
where u(t) and v(t) denote, respectively, the concentration of isolated and attached bacteria of the first species at time t. Functions f (•) and g(•) represent, respectively, the per-capita growth rates of the isolated and attached bacteria. The variable x 2 (t) denotes the concentration of the second species and f 2 (•) represents the per-capita growth rate of the second species. We assume that two isolated bacteria can stick together to form a new floc with rate a, and that a floc can split and liberate isolated bacteria with rate b. 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.

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). One has the following property. Proposition 2.1. For any non-negative initial condition, the forward solution of (2) 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 (2).

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

Study of the one species model

In order to better understand the qualitative behaviors of trajectories of the flocculation model [START_REF] Fekih-Salem | Extensions of the chemostat model with flocculation[END_REF], 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. (3)
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.

Existence of equilibria

The equilibria of (3) 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. (4)
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.
The sum of the second and the third equation of (4) gives the following equation

ϕ(S )u + ψ(S )v = 0. ( 5
)
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). So, we shall consider only the two first possible cases in the following. Equation ( 5) 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 1. 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 ), +∞[.

(a) Case I J Then we have to seek solutions (S , u, v) of (4) with S ∈ I ∪ J. As one has ψ(S ) 0, equation ( 5) can be rewritten as

g f D I J ψ ϕ b I J b λ 0 λ 1 λ b µ 0 S (b) g f D I J ψ ϕ b I J b λ 0 µ 0 λ 1 λ b S
λ 1 < µ 0 ]λ 0 , λ 1 [ ]µ 0 , +∞[ µ 0 < λ 1 ]λ 0 , µ 0 [ ]λ 1 , +∞[
v = - ϕ(S ) ψ(S ) u. (7)
Replacing v by its expression [START_REF] Haegeman | Modeling bacteria flocculation as density-dependent growth[END_REF] in the second equation of (4), we obtain (for positive u, v)

u = U(S ) with U(S ) := ϕ(S ) aψ(S ) ψ(S ) -b . (8)
By replacing u by ( 8) in [START_REF] Haegeman | Modeling bacteria flocculation as density-dependent growth[END_REF], we obtain

v = V(S ) with V(S ) := - ϕ 2 (S ) aψ 2 (S ) ψ(S ) -b . (9)
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 ( 8) and ( 9), respectively, are positive if and only if

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 (4), we deduce that S in -S = u + v. Replacing u and v by ( 8) and ( 9), we obtain S in -S = H(S ) where

H(S ) := ϕ(S ) aψ 2 (S ) ψ(S ) -b ψ(S ) -ϕ(S ) . ( 10
)
We can then state the following result.

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

1. the washout equilibrium E 0 = (S in , 0, 0), that always exists. A straightforward calculation gives the following expression of the derivative of H(•).

H = f (ψ -b)(ψ -2ϕ) aψ 2 + g ϕ -ϕψ + 2ϕ(ψ -b) + bψ aψ 3 (11)
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.

Fig. 2 illustrates the number of positive equilibria in the case λ 1 < µ 0 , depending on S in , where the function H(•) is increasing (a) or non-monotonic (b) on I. x 1 = u + v is the total biomass of the first species.

To illustrate the existence of three positive equilibria on I (see Fig. 2(c)), we have considered the parameter values given in Table B.4 with the growth rates f (•) of Haldane-type and g(•) of Monod-type:

f (S ) = m 11 S K 11 + S + S 2 K i and g(S ) = m 12 S K 12 + S .
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) (a)

x 1 H S in S in S in E 1 1 E 2 1 E 3 1 E 0 λ 0 I λ 1 µ 0 J b λ b S in S (b) x 1 S in S in E 1 1 E 2 1 E 3 1 E 4 1 λ 0 λ 1 µ 0 λ b S (c) x 1 S in H E 1 1 E 2 1 E 3 1 λ 0 λ 1 S
x 1 H S in S in S in E 1 1 E 2 1 E 3 1 E 0 λ 0 I µ 0 λ 1 J b λ b S in S (b) x 1 S in S in E 1 1 E 2 1 E 3 1 E 4 1 E 0 λ 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 (3). When this condition is not fulfilled, dynamics (3) admits multi-equilibria and we focus then on the study of their local asymptotic stability. Proposition 3.3. If S in < λ 0 , then the washout equilibrium E 0 is globally asymptotically stable for (3) with respect to initial conditions in R 3 + . To study the local stability of each equilibrium point of dynamics (3), 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 (3) possesses a cascade structure 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) .
Thus, the three order system (3) can be reduced 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) (12)
with z = S in the equilibrium of the first dynamics. Let A 1 denote the Jacobian matrix of reduced system (12) at the equilibrium E 1 = ( S , ū) corresponding to the equilibrium E 1 of (3). We define

D 1 ( S ) := det A 1 and T 1 ( S ) := tr A 1 .
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 .

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

H ( S ) > -1 if S ∈ I H ( S ) < -1 and T 1 ( S ) < 0 if S ∈ J b . (13)
We summarize these results in Table 2 Equilibria Existence condition Stability condition

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

Sub-critical Hopf bifurcation

When S ∈ J b , the equilibrium E 1 can change its stability through a Hopf bifurcation. More precisely, when H ( S ) < -1, T 1 ( S ) can change its sign from positive to negative as S in increases because f (•) is a non-monotonic function. Denote S in 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 . If S in ∈]λ b , S in [, 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 while the positive equilibrium E 3 1 can change its stability (see Fig. 2(a)). In this case, E 3 1 satisfies the condition H ( S ) < -1 whereas T 1 ( S ) can change its sign as S in increases. Fig. 4(a) illustrates this change of sign of T 1 ( S ) where D 1 ( S ) is positive or, equivalently, H ( S ) < -1 from expression (A.6). In fact, the Jacobian matrix of reduced system (12) at E 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. One can check (numerically) that the following inequality is fulfilled

dα dS in (S c in ) > 0 .
With the values of the parameters given in Table B.4, we have first considered the two following values of the parameter S in . i. For S in = 9.8, E 3 1 is a stable focus since the Jacobian matrix has one pair of complex-conjugate eigenvalues with negative real part λ1,2 -0.003 ± 0.259 i .

ii. For S in = 9.95, E 3 1 changes behavior and becomes an unstable focus whose eigenvalues are given by λ1,2 0.001 ± 0.261 i .

To better 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. (a)

T 1 ( S ) We conclude that the one species model (3) 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. In the next section, we study the effects on the asymptotic behavior of system (2) when another species that does not aggregate and with monotonic growth is present.

D 1 ( S ) µ 0 λ b S (b) β(S in ) λ1 λ2 α(S in ) ¡ ¡ ¡ ! e e e (c) u S in E 1 1 E 2 1 E 3 1 E 0 • S
(c) u S in E 1 1 E 2 1 E 3 1 E 0 • S

Study of the two species model

In this section, we extend the previous analysis to the full model (2) 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 (2) 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 ( 8), ( 9) and [START_REF] Hardin | The competitive exclusion principale[END_REF].

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 (2) 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 (2) 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 (14)
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.2.
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 ( 13) is satisfied.

Thus we find the stability condition of the one species model (3) 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 sufficient conditions for the stability or instability of the positive equilibrium F * = (λ 2 , u * , v * , x * 2 ), depending on the break-even concentration λ 2 .

Proposition 4.3. 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 3. The result of Prop. 4.3 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 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 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).

Equilibria Existence condition

Stability condition F 0 always exists ϕ(S in ) < 0 and S in < min(λ b , λ 2 ) (a)

F 2 S in > λ 2 ϕ(λ 2 ) < 0 and λ 2 < λ b F 1 H(S ) = S in
x 1 H S in F 1 1 F * F 2 F 2 1 F 3 1 F 0 λ 0 λ 2 λ 1 µ 0 λ b S in S (b) x 1 H S in F 1 1 F * F 2 F 2 1 F 3 1 F 0 λ 0 λ 2 λ 1 µ 0 λ b S in S Figure 6: When λ 1 < µ 0 : Case a: λ 2 ∈ I; Case b: λ 2 ∈ J b .

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 m 2 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.5. Recall from the previous section that S in 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 , S in [.
Depending on the value λ 2 , the break-even concentration of the second species, the following change of stability occur -For λ 2 ∈]λ 0 , S 1 [, 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 = S 1 , F * coalesces with F 1 1 . -For λ 2 ∈] S 1 , S 2 [, F * disappears in a saddle-node bifurcation and transfers stability to

F 1 1 . -For λ 2 = S 2 , F 2 1 coalesces with F * . -For λ 2 ∈] S 2 , λ c 1 2 [, F * is unstable. -For λ 2 ∈]λ c 1 2 , µ 0 [, F * is LES.
The Jacobian matrix of reduced system (14) at the equilibrium F * = (λ 2 , u * , x * 2 ) has one negative eigenvalue and one pair of complex-conjugate eigenvalues λ j (m 2 ) = α(m 2 ) ± iβ(m 2 ), j = 1, 2 which becomes purely imaginary for a particular value m 2 = m c 1 2 such that α(m c 1 2 ) = 0 with β(m c 1 2 ) 0. One can check (numerically) that the following inequality is fulfilled

dα dm 2 (m c 1 2 ) > 0 .
We have first considered numerically the two following cases.

-For m 2 = 1.78 (then one has λ 2 2.045) and S in < S in 3.61, 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.875 and λ1,2 0.0012 ± 0.0556 i .

-For m 2 = 1.7 (then one has λ 2 = 2.25), F * changes its stability and becomes a stable focus. The eigenvalues are given by λ3 -4.015 and λ1,2 -0.0251 ± 0.0552 i .

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

(a) S in x 1

F 1 1 H F 2 1 F * F 2 F * F 2 F 0 λ 0 I S 1 S 2 λ c 1 2 µ 0 S in λ 1 S (b) β(m 2 ) λ1 λ2 α(m 2 ) 0 t t (c) α(m 2 ) m 2 m c 1 2 •
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 m 2 increases. Fig. 8 shows the ω-limit set projected in coordinates S , x 1 and x 2 depending on the value of the parameter m 2 , which shows the existence of a limit cycle for a certain range of the values of m 2 .

Finally, one has the following possible pictures. 1. For m 2 < m c 1 2 , the trajectories converge to the stable focus F * or to the stable node F 1 1 : Fig. 9(a) shows the bi-stability in three dimensions for m 2 = 1.7, with either convergence to F * (2.25, 0.492, 0.225, 0.533) or

F 1 1 (1.184, 1.163, 1.152, 0). 2. For m c 1 2 < m 2 < m c 2 2
1.790, F * changes its stability and becomes a saddle-focus with a super-critical Hopf bifurcation: Fig. 9(b) shows the bi-stability with two basins of attraction for m 2 = 1.78, one to the limit cycle and the other one to the stable node F 1 1 . Thus, the trajectories converge to a limit cycle whose radius increases with respect to m 2 , until the critical value m c 2 2 when the limit cycle crosses the saddle point F 2 1 with a Andronov-Leontovich or homoclinic bifurcation [START_REF] Yu | Elements of Applied Bifurcation Theory[END_REF]. 3. For m 2 > m c 2 2 , all trajectories converge to the stable node F 1 1 where there is competitive exclusion of the second species: Fig. 9(c) shows for m 2 = 1.792 the convergence to the equilibrium F 1 1 from any positive initial condition. Figure 8: Projections of the ω-limit set in variables S , x 1 and x 2 , as a function of m 2 , reveal the occurrence of limit cycles through a super-critical Hopf bifurcation and disappearance through a homoclinic bifurcation.

(a)

x 2 x 1 S F 2 • F 0 • F * • F 2 1 • F 1 1 • (b) x 2 x 1 S F 2 • F 0 • F * • F 2 1 • F 1 1 • (c) x 2 x 1 S F 2 • F 0 • F * • F 2 1 • F 1 1 • Figure 9
: Super-critical Hopf bifurcation and homoclinic bifurcation: bi-stability, coexistence and limit cycle.

Pictures when S in > S in .

Depending on the value of λ 2 , one has the following change of stability.

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

λ 2 ∈]λ c 2 2 , λ c 1 2 [, F * is unstable. -For λ 2 ∈]λ c 1
2 , µ 0 [, F * is LES. The Jacobian matrix of reduced system (14) 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 c 1 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 c 2 2 2.085 (see Fig. 10(b)). Fig. 10(c) illustrates the corresponding real part of the complex-conjugate eigenvalues, as a function of m 2 .

Fig. 11 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.

1. For m 2 < m c 1 2 , all solutions of (2) converge to the stable focus F * . Fig. 12(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 c 1 2 , m c 2 2 [
, there is a super-critical Hopf bifurcation with appearance of limit cycles where F * becomes a saddle-focus. Fig. 12(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. 12(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. 3. For m 2 > m c 2 2 , the equilibrium F * returns to be a stable focus and one has the global convergence towards the coexistence of two species from any positive initial condition.

a) S in x 1 H F * F 2 F * F 2 F * F 2 λ 0 I λ c 2 2 λ c 1 2 µ 0 λ 1 F 0 S (b) β(m 2 ) λ1 λ2 α(m 2 ) ¢ ¢ f f x (c) α(m 2 ) m 2 m c 1 2 m c 2 2
x 2 x 1 S F 2 • F 0 • F * • (b) x 2 x 1 S F 2 • F 0 • F * • (c) S x 1 x 2

T ime

Figure 12: Coexistence: global convergence to the positive equilibrium F * or to the limit cycle. Fig. 13 shows the appearance and disappearance of limit cycles in the three-dimensional space (S , x 1 , x 2 ) through a double super-critical Hopf bifurcation 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 (2) 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

λ 0 < λ 1 < µ 0 < µ 1 or λ 0 < µ 0 < λ 1 < µ 1 .
J b = J ∩ (]0, λ b [∪]µ b , +∞[) J =] max(µ 0 , λ 1 ), µ 1 [.
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. 14(a) illustrates the interval of existence of equilibria F 1 and F * for the case λ 1 < µ 0 . Fig. 14(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. As p increases to 1, the curves f (•) and g(•) coalesce and the graph of the curve H on I and J b converges to the graphs of S = λ 0 and S = µ 0 , respectively. Moreover, the positive equilibrium F * disappears with a saddle-node bifurcation with F 1 1 , which becomes LES when F 3 1 and F 4 1 disappear. For p = 1, (a) unstable if λ 2 ∈]λ 0 , µ 0 [. Then, the system exhibits a bi-stability of F 1 1 and F 2 when λ 2 > µ 0 . Finally, 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.

F 1 1 is LES while F 2 1 , F 2 and F 0 are 14 (a) g f D ψ ϕ b I J b λ 0 λ 1 µ 0 λ b µ b µ 1 S (b) x 1 S in H F * F 2 F 1 1 F 2 1 F 3 1 F 4 1 F 0 λ 0 λ 1 µ 0 µ b λ b µ 1 S
g f D ψ ϕ b I J b λ 0 µ 0 λ 1 λ b µ b µ 1 S (b) x 1 S in H F * F 2 F 1 1 F 2 1 F 3 1 F 4 1 F 0 λ 0 µ 0 λ 1 µ b λ b µ 1 S
Fig. 15(a) illustrates the interval of existence of equilibria F 1 and F * in the case µ 0 < λ 1 . 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 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. 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. 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.

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 maximum growth rate of the second species 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. with a 11 = -ψ(S ) + f (S )u + g (S )(S in -Su), a 12 = f (S )g(S ), a 21 = f (S )ub, m 22 = -ϕ(S ) + 2au + b.

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.

Thus, it suffices to prove that E 0 = (S in , 0) is a LES equilibrium of the reduced system [START_REF] Hsu | Limiting behavior for competing species[END_REF]. From (A.2), the Jacobian matrix of ( 12) 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] Hsu | Limiting behavior for competing species[END_REF]. From expression (A. Replacing u by its expression (8), 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ϕ) ψ .
According to [START_REF] Haegeman | Modeling bacteria flocculation as density-dependent growth[END_REF], we can replace ϕ/ψ by -v/u in the function H (•) given by ( 11) and so we obtain

H = -f a 12 + m 22 aψ -g v m 22
auψ .

Multiplying (H + 1) by auψ finally yields the following expression

D 1 ( S ) = -aūψ( S ) 1 + H ( S ) . (A.6)
From the expression (A.5) of m 22 and the equations (7-8), a straightforward calculation of T 1 ( S ) := tr A 1 gives

T 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. 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.

Proof of Prop. 4.1. The equilibria of (2) are solutions of the following system of equations

                 D (S 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 .
(A.8)

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 ).

For convenience, we denote ϕ 2 (S ) = f 2 (S )-D for the second species. Let J be the Jacobian matrix of ( 14) at (S , u, x 2 ), that we write with the following notation

J =             -a 11 -a 12 -a 13 a 21 -m 22 -b m 31 0 a 33             (A.9)
where

a 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.2.

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] Jones | The Freter model: A simple model of biofilm formation[END_REF]. From (A.9), the Jacobian matrix of ( 14) at F 0 is given by

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 ( 12) 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] Jones | The Freter model: A simple model of biofilm formation[END_REF]. From (A.9), the Jacobian matrix of ( 14) at F 2 , is given by

J 2 =             ψ(λ 2 ) -f 2 (λ 2 ) x2 ψ(λ 2 ) -ϕ(λ 2 ) ψ(λ 2 ) -b ϕ(λ 2 ) -b -b f 2 (λ 2 ) x2 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 ) x2 -λ ψ(λ 2 ) -ϕ(λ 2 ) ψ(λ 2 ) 0 ϕ(λ 2 ) -b -λ -b 0 ψ(λ 2 ) -ϕ(λ 2 ) ψ(λ 2 ) -λ                   .
Moreover, the replacements of C 1 by C 1 -C 2 and L 2 by L 1 + L 2 lead to

P(λ) = -f 2 (λ 2 ) x2 -λ ϕ(λ 2 ) -λ -b 0 ψ(λ 2 ) -b -λ .
Therefore, f 2 (λ 2 ) x2 , ϕ(λ 2 ) and ψ(λ 2 )b 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] Jones | The Freter model: A simple model of biofilm formation[END_REF]. From (A.9), the Jacobian matrix of ( 14) at F 1 , is given by

J 1 =              A 1 g( S ) -f 2 ( S ) -b 0 0 ϕ 2 ( S )             
.

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 (13) holds.

Proof of Prop. 4.3. 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] Jones | The Freter model: A simple model of biofilm formation[END_REF]. Let J * denote the Jacobian matrix of reduced system [START_REF] Jones | The Freter model: A simple model of biofilm formation[END_REF] 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 (λ 2 ) > -1. From (A.6) and (A.

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.

2 .

 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 .

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.

Figure 3 :

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

1 ( 3 .

 13 4(b)) where the pair of complex-conjugate eigenvalues crosses the imaginary axis at S c in 9.911 from negative half plane to positive half plane. Note that the equilibria E 0 and E 2 i. For S in = 9.8, Fig.5(c) shows the bi-stability in the plan (S , u) with convergence to stable focus E 3 624, 3.565) or to stable node E 1 1 (0.561, 8.669): see on Fig.5(a) the convergence (after a long transient of oscillations) to E 3 1 (for the initial condition S (0) = 3.5 and u(0) = 3.5) and on Fig.5(b) the convergence to E 1 1 (for the initial condition S (0) = 3.4 and u(0) = 3.6). ii. For S in = 9.95, Fig.4(c)shows the converge to the stable node E 1 1 (0.562, 8.8) for several positive initial conditions and E3 1 becomes an unstable focus.

Figure 4 :

 4 Figure 4: Sub-critical Hopf bifurcation: (a) change of sign of T 1 ( S ) where D 1 ( S ) > 0. (b) Variation of a pair of complex-conjugate eigenvalues as S in increases. (c) The phase portrait of (12) for S in = 9.95 .

Figure 5 :

 5 Figure 5: Bi-stability of E 3 1 and E 1 1 for S in = 9.8 .

F * λ 2 ∈

 2 -S has solution S ∈ I ∪ J b S < λ 2 and condition (13) I ∪ J b and S in -λ 2 > H(λ 2 ) λ 2 ∈ I and Routh-Hurwitz condition (A.10)

(

Figure 10 : 2 Figure 11 :

 10211 Figure 10: (a) Change of asymptotic behavior for S in > S in and λ 2 ∈ I. Double super-critical Hopf bifurcation: Variation of a pair of complexconjugate eigenvalues (b) and the corresponding real part (c) as m 2 increases.S

m 2 = 1 1 SFigure 13 :

 21113 Figure 13: Occurrence and disappearance of limit cycles in the three-dimensional space (S , x 1 , x 2 ) for different values of m 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 and define

Figure 14 :

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

Figure 15 :

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

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 12ba 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 22aūψ( S ). (A.4)

  at the equilibrium F * . The characteristic polynomial is given by P(λ) = det(J * -λI) = c 0 λ 3 + c 1 λ 2 + c 2 λ + c 3 , where c 0 = -1, c 1 = -(a 11 + m 22 ), c 2 = -(m 22a 11 + a 21 a 12 + a 13 m 31 ) and c 3 = m 31 (ba 12m 22 a 13 ).According to Routh-Hurwitz criterion, F * is LES if and only ifc i < 0, i = 0, . . . , 3 c 1 c 2c 0 c 3 > 0.(A.10)From equation (A.3), we obtain c 3 = m 31 au * ψ(λ 2 )

Table 1 :

 1 Intervals of existence of positive equilibria according to the case.

Table 3 :

 3 Existence and local stability of equilibria of system (2).

  [START_REF] Grover | Resource Competition[END_REF], 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 ). 12 + m 22 = aū + P with P and -ψ( S ) are positive. Thus a 11 > 0 and so one hasT 1 ( S ) = -(a 11 + m 22 ) < 0. 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.

					(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

Table B .

 B From (A.3), we obtain m 22 a 11 + a 21 a 12 = m 22 f (λ 2 )u * + g (λ 2 )v * + m 31au * ψ(λ 2 ) + a 12 f (λ 2 )u * which is positive from(A.11). Consequently one has c 2 < 0 because a 13 m 31 is positive. Since a 12 > 0 in this case, one can divide inequality (A.11) by a 12 + m 22 and obtainf (λ 2 )u * < -au * a 12 + m 22 ψ(λ 2 ) + g (λ 2 )v * m 22 a 12 + m 22 < -ψ(λ 2 ) + g (λ 2 )v *because a 12 + m 22 = au * + P with P and -ψ(λ 2 ) are positive. Thus, a 11 > 0 and so c 1 < 0. Then, a straightforward calculation gives c 1 c 2c 0 c 3 = (a 11 + m 22 )(m 22 a 11 + a 21 a 12 ) + a 11 a 13 m 31 + ba 12 m 31 which is positive. Therefore all the conditions of the Routh-Hurwitz criterion are satisfied and thus F * is LES. Appendix B. Parameters used in numerical simulations All the values of the parameters values, for the one and two species models, respectively, are given in Tables B.4 and B.5. 4: Parameter values and corresponding values λ 0 , µ 0 , λ 1 and λ b .

	Parameter	m 11	K 11	K i	m 12	K 12	D	a	b	S in	λ 0	µ 0	λ 1	λ b
	Fig. 2(a)	2.7	1	1.2	1.5	1		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(a-b)									Variable
	Fig. 5	3	1	1.2	1.5	1	0.9	0.005	0.3	9.8	0.528	2.272	1.5	4
	Fig. 4(c)									9.95
		Parameter		m 11	K 11	K i	m 12	K 12	D	a	b	S in	m 2	K 2
		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 and 9 Figs. 10, 11, 12 and 13	3	0.5	1.2	1.25	1.5	0.9	2	2	3.5 3.65	Variable	2

[START_REF] Haegeman | Modeling bacteria flocculation as density-dependent growth[END_REF]

, we deduce that

-au * ψ(λ 2 ) + f (λ 2 )u * (a 12 + m 22 ) + g (λ 2 )v * m 22 > 0 . (A

.11)

Table B .

 B 5: Parameter values for the two species model.

are a saddle points in the Figs. 4(c) and 5(c). As the parameter S in approaches its critical value, the region of attraction of E[START_REF] Fekih-Salem | Sur un modèle de compétition et de coexistence dans le chemostat[END_REF] 1 is bounded by the unstable cycle, which shrinks and disappears in a sub-critical Hopf bifurcation.

Acknowledgments

The authors wish to thank the financial support of TREASURE euro-Mediterranean research network (https:// project.inria.fr/treasure/). This work has been achieved during the preparation of the PhD 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

and then u(t) 0 for all t 0.

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

One can write Z(t) max(S in , Z(0)) for all t 0. the solution of ( 2) is positively bounded and is defined for all t 0. From (A.1), it can deduced that the set Ω is positively invariant and is a global attractor for (2).

Proof Prop. 3.3. If (u + v)(0) > 0, then from the first equation of (3), 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 )). We recall that λ 0 < λ 1 . For S in < λ 0 , 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 ( 12) at (S , u), that we write with the following notation A = -a 11 -a 12 a 21 -m 22 (A.2)


Pierre Chatelain

Alexandre Krupa

Nassir Navab

Optimization of ultrasound image quality via visual servoing

In this paper we propose a new ultrasoundbased visual servoing framework, for the optimization of the positioning of an ultrasound probe manipulated by a robotic arm, in order to improve the quality of the acquired ultrasound images. To this end, we use the recent framework of ultrasound confidence map, which aims at estimating the perpixel quality of the ultrasound signal based on a model of sound propagation in soft tissues. More specifically, we treat the ultrasound confidence maps as a new modality to design a visual servoing control law for image quality optimization. The proposed framework aims at improving ultrasound imaging techniques, such as robotic tele-echography, target tracking or volume reconstruction. Here we illustrate our approach with the application of robotic tele-echography. Experiments are performed on both an ultrasound examination training phantom and ex vivo tissue samples.

I. INTRODUCTION

Ultrasound is a widely used modality for diagnostic and intraoperative imaging. Despite its low signal-to-noise ratio compared to other imaging techniques such as computed tomography or X-rays, ultrasound imaging is attractive because of its non-invasiveness, low cost, and high frame rate. However, the quality of ultrasound images depends on both the intrinsic properties of the device (the transducer design and various imaging parameters) and extraneous factors such as the placement of the probe, the contact force with the body, and the distribution of ultrasound gel.

The impact of the transducer design on the image quality has been extensively studied [START_REF] Goldberg | Multilayer piezoelectric ceramics for two-dimensional array transducers[END_REF], [START_REF] Kofler | Optimization of medical ultrasound transducers with simulation and genetic algorithms[END_REF], leading to significant improvements of the hardware over the last decades. On the other hand, a wide community of researchers has focused on the enhancement of ultrasound image quality via postprocessing algorithms such as denoising or speckle reduction [START_REF] Loizou | Quality evaluation of ultrasound imaging in the carotid artery based on normalization and speckle reduction filtering[END_REF], [START_REF] Tsantis | Multiresolution edge detection using enhanced fuzzy c-means clustering for ultrasound image speckle reduction[END_REF], or the use of contrast-enhanced ultrasound [START_REF] Goldberg | Ultrasound contrast agents: A review[END_REF], where ultrasound contrast agents consisting of gas-filled microbubbles are injected into the circulatory system in order to increase the contrast of the ultrasound images. A review of pre-and post-processing ultrasound image enhancement techniques can be found in [START_REF] Ortiz | Ultrasound image enhancement: A review[END_REF].

More recently, El-Zehiry et al. [START_REF] El-Zehiry | Learning the manifold of quality ultrasound acquisition[END_REF] have proposed an automatic optimization of the ultrasound imaging parameters to improve the quality of the images. The main advantage of this approach compared to conventional image enhancement techniques is that the optimization is performed on the settings of the imaging device rather than by post-processing, therefore actually improving the quality of the raw signal.

Although the placement of the probe is known to have a significant impact on the quality of the image content (see e.g. [START_REF] Klimstra | The effect of ultrasound probe orientation on muscle architecture measurement[END_REF]), very little work has been published on this topic. In this paper we investigate the relation between the probe positioning and the image quality, and design a new visual servoing technique to guide a probe attached to the end-effector of a robotic arm, in order to improve the image quality. The proposed approach is intended to be combined with other ultrasound-based visual servoing methods. It could for instance be used in the context of tracking, where improving the image quality would in turn ease the detection of the target.

Ultrasound-based visual servoing is a visual servo-control scheme where the control law is a function of features extracted from the ultrasound images. This framework has been used in different applications of robot-assisted ultrasound imaging. For instance, in [START_REF] Abolmaesumi | Visual servoing for robot-assisted diagnostic ultrasound[END_REF], Abolmaesumi et al. control the in-plane motion of a 2D probe to track a section of the carotid artery. In [START_REF] Nakadate | Out-of-plane visual servoing method for tracking the carotid artery with a robot-assisted ultrasound diagnostic system[END_REF] an intensity-based method is proposed to track the out-of-plane translation of the carotid artery. In [START_REF] Nadeau | Intensity-based ultrasound visual servoing: modeling and validation with 2D and 3D probes[END_REF] intensity-based visual servoing is used to compensate for organ motion. A method to maintain visibility constraints during tele-echography is proposed in [START_REF] Krupa | Robotized teleechography: an assisting visibility tool to support expert diagnostic[END_REF]. The main objective of robot-assisted ultrasound imaging is to ease or automatize the manipulation of the probe, which can be exhausting during prolonged examination. In this spirit, we propose an example application of the new visual servoing method designed in this paper, where the goal is to robotically assist the operator by automating the in-plane rotation of the probe in a way that improves the image quality. To this end, our visual servoing scheme relies on ultrasound quality features.

Evaluating the quality of ultrasound image is a challenging issue, because quality is a subjective concept, and its definition can vary a lot from an observer to another, or depending on the nature of the examination. Few quality measures for ultrasound have been proposed in the literature, and most of them are based on general image analysis techniques [START_REF] Loizou | Quality evaluation of ultrasound imaging in the carotid artery based on normalization and speckle reduction filtering[END_REF], [START_REF] Swathi Rangaraju | Review paper on quantitative image quality assessment -medical ultrasound images[END_REF], [START_REF] Zhao | Quality evaluation of ultrasound imaging using a MATLAB test-bench[END_REF]. In [START_REF] El-Zehiry | Learning the manifold of quality ultrasound acquisition[END_REF] the authors proposed a new quality evaluation method based on statistical learning from expert ratings. However all these methods provide only a per-image quality measure, from which it is difficult to design a control law for quality optimization. In [START_REF] Karamalis | Ultrasound confidence maps using random walks[END_REF] Karamalis et al. designed a new framework for ultrasound signal loss estimation, providing a per-pixel measure of the confidence in the ultrasound signal. The computation of Karamalis' confidence maps is based on a simplified model of sound propagation in soft tissues, and therefore provides an objective measure of quality. This framework has been successfully applied to different applications such as shadow detection and ultrasound reconstruction [START_REF] Karamalis | Ultrasound confidence maps using random walks[END_REF], or ultrasound image compounding [START_REF] Zu Berge | Orientation-driven ultrasound compounding using uncertainty information[END_REF]. Hennersperger et al. [START_REF] Hennersperger | A quadratic energy minimization framework for signal loss estimation from arbitrarily sampled ultrasound data[END_REF] have recently generalized this signal loss estimation method to arbitrarily sampled 2D and 3D data. We propose here to use the confidence maps as a new signal for servo-control, and to this end we investigate the relation between the probe motion and the variations of the confidence map.

After recalling the definition of the confidence map (section II-A), we describe in section II-B our confidence-based control law and the tele-echography framework. We present in section III the results of experiments performed on an ultrasound training phantom and ex vivo tissue samples.

II. METHODS

A. Ultrasound Confidence maps

In medical ultrasonography, focused ultrasound waves are generated and sent into the body by a piezoelectric transducer array. The waves propagating through the tissues undergo acoustic attenuation due to the medium viscosity and reflections at the boundaries between tissues of different densities. The echoes returning to the transducer are converted back into electric pulses and processed in order to form a 2D image. The resulting images have a typically low signal-tonoise ratio compared to other imaging modalities, because of the important attenuation of the ultrasound waves and the presence of artifacts due to reverberation. Therefore the quality of the signal is usually non-uniform in the image.

Ultrasound confidence maps were introduced by Karamalis et al. [START_REF] Karamalis | Ultrasound confidence maps using random walks[END_REF] as a per-pixel measure of the confidence for ultrasound images. This signal loss estimation method relies on a graphical representation of the image, where the nodes represent the pixels, and the edges model the sound propagation constraints (attenuation, reflection, refraction). The graph neighborhood can be either 4-or 8-connected. In our setup we consider a 4-connected graph. Indeed, the computation of the confidence map on a 4-connected graph provides a precision which is sufficient for our application, while complying with the real-time constraint. The computation with a 8-connected graph is more time-consuming and provides no noticeable improvement on the servoing behavior for the considered application.

Karamalis' confidence is defined as the probability of a random walk starting from a pixel to reach each of the virtual transducer elements under ultrasound-specific constraints. Thus it provides a per-pixel measure of uncertainty of the ultrasound signal, based on a model of sound wave propagation in soft tissues.

Given an ultrasound image I : Ω -→ R defined over the field of view Ω, the confidence map C : Ω -→ [0, 1] provides for each pixel (x, y) ∈ Ω the confidence C(x, y) related to the signal component I(x, y) (a confidence of 1 corresponds to a perfect signal confidence). In practice, it is computed using the random walks algorithm [START_REF] Grady | Random walks for image segmentation[END_REF], with three constraints: • The confidence at the top of the image (next to the transducer) is 1.

• The confidence at the bottom of the image is 0.

• The signal propagation obeys ultrasound-specific constraints. The reader may refer to [START_REF] Karamalis | Ultrasound confidence maps using random walks[END_REF] for the details of the propagation model. We show in Fig. 1 an example of confidence map computed from an ultrasound image acquired on an ex vivo porcine tissue sample.

B. Control law

We propose to use the confidence map as an additional modality to design a visual servoing framework for optimizing the quality of ultrasound imaging. We consider a convex ultrasound probe attached to the end-effector of a 6-DOF robotic arm. The field of view of the convex probe can be defined in polar coordinates as

Ω = [r min , r max ] × [θ min , θ max ],
where r min is the probe radius, r max is the imaging depth, and [θ min , θ max ] is the imaging sector. The confidence map C is computed in real-time using the method described in [START_REF] Karamalis | Ultrasound confidence maps using random walks[END_REF]. Our goal is to maximize the overall confidence in the image during teleoperation.

To this end, we control the in-plane rotation of the probe automatically based on the confidence map (section II-B.1). For security reasons and to maintain contact with the body, we use in addition a force control scheme (section II-B.2) for the longitudinal translation (in the direction of the image depth). The remaining 4 degrees of freedom are teleoperated by the user from a control station.

The different Cartesian frames used in our framework are represented in Fig. 2. The force sensor frame F s is attached to a 6-axes force/torque sensor, the probe frame F p to the imaging center of the convex probe, and the probe contact frame F pc to the contact point between the probe and the object's surface. These three frames are aligned along the y-axis of the probe. The probe gravity frame F g is centered on the probe's mass center and its z-axis is aligned with the direction of the gravity force.

To design our control law we refer to the visual servoing framework [START_REF] Chaumette | Visual servo control, part I. Basic approaches[END_REF], where the control screw v of the visual sensor is computed from a visual feature vector s as:

v = -λ L s + (s -s * ) (1)
where λ > 0 is the control gain, s * is the desired value of the feature vector, and L s + is the pseudo-inverse of an approximation of the interaction matrix L s relating the variation of the features to the velocity screw, such as:

ṡ = L s v (2)
1) Confidence control: We propose to define the feature for the confidence-based control law as the angular coordinate of the confidence-weighted barycenter of the image. Considering the confidence map C defined over the polar space Ω = [r min , r max ] × [θ min , θ max ] (see Fig. 2(a)), this feature can be computed from the image moments as:

s c = 1 C (θ ,r)∈Ω θC(r, θ )rdrdθ (3)
where C = (θ ,r)∈Ω C(r, θ )rdrdθ is the mean confidence of the current image, and θ = 0 • corresponds to the central scan-line in the image.

The geometric feature defined in (3) is used to compute the in-plane angular velocity ω Z p to apply to the probe. In our setup, this angular velocity corresponds to the sixth component of the velocity screw in the probe frame F p : v c = [0 0 0 0 0 ω Z p ] (see Fig. 2). As we choose to define the control law in the probe frame, the polar image coordinate θ corresponds to the rotation around Z p , so that we can assume the interaction matrix to be equal to the identity. Thus we can write the confidence control law as:

ω Z p = -λ c a c (C)(s c -s * c ) (4)
where λ c > 0 is the confidence control gain, a c : R -→ [0, 1] is an activation function, and s * c = 0 • is the desired value of the confidence-weighted angular feature. As θ = 0 • corresponds to the central scan-line, the effect of choosing s * c = 0 • is intuitively to balance the distribution of the confidence in the field of view, bringing high-confidence areas to the center of the image.

The behavior induced by this control law is twofold:

• While scanning an homogeneous region, the probe will be oriented such as the contact along the transducer is as uniform as possible (see Fig. 3). As a result, the probe will remain orthogonal to the object's surface during the navigation. • While scanning a region which contains a strongly shadowing object, such as a bone, the probe will be oriented in a way that the ultrasound beam does not go through this object (see Fig. 4). This provides a kind of obstacle-avoidance behavior during the navigation. Note that the feature s c is only consistent when a sufficient contact force is applied, i.e. when the mean confidence is high enough. Therefore we use the activation function a c to apply the confidence control only when the mean confidence is above a certain threshold. In order to ensure smoothness of the control, we use a sigmoidal activation function to activate progressively the confidence control when the mean confidence is between c -and c + :

a c (c) =      0 if c ≤ c - 1 1+exp c + -c - c-c --c + -c - c + -c if c -< c < c + 1 if c ≥ c + (5)
2) Force control: In addition to the confidence control, we apply a force control law in order to guarantee a constant contact force between the probe and the object being scanned. We recall here briefly the force control scheme as described in [START_REF] Nadeau | Intensity-based ultrasound visual servoing: modeling and validation with 2D and 3D probes[END_REF].

The robot is equipped with a 6-axes force/torque sensor providing a force tensor measurement s H s , expressed in the sensor frame F s . We consider the transformation matrices pc F s and s F g , which allow expressing in the frame F pc (resp. F s ) a force tensor known in the frame F s (resp. F g ). These transformation matrices depend on the relative position of the different frames, and are estimated experimentally. The gravity force tensor g H g is expressed in the gravity frame F g as g H g = [0 0 9.81m p 0 0 0] , where m p is the mass of the probe. Knowing these calibration parameters, the force tensor pc H pc in the contact frame F pc can be computed as:

pc H pc = pc F s ( s H s -s F g g H g ) (6)
As we aim to apply the force control only along the yaxis of the probe frame, this force tensor is multiplied by a selection matrix M s = diag(0, 1, 0, 0, 0, 0), and the desired contact force is defined as pc H * pc = [0 f * y 0 0 0 0] The resulting control law is defined in the probe frame as:

v s = -λ s M s pc H pc -pc H * pc k ( 7
)
where k is an estimation of the contact stiffness, and λ s is the force control gain.

3) Control fusion: Finally, the teleoperator controls the remaining 4 degrees of freedom (DOF) of the robotic system. The teleoperation control outputs v X p , v Z p , ω X p and ω Y p sent by the user are expressed in the probe frame F p , so that the operator controls the lateral translation and the out-of-plane motion of the probe:

v op = [v X p 0 v Z p ω X p ω Y p 0] ( 8
)
This teleoperation control is combined with the confidence and force controls, yielding the following control law for robot's end-effector:

v e = e V p (v s + v c + v op ) (9)
where e V p is the velocity transformation matrix from the probe frame F p to the end-effector frame F e . Note that the force, confidence and manual controls are decoupled, since they are performed on different axes of the probe frame.

III. RESULTS

A. Experimental setup

For all experiments we used a Viper s650 robot (Adept Technology Inc., USA), equipped with a 6-axes force/torque sensor and a SonoSite C60 5-2 MHz convex ultrasound probe (FUJIFILM SonoSite Inc., USA) rigidly attached to its end-effector. The ultrasound frames were grabbed from the SonoSite 180PLUS ultrasound system and processed on a workstation (Intel Xeon CPU @2.1 GHz), which provided user interaction for the teleoperation, and communication with the robot. The images were grabbed at a frequency of 25 Hz and processed in less than 40 ms. While the confidence and teleoperation control velocities v c and v op were updated every 40 ms, the force control was performed with a higher frequency (200 Hz), and the control law ( 9) was updated with a period of 5 ms. The experiments were performed on an ABDFAN ultrasound examination training model (Kyoto Kagaku Co., Japan) simulating a patient's abdomen, and on ex vivo porcine samples.

B. Confidence vs. contact force

In our first experiment, we studied the variation of the mean confidence value with the force contact. Although this experiment is not directly related to our control framework, we include it here as we believe it provides a first insight on the dynamic behavior of the confidence maps, and validates its use as the measure of quality. The probe was progressively put in contact with the phantom, the contact force varying from 0 N (no contact) to 5 N. The corresponding confidence values are displayed in Fig. 5. We observe that the mean confidence value increases rapidly while contact is made (between 0 and 1 N), and then varies only slightly. This is consistent with empirical observation, as we know that a sufficient contact force is necessary to obtain a good quality image.

C. Convergence

The second experiment, performed on the ABDFAN phantom, demonstrates the fitness of the chosen feature. The probe was initially positioned above the phantom, without contact, and inclined with respect to the phantom's surface with an angle of approximately 10 • . Then the control law ( 9) was used with a desired force f * y = 3 N, a desired confidencebased angular feature s * c = 0 • and no manual interaction on the 4 teleoperated DOF, so that the probe was progressively put in contact with the phantom and converged towards an optimal placement. The control gains where set empirically to λ s k = 0.001 for the force control and λ c = 0.5 for the confidence-driven control. The activation thresholds for the confidence control, defined in [START_REF] Goldberg | Ultrasound contrast agents: A review[END_REF], were set empirically to c -= 0.1 and c + = 0.5.

We present in Fig. 6(a) the evolution of the confidence feature s c and in Fig. 6(b) the control velocities in the probe frame. The system successfully converged, with a final feature error of 0.078 • . The mean image confidence and the contact force value are represented in Fig. 6(c) and (d) respectively. This experiment also shows that the overall quality of the image was improved by our confidence-based servoing method, as the mean confidence increased from 0.25 (at t = 3.4 s) to 0.58 (at t = 14 s). The measured contact force at t = 3.4 s was already 1.5 N, hence the behavior depicted in Fig. 5 suggests that this increase in confidence is mainly due to the rotation of the probe rather than to the force control.

We also performed the same experiment on an ex vivo porcine sample, in order to validate our approach on biological tissues. For this experiment the same control gains were used, but we allowed a higher force contact f * y = 4 N. We also adjusted the thresholds to the tissue properties, setting c -= 0.1 and c + = 0.4. As illustrated in Fig. 7, the system also converged successfully in this experiment, with a final feature error of 0.050 • .

D. Tele-echography

Finally we illustrate here the applicability of the proposed method in the context of teleoperation. The objective is to ease the manipulation of the probe by optimizing its orientation during the teleoperated examination. The probe was controlled according to [START_REF] Abolmaesumi | Visual servoing for robot-assisted diagnostic ultrasound[END_REF], so that 4 degrees of freedom were controlled manually by the operator. The desired contact force was set to f * y = 5 N. During the navigation, the probe orientation adapted automatically an orientation providing a good image quality, even when scanning shadowed areas such as behind the ribs. We display in Fig. the evolution of system (confidence feature value, velocity, mean confidence and contact force) during a selected period of a teleoperation experiment. We invite the reader to refer to the corresponding video supplied in supplementary material for a full demonstration of the method. A selection of frames extracted from this video is reproduced in Fig. 9. Between the frames presented in Fig. 9(c) and (d), the probe is guided by the operator towards the side of the phantom, over the ribs. The peak in feature error observed around t = 40 s in Fig. 8(a) corresponds to the instant when the probe reaches the rib cage. A peak in the contact force, due to the salient ribs, occurs at the same instant. The feature error then quickly decreases thanks to our confidence-based control scheme. As a result, although the ribs are strongly shadowing, our confidence-based visual servoing method allows the probe to be automatically re-oriented in a way that maintains the visibility of the inner tissues.

IV. CONCLUSION

We proposed a new approach for optimizing the quality of robot-assisted ultrasound imaging, based on confidencedriven visual servoing. Although we presented here an application to tele-echography, our method is also intended to be integrated into other ultrasound-based visual servoing frameworks, in order to improve the efficiency of servoing tasks. Potential applications are for example tracking a soft tissue target, or reconstructing a 3D region from 2D frames, where our method could be used as a secondary task to improve the image quality.

In this article we considered the control of the in-plane rotation. The results show that the confidence map can be used in control tasks to improve the image quality. In future work we will consider the confidence-driven servoing of other degrees of freedom, in particular out-of-plane motions, which could provide finer control of the image quality.

Fig. 1 .

 1 Fig. 1. (a) Ultrasound image of an ex vivo porcine tissue sample and (b) the corresponding confidence map.

Fig. 2 .

 2 Fig. 2. Illustration of the different coordinate systems. (a) Probe frame and image polar coordinates. (b) The Viper s650 robot's end-effector and the corresponding force sensor frame (F s ), probe frame (F p ), probe contact frame (F pc ), and gravity frame (F g ).

Fig. 3 .Fig. 4 .

 34 Fig. 3. Illustration of the system's behavior without obstacle. (a) Initial state: the probe orientation induces a low confidence (LC) in the left part of the image due to weaker contact force, while the right part has a high confidence (HC). (b) Final state: the probe is orthogonal to the object's surface, so that the confidence is balanced in the image.

Fig. 5 .

 5 Fig. 5. Relation between the mean image confidence and the contact force. The confidence measures are grouped in force bins spaced of 0.16 N. The boxes represent the first and last quartiles of the measured mean confidence for each force bin, while the whiskers represent the minimum and maximum values. Each bin contains approximately 20 measures of confidence.

Fig. 6 .

 6 Fig. 6. Phantom experiment: Evolution of the system while establishing contact with the phantom. (a) Confidence-based feature s c in degrees. (b) Control velocities: v Yp in mm/s (-) and ω Zp in • /s (---). (c) Mean image confidence C. (d) Contact force f y in N. The force control was activated from t = 0 s. The confidence control started at t = 3.1 s (C > 0.1), and was fully activated from t = 6.0 s (C > 0.5) until convergence.

8

 8

Fig. 7 .

 7 Fig. 7. Ex vivo experiment: Evolution of the system while establishing contact with the ex vivo porcine tissue sample. (a) Confidence-based feature s c in degrees. (b) Control velocities: v Yp in mm/s (-) and ω Zp in • /s (---). (c) Mean image confidence C. (d) Contact force f y in N. The force control was activated from t = 0 s. The confidence control started at t = 2.0 s (C > 0.1), and was fully activated from t = 3.5 s (C > 0.4) until convergence.

Fig. 8 .Fig. 9 .

 89 Fig. 8. Teleoperation experiment: (a) Confidence-based feature s c in degrees. (b) Evolution of the control velocity v Yp in mm/s (-), v Xp in mm/s (• • • ) and ω Zp in • /s (---). (c) Mean image confidence C. (d) Contact force f y in N.

Pierre Chatelain is with Université de Rennes 1, IRISA, France and Computer Aided Medical Procedures, Technische Universität München, Germany pierre.chatelain@irisa.fr

Alexandre Krupa is with Inria Rennes -Bretagne Atlantique, IRISA, France alexandre.krupa@inria.fr

Nassir Navab is with Computer Aided Medical Procedures, Technische Universität München, Germany and the Johns Hopkins University, MD, USA nassir.navab@in.tum.de

ACKNOWLEDGMENT

This work is partially supported by the BFS funding for project RoBildOR.