# A flocculation model in the chemostat with distinct removal rates and yields

Radhouane Fekih-Salem, Tewfik Sari

► To cite this version:

Radhouane Fekih-Salem, Tewfik Sari. A flocculation model in the chemostat with distinct removal rates and yields. Journal of Biological Systems, 2025, 10.1142/S0218339025500342. hal-04327081v3

HAL Id: hal-04327081

https://hal.science/hal-04327081v3

Submitted on 14 May 2025

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
A FLOCCULATION MODEL IN THE CHEMOSTAT WITH
DISTINCT REMOVAL RATES AND YIELDS

RADHOUANE FEKIH-SALEM

University of Tunis El Manar, National Engineering School of Tunis, LAMSIN
BP 37, Le Belvédère 1002, Tunis, Tunisia
radhouane.fekih-salem@enit.utm.tn

TEWFIK SARI

ITAP, Univ Montpellier, INRAE, Institut Agro
361 rue Jean-François Breton, 34196, Montpellier, France
tewfik.sari@inrae.fr

Received (03 July 2024)

Revised (22 April 2025)

Accepted (06 May 2025)

In this paper, a flocculation model is considered in a chemostat with a single resource and one species in two forms: planktonic and aggregated bacteria. The removal rates are distinct and include the specific death rates. Considering distinct yield coefficients with a large class of growth rates, we present a mathematical analysis of the model by establishing the necessary and sufficient conditions of existence and local asymptotic stability of all equilibria. From these conditions, we first theoretically determine the operating diagram describing the asymptotic behavior of the system according to two operating parameters. The bifurcation analysis shows a rich behavior: transcritical, saddle-node, Hopf, and homoclinic bifurcations. Using the numerical method with MATCONT software, we detect fold bifurcation points of two stable and unstable limit cycles that emerge through supercritical and subcritical Hopf bifurcations, as well as two-parameter bifurcations such as Bogdanov-Takens, Cusp, and Generalized Hopf bifurcations. Our analysis shows that the flocculation process promotes the coexistence of isolated and attached bacteria around a stable limit cycle or positive equilibrium. These theoretical messages explain and reinforce the experimental observations about the positive effect of flocculation on the maintenance of species coexistence and the protection of the most relevant species among microbial ecosystems.

Keywords: Bifurcations theory; Coexistence; Flocculation; Hopf bifurcation; Limit cycle

# 1. Introduction

The chemostat is an important laboratory apparatus used for experiments on the con-
trolled growth of microorganisms in microbiology and ecology. It has played an important
role in many fields, such as the wastewater treatment process, biomass energy recovery,
and biotechnologies in a broad sense. Mathematical models of competition on a single
limiting nutrient in a chemostat have played a central role in microbial ecology, micro-
biology, and evolutionary and applied biology. The mathematical study of the classical
chemostat model of several species competing on the same limiting resource can be
found in the monograph by Smith and Waltman.¹ They have shown that only the most
2 Fekih-Salem and Sari

competitive species that consumes less substrate to reach its steady state (or that has the lowest break-even concentration) survives the competition of several species on a single nutrient, while all other species are excluded. This result is well known as the Competitive Exclusion Principle (CEP), which states that two species competing for identical limited resources cannot coexist indefinitely. However, the CEP contradicts the biodiversity observed in nature and microbial ecosystems.

In order to reconcile the mathematical results of the classical chemostat model as-
serting the CEP and the experimental results and nature showing the biodiversity of
microbial species, various recent studies have revised the mathematical modeling of the
competition of several microbial species competing on a single resource. More specifi-
cally, a lot of research has tried to understand and explain the biodiversity in microbial
ecosystems by analyzing the various types of interactions favoring the coexistence of
microbial species. In 2, the constant input of some species in a chemostat of n species
competing on a single nutrient can lead to coexistence. In the literature, we can cite
these various mechanisms of coexistence: flocculation 3–7, intra- and interspecific inter-
ference 8,9, density-dependence 10–15, presence of internal or external inhibitors 16–20,
predator-prey interaction 21, 22, simple or complex food web 23–25, and the references
therein.

In this paper, we consider a flocculation model of one microbial species that is
decomposed into isolated (or planktonic) bacteria and attached (or aggregate) bacteria
with a single nutrient *S* in a chemostat. This distinction allows us to take into account
different growth and death characteristics according to whether the microorganisms are
attached or not. Moreover, isolated bacteria of density *u* can aggregate with isolated
bacteria or flocs of density *v* to form new flocs, with a rate *a*(*u* + *v*)*u*, while flocs can
split and liberate isolated bacteria, with a rate *bv*.

The phenomena of bacterial attachment and detachment, in the form of biofilms on
a support or in the form of aggregates or flocs, are well known and frequently observed in
bacterial growth in a chemostat.<sup>26-28</sup> Freter's model <sup>29,30</sup> of the gut bacterial ecosystem
explicitly distinguishes planktonic biomass from attached biomass and takes into account
specific attachment and detachment terms. Several works in the literature consider ex-
tensions of the chemostat model with attachment to a wall by <sup>31-33</sup>, and the formation
of flocs <sup>28,34</sup>, as well as simplified models that distinguish only two biomass compart-
ments, that which is planktonic and that which is aggregated in the form of flocs.<sup>3,7,35</sup>
The purpose of the present work is to generalize the existing results concerning these
simplified models.

In 7, the mass densities $u = u_1 m$ and $v = 2u_2 m$, where $m$ is the mass of bacteria, $u_1$ is the density of planktonic bacteria, and $u_2$ is the density of flocs of size 2 (bacteria). Hence, the specific velocity of attachment of free biomass is given by $au^2$ and corresponds to the aggregation of two planktonic bacteria to form a floc of size 2. The specific detachment velocity is given by $bv$ and corresponds to the breakage of a floc into planktonic bacteria. The case of flocs with an arbitrary number of bacteria was studied in 34 using Smoluchowski's modeling of flocculation.<sup>28</sup> In this work, we will not consider the size or the number of bacteria in flocs. Indeed, we simply distinguish the total biomass $v = (2u_2 + 3u_3 + \dots) m$ in flocs and the planktonic biomass $u = u_1 m$, where $u_i, i \ge 2$ is the density of flocs of sizes $i$ (bacteria). For simplicity, we assume that the coefficient of proportionality is also equal to $a$ and we take the specific velocity of attachment equal to $a(u+v)$ instead of the more general choice $au_u + av_v$. This choice of a linear attachment term was first proposed in 3 and was also considered in.<sup>35,36</sup> For more details on the modeling of this flocculation mechanism, see.<sup>5</sup> The model is given by the following three-dimensional system of ordinary differential equations

$$
\begin{equation}
\begin{cases}
\dot{S} = D(S_{in} - S) - \frac{1}{y_u} f(S)u - \frac{1}{y_v} g(S)v \\
\dot{u} = [f(S) - D_u]u - a(u + v)u + bv \\
\dot{v} = [g(S) - D_v]v + a(u + v)u - bv
\end{cases}
\tag{1.1}
\end{equation}
$$

where *S*(*t*) is the concentration of the substrate at time *t*; *f*(*S*) and *g*(*S*) represent, respectively, the growth rates of isolated and attached bacteria; *D* and *S*<sub>in</sub> are, respec-
A flocculation model in the chemostat with distinct removal rates and yields 3

tively, the dilution rate and the concentration of the substrate in the feed device; $D_u$ and $D_v$ represent, respectively, the disappearance rates of planktonic and attached bacteria; $y_u$ and $y_v$ are “yield” constants reflecting the conversion of nutrient to planktonic and aggregated bacteria, respectively. This model was studied in 35, 36 in the case where the yields are equal and the removal rates are equal to the dilution rate. In this case, the yields can be normalized to 1, and the system can be reduced to a planar system. Moreover, because of the structure of model (1.1), when the yields are distinct, they cannot be normalized to 1 by the usual change of variable where $u$ and $v$ are replaced by $u/y_u$ and $v/y_v$, respectively. The model (1.1) was also studied in 3, 5, 6 in the case where the yields are equal, and are normalized to 1. However, the consideration of these yield coefficients is very important in the mathematical models of the chemostat to model reproduction by nutrient uptake as mentioned in.¹, ³⁶ Consequently, the present work fills this gap by examining the general case, which is the most relevant from a biological point of view. In this work, we study model (1.1) where $D_u$ and $D_v$ can be modeled as in 37 by:

$$D_u = \alpha D + m_u, \quad D_v = \beta D + m_v$$

where the non-negative parameters $m_u$ and $m_v$ representing the mortality rate are taken into consideration.

In 5, we have determined the existence and local stability of all steady states of system (1.1) with the same yield coefficients $y_u$ and $y_v$. The model presents a multiplicity of positive steady states that can only appear or disappear through saddle-node or transcritical bifurcations. Under the joint effect of flocculation and mortality, the coexistence steady states may destabilize via a supercritical Hopf bifurcation with the emergence of a stable limit cycle that can disappear through a homoclinic bifurcation. However, the study of bifurcations is limited to the one-parameter diagrams by fixing $D$ and varying $S_{in}$. In 6, the theoretical study of the operating diagram of model (1.1) with the same yields coefficients shows that the system can exhibit bistability between the washout steady state $E_0$ and the coexistence steady state $E_1$. There may also only be coexistence around the positive steady state $E_1$. The construction of the operating diagram of model (1.1) in 6 has omitted the existence of the region of destabilization of a positive steady state where there can be the emergence of a stable limit cycle via a Hopf bifurcation for very small values of $D$ as demonstrated in Appendix E.

Indeed, the operating diagram is a very useful tool to visualize and summarize the asymptotic behavior of a process according to the operating parameters, which are the most easily manipulated parameters in a chemostat as explained in.¹, ³⁶ In the existing literature, the study of the operating diagram can be purely numerical. By exploring the set of operating parameters with a certain discretization step, a significant steady state (i.e. with nonnegative components) is determined by solving numerically algebraic equations giving steady states. Their asymptotic behaviors are established by solving the characteristic polynomial and the sign of its roots.³⁸ This method can be applied to complex processes with a large number of state variables and parameters.³⁹⁻⁴³

Another numerical alternative consists in constructing the boundaries of the various regions of the operating diagram using a continuation and correction algorithm. Various software packages have been developed to determine the values of the critical parameters corresponding to the different types of bifurcations for autonomous dynamic systems. The most used software packages are MATCONT, CONTENT, AUTO, and XPPAUT (see 44 and the reference therein).

However, the theoretical determination of the operating diagram consists in constructing the boundaries of the different regions from the theoretical analysis of the dynamic system. More precisely, using a scientific numerical platform (like MAPLE 45), these boundaries are drawn from the conditions of existence and stability of all steady states according to the operating parameters when all other parameters are fixed.⁶, ⁸, ¹⁰, ¹⁴⁻²⁰, ⁴⁴, ⁴⁶⁻⁵⁵ Note that the single-parameter or two-parameter bifurcation diagrams obtained with MATCONT 56 allow additional phenomena to be detected (such as homoclinic, Cusp, and Bogdanov-Takens bifurcations) compared with those obtained
4 Fekih-Salem and Sari

theoretically from the existence and stability conditions.

Our main objective in this paper is to establish an in-depth mathematical study of model (1.1) by considering distinct yields and removal rates, including the specific death rates. Moreover, we establish theoretically and numerically the operating diagram of model (1.1) to describe the asymptotic behavior of the process according to the two operating parameters $S_{in}$ and $D$. The operating diagrams show the emergence of a region of destabilization of the interior steady state via a Hopf bifurcation with coexistence around a stable limit cycle. The one- and two-parameter diagrams are also obtained by the numerical continuation method using MATCONT software 56, which allowed us to detect other types of bifurcations according to two parameters. However, a synthetic biologist would manipulate exactly the species parameters, e.g. the growth rates, or the adhesion and dissociation rates $a$ and $b$, and it would also be interesting to study the operating diagrams of these parameters. The study of the operating diagram in $(a, b)$-plan when $(S_{in}, D)$ is fixed shows that the process presents new phenomena such as the emergence of unstable limit cycles through subcritical Hopf bifurcations. They can disappear by a Limit Point of Cycles (LPC) or saddle-node bifurcation point of two stable and unstable cycles. The effect of attachment and detachment on the operating diagram shows the importance of considering the phenomenon of flocculation as a coexistence mechanism in the classic chemostat model.

This paper is organized as follows. First, we present in Section 2 a general hypothesis about the growth functions of the flocculation model (1.1). Then we determine the existence and local stability conditions of all steady states according to operating parameters. In Section 3, we theoretically analyze the operating diagram according to the input concentration of the substrate $S_{in}$ and the dilution rate $D$. First, in Subsection 3.1, a simple case is considered where there is only a Branch Point (BP) and no Limit Point (LP) or Hopf bifurcation. In Subsection 3.2, a case with LP and Hopf bifurcations is considered. In Subsection 3.3, another case with LP and Hopf bifurcations is considered, but a new region of instability of two positive steady states emerges in the operating diagram. In Section 4, we study numerically the one- and two-parameter bifurcation diagrams using MATCONT software for the two cases in Subsections 3.2 and 3.3. In Subsection 4.1, we analyze the operating diagram in the $(S_{in}, D)$-plane of Section 3.2 as well as its corresponding one-parameter bifurcation diagram in $S_{in}$ when $D = 0.1$ fixed. In Subsection 4.2, we analyze the operating diagram in the $(S_{in}, D)$-plane of Section 3.3. In Subsection 4.3, we analyze the one-parameter bifurcation diagram in $S_{in}$ when $D = 0.1$ fixed in the operating diagram of Section 3.3. In Subsection 4.4, we analyze the one-parameter bifurcation diagram in $D$ when $S_{in} = 3.845$ fixed in the operating diagram of Section 3.3.

In Section 5, we study the effect of flocculation and deflocculation on the operating diagrams for the set of parameters considered in Sections 3.3 but where $a$ and $b$ are variable. In Subsection 5.1, we analyze the operating diagram in the $(S_{in}, D)$-plane for different values of $a$ and $b$, which tend to zero. In Subsection 5.2, we analyze the operating diagram in the $(a, b)$-plane and the one-parameter bifurcation diagrams in $a$ or $b$ when all other parameters are those in the operating diagram of Subsection 5.1. Finally, conclusions are drawn in the last Section 6. The proof of the existence and local stability conditions of all steady states is reported in Appendix A. In Appendix B, we show that a stability condition of the positive steady state holds for the set of parameters considered in Section 3.1. In Appendix C, we illustrate that this stability condition does not hold for the set of parameters considered in Section 3.2. In Appendix D, we show the destabilization of a positive steady state and then illustrate the stable limit cycles in the three-dimensional space $(S, u, v)$ and their disappearance via a homoclinic bifurcation for the set of parameters considered in Section 3.3. In Appendix E, we show that the region of destabilization of the positive steady state is omitted in the construction of the operating diagram in.⁶ All the values of parameters used throughout this paper are provided in Appendix F.
## 2. Hypothesis and model analysis

In this paper, we make the following general assumption on the growth functions $f(S)$ and $g(S)$, which are continuously differentiable ($C^1$).

(H1) $f(0) = g(0) = 0$ and $f'(S) > 0$ and $g'(S) > 0$ for all $S > 0$.

Assumption (H1) means that no growth can occur for isolated bacteria u and attached bacteria v without the presence of the substrate S. Moreover, the growth rates of isolated and attached bacteria increase with the concentration of the substrate S.

In this section, we summarize the main results of the existence and stability of all steady states of system (1.1). A steady state exists if and only if all its components are nonnegative. This predicts two types of steady states, labeled as follows:

* $E_0$ ($u = 0, v = 0$): the washout of planktonic and attached bacteria.

* $E_1$ ($u > 0, v > 0$): both planktonic and attached bacteria are present.

To determine these steady states, we define the following auxiliary functions

$$H(S) := \frac{1}{y_u} f(S)U(S) + \frac{1}{y_v} g(S)V(S) \quad (2.1)$$

where

$$U(S) := \frac{\varphi(S)(\psi(S) - b)}{a[\psi(S) - \varphi(S)]} \quad \text{and} \quad V(S) := -\frac{\varphi^2(S)(\psi(S) - b)}{a[\psi(S) - \varphi(S)]\psi(S)} \quad (2.2)$$

and

$$\varphi(S) := f(S) - D_u \quad \text{and} \quad \psi(S) := g(S) - D_v. \quad (2.3)$$

In addition, we need to define the following interval of existence of the positive steady states:

$$I =]\lambda_u, \lambda_v[ \quad \text{if} \quad \lambda_u < \lambda_v, \quad \text{else} \quad I =]\lambda_v, \min(\lambda_u, \lambda_b)[ \quad (2.4)$$

where $\lambda_u, \lambda_v$ and $\lambda_b$ are defined in Table 3. For convenience, we shall use the abbreviation LES for Locally Exponentially Stable. Any reference to steady state stability should be considered as local exponential stability, that is to say, the real parts of the eigenvalues of the Jacobian matrix are negative. To determine the stability of the positive steady state $E_1 = (S^*, u^*, v^*)$, we define the Routh-Hurwitz coefficients by

$$\begin{align}
c_1 &= m_{11} + m_{22} + m_{33}, \nonumber \\
c_2 &= m_{12}m_{21} + m_{13}m_{31} - m_{32}a_{23} + m_{11}m_{22} + m_{11}m_{33} + m_{22}m_{33}, \nonumber \\
c_3 &= m_{11}(m_{22}m_{33} - m_{32}a_{23}) + m_{21}(m_{12}m_{33} + m_{32}m_{13}) + m_{31}(m_{12}a_{23} + m_{13}m_{22}) \tag{2.5} \\
c_4 &= c_1c_2 - c_3. \nonumber
\end{align}$$

where

$$\left\{
\begin{array}{l}
m_{11} = D + \frac{1}{y_u} f'(S^*) u^* + \frac{1}{y_v} g'(S^*) v^*, \quad m_{12} = \frac{1}{y_u} f(S^*), \quad m_{13} = \frac{1}{y_v} g(S^*), \\
m_{21} = f'(S^*) u^*, \quad m_{22} = a(2u^* + v^*) - \varphi(S^*), \quad a_{23} = b - au^*, \\
m_{31} = g'(S^*) v^*, \quad m_{32} = a(2u^* + v^*) \quad \text{and} \quad m_{33} = b - au^* - \psi(S^*).
\end{array}
\right.
\quad (2.6)$$

Now, we can state the main result, which establishes the components of all steady states of (1.1) and their existence and local asymptotic stability conditions.

**Theorem 2.1** Assume that Hypothesis (H1) holds. The steady states of (1.1) and the necessary and sufficient conditions of existence and local asymptotic stability are given in Tables 1 and 2, respectively.

The proof is given in Appendix A. It was shown in 36 (see also 35) that when $D_u = D_v = D$ and $Y_u = Y_v$, then the positive steady state $E_1$ exists and is unique and LES if and only if $S_{in} > \lambda_u$.
6 Fekih-Salem and Sari

<caption>Table 1. Steady states of (1.1). The functions $H(S)$, $U(S)$ and $V(S)$ are defined by (2.1) and (2.2).</caption>

$$ \begin{array}{l} \text{S, u, v components} \\ E_0 \quad S = S_{in}, u = 0, v = 0 \\ E_1 \quad S^* \text{ solution of equation } D(S_{in} - S) = H(S), u^* = U(S^*) \text{ and } v^* = V(S^*) \end{array} $$

<table><caption>Table 2. Necessary and sufficient existence and local stability conditions of steady states of (1.1) where $c_4$ is defined by (2.5).</caption><thead><tr><th></th><th>Existence conditions</th><th>Stability conditions</th></tr></thead><tbody><tr><td>$E_0$</td><td>always exists<br>equation $D(S_{in} - S) = H(S)$</td><td>$S_{in} < \min(\lambda_u, \lambda_b)$<br>$c_3 = \varphi(S^*)(b - \psi(S^*)) (D + H'(S^*)) > 0$</td></tr><tr><td>$E_1$</td><td>has a solution $S^* \in I$</td><td>and $c_4 > 0$</td></tr></tbody></table>

### 3. Operating diagrams

In this section, we theoretically study the operating diagrams of model (1.1) to determine the various qualitative asymptotic behaviors of the process according to the operating parameters, which are the concentration of substrate in the feed bottle $S_{in}$ and the dilution rate $D$. Each region of the diagram is characterized by a different color according to the number of existing steady states and their various asymptotic behaviors.

From definition of $\lambda_v(D)$ and $\lambda_b(D)$ in Table 3, we have $\lambda_v(D) < \lambda_b(D)$ for all $D \in [0, \bar{D}_v]$. To construct theoretically the operating diagram of system (1.1) by determining the various curves, we define the auxiliary functions according to the dilution rate $D$ in Table 3 and the set of curves $\Gamma_i, i = \{u, v, \text{BP}, \text{LP}, \text{H}\}$ in Table 4, which are the boundaries of different regions of the $(S_{in}, D)$-plane. As in this work, the construction of the operating diagram will be done with the specific growth rates of Monod-type (3.1) satisfying hypothesis (H1) and we know that in this case the function $S \mapsto H(S)$ is convex, we can define the functions $D \mapsto S_{LP}(D)$ and $D \mapsto \lambda_{LP}(D)$ in Table 3 by following.⁶

<table><caption>Table 3. Notations, auxiliary functions, and their domains of definition.</caption><thead><tr><th></th><th>Definition</th></tr></thead><tbody><tr><td>$\lambda_u(D)$</td><td>$\lambda_u(D) = f^{-1}(\alpha D + m_u)$. It is defined for $0 \le D < \bar{D}_u := (f(+\infty) - m_u)/\alpha$.</td></tr><tr><td>$\lambda_v(D)$</td><td>$\lambda_v(D) = g^{-1}(\beta D + m_v)$. It is defined for $0 \le D < \bar{D}_v := (g(+\infty) - m_v)/\beta$.</td></tr><tr><td>$\lambda_b(D)$</td><td>$\lambda_b(D) = g^{-1}(\beta D + m_v + b)$. It is defined for $0 \le D < \bar{D}_b := (g(+\infty) - m_v - b)/\beta$.</td></tr><tr><td>$\lambda_{\text{BP}}(D)$</td><td>$\lambda_{\text{BP}}(D) = \min(\lambda_u(D), \lambda_b(D))$. It is defined for $0 \le D < \max(\bar{D}_u, \bar{D}_b)$.</td></tr><tr><td>$S = S_{\text{LP}}(D)$</td><td>$S = S_{\text{LP}}(D)$ is the unique solution of equation $H'(S) = -D$ on $]\lambda_v(D), \lambda_{\text{BP}}(D)]$. It is defined for $\bar{D} \le D < \bar{D}_v$ where $\bar{D}$ is the unique solution in $]0, \max(\bar{D}_u, \bar{D}_b)[$ of equation $H'(\lambda_{\text{BP}}(D)) = -D$.</td></tr><tr><td>$\lambda_{\text{LP}}(D)$</td><td>$\lambda_{\text{LP}}(D) = H(S_{\text{LP}}(D))/D + S_{\text{LP}}(D)$. It is defined for $\bar{D} \le D < \bar{D}_v$.</td></tr></tbody></table>

The passage through the $\Gamma_{\text{BP}}$ curve corresponds to a transcritical bifurcation or Branch Point (BP) between $E_0$ and $E_1^1$ or between $E_0$ and $E_1^2$ as we will see in the
<table><caption>Table 4. Definitions of the curves $\Gamma_i$, $i = \{u, v, BP, LP, H\}$, and the corresponding colors where all functions $\lambda_i$ and $c_4$ are defined in Table 3 and (2.5), resp. The abbreviations BP, LP, and H mean a Branch Point, Limit Point, and Hopf bifurcations, respectively.</caption><thead><tr><th>Curves</th><th>Color</th><th>Bifurcation</th></tr></thead><tbody><tr><td>$\Gamma_u = \{(S_{in}, D) : S_{in} = \lambda_u(D)\}$</td><td>Red</td><td>BP</td></tr><tr><td>$\Gamma_b = \{(S_{in}, D) : S_{in} = \lambda_b(D)\}$</td><td>Blue</td><td>BP</td></tr><tr><td>$\Gamma_{BP} = \{(S_{in}, D) : S_{in} = \lambda_{BP}(D)\}$</td><td>Red or Blue</td><td>BP</td></tr><tr><td>$\Gamma_{LP} = \{(S_{in}, D) : S_{in} = \lambda_{LP}(D)\}$</td><td>Green</td><td>LP</td></tr><tr><td>$\Gamma_H = \{(S_{in}, D) : c_4(S_{in}, D) = 0\}$</td><td>Magenta</td><td>H</td></tr></tbody></table>

following section. As we shall see later, passing through the curve $\Gamma_{LP}$ in the operating plan $(S_{in}, D)$ gives rise to the two positive steady states $E_1^1$ and $E_1^2$ via a Limit Point (LP) or saddle-node bifurcation. In addition, the passage through the curve $\Gamma_H$ corresponds to a Hopf bifurcation with the appearance or disappearance of a stable limit cycle.

To illustrate the operating diagram of model (1.1), we choose the following specific growth rates of Monod-type satisfying hypothesis (H1):

$$f(S) = \frac{m_1 S}{k_1 + S} \quad \text{and} \quad g(S) = \frac{m_2 S}{k_2 + S}, \qquad (3.1)$$

where $m_i$ is the maximum growth rate and $k_i$ is the Michaelis-Menten constant, $i = 1, 2$. The values of these biological parameters are provided in Table 20. In addition, the construction of the operating diagram is similar to any other specific growth rate satisfying hypothesis (H1).

In the next section, we start with the simplest case where the operating diagram does not present the regions of destabilization of the positive steady state and the emergence of two positive steady states. Then, we study the general case with the emergence of the limit cycle and BP bifurcation. Then, we find these results using the numerical continuation method with the MATCONT software. Finally, we determine the effect of flocculation on the appearance and disappearance of various regions.

### 3.1. A case where the positive steady state is unique and stable if it exists

In this section, we consider a case where there is only BP of two steady states and no LP or Hopf bifurcation. Therefore, the positive steady state $E_1$ is unique and stable if it exists and can only bifurcate with the washout steady state $E_0$. For this purpose, we consider the biological parameter values that were used in [5, Fig. 13] (see Table 20, line 1). However, $S_{in}$ and $D$ are variable and not fixed as in 5 where $S_{in} = 5$ and $D = 3.5$. The next result describes the operating diagram in the case of Table 20 (line 1).

**Proposition 3.1** For the specific growth rates $f$ and $g$ defined in (3.1) and the set of biological parameter values in Table 20 (line 1), we have $\lambda_u(D) < \lambda_v(D)$ for all $D \in [0, \bar{D}_u[$. In addition, the existence and local stability of the steady states $E_0$ and $E_1$ of model (1.1) in the two regions $I_0$ and $I_1$ of the operating diagram shown in Fig. 1(a) are described in Table 5.

**Proof.** With this set of parameters, we have $\lambda_u(D) < \lambda_v(D)$ for all $D$ in their definition domain so that the function $S \mapsto H(S)$ is defined and increasing on $[\lambda_u(D), \lambda_v(D))$ as shown in Fig. 13(a). Hence, the positive steady state $E_1$ is unique if it exists. Since $\varphi(S^*)$ is positive, $\psi(S^*)$ is negative, and $H'(S^*) > -D$, it follows that $c_3$ is positive from the expression of $c_3$ defined in Table 2. Moreover, Appendix B shows that the stability condition $c_4 > 0$ holds for all $S_{in}$ and $D$ in the existence domain of $E_1$. Therefore, the curve $\Gamma_u$, defined in Table 4 divides the $(S_{in}, D)$-plane of the operating parameters in
<figure><img src="image_2.png" /><figcaption>Fig. 1. Operating diagram of (1.1). (a) The case considered in Section 3.1. (b) The case considered in Section 3.2. (c) Magnification of (b) showing the curve $\Gamma_H$.</figcaption></figure>

<table><caption>Table 5. Existence and local stability of steady states according to the regions in the operating diagram of Fig. 1(a). The letter S [resp. U] means stable [resp. unstable]. No letter means that the steady state does not exist.</caption><thead><tr><th>Condition</th><th>Region</th><th>Color</th><th>E<sub>0</sub></th><th>E<sub>1</sub></th></tr></thead><tbody><tr><td>S<sub>in</sub> &lt; &lambda;<sub>u</sub>(D)</td><td>I<sub>0</sub></td><td>Cyan</td><td>S</td><td></td></tr><tr><td>S<sub>in</sub> &gt; &lambda;<sub>u</sub>(D)</td><td>I<sub>1</sub></td><td>Red</td><td>U</td><td>S</td></tr></tbody></table>

two regions $I_0$ and $I_1$ shown in Fig. 1(a). From Table 2, the existence and local stability conditions of the steady states are shown in Table 5. □

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

Note that the cyan region $I_0$ in the operating diagram in Fig. 1(a) corresponds to the washout of isolated and attached bacteria while the red region $I_1$ corresponds to the coexistence of both species around a steady state. The theoretical study of the operating diagram determines the asymptotic behavior of the solutions for the set of biological parameters in [5, Fig. 13]. Note that this figure illustrates the trajectories over time that converge towards the positive steady state $E_1$ for various initial conditions when the point of coordinates (5, 3.5) in the $(S_{in}, D)$-plane belongs to region $I_1$. Thus, this operating diagram presents a global vision of the behavior of the solutions according to the two operating parameters.

## 3.2. A case with multiplicity of positive steady states and Hopf bifurcations

In this section, we consider a case where there are LP and Hopf bifurcations. In this example, we can have two positive steady states, one being stable and the other unstable, and, in addition, the stable one can be destabilized through a Hopf bifurcation. For this purpose, we consider the biological parameter values that were used in [5, Fig. 12] (see Table 20, line 2). However, $S_{in}$ and $D$ are variable and not fixed as in 5 where $D = 0.1$. The next result describes the operating diagram in the case of Table 20 (line 2).

**Proposition 3.2** For the specific growth rates $f$ and $g$ defined in (3.1) and the set of biological parameter values in Table 20 (line 2), we have $\lambda_v(D) < \lambda_{BP}(D)$ for all $D \in [0, \bar{D}_v[$. In addition, the existence and local stability of the steady states $E_0$, $E_1^1$ and $E_1^2$ of model (1.1) in the four regions $I_i$, $i = 0, \dots, 3$ of the operating diagram shown in Fig. 1(b-c) are described in Table 6.

**Proof.** With this set of parameters, we have $\lambda_v(D) < \lambda_{BP}(D) = \lambda_b(D)$ for all $D \in [0, \bar{D}_v[$ so that the function $S \mapsto H(S)$ is defined and decreasing on $]\lambda_v(D), \lambda_b(D)]$ as shown in Fig. 14. Hence, there can be up to two positive steady states denoted by $E_1^1$ and
A flocculation model in the chemostat with distinct removal rates and yields 9

<table><caption>Table 6. Existence and local stability of steady states according to the regions in the operating diagram of Fig. 1(b-c).</caption><thead><tr><th>Condition 1</th><th>Condition 2</th><th>Region</th><th>Color</th><th>E<sub>0</sub></th><th>E<sub>1</sub><sup>1</sup></th><th>E<sub>1</sub><sup>2</sup></th></tr></thead><tbody><tr><td>S<sub>in</sub> &lt; min(λ<sub>LP</sub>(D), λ<sub>BP</sub>(D))</td><td></td><td>I<sub>0</sub></td><td>Cyan</td><td>S</td><td></td><td></td></tr><tr><td>λ<sub>BP</sub>(D) &lt; S<sub>in</sub></td><td>c<sub>4</sub>(S<sub>in</sub>, D) &gt; 0</td><td>I<sub>1</sub></td><td>Red</td><td>U</td><td>S</td><td></td></tr><tr><td>λ<sub>LP</sub>(D) &lt; S<sub>in</sub> &lt; λ<sub>BP</sub>(D)</td><td>c<sub>4</sub>(S<sub>in</sub>, D) &gt; 0</td><td>I<sub>2</sub></td><td>Green</td><td>S</td><td>S</td><td>U</td></tr><tr><td>λ<sub>BP</sub>(D) &lt; S<sub>in</sub></td><td>c<sub>4</sub>(S<sub>in</sub>, D) &lt; 0</td><td>I<sub>3</sub></td><td>Blue</td><td>U</td><td>U</td><td></td></tr></tbody></table>

$E_1^2$ so that $H'(S^*) < -D$ and $H'(S^*) > -D$, respectively. They emerge in the positive quadrant via a Limit Point (LP) bifurcation. Thus, the corresponding curve $\Gamma_{LP}$ exists and is unbounded (see Fig. 1(b-c)). In this case $\lambda_v(D) < \lambda_{BP}(D) = \lambda_b(D)$, $\varphi(S^*) < 0$ and $0 < \psi(S^*) < b$. From Table 2, it follows that $c_3 > 0$ [resp. $c_3 < 0$] for $E_1^1$, [resp. $E_1^2$]. Thus, $E_1^2$ is unstable whenever it exists. In Appendix C and Fig. 14, we show that the stability condition $c_4 > 0$ of $E_1^1$ does not always hold for all $S_{in}$ and $D$ in the domain of existence of $E_1^1$ so that the curve $\Gamma_H$ exists and is unbounded. Note that $c_4$ is a function of $S$, that is, $c_4 = c_4(S)$ because it depends on the three state variables $S$, $u = U(S)$ and $v = V(S)$ defined in (2.2). Moreover, $c_4$ is a function of $(S_{in}, D)$, that is, $c_4 = c_4(S_{in}, D)$ because we can determine $S$ from the equation $D(S_{in} - S) = H(S)$. With the set of parameters in Table 20 (line 2), we provide numerical evidence in Appendix C of the change in sign of the function $c_4(S)$ on the existence interval of the positive steady state $E_1^1$ according to $D$. Indeed, for $D < D_H^{max} \approx 0.165$ and fixed, the function $c_4(S)$ changes sign in $]\lambda_v(D), \lambda_{BP}(D)]$ so that the equation $c_4(S) = 0$ has two solutions noted by

$$S_H^2(D) < S_H^1(D) < \lambda_{BP}(D).$$

For all $D < D_H^{max}$, we define the two solutions of the equation $c_4(S_{in}, D) = 0$ by the following critical values of $S_{in}$ which corresponds to a Hopf bifurcation

$$S_{in} = S_{in}^{H}(D) = \frac{1}{D} H(S_{H}^{i}(D)) + S_{H}^{i}(D), \quad i = 1, 2.$$

Consequently, the $\Gamma_H$ curve of equation $c_4(S_{in}, D) = 0$ is given by the union of the two curves $\Gamma_H^1$ (on the left of the maximum) of equation $S_{in} = S_{in}^{H1}(D)$ and $\Gamma_H^2$ (on the right of the maximum) of equation $S_{in} = S_{in}^{H2}(D)$, see Fig. 1(c).

To construct the operating diagram in Fig. 1(b-c), we plot the various curves $\Gamma_i$, $i = \{\text{BP, LP, H}\}$ defined in Tables 4 and corresponding to the existence and stability conditions of all steady states provided in Table 2. These curves divide the $(S_{in}, D)$-plane of operating parameters in four regions $I_0$ to $I_3$ as shown in Fig. 1(b-c). From Table 2, the existence and local stability conditions of the steady states are shown in Table 5. □

Note that the green region $I_2$ in the operating diagram of Fig. 1(b-c) corresponds to the bistability with either the coexistence around a steady state or the washout of the isolated and attached bacteria according to the initial condition. The blue region $I_3$ corresponds to the instability of the positive steady state $E_1^1$, where there can be coexistence around a stable limit cycle.

**Remark 3.1** In Appendix E, we establish the operating diagram with the parameter set in.⁶ It is similar to that in Fig. 1(b-c) where we find the same regions in Table 6 (just $\lambda_{BP}(D)$ is equal to $\lambda_u(D)$ instead of $\lambda_b(D)$). However, the region $I_3$ of destabilization of the positive steady state $E_1^1$ with the appearance of a stable limit cycle was not detected in 6 because of the order of magnitude of $D_H^{max}$ as demonstrated in Appendix E.
### 3.3. Another case with multiplicity of positive steady states and Hopf bifurcations

In the operating diagram of the case considered in Section 3.2, we do not have a region where the two positive steady states are both unstable, see Table 6. Thus, this section aims to provide an example where there is a new region (labeled $I_4$, see Table 7) of instability of the two positive steady states. For this purpose, we consider the biological parameter values that were used in [5, Fig. 6] (see Table 20, line 3). However, $S_{in}$ and $D$ are variable and not fixed as in 5 where $D = 0.1$. Indeed, the study in 5 was limited to one-parameter bifurcation diagrams according to $S_{in}$. The next result describes the operating diagram in the case of Table 20 (line 3).

**Proposition 3.3** For the specific growth rates $f$ and $g$ defined in (3.1) and the set of biological parameter values in Table 20 (line 3), we have $\lambda_v(D) < \lambda_{BP}(D)$ for all $D \in [0, \bar{D}_v[$. In addition, the existence and local stability of the steady states $E_0, E_1^1$, and $E_1^2$ of model (1.1) in the five regions $I_k$, $k = 0, \dots, 4$ of the operating diagram shown in Fig. 2 are described in Table 7.

<figure><img src="image_4.png" /><figcaption>Fig. 2. Operating diagram of (1.1). (a) The case considered in Section 3.3. (b-c) Magnifications of (a) showing the curve $\Gamma_H$.</figcaption></figure>

<table><caption>Table 7. Existence and local stability of steady states according to the regions in the operating diagram of Fig. 2.</caption><thead><tr><th>Condition 1</th><th>Condition 2</th><th>Region</th><th>Color</th><th>$E_0$</th><th>$E_1^1$</th><th>$E_1^2$</th></tr></thead><tbody><tr><td>$S_{in} < \min(\lambda_{LP}(D), \lambda_{BP}(D))$</td><td></td><td>$I_0$</td><td>Cyan</td><td>S</td><td></td><td></td></tr><tr><td>$\lambda_{BP}(D) < S_{in}$</td><td>$c_4(S_{in}, D) > 0$</td><td>$I_1$</td><td>Red</td><td>U</td><td>S</td><td></td></tr><tr><td>$\lambda_{LP}(D) < S_{in} < \lambda_{BP}(D)$</td><td>$c_4(S_{in}, D) > 0$</td><td>$I_2$</td><td>Green</td><td>S</td><td>S</td><td>U</td></tr><tr><td>$\lambda_{BP}(D) < S_{in}$</td><td>$c_4(S_{in}, D) < 0$</td><td>$I_3$</td><td>Blue</td><td>U</td><td>U</td><td></td></tr><tr><td>$\lambda_{LP}(D) < S_{in} < \lambda_{BP}(D)$</td><td>$c_4(S_{in}, D) < 0$</td><td>$I_4$</td><td>Yellow</td><td>S</td><td>U</td><td>U</td></tr></tbody></table>

**Proof.** With this set of parameters, we have $\lambda_v(D) < \lambda_{BP}(D)$ so that the function $S \mapsto H(S)$ is defined and decreasing on $(\lambda_v(D), \lambda_{BP}(D)]$ as shown in Fig. 15. Hence, there can be up to two positive steady states denoted by $E_1^1$ and $E_1^2$ so that $H'(S^*) < -D$ and $H'(S^*) > -D$, respectively. They emerge in the positive quadrant via a Limit Point (LP) bifurcation. Thus, the corresponding curve $\Gamma_{LP}$ exists and is unbounded (see Fig. 2). In addition, the two curves $\Gamma_u$ and $\Gamma_b$ intersect when $(S_{in}, D) = (14.588, 1.147)$ so that $\Gamma_{BP} = \Gamma_u$ for all $D \in [0, 1.147]$ and $\Gamma_{BP} = \Gamma_b$ for all $D \in [1.147, \bar{D}_b[$. In this case
$\lambda_v(D) < \lambda_{BP}(D) = \min(\lambda_u(D), \lambda_b(D))$, $\varphi(S^*) < 0$ and $0 < \psi(S^*) < b$. From Table 2, it follows that $c_3 > 0$ [resp. $c_3 < 0$] for $E_1^1$, [resp. $E_1^2$]. Thus, $E_1^2$ is unstable whenever it exists. In Appendix D and Fig. 15, we show that the stability condition $c_4 > 0$ of $E_1^1$ does not always hold for all $S_{in}$ and $D$ in the domain of existence of $E_1^1$ so that the curve $\Gamma_H$ exists and is unbounded. Similarly to the previous cases, we show that the $\Gamma_H$ curve of equation $c_4(S_{in}, D) = 0$ is given by the union of the two curves $\Gamma_H^1$ (on the left of the maximum) of equation $S_{in} = S_{in}^{H1}(D)$ and $\Gamma_H^2$ (on the right of the maximum) of equation $S_{in} = S_{in}^{H2}(D)$, see Fig. 2. To construct the operating diagram in Fig. 2, we plot the various curves $\Gamma_i$, $i = \{BP, LP, H\}$ defined in Tables 4 and corresponding to the existence and stability conditions of all steady states provided in Table 2. These curves divide the $(S_{in}, D)$-plane of operating parameters in five regions $I_0$ to $I_4$ as shown in Fig. 2. From Table 2, the existence and local stability conditions of the steady states are shown in Table 7.

Fig. 2(a) illustrates the operating diagram of model (1.1) while Figs. 2(b-c) illustrate magnifications of regions $I_1, I_3$ and $I_4$. The blue region $I_3$ corresponds to the instability of the positive steady state. The yellow region $I_4$ corresponds to the instability of the two positive steady states $E_1^1$ and $E_1^2$ where the system can exhibit bistability with either coexistence around a stable limit cycle or the washout of the isolated and attached bacteria.

## 4. One- and two-parameter bifurcation diagrams with MATCONT

In this section, we use MATCONT 56 to numerically analyze the one- and two-parameter diagrams of model (1.1) and to detect two-parameter bifurcations that cannot be established theoretically. It also allows us to validate our theoretical results. In fact, MATCONT is a MATLAB numerical continuation package used to analyze the different types of bifurcations of the continuous and discrete parameterized systems of ODEs. It allows to trace the trajectories over time according to the initial condition and the bifurcation diagrams with a single parameter or two parameters. More precisely, it allows one to visualize the curves of steady states according to a parameter by determining their local asymptotic behavior thanks to the calculation of the eigenvalues of the Jacobian matrix evaluated at the steady state. Moreover, it also allows one to determine the stable or unstable limit cycles by calculation of the sign of the First Lyapunov coefficient. Thanks to test functions, MATCONT detects all types of bifurcations such as the transcritical bifurcation or Branch Points (BP), saddle-node or Limit Point (LP) bifurcation, Cusp (CP) bifurcation, Hopf (H) bifurcation, Limit Point of Cycles (LPC) or fold bifurcation point of limit cycles, period-doubling bifurcation point of limit cycles.

From these critical bifurcation points, MATCONT can determine the various curves in the operating diagram according to two parameters by numerical continuation. These curves of objects of a given type (e.g. steady states, limit cycles, Hopf bifurcation points, homoclinic orbits, etc.) are calculated under variation of one or more system parameters. The reader is addressed to the relevant paper of Dhooge et al. 44 for more on this interesting subject.

### 4.1. One- and two-parameter bifurcation diagrams of Section 3.2

In this section, we determine the one and two-parameter bifurcation diagrams in Fig. 3 using MATCONT for the set of biological parameter values in [5, Fig. 12] where the one-parameter bifurcation diagram is obtained using SCILAB.⁵⁸ The corresponding set of parameters is provided in Table 20 (line 2). The intersection point between $\Gamma_b$ and $\Gamma_{LP}$ is a two-parameter bifurcation of type Cusp (CP) while the intersection points between $\Gamma_H$ and $\Gamma_b$ with the $D = 0$ axis are of type Bogdanov-Takens (BT). These types of bifurcation are not detected in the theoretical study of the operating diagram obtained in Section 3.2. Table 8 summarizes the critical operating parameters, the state, and the
12 Fekih-Salem and Sari

normal form coefficient for BT and CP bifurcations (for more details and additional information on the normal form coefficients, see [57, Section 3]).

<table><caption>Table 8. Operating parameters, state, normal form coefficient values ($\theta_1, \theta_2$) [resp. $\theta_3$] for BT [resp. CP], at the bifurcation points in Fig. 3. The abbreviation BT [resp. CP] means a Bogdanov-Takens point [resp. Cusp] bifurcation.</caption><thead><tr><th>Bifurcation</th><th>Parameter (S<sub>in</sub>, D)</th><th>State (S, u, v)</th><th>Normal form coefficient</th></tr></thead><tbody><tr><td>BT</td><td>(2,0)</td><td>(2,0,0)</td><td>(&theta;<sub>1</sub>, &theta;<sub>2</sub>) = (1.03 &times; 10<sup>-6</sup>, -0.6)</td></tr><tr><td>CP</td><td>(2.204,0.130)</td><td>(2.204,0,0)</td><td>&theta;<sub>3</sub> = -0.523</td></tr><tr><td>BT</td><td>(2.236,0)</td><td>(2,0,0)</td><td>(&theta;<sub>1</sub>, &theta;<sub>2</sub>) impossible</td></tr></tbody></table>

Fig. 3(c) illustrates the one-parameter bifurcation diagram in variable S when D is fixed at D = 0.1 in the case considered in Section 3.2. It reveals the appearance and the disappearance of stable limit cycles via two Hopf bifurcations.

<figure><img src="image_5.png" /><figcaption>Fig. 3. MATCONT: (a) operating diagram of (1.1) in the case considered in Section 3.2. (b) Magnification of (a) showing the curve $\Gamma_H$. (c) The corresponding one-parameter bifurcation diagram in variable S when D = 0.1.</figcaption></figure>

**Remark 4.1** In all the figures showing bifurcation diagrams, stable [resp. unstable] equilibria are plotted in Red [resp. Blue] and stable [resp. unstable] limit cycles are plotted in Magenta [resp. Cyan]. In addition, H and LPC denote Hopf bifurcations and Limit Point of Cycles, respectively.

## 4.2. Operating diagram in the (S<sub>in</sub>, D)-plane of Section 3.3

Fig. 4 illustrates the operating diagram obtained numerically using MATCONT. It is identical to the operating diagram obtained theoretically in Fig. 2. However, MATCONT detects the nature of bifurcations at the intersection points between the curves $\Gamma_u$ and $\Gamma_{LP}$, which is of Cusp (CP) type and between the curves $\Gamma_u$ and $\Gamma_H$ with the S<sub>in</sub>-axis which are Bogdanov-Takens points (BT). The critical operating parameters, the state, and the normal form coefficient for BT and CP bifurcations are summarized in Table 9.

## 4.3. One-parameter bifurcation diagram in S<sub>in</sub> when D = 0.1 fixed in the operating diagram of Section 3.3

In what follows, we will analyze the various types of bifurcation by crossing one region to another in the operating diagram of Fig. 2 (or Fig. 4). The nature of all the bifurcations by passing through the various curves $\Gamma_i$ defined in Table 4 is described in the following result.
<figure><img src="image_6.png" /><figcaption>Fig. 4. MATCONT: (a) operating diagram of (1.1) in the case considered in Section 3.3. (b-c) Magnifications of (a) showing the curve $\Gamma_H$.</figcaption></figure>

<table><caption>Table 9. Operating parameters, state, normal form coefficient values ($\theta_1, \theta_2$) [resp. $\theta_3$] for BT [resp. CP], at the bifurcation points in Fig. 4. The abbreviation BT [resp. CP] means a Bogdanov-Takens point [resp. Cusp] bifurcation.</caption><thead><tr><th>Bifurcation</th><th>Parameter ($S_{in}, D$)</th><th>State ($S, u, v$)</th><th>Normal form coefficient</th></tr></thead><tbody><tr><td>BT</td><td>(3.714,0)</td><td>(3.714,0,0)</td><td>($\theta_1, \theta_2$) = (3.23×10<sup>-8</sup>, -0.306)</td></tr><tr><td>CP</td><td>(3.819,0.032)</td><td>(3.819,0,0)</td><td>$\theta_3$ = -0.483</td></tr><tr><td>BT</td><td>(3.823,0)</td><td>(3.714,0,0)</td><td>($\theta_1, \theta_2$) = (2.38×10<sup>-6</sup>, -0.306)</td></tr><tr><td>BT</td><td>(14.588,1.147)</td><td>(14.588,0,0)</td><td>($\theta_1, \theta_2$) = (2.019, -1.127)</td></tr></tbody></table>

**Proposition 4.1** Let $f$ and $g$ be the specific growth rates defined in (3.1). Let the set of the biological parameter values be in Table 20 (line 3). The nature of all the bifurcations of model (1.1) by crossing the different regions of the operating diagram in Fig. 2 is provided in Table 10.

<table><caption>Table 10. Nature of all the bifurcations of system (1.1) by passing the different curves $\Gamma_i$, $i = \{\text{BP, LP, H}\}$ defined in Table 4.</caption><thead><tr><th>Transition</th><th>Curve</th><th>Bifurcation</th><th>Steady states</th></tr></thead><tbody><tr><td>$\mathcal{I}_0$ to $\mathcal{I}_2$</td><td>$\Gamma_{LP}$</td><td>LP</td><td>$E_1^1 = E_1^2$</td></tr><tr><td>$\mathcal{I}_0$ to $\mathcal{I}_1$</td><td>$\Gamma_u$</td><td>BP</td><td>$E_0 = E_1^1$</td></tr><tr><td>$\mathcal{I}_2$ to $\mathcal{I}_1$</td><td>$\Gamma_{BP}$</td><td>BP</td><td>$E_1^2 = E_0$</td></tr><tr><td>$\mathcal{I}_2$ to $\mathcal{I}_4$</td><td>$\Gamma_H$</td><td>H</td><td>$E_1^1$</td></tr><tr><td>$\mathcal{I}_1$ to $\mathcal{I}_3$</td><td>$\Gamma_H$</td><td>H</td><td>$E_1^1$</td></tr><tr><td>$\mathcal{I}_4$ to $\mathcal{I}_3$</td><td>$\Gamma_u$</td><td>BP</td><td>$E_1^2 = E_0$</td></tr></tbody></table>

**Proof.** The transition from region $\mathcal{I}_0$ to $\mathcal{I}_2$ via the $\Gamma_{LP}$ curve corresponds to an LP bifurcation between $E_1^1$ and $E_1^2$ since we have $S_{in} = \lambda_{LP}(D) := H(S_{LP}(D))/D + S_{LP}(D)$, that is,

$$D(S_{in} - S_{LP}(D)) = H(S_{LP}(D)), \text{ with } S_{LP}(D) \text{ is solution of equation } H'(S_{LP}(D)) = -D.$$

The transition from region $\mathcal{I}_0$ to $\mathcal{I}_1$ via the $\Gamma_{BP} = \Gamma_u$ curve corresponds to a BP bifurcation between $E_0$ and $E_1^1$ since $S_{in} = \lambda_{BP}(D) = \min(\lambda_u(D), \lambda_b(D)) = \lambda_u(D)$. If $S_{in} = \lambda_u(D)$, then $E_0 = E_1^1 = (S_{in}, 0, 0)$. If $S_{in} < \lambda_{BP}(D)$, then $E_0$ is stable and $E_1^1$ does not exist. If $S_{in} > \lambda_{BP}(D)$, then $E_0$ is unstable and $E_1^1$ emerges LES. The transition from region $\mathcal{I}_2$ to $\mathcal{I}_4$ (or $\mathcal{I}_1$ to $\mathcal{I}_3$) via the $\Gamma_H$ curve corresponds to a Hopf bifurcation at $E_1^1$ since $c_4(S_{in})$ changes sign with the emergence of a periodic solution.
14 Fekih-Salem and Sari

Similarly to the previous BP bifurcation, the transition from region $\mathcal{I}_2$ to $\mathcal{I}_1$ (or $\mathcal{I}_4$ to $\mathcal{I}_3$) via the $\Gamma_{BP} = \Gamma_u \cup \Gamma_b$ curve corresponds to a BP bifurcation between $E_0$ and $E_1^2$. $\square$

Let $D$ be fixed at $D = D^* = 0.1$. Next, we analyze the one-parameter bifurcation diagram with respect to $S_{in}$ as the bifurcating parameter to show the nature of bifurcations by crossing various boundaries between the different regions in the operating diagram. Note that the one-parameter bifurcation diagram in $D$ can be obtained in the same way. Using MATCONT, we illustrate in Fig. 5 the one-parameter bifurcation diagram in $S_{in}$, with $S$ on the $y$-axis. Similarly, we can obtain the one-parameter bifurcation diagram for the concentrations of isolated and attached bacteria, $u$ and $v$, respectively. Note that the two-parameter bifurcation diagram does not show the disappearance of the limit cycle like the one-parameter bifurcation diagram.

<figure><img src="image_7.png" /><figcaption>Fig. 5. MATCONT: one-parameter bifurcation diagram of (1.1) in variable $S$ in the case considered in Section 3.3. (b) Magnification of the two homoclinic bifurcations; (c) magnification of the first supercritical Hopf bifurcation.</figcaption></figure>

The following result establishes the one-parameter bifurcation diagram in $S_{in}$ from the operating diagram in Fig. 2 when $D = 0.1$ is fixed.

**Proposition 4.2** For the specific growth rates $f$ and $g$ defined in (3.1) and the set of biological parameter values in Table 20 (see line 3), the existence and local stability of all steady states of (1.1) according to $S_{in}$ are described in Table 11 when $D = 0.1$ is fixed. The critical values $\sigma_i$, $i = 1, \dots, 6$ of different bifurcations according to the parameter $S_{in}$ and the corresponding nature are defined in Table 12.

<table><caption>Table 11. Existence and stability of all steady states of (1.1) according to $S_{in}$ for the set of parameters in Fig. 2 when $D = 0.1$. The critical values $\sigma_i$, $i = 1, \dots, 6$ are defined in Table 12.</caption><thead><tr><th>Interval of S<sub>in</sub></th><th>E<sub>0</sub></th><th>E<sub>1</sub><sup>1</sup></th><th>E<sub>1</sub><sup>2</sup></th></tr></thead><tbody><tr><td>(0, σ<sub>1</sub>)</td><td>S</td><td></td><td></td></tr><tr><td>(σ<sub>1</sub>, σ<sub>2</sub>)</td><td>S</td><td>S</td><td>U</td></tr><tr><td>(σ<sub>2</sub>, σ<sub>3</sub>)</td><td>S</td><td>U</td><td>U</td></tr><tr><td>(σ<sub>3</sub>, σ<sub>4</sub>)</td><td>S</td><td>U</td><td>U</td></tr><tr><td>(σ<sub>4</sub>, σ<sub>5</sub>)</td><td>S</td><td>U</td><td>U</td></tr><tr><td>(σ<sub>5</sub>, σ<sub>6</sub>)</td><td>U</td><td>U</td><td></td></tr><tr><td>(σ<sub>6</sub>, +∞)</td><td>U</td><td>S</td><td></td></tr></tbody></table>
<table><caption>Table 12. Definitions of the critical values $\sigma_i$, $i = 1, \dots, 6$ of $D$ and their corresponding nature of bifurcation when $D = 0.1$ is fixed. The abbreviations Hom and T mean homoclinic and T period of solutions, respectively.</caption><thead><tr><th>Definition</th><th>Value</th><th>Bifurcation</th></tr></thead><tbody><tr><td>$\sigma_1 = \lambda_{LP}(D)$</td><td>3.837</td><td>LP</td></tr><tr><td>$\sigma_2$ is the first solution of equation $c_4(S_{in}) = 0$</td><td>3.842</td><td>H</td></tr><tr><td>$\sigma_3$: $\lim_{S_{in} \to \sigma_3} T(S_{in}) = +\infty$</td><td>3.84770</td><td>Hom</td></tr><tr><td>$\sigma_4$: $\lim_{S_{in} \to \sigma_4} T(S_{in}) = +\infty$</td><td>4.03468</td><td>Hom</td></tr><tr><td>$\sigma_5 = \lambda_{BP}(D)$</td><td>4.061</td><td>BP</td></tr><tr><td>$\sigma_6$ is the second solution of equation $c_4(S_{in}) = 0$</td><td>8.179</td><td>H</td></tr></tbody></table>

**Proof.** Increasing $S_{in}$ from zero, the bifurcation diagram in Fig. 5(a) illustrates the BP bifurcation occurring at $S_{in} = \sigma_5 \approx 4.061$ between $E_0$ and $E_1^2$. Increasing $S_{in}$ further, the washout steady state $E_0$ changes stability and becomes unstable (see Fig. 5(a-b)). Continuing this BP bifurcation backward with decreasing $S_{in}$, $E_1^2$ emerges at $S_{in} = \sigma_5$ unstable. A first close-up is illustrated in Fig. 5(b) and a second close-up is illustrated in Fig. 5(c). Next, there is an LP bifurcation between $E_1^1$ and $E_1^2$ at $S_{in} = \sigma_1 \approx 3.837$ when $S \approx 3.492$, $u \approx 996 \times 10^{-5}$ and $v \approx 107 \times 10^{-5}$ so that these two interior steady states disappear by decreasing $S_{in}$ further. Inversely, increasing $S_{in}$ from LP, $E_1^1$ and $E_1^2$ appear stable and unstable, respectively. After that, increasing $S_{in}$ further, a Hopf bifurcation occurs at $E_1^1$ when $S_{in} = \sigma_2 \approx 3.842$, $S \approx 3.422$, $u \approx 0.012$ and $v \approx 1.5 \times 10^{-3}$. A stable limit cycle emerges through a supercritical Hopf bifurcation where the first Lyapunov coefficient is given by $-0.430$. Moreover, $E_1^1$ changes stability and becomes unstable. Increasing $S_{in}$ further, once again a Hopf bifurcation occurs at $E_1^1$ when $S_{in} = \sigma_6 \approx 8.179$, $S \approx 1.963$, $u \approx 0.140$ and $v \approx 0.139$. A stable limit cycle disappears through a supercritical Hopf bifurcation where the first Lyapunov coefficient is given by $-34 \times 10^{-3}$. Moreover, $E_1^1$ changes stability and becomes LES. Fig. 16 shows the stable limit cycles in the three-dimensional space $(S, u, v)$ for different values of $S_{in}$ between $\sigma_4$ and $\sigma_6$. Starting from the first Hopf bifurcation at $\sigma_2$ and increasing $S_{in}$, the radius of the stable limit cycle increases until its disappearance through a homoclinic bifurcation when $S_{in} = \sigma_3 \approx 3.8477$. Fig. 17(a) shows that the period of the cycle tends to infinity when $S_{in}$ tends to $\sigma_3$. Starting from the second Hopf bifurcation at $\sigma_6$ and decreasing $S_{in}$, the radius of the stable limit cycle first increases and then decreases until its disappearance through a homoclinic bifurcation when $S_{in} = \sigma_4 \approx 4.03468$. Fig. 17(b) shows that the period of the cycle tends to infinity when $S_{in}$ tends to $\sigma_4$, confirming the homoclinic bifurcation. □

### 4.4. One-parameter bifurcation diagram in $D$ when $S_{in} = 3.845$ fixed in the operating diagram of Section 3.3

Next, we study the one-parameter bifurcation diagram with respect to the dilution rate $D$ when $S_{in} = 3.845$ and all other parameter values are those used in the operating diagram of Fig. 2 (or Fig. 4) considered in Section 3.3.

**Proposition 4.3** For the specific growth rates $f$ and $g$ defined in (3.1) and the set of biological parameter values in Table 20 (see line 3), the existence and local stability of the limit cycle and equilibria of (1.1) according to $D$ are shown in Fig. 6 and described in Table 13 when $S_{in} = 3.845$ is fixed. The critical values $\sigma_i$, $i = 1, \dots, 7$ of different bifurcations according to the parameter $D$ and their corresponding natures are defined in Table 14. The corresponding transitions in the operating diagram regions of Fig. 2 (or Fig. 4) for fixed $S_{in}$ are provided in Table 15.
<figure><img src="image_8.png" /><figcaption>Fig. 6. MATCONT: (a) one-parameter bifurcation diagram in the parameter *D* when *S*<sub>in</sub> = 3.845 fixed in the operating diagram of Fig. 2 (or Fig. 4). (b) Enlargement showing the leftmost limit cycle.</figcaption></figure>

<table><caption>Table 13. Existence and stability of the limit cycles and equilibria shown in Fig. 6 according to *D* for the set of parameters in 2 (or Fig. 4) when *S*<sub>in</sub> = 3.845. The critical values *σ*<sub>*i*</sub>, *i* = 1, . . . , 7 are defined in Table 14.</caption><thead><tr><th>Interval of D</th><th>E<sub>0</sub></th><th>E<sub>1</sub><sup>1</sup></th><th>E<sub>1</sub><sup>2</sup></th><th>C<sub>1</sub></th></tr></thead><tbody><tr><td>(0, σ<sub>1</sub>)</td><td>U</td><td>U</td><td></td><td>S</td></tr><tr><td>(σ<sub>1</sub>, σ<sub>2</sub>)</td><td>U</td><td>S</td><td></td><td></td></tr><tr><td>(σ<sub>2</sub>, σ<sub>3</sub>)</td><td>S</td><td>S</td><td>U</td><td></td></tr><tr><td>(σ<sub>3</sub>, σ<sub>4</sub>)</td><td>S</td><td></td><td></td><td></td></tr><tr><td>(σ<sub>4</sub>, σ<sub>5</sub>)</td><td>S</td><td>S</td><td>U</td><td></td></tr><tr><td>(σ<sub>5</sub>, σ<sub>6</sub>)</td><td>S</td><td>U</td><td>U</td><td>S</td></tr><tr><td>(σ<sub>6</sub>, σ<sub>7</sub>)</td><td>U</td><td>S</td><td>U</td><td></td></tr><tr><td>(σ<sub>7</sub>, +∞)</td><td>S</td><td></td><td></td><td></td></tr></tbody></table>

**Proof.** The one-parameter bifurcation diagram in the variables *S*, *u*, and *v* corresponds to the vertical line *S*<sub>in</sub> = 3.845 in the operating diagram shown in Fig. 2 (or Fig. 4). Hence, there exist three critical values *σ*<sub>1</sub>, *σ*<sub>5</sub>, and *σ*<sub>6</sub> [resp. *σ*<sub>3</sub>, *σ*<sub>4</sub> and *σ*<sub>7</sub>] of the parameter *D*, corresponding to the intersection of the vertical line *S*<sub>in</sub> = 3.845 with the Γ<sub>H</sub> curve [resp. Γ<sub>LP</sub> curve]. Another critical value *σ*<sub>2</sub> corresponds to the intersection of the vertical line *S*<sub>in</sub> = 3.845 with the Γ<sub>u</sub> curve. Using MATCONT, we obtain the critical values *σ*<sub>*i*</sub>, *i* = 1, . . . , 7 of *D* in Table 14 and their corresponding bifurcation types.

Increasing *D* from zero, the one-parameter bifurcation diagram in Fig. 6 shows the BP bifurcation occurring at *D* = *σ*<sub>2</sub> ≈ 0.039 between *E*<sub>0</sub> and *E*<sub>1</sub><sup>2</sup>, which appears unstable while *E*<sub>0</sub> becomes stable. This bifurcation corresponds to the passage from the region *I*<sub>1</sub> to *I*<sub>2</sub> through the Γ<sub>BP</sub> curve. Next, there is an LP bifurcation between *E*<sub>1</sub><sup>1</sup> and *E*<sub>1</sub><sup>2</sup> at *D* = *σ*<sub>3</sub> ≈ 0.045 so that these two interior equilibria disappear by increasing *D* further. This bifurcation corresponds to the passage from the region *I*<sub>2</sub> to *I*<sub>0</sub> through the Γ<sub>LP</sub> curve. Inversely, decreasing *D* from LP, *E*<sub>1</sub><sup>1</sup> appears stable and then it destabilizes at *D* = *σ*<sub>1</sub> ≈ 0.018 with the emergence of a stable limit cycle through a supercritical Hopf
<table><caption>Table 14. Definitions of the critical values $\sigma_i$, $i = 1, \dots, 7$ of $D$ and their corresponding bifurcation types when $S_{in} = 3.845$ is fixed, as well as the first Lyapunov coefficient $l_1$ at Hopf points.</caption><thead><tr><th>Definition</th><th>Value</th><th>Bifurcation</th></tr></thead><tbody><tr><td>$\sigma_1$ is the first solution of equation $c_4(D) = 0$</td><td>0.018</td><td>H ($l_1 \approx -0.216$)</td></tr><tr><td>$\sigma_2 = \lambda_u^{-1}(D) = f(\alpha D + m_u)$</td><td>0.039</td><td>BP</td></tr><tr><td>$\sigma_3$ is the first solution of equation $\lambda_{LP}(D) = S_{in}$</td><td>0.045</td><td>LP</td></tr><tr><td>$\sigma_4$ is the second solution of equation $\lambda_{LP}(D) = S_{in}$</td><td>0.081</td><td>LP</td></tr><tr><td>$\sigma_5$ is the second solution of equation $c_4(D) = 0$</td><td>0.088</td><td>H ($l_1 \approx -0.598$)</td></tr><tr><td>$\sigma_6$ is the third solution of equation $c_4(D) = 0$</td><td>0.127</td><td>H ($l_1 \approx -0.180$)</td></tr><tr><td>$\sigma_7$ is the third solution of equation $\lambda_{LP}(D) = S_{in}$</td><td>0.288</td><td>LP</td></tr></tbody></table>

<table><caption>Table 15. Type of bifurcations according to the various transitions in the operating diagram regions of Fig. 2 (or Fig. 4) through the different curves $\Gamma_i$, $i = \{BP, LP, H\}$ defined in Table 4.</caption><thead><tr><th>Transition</th><th>Curve</th><th>Bifurcation</th><th>Equilibria & limit cycle</th></tr></thead><tbody><tr><td>$\mathcal{I}_3$ to $\mathcal{I}_1$</td><td>$\Gamma_H$</td><td>H</td><td>$E_1^1 = C_1$</td></tr><tr><td>$\mathcal{I}_1$ to $\mathcal{I}_2$</td><td>$\Gamma_{BP}$</td><td>BP</td><td>$E_0 = E_1^2$</td></tr><tr><td>$\mathcal{I}_0$ to $\mathcal{I}_2$</td><td>$\Gamma_{LP}$</td><td>LP</td><td>$E_1^1 = E_1^2$</td></tr><tr><td>$\mathcal{I}_2$ to $\mathcal{I}_4$</td><td>$\Gamma_H$</td><td>H</td><td>$E_1^1 = C_1$</td></tr></tbody></table>

bifurcation since the first Lyapunov coefficient is negative (see Table 14). This bifurcation corresponds to the passage from the region $\mathcal{I}_1$ to $\mathcal{I}_3$ through the $\Gamma_H$ curve.

Now, by increasing $D$ further, there is a second LP bifurcation between $E_1^1$ and $E_1^2$ at $D = \sigma_4 \approx 0.081$, which appear stable and unstable, respectively. This bifurcation corresponds to the passage from the region $\mathcal{I}_0$ to $\mathcal{I}_2$ through the $\Gamma_{LP}$ curve. Then, $E_1^1$ destabilizes at $D = \sigma_5 \approx 0.088$ with the emergence of a stable limit cycle through a supercritical Hopf bifurcation since the first Lyapunov coefficient is negative (see Table 14). This bifurcation corresponds to the passage from the region $\mathcal{I}_2$ to $\mathcal{I}_4$ through the $\Gamma_H$ curve. As $D$ increases, the radius of the limit cycle $C_1$ increases and then decreases until it disappears at $D = \sigma_6 \approx 0.127$ through a supercritical Hopf bifurcation since the first Lyapunov coefficient is negative (see Table 14). This bifurcation corresponds to the passage from the region $\mathcal{I}_4$ to $\mathcal{I}_2$ through the $\Gamma_H$ curve. Finally, there is a third LP bifurcation at $D = \sigma_7 \approx 0.288$ with the coalescence and disappearance of a node and a saddle point ($E_1^1$ and $E_1^2$, respectively). This bifurcation corresponds to the passage from the region $\mathcal{I}_2$ to $\mathcal{I}_0$ through the $\Gamma_{LP}$ curve. $\square$

**Remark 4.2** Note that MATCONT could not continue tracing limit cycles in the one-parameter bifurcation diagram of Figure 6 when $D$ is quite close to zero because of numerical problems related to precision.

## 5. Effect of flocculation on the operating diagram

### 5.1. Operating diagram in the $(S_{in}, D)$-plan for several values of a and b

To see the effect of the attachment and detachment rates on the asymptotic behavior of the process, we study the operating diagram in the $(S_{in}, D)$-plan by varying the parameters a and b to zero. We consider the same parameter values as in Section 3.3 (or in [5, Fig. 6]). Fig. 7(a-c) illustrates the reduction in the size of the coexistence region $\mathcal{I}_3$ by decreasing the rates of attachment a and/or detachment b.

In Fig. 7(d-f), region $\mathcal{I}_3$ has disappeared and region $\mathcal{I}_2$ is reduced to disappearance in the limiting case $a = b = 0$, where we obtain the operating diagram of the classic
<figure><img src="image_9.png" /><figcaption>Fig. 7. Operating diagram in the case considered in Section 3.3 when (a) $a = 0.5$ and $b = 2$; (b) $a = 0.1$ and $b = 2$; (c) $a = 0.5$ and $b = 1.05$; (d) $a = 0.01$ and $b = 0.2$; (e) $a = 0.01$ and $b = 0.01$; (f) $a = 0$ and $b = 0$.</figcaption></figure>

chemostat model. Fig. 8 illustrates the operating diagrams with the various colors of regions by decreasing the rates of attachment and detachment, where there is a reduction in the size of regions $I_2$ and $I_3$ until their disappearance.

<figure><img src="image_10.png" /><figcaption>Fig. 8. Operating diagram in the case considered in Section 3.3 when (a) $a = 0.5$ and $b = 2$; (b) $a = 0.5$ and $b = 1.05$; (c) $a = 0.01$ and $b = 0.01$; (d) $a = 0$ and $b = 0$.</figcaption></figure>

## 5.2. Operating diagram in the (a, b)-plan

In the following, we study the operating diagram with respect to the parameters $a$ and $b$, by fixing $(S_{in}, D)$. To this end, we use MATCONT to construct the operating diagram in the (a, b)-plane when $(S_{in}, D) = (5, 0.08)$. It corresponds to a point in region $I_3$ of the operating diagram of Fig. 8(a) when $(a, b) = (0.5, 2)$. All other parameter values are identical to those in Fig. 8 and are provided in Table 20 (line 4).

In Fig. 9, we obtain the $\Gamma_H$ curve (in red) denoted $\Gamma_H$ which corresponds to the Hopf bifurcation and the LPC curve (in green) denoted $\Gamma_{LPC}$, which corresponds to
<figure><img src="image_11.png" /><figcaption>Fig. 9. MATCONT: (a) operating diagram of (1.1) in the (a, b)-plane in case of Table 20 (line 4) when ($S_{in}$, $D$) = (5, 0.08). (b) Magnification near the GH bifurcation. The asymptotic behavior of (1.1) in the regions $\mathcal{I}_1$, $\mathcal{I}_1^b$, and $\mathcal{I}_3$ is given in Table 16. $\Gamma_{LPC}$ is the LPC curve.</figcaption></figure>

<table><caption>Table 16. Existence and local stability of steady states and limit cycles in the regions of the operating diagram in Fig. 9. The $\mathcal{I}_1^b$ region exhibits bistability between one cycle and $E_1$.</caption><thead><tr><th>Region</th><th>$E_0$</th><th>$E_1$</th><th>Cycle</th></tr></thead><tbody><tr><td>$\mathcal{I}_1$</td><td>U</td><td>S</td><td></td></tr><tr><td>$\mathcal{I}_1^b$</td><td>U</td><td>S</td><td>Two cycles (one stable and one unstable)</td></tr><tr><td>$\mathcal{I}_3$</td><td>U</td><td>U</td><td>One stable cycle</td></tr></tbody></table>

Limit Point of Cycles (LPC). The LPC curve defines a sub-region, denoted $\mathcal{I}_1^b$ of region $\mathcal{I}_1$ in which the system has a multiplicity of limit cycles and hence, presents a bistability behavior. In the region $\mathcal{I}_1^b$, a limit cycle and $E_1$ are both stable. The asymptotic behavior of system (1.1) according to the various regions of the operating diagram is summarized in Table 16.

In addition to the LPC curve, MATCONT detects a new two-parameter bifurcation in the operating diagram of Fig. 9, corresponding to the intersection of H and LPC curves, which is a Generalized Hopf (GH) bifurcation point occurring at coordinates $(a, b) \approx (24.685, 1.735)$. Crossing the $\Gamma_H$ curve on the right side of the GH point produces a subcritical Hopf bifurcation with the appearance of an unstable limit cycle while crossing it on the left side of the GH point or the branch at the top of $\Gamma_H$ curve gives a supercritical Hopf bifurcation with the appearance of a stable limit cycle.

These bifurcations and phenomena are new and have not been detected in the previous operating diagram according to the parameters ($S_{in}$, $D$). Note that with this parameter set, when ($S_{in}$, $D$) = (5, 0.08), the LP curve of the emergence of two positive equilibria and the BP curve of the stabilization of the washout equilibrium do not exist. Thus, $E_0$ always remains unstable and the coexistence equilibrium $E_1$ always exists and is unique for any values of attachment and detachment rates.

Next, we study the asymptotic behavior in the various regions of the operating diagram depicted in Fig. 9 with the one-parameter bifurcation diagrams by fixing either a or b.

### 5.3. One-parameter bifurcation diagram in the attachment rate a when b = 2

In this section, we study the one-parameter bifurcation diagram in the attachment rate a when b = 2 and all other parameter values are those used in the operating diagram of Fig. 9 such that ($S_{in}$, $D$) = (5, 0.08).

Therefore, the asymptotic behavior of system (1.1) corresponds to the horizontal line b = 2 of the operating diagram shown in Fig. 9. Hence, there exists a critical value $\sigma_1$ of
<figure><img src="image_12.png" /><figcaption>Fig. 10. MATCONT: (a) one-parameter bifurcation diagram in the parameter *a* when *b* = 2 fixed in the operating diagram of Fig. 9 such that (*S*<sub>in</sub>, *D*) = (5, 0.08). (b) Enlargement showing the second Hopf (subcritical) bifurcation.</figcaption></figure>

<table><caption>Table 17. Existence and stability of the limit cycles and equilibria shown in Fig. 10. The last column gives the corresponding region of Fig. 9 and the type of the bifurcations occurring at the critical values σᵢ, as well as the first Lyapunov coefficient l₁ at Hopf points.</caption><thead><tr><th>Interval of a</th><th>E<sub>0</sub></th><th>E<sub>1</sub></th><th>C<sub>1</sub></th><th>C<sub>2</sub></th><th>Region &amp; Bifurcation</th></tr></thead><tbody><tr><td>(0, σ<sub>1</sub>)</td><td>U</td><td>S</td><td></td><td></td><td>I<sub>1</sub></td></tr><tr><td>a = σ<sub>1</sub> ≈ 1.173</td><td></td><td></td><td></td><td></td><td>Hopf (l<sub>1</sub> ≈ -2.305 × 10<sup>-2</sup>)</td></tr><tr><td>(σ<sub>1</sub>, σ<sub>2</sub>)</td><td>U</td><td>U</td><td>S</td><td></td><td>I<sub>3</sub></td></tr><tr><td>a = σ<sub>2</sub> ≈ 30.898</td><td></td><td></td><td></td><td></td><td>Hopf (l<sub>1</sub> ≈ 1.754 × 10<sup>-2</sup>)</td></tr><tr><td>(σ<sub>2</sub>, σ<sub>3</sub>)</td><td>U</td><td>S</td><td>S</td><td>U</td><td>I<sub>1</sub><sup>b</sup></td></tr><tr><td>a = σ<sub>3</sub> ≈ 31.465</td><td></td><td></td><td></td><td></td><td>LPC (C<sub>1</sub> = C<sub>2</sub>)</td></tr><tr><td>(σ<sub>3</sub>, +∞)</td><td>U</td><td>S</td><td></td><td></td><td>I<sub>1</sub></td></tr></tbody></table>

the parameter *a*, corresponding to the intersection of the horizontal line *b* = 2 with the Γ<sub>H</sub> curve, and two other critical values σ₂ and σ₃ of the parameter *a*, corresponding to the intersection of the horizontal line *b* = 2 with the H and LPC curves, respectively.

Fig. 10 shows that a supercritical Hopf bifurcation occurs at *a* = σ₁ by increasing the attachment rate *a* from zero, where *E*₁ destabilizes and the limit cycle *C*₁ appears stable since the first Lyapunov coefficient is negative (see Table 17). This bifurcation corresponds to the passage from the region *I*₁ to *I*₃ through the Γ<sub>H</sub> curve. As *a* increases, the radius of the limit cycle *C*₁ increases and then decreases until it disappears at *a* = σ₃ through an LPC bifurcation. This bifurcation corresponds to the passage from the region *I*₁ᵇ to *I*₁ through the LPC curve. As *a* decreases from this LPC bifurcation, the limit cycle *C*₂ appears unstable (see Fig. 10). Then, the radius of this cycle limit decreases until its collapse to the positive equilibrium *E*₁ at *a* = σ₂ by a subcritical Hopf bifurcation since the first Lyapunov coefficient is positive (see Table 17). This bifurcation corresponds to the passage from the region *I*₁ᵇ to *I*₃ through the Γ<sub>H</sub> curve. When *a* > σ₂, the positive equilibrium *E*₁ recovers its stability. The type of the bifurcations with their
critical bifurcation values, as well as the Lyapunov coefficient and asymptotic behavior of system (1.1) as a function of *a* are provided in Table 17.

5.4. One-parameter bifurcation diagram in the detachment rate *b* when *a* = 20

In this section, we study the one-parameter bifurcation diagram in the detachment rate *b* when *a* = 20 and all other parameter values are those used in the operating diagram of Fig. 9 such that (*S*<sub>in</sub>, *D*) = (5, 0.08). Therefore, the asymptotic behavior of system (1.1) corresponds to the vertical line *a* = 20 of the operating diagram shown in Fig. 9. Hence, there exist only two critical values *σ*<sub>1</sub> and *σ*<sub>2</sub> of the parameter *b*, corresponding to the intersection of the vertical line *a* = 20 with the two branches of the Γ<sub>H</sub> curve.

<figure><img src="image_13.png" /><figcaption>Fig. 11. MATCONT: one-parameter bifurcation diagram in *b* when *a* = 20 fixed in the operating diagram of Fig. 9 such that (*S*<sub>in</sub>, *D*) = (5, 0.08).</figcaption></figure>

<table><caption>Table 18. Existence and stability of the limit cycle and equilibria shown in Fig. 11. The last column gives the corresponding region of Fig. 9 and the type of the bifurcations occurring at the critical values σᵢ, as well as the first Lyapunov coefficient l₁ at Hopf points.</caption><thead><tr><th>Interval of b</th><th>E<sub>0</sub></th><th>E<sub>1</sub></th><th>C<sub>1</sub></th><th>Region &amp; Bifurcation</th></tr></thead><tbody><tr><td>(0, σ<sub>1</sub>)</td><td>U</td><td>S</td><td></td><td>I<sub>1</sub></td></tr><tr><td>b = σ<sub>1</sub> ≈ 1.541</td><td></td><td></td><td></td><td>Hopf (l<sub>1</sub> ≈ -1.571 × 10<sup>-2</sup>)</td></tr><tr><td>(σ<sub>1</sub>, σ<sub>2</sub>)</td><td>U</td><td>U</td><td>S</td><td>I<sub>3</sub></td></tr><tr><td>b = σ<sub>2</sub> ≈ 13.139</td><td></td><td></td><td></td><td>Hopf (l<sub>1</sub> ≈ -2.043 × 10<sup>-2</sup>)</td></tr><tr><td>(σ<sub>2</sub>, +∞)</td><td>U</td><td>S</td><td></td><td>I<sub>1</sub></td></tr></tbody></table>

Fig. 11 shows that a supercritical Hopf bifurcation occurs at *b* = *σ*₁ by increasing *b* from zero, where *E*₁ destabilizes and the limit cycle *C*₁ appears stable since the first Lyapunov coefficient is negative (see Table 18). This bifurcation corresponds to the passage from the region *I*₁ to *I*₃ through the Γ<sub>H</sub> curve. As *b* increases, the radius of the limit cycle *C*₁ increases and then decreases until its collapse to the positive equilibrium *E*₁ at *b* = *σ*₂ by a supercritical Hopf bifurcation since the first Lyapunov coefficient is negative (see Table 18). This bifurcation corresponds to the passage from the region *I*₃ to *I*₁ through the Γ<sub>H</sub> curve. When *b* > *σ*₂, the positive equilibrium *E*₁ recovers its stability. The type of the bifurcations with their critical bifurcation values, as well as the Lyapunov coefficient and asymptotic behavior of system (1.1) as a function of *b* are provided in Table 18.

5.5. One-parameter bifurcation diagram in *b* when *a* = 30

In this section, we study the one-parameter bifurcation diagram in the detachment rate *b* when *a* = 30 and all other parameter values are those used in the operating diagram of
Fig. 9 such that $(S_{in}, D) = (5, 0.08)$. Therefore, the asymptotic behavior of system (1.1) corresponds to the vertical line $a = 30$ to the right of point GH in the operating diagram of Fig. 9. Hence, there exists a critical value $\sigma_1$ of the parameter $b$, corresponding to the intersection of the vertical line $a = 30$ with the LPC curve, and two other critical values $\sigma_2$ and $\sigma_3$ of the parameter $b$, corresponding to the intersection of this vertical line with the two branches of $\Gamma_H$ curve.

<figure><img src="image_14.png" /><figcaption>Fig. 12. MATCONT: (a) one-parameter bifurcation diagram in $b$ when $a = 30$ fixed in the operating diagram of Fig. 9 such that $(S_{in}, D) = (5, 0.08)$. (b) Enlargement showing the second Hopf (subcritical) bifurcation.</figcaption></figure>

<table><caption>Table 19. Existence and stability of the limit cycles and equilibria shown in Fig. 11. The last column gives the corresponding region of Fig. 9 and the type of the bifurcations occurring at the critical values $\sigma_i$, as well as the first Lyapunov coefficient $l_1$ at Hopf points.</caption><thead><tr><th>Interval of b</th><th>E<sub>0</sub></th><th>E<sub>1</sub></th><th>C<sub>1</sub></th><th>C<sub>2</sub></th><th>Region &amp; Bifurcation</th></tr></thead><tbody><tr><td>(0, σ<sub>1</sub>)</td><td>U</td><td>S</td><td></td><td></td><td>I<sub>1</sub></td></tr><tr><td>b = σ<sub>1</sub> ≈ 1.946</td><td></td><td></td><td></td><td></td><td>LPC (C<sub>1</sub> = C<sub>2</sub>)</td></tr><tr><td>(σ<sub>1</sub>, σ<sub>2</sub>)</td><td>U</td><td>S</td><td>S</td><td>U</td><td>I<sub>1</sub><sup>b</sup></td></tr><tr><td>b = σ<sub>2</sub> ≈ 1.961</td><td></td><td></td><td></td><td></td><td>Hopf (l<sub>1</sub> ≈ 1.522 × 10<sup>-2</sup>)</td></tr><tr><td>(σ<sub>2</sub>, σ<sub>3</sub>)</td><td>U</td><td>U</td><td>S</td><td></td><td>I<sub>3</sub></td></tr><tr><td>b = σ<sub>3</sub> ≈ 18.957</td><td></td><td></td><td></td><td></td><td>Hopf (l<sub>1</sub> ≈ -1.974 × 10<sup>-2</sup>)</td></tr><tr><td>(σ<sub>3</sub>, +∞)</td><td>U</td><td>S</td><td></td><td></td><td>I<sub>1</sub></td></tr></tbody></table>

Fig. 12 shows that a subcritical Hopf bifurcation occurs at $b = \sigma_2$ by increasing $b$ from zero where $E_1$ destabilizes. By decreasing $b$ from $\sigma_2$, the limit cycle $C_2$ appears unstable since the first Lyapunov coefficient is positive (see Table 19). This bifurcation corresponds to the passage from the region $I_3$ to $I_1^b$ through the $\Gamma_H$ curve.

As $b$ decreases, the radius of the limit cycle $C_1$ increases until it disappears at $b = \sigma_1$ through an LPC bifurcation. This bifurcation corresponds to the passage from the region
$\mathcal{I}_1^b$ to $\mathcal{I}_1$ through the LPC curve. As $b$ increases from this LPC bifurcation, the limit cycle $C_1$ appears stable (see Fig. 12). Then, the radius of this cycle limit increases and then decreases until its collapse to the positive equilibrium $E_1$ at $b = \sigma_3$ by a supercritical Hopf bifurcation since the first Lyapunov coefficient is negative (see Table 19). This bifurcation corresponds to the passage from the region $\mathcal{I}_3$ to $\mathcal{I}_1$ through the $\Gamma_H$ curve. When $b > \sigma_3$, the positive equilibrium $E_1$ recovers its stability. The type of the bifurcations with their critical bifurcation values, as well as the Lyapunov coefficient and asymptotic behavior of system (1.1) as a function of $b$ are provided in Table 19.

## 6. Conclusion

In this work, we have extended our mathematical study in 5 by considering distinct yields in the flocculation model (1.1) involving the attachment and detachment dynamics of isolated and aggregated bacteria in the presence of a single resource in a chemostat. Considering distinct removal rates and without ignoring the yield coefficients, we have provided a complete analysis of the existence and local asymptotic stability of all steady states for general monotonic growth rates. Using the necessary and sufficient conditions of existence and stability, we have studied theoretically and numerically the operating diagrams of the flocculation model (1.1) according to the operating parameters, which are the dilution rate $D$ and the input concentration of the substrate $S_{in}$ or the rates of attachment $a$ and detachment $b$.

To have a better understanding of the theoretical study of the operating diagram in the $(S_{in}, D)$-plane, we start with a simple case where the positive steady state is unique and stable if it exists (see Section 3.1). It is revealed that there can only be two regions: the region $\mathcal{I}_0$ of the washout ($E_0$ is the only steady state) or the region $\mathcal{I}_1$ of coexistence of isolated and attached bacteria around the positive steady state (the only steady states are $E_0$ which is unstable and $E_1$ which is LES).

Next, we have considered a case with the emergence of two positive steady states and the destabilization of one positive steady state via a Hopf bifurcation (see Section 3.2). The operating diagram in the $(S_{in}, D)$-plane shows the emergence of the green region $\mathcal{I}_2$ corresponding to the bistability between $E_0$ and $E_1^1$, and the blue region $\mathcal{I}_3$ corresponding to the destabilization of the positive steady state $E_1^1$, where there can be coexistence around a stable limit cycle. In Appendix E, we have considered a similar case for the set of parameter values in.⁶ Indeed, we obtain the operating diagram in Fig. 1(b-c), which is similar to one in Fig. 18. However, the region $\mathcal{I}_3$ has not been detected numerically in 6 because of its size, where the maximum value of $D$ is around $10^{-5}$. To detect it with good accuracy, we have changed the default value of “Digits” in MAPLE to 20.

Then, we have considered another case with the emergence of two positive steady states and the destabilization of one positive steady state via a Hopf bifurcation (see Section 3.3). The operating diagram in the $(S_{in}, D)$-plane is divided into five regions where there can be one more behavior (yellow region $\mathcal{I}_4$): the bistability with either coexistence around a stable limit cycle or the washout of the isolated and attached bacteria according to the initial condition.

Using the software MATCONT 56, we found numerically the operating diagram obtained theoretically in the case considered in Sections 3.2 and 3.3. However, we have also detected new bifurcations with two parameters like those of type Bogdanov-Takens (BT) or Cusp (CP). Moreover, the one bifurcation diagram shows the various types of bifurcations by crossing the different regions in the two-dimensional $(S_{in}, D)$-plane. Then, the study of the operating diagram with the two control parameters $(S_{in}, D)$ using MATCONT gives a more general vision of the asymptotic behavior of the system compared to the study of the bifurcation diagram as a function of $S_{in}$.

Finally, we analyze the effect of flocculation and deflocculation on the size and shape of various regions of the operating diagram in the $(S_{in}, D)$-plane. Decreasing the rates of attachment and/or detachment, the regions $\mathcal{I}_2$ and $\mathcal{I}_3$ are reduced until their disap-
24 Fekih-Salem and Sari

pearance. In the limiting case $a = b = 0$, we obtain the operating diagram of the classic chemostat model, where the CEP asserts that generically at most one species can survive the competition.

The operating diagram in the (a, b)-plan shows the emergence of a novel region $I_1^b$ where the system exhibits bistability between one cycle limit and the coexistence steady state $E_1$. In addition, MATCONT detects a new two-parameter bifurcation of type Generalized Hopf with the destabilization of $E_1$ via a subcritical Hopf bifurcation and the appearance of an unstable limit cycle that can disappear through a bifurcation of type Limit Point of Cycles.

One concludes that the flocculation process promotes the coexistence of isolated and attached bacteria of a microbial species around a limit cycle or positive steady state. This flocculation mechanism also favors bistability, where the asymptotic behavior of the solutions depends on the initial condition.

Studying the flocculation model with multiple species is quite complex and warrants further research. More precisely, we assume that $n$ species are competing on the same limiting resource, and that each species is present in two forms: isolated bacteria, of density $u_i$, and bacteria in flocks, of density $v_i$, for $i = 1, \dots, n$. We assume that isolated bacteria can stick with isolated bacteria and flocks to form new flocks, with rate $\alpha_i(\cdot)u_i$. In addition, flocks can split and liberate isolated bacteria with a rate $\beta_i(\cdot)v_i$. Thus, the flocculation model with $n$ species is given by the following system:

$$ egin{cases} 
\dot{S} = D(S_{in} - S) - \sum_{i=1}^{n} \left( \frac{1}{y_{ui}} f_i(S)u_i + \frac{1}{y_{vi}} g_i(S)v_i \right) \\ 
\dot{u}_i = (f_i(S) - D_{ui})u_i - \alpha_i(\cdot)u_i + \beta_i(\cdot)v_i, \quad 1 \le i \le n \\ 
\dot{v}_i = (g_i(S) - D_{vi})v_i + \alpha_i(\cdot)u_i - \beta_i(\cdot)v_i 
\end{cases} 
\quad (6.1) $$

where $f_i(S)$ and $g_i(S)$ are, respectively, the growth rates of isolated and attached bacteria of each species $i$; $D_{ui}$ and $D_{vi}$ are, respectively, the disappearance rates of planktonic and attached bacteria of each species $i$; $y_{ui}$ and $y_{vi}$ are, respectively, “yield” constants reflecting the conversion of nutrient to planktonic and aggregated bacteria of each species $i$.

This model with 2 species in which one can flocculate (the one which would be the 'stronger' competitor if it would not flocculate) and the other, weaker one, does not flocculate has been studied in.⁴ However, the behavior of the process in the various regions of the operating diagram of model (6.1) with $n$ species including the mechanism of flocculation is a question of major interest and importance from the biological and ecological point of view. This question deserves further attention and will be the object of future work.

Acknowledgments

The first author thanks the financial support of Cimpa "research in pairs", the European Mathematical Society EMS-Simons for Africa, and also the I-SITE Excellence Program of the University of Montpellier for projects "Support for international mobility - EX-PLORE#4". We thank the Euro-Mediterranean research network TREASURE.

Appendix A. Proof of Theorem 2.1

To determine the existence of the steady states of (1.1), we set the right-hand sides of equations of this system equal to zero:

$$ 0 = D(S_{in} - S) - f(S)u/y_u - g(S)v/y_v, \quad (A.1) $$

$$ 0 = [f(S) - D_u]u - a(u + v)u + bv, \quad (A.2) $$

$$ 0 = [g(S) - D_v]v + a(u + v)u - bv. \quad (A.3) $$

The equation (A.2) implies $v = 0$ when $u = 0$ and inversely, $u = 0$ when $v = 0$ from equation (A.3). Thus, a steady state of exclusion of isolated or aggregated bacteria does not exist.
A flocculation model in the chemostat with distinct removal rates and yields 25

* For $E_0$, one has $u = v = 0$. Hence, (A.1) implies that $S = S_{in}$. Thus, $E_0 = (S_{in}, 0, 0)$ always exists.

* For $E_1$, one has $u > 0$, and $v > 0$. The sum of equations (A.2) and (A.3) results in

$$\varphi(S)u + \psi(S)v = 0,$$

where $\varphi$ and $\psi$ are defined in (2.3). Therefore, $u$ and $v$ are positive if and only if $\varphi(S)$ and $\psi(S)$ have opposite signs or equivalently $S \in (\lambda_u, \lambda_v)$ when $\lambda_u < \lambda_v$ and $S \in (\lambda_v, \lambda_u)$ otherwise. Consequently, one can write $v$ as a function of $S$ and $u$ such that

$$v = -\frac{\varphi(S)}{\psi(S)}u, \qquad (A.4)$$

where $\psi(S) \neq 0$ or equivalently $S \neq \lambda_v$. Using (A.4) to replace $u$ in (A.2), leads to

$$u = U(S), \qquad (A.5)$$

where the function $S \mapsto U$ is defined in (2.2). Hence, $u$ is positive if and only if $S \in I$, which is defined by (2.4). Using (A.4) to replace $u$ by (A.5), it follows that

$$v = V(S), \qquad (A.6)$$

where $v$ is defined in (2.2). Replacing $u$ by (A.5) and $v$ by (A.6) in (A.1) implies that $S$ is a solution of equation $D(S_{in} - S) = H(S)$ where the function $S \mapsto H(S)$ is defined by (2.1).

Now, to establish the local asymptotic stability of the two steady states of system (1.1), we determine the the Jacobian matrix $J$ at $(S, u, v)$:

$$J = \begin{bmatrix} -D - f'(S)u/y_u - g'(S)v/y_v & -f(S)/y_u & -g(S)/y_v \\ f'(S)u & \varphi(S) - a(2u + v) & -au + b \\ g'(S)v & a(2u + v) & \psi(S) + au - b \end{bmatrix}. \qquad (A.7)$$

* For $E_0$, the three eigenvalues are on the diagonal of $J$ evaluated at $(S_{in}, 0, 0)$:

$$-D, \quad \varphi(S_{in}), \quad \psi(S_{in}) - b.$$

Hence, $E_0$ is LES if and only if the three eigenvalues are negative, that is, if and only if $S_{in} < \min(\lambda_u, \lambda_b)$.

* For $E_1$, the Jacobian matrix $J$ evaluated at $E_1 = (S^*, u^*, v^*)$ is defined by

$$J_1 = \begin{bmatrix} -m_{11} & -m_{12} & -m_{13} \\ m_{21} & -m_{22} & a_{23} \\ m_{31} & m_{32} & -m_{33} \end{bmatrix}, \qquad (A.8)$$

where the coefficients of $J$ are given by (2.6) so that all $m_{ij}$ are positive for all $i, j = 1, \dots, 3$ with $(i, j) \neq (2, 3)$. Indeed, from (A.2) and (A.3), a straightforward calculation shows that

$$m_{22} = au^* + b\frac{v^*}{u^*} > 0, \quad m_{33} = a\frac{(u^*)^2}{v^*} > 0.$$

The characteristic polynomial is

$$P(\lambda) = \lambda^3 + c_1\lambda^2 + c_2\lambda + c_3,$$

where their coefficients are given by (2.5). Using [5, Lemma A.1], a direct calculation gives

$$c_3 = -\det(J_1) = -\varphi(S^*)(\psi(S^*) - b)(D + H'(S^*)).$$

From the Routh-Hurwitz criterion, we conclude that $E_1$ is LES if and only if

$$c_3 > 0 \quad \text{and} \quad c_4 = c_1c_2 - c_3 > 0,$$

since $c_1$ is positive. This completes the proof.
26 Fekih-Salem and Sari

<figure><img src="image_15.png" /><figcaption>Fig. 13. Case of Section 3.1: (a) curves of the function $H(S)$ and the corresponding curves of $c_4(S)$ in their existence domain when $D \in \{0.01, 0.1, 0.5, 1, 2, 3\}$ corresponding to red, blue, green, magenta, cyan, and grey colors, respectively. (b) Only the curves of the function $c_4(S)$.</figcaption></figure>

## Appendix B. Case of Section 3.1: positivity of the stability condition $c_4$

In the following, we show that the stability condition $c_4 > 0$ holds for the positive steady state $E_1$ in the case considered in Section 3.1 so that the curve $\Gamma_H$ corresponding to $c_4 = 0$ does not exist in the operating diagram of Fig. 1(a) for model (1.1). Fig. 13 illustrates the positivity of the function $c_4(S)$ for several values of $D$ and the corresponding curves of $H(S)$.

## Appendix C. Case of Section 3.2: sign of $c_4$

In the section, we show that the stability condition $c_4 > 0$ of the positive steady state $E_1$ is not always verified, so that $c_4(S)$ changes sign in the interval $I(D) = ]\lambda_v(D), \lambda_{BP}(D)]$ of the existence of $E_1$. Fig. 14(a) shows that the equation $c_4 = 0$ has two roots $S_H^1(D)$ and $S_H^2(D)$ as defined in Section 3.2 for all $D < D_H^{max}$. Let $D$ be fixed at $D = D^* = 0.142$ (the green curve in 14(a)). Fig. 14(b) shows the curve of the function $S \mapsto H(S)$ in red [resp. in blue] when the function $S \mapsto c_4(S)$ is positive [resp. negative]. More precisely, $c_4(S)$ is positive for all $S \in ]\lambda_v, S_H^2[ \cup ]S_H^1, \lambda_{BP}]$ and negative for all $S \in ]S_H^2, S_H^1[$ where $\lambda_v(D^*) \approx 0.078, \lambda_{BP}(D^*) = \lambda_b \approx 2.222$ and the critical values according to $S$ and corresponding to Hopf bifurcation are given by

$$S_H^2(D^*) \approx 1.284, \quad S_H^1(D^*) \approx 1.748.$$

These critical values are equivalent to the following critical values according to $S_{in}$,

$$S_{in}^{H2}(D^*) \approx 3.674, \quad S_{in}^{H1}(D^*) \approx 2.640, \text{ respectively.}$$

By increasing the value of $S_{in}$ from zero to $\lambda_{BP}(D^*)$, $E_1$ emerges LES via a Branch Point (BP) with $E_0$ when $S = S_{in} = \lambda_{BP}(D^*)$. Increasing $S_{in}$ further, $E_1$ destabilizes through the first Hopf bifurcation at $S_{in}^{H1}$ and remains unstable up to the value of $S_{in}^{H2}$. Finally, $E_1$ returns LES for all $S_{in} > S_{in}^{H2}$ via a second Hopf bifurcation.

## Appendix D. Case of Section 3.3

In this appendix, we give numerical evidence of the change of sign of $c_4(S)$ and the appearance of a stable limit cycle as well as its disappearance by a homoclinic bifurcation for fixed $D$. Fig. 15 shows the curve of the function $S \mapsto H(S)$ in red [resp. in blue] when the function $S \mapsto c_4(S)$ is positive [resp. negative] and $D$ is fixed at $D^* = 0.1$. More precisely, $c_4(S)$ is positive for all $S \in ]\lambda_v, S_H^2[ \cup ]S_H^1, S_{LP}]$ and negative for all $S \in ]S_H^2, S_H^1[$ where $\lambda_v(D^*) \approx 0.846$, and the critical values according to $S$ and corresponding to Hopf bifurcation are given by

$$S_H^2(D^*) \approx 1.963, \quad S_H^1(D^*) \approx 3.422, \quad S_{LP}(D^*) \approx 3.492.$$

These critical values are equivalent to the following critical values according to $S_{in}$,

$$S_{in}^{H2}(D^*) \approx 8.179, \quad S_{in}^{H1}(D^*) \approx 3.842, \quad \lambda_{LP}(D^*) \approx 3.837, \text{ respectively.}$$
<figure><img src="image_16.png" /><figcaption>Fig. 14. Case of Section 3.2: (a) the functions $H(S)$ and $c_4(S)$ on the existence domain of $E_1$ when $D \in \{0.09, 0.12, 0.142, D_H^{max}, 0.21\}$ ($D_H^{max} \approx 0.165$) corresponding to red, blue, green, magenta, and cyan colors, respectively. (b) The function $H(S)$ in red [resp. in blue] when the function $c_4(S)$ is positive [resp. negative], for $D = D^* = 0.142$.</figcaption></figure>

By increasing the value of $S_{in}$ from zero to $\lambda_{LP}(D^*)$, the two positive steady states $E_1^1$ and $E_1^2$ emerge LES and unstable, respectively, via a Limit Point (LP) bifurcation where $S = S_{LP}(D^*)$. Increasing $S_{in}$ further, $E_1^1$ destabilizes trough the first Hopf bifurcation at $S_{in}^{H1}$ and remains unstable up to the value of $S_{in}^{H2}$. Finally, $E_1^1$ returns LES for all $S_{in} > S_{in}^{H2}$ via a second Hopf bifurcation. In Fig. 15(b), we have chosen the red color for LES steady states and the blue color for unstable steady states.

<figure><img src="image_17.png" /><figcaption>Fig. 15. Case of Section 3.3: (a) the function $H(S)$ and $c_4(S)$ when $D = D^* = 0.1$ showing the changes of the sign of $c_4(S)$; (b) a magnification for $3.25 < S < \lambda_u = 4.061$ where $S_{in} = 3.86$.</figcaption></figure>

Fig. 16 shows the limit cycles in the three-dimensional phase plot $(S, u, v)$ for various values of $S_{in}$ between $\sigma_4$ and $\sigma_6$ (defined in Table 12) until their disappear by homoclinic bifurcation at $\sigma_4$. Decreasing $S_{in}$ from the Hopf bifurcation at $\sigma_6$, the radius of the limit cycle increases (see Fig. 16(a)). Decreasing again $S_{in}$, the radius of the limit cycle decreases until its disappearance by approaching a homoclinic orbit when $S_{in} = \sigma_4 \approx 4.03468$ (see Fig. 16(b-c)).

Fig. 17 reveals the homoclinic bifurcations at $S_{in} = \sigma_3$ and $S_{in} = \sigma_4$ where the time period $T$ of the limit cycle solutions of model (1.1) tends to $+\infty$ as $S_{in}$ tends to these critical values.

## Appendix E. Case of parameter set in Table 20 (line 8)

The main purpose of this appendix is to show that the region of destabilization of the positive steady state is omitted in the construction of the operating diagram in.⁶ With
<figure><img src="image_18.png" /><figcaption>Fig. 16. Case of Section 3.3: the three-dimensional space (S, u, v) in MATCONT when D = 0.1: (a) a family of limit cycles, starting from a Hopf point (H) at σ₆ = 8.179 and decreasing Sᵢₙ until the maximum radius. (b) A family of limit cycles, starting from the maximum radius and approaching a homoclinic orbit by decreasing Sᵢₙ. (c) A homoclinic orbit for Sᵢₙ = σ₄ ≈ 4.03468.</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>Fig. 17. Case of Section 3.3: a plot in MATCONT of time period T of the limit cycle solutions of model (1.1) for D = 0.1 starting from the Hopf bifurcation at (a) σ₂ [(b) resp. σ₆]; homoclinic bifurcation at Sᵢₙ = σ₃ ≈ 3.8477 [resp. Sᵢₙ = σ₄ ≈ 4.03468].</figcaption></figure>

the same set of parameters in 6, see Table 20 (line 8), we find the operating diagram in Fig. 18, which is similar to that in Fig. 1(b-c) but where λ<sub>BP</sub>(D) = λ<sub>u</sub>(D). Thus, the existence and local stability of all steady states of model (1.1) in the four regions I<sub>1</sub>, i = 0, . . . , 3 of the operating diagram shown in Fig. 18 can be obtained from Table 6. Note that the numbering of the I<sub>1</sub> and I<sub>2</sub> regions is reversed in.⁶ Similarly, for the regions I<sub>3</sub> and I<sub>4</sub>. Next, we will show that the region I<sub>3</sub> corresponds to the emergence of the stable limit cycle via Hopf bifurcations.

Since the order of D can reach 10⁻¹⁰ to plot the curve Γ<sub>H</sub>, we modified the value of “Digits” in MAPLE to 20 instead of the default value 10 to avoid the introduction of round-off error. This allows, in particular, to have precision in the tracing of the Γ<sub>H</sub> curve and that of the function c₄(S) for fairly small D.

To give numerical evidence of the Hopf bifurcation occurring through the transition through the curve Γ<sub>H</sub> from region I<sub>1</sub> to region I<sub>3</sub>, we determine numerically the eigenvalues of the Jacobian matrix of system (1.1) at E<sub>1</sub><sup>1</sup> by computing the roots of the characteristic polynomial as we vary the parameter S<sub>in</sub>. Let D be fixed such that D = D* = 2.5 × 10<sup>-5</sup>. Fig. 19(a) shows that one eigenvalue denoted by λ<sub>1</sub>(S<sub>in</sub>) is real and remains negative for all S<sub>in</sub> ∈ [η*, η<sub>f</sub>] where η* = λ<sub>u</sub>(D*) ≈ 0.625 denotes the value of S<sub>in</sub> at which the positive steady state appears and η<sub>f</sub> = 100 denotes the final value of the variation of S<sub>in</sub>. Fig. 19(b) shows that the two other eigenvalues λ<sub>2</sub>(S<sub>in</sub>) and λ<sub>3</sub>(S<sub>in</sub>) defined by

$$
\lambda_{2,3}(S_{in}) = \alpha_{2,3}(S_{in}) \pm i\beta_{2,3}(S_{in}), \quad \text{for all } S_{in} \in [\eta^*, \eta_f]
$$
<figure><img src="image_20.png" /><figcaption>Fig. 18. MAPLE: (a) operating diagram of (1.1) in case of Table 20 (line 8). (b) Magnification on the region $\mathcal{I}_3$ and the curve $\Gamma_H = \Gamma_H^1 \cup \Gamma_H^2$ when $D \in [0, 0.000035]$.</figcaption></figure>

are complex-conjugate so that the real part $\alpha_{2,3}(S_{in})$ is negative for all $S_{in} \in [\eta^*, \eta_1) \cup (\eta_2, \eta_f]$ and positive for all $S_{in} \in (\eta_1, \eta_2)$. When $S_{in} = \eta_i$, $i = 1, 2$, the pair $\lambda_{2,3}(\eta_i)$ is purely imaginary such that $\alpha_{2,3}(\eta_i) = 0$, with $\beta_{2,3}(\eta_i) \neq 0$. Moreover, the following transversality condition is checked numerically

$$ \frac{d\alpha_{2,3}}{dS_{in}}(\eta_1) > 0 \quad \text{and} \quad \frac{d\alpha_{2,3}}{dS_{in}}(\eta_2) < 0. \qquad (E.1) $$

That is, the two complex-conjugate eigenvalues cross the imaginary axis with non-zero speed. Thus, the positive steady state $E_1$ is destabilized via two Hopf bifurcations with the occurrence or disappearance of a stable limit cycle when $S_{in}$ increases and crosses the critical values $\eta_1$ and $\eta_2$. This result is consistent with the numerical simulation in Fig. 20(b) showing the emergence of a stable limit cycle where the oscillations are sustained.

<figure><img src="image_21.png" /><figcaption>Fig. 19. Case of the parameter set in Table 20 (line 8): variation of $S_{in}$ from $\eta^*$ to $\eta_f$ when $D = D^*$; (a) the real eigenvalue $\lambda_1(S_{in})$. (b) The pair of complex-conjugate eigenvalues $\lambda_{2,3}(S_{in})$.</figcaption></figure>

Recall that $D$ is fixed at $D = D^* = 2.5 \times 10^{-5}$. Fig. 20(a) illustrates the convergence towards $E_1$ in the three-dimensional phase space $(S, u, v)$ when $S_{in} = 1$ where the pair of complex-conjugate eigenvalues has negative real parts. In this case, the point $(S_{in}, D)$ belongs to region $\mathcal{I}_1$ where there are only two steady states: $E_0$ is unstable while $E_1$ is LES.

Fig. 20(b) illustrates the convergence towards a stable limit cycle when $S_{in} = 48$ where the pair of complex-conjugate eigenvalues have positive real parts. In this case, the point $(S_{in}, D)$ belongs to region $\mathcal{I}_3$ where there are only two steady states $E_0$ and
<figure><img src="image_22.png" /><figcaption>Fig. 20. Case of the parameter set in Table 20 (line 8): the three-dimensional space (S, u, v) in MATCONT; (a) convergence to the positive steady state $E_1$ when $S_{in} = 1$; (b) convergence to the stable limit cycle when $S_{in} = 48$.</figcaption></figure>

$E_1$ which are unstable. To solve the problem of the calculation time of the solution of (1.1) until convergence to the limit cycle where $D$ is small enough, we have changed the default solver “ode45” to “ode23” in MATCONT.

<figure><img src="image_23.png" /><figcaption>Fig. 21. Case of the parameter set in Table 20 (line 8): curve of the function $H(S)$ where $c_4(S)$ is positive when $S \in [0.622, S_H^2 \cup [S_H^1, \lambda_u]]$ and negative when $S \in [S_H^2, S_H^1[$. Magnifications when (b) $(S, H(S)) \in [0.622, \lambda_u] \times [-5 \times 10^{-6}, 2.1 \times 10^{-3}]$ and (c) $(S, H(S)) \in [0.620, \lambda_u] \times [-2 \times 10^{-6}, 4 \times 10^{-5}]$.</figcaption></figure>

Fig. 21 illustrates the curve of the function $S \mapsto H(S)$ in red [resp. in blue] when the function $S \mapsto c_4(S)$ is positive [resp. negative] and $D$ is fixed at $D^* = 2.5 \times 10^{-5} \in ]0, D_H^{max}]$. The solutions $S_H^1$ and $S_H^2$ of the equation $c_4(S) = 0$ correspond to the critical values $S_{in}^{H1}$ and $S_{in}^{H2}$. They are the intersections of the horizontal line of equation $D = D^*$ with the curves $\Gamma_H^1$ and $\Gamma_H^2$, respectively, in the $(S_{in}, D)$-plane of the operating diagram in Fig. 18(b). By increasing the value of $S_{in}$ from zero to $\eta^* = \lambda_u(D^*) \approx 0.625$, $E_0$ becomes unstable by a Branch Point (BP) with $E_1$ that appears stable until the first Hopf bifurcation at $S_{in}^{H1} \approx 28.990$ (or equivalently $S = S_H^1 \approx 0.62398$). Then, $E_1$ remains unstable up to the value of $S_{in}^{H2} \approx 64.878$ (or equivalently $S = S_H^2 \approx 0.62267$), that is, for all $S \in [S_H^2, S_H^1[$. Finally, for $S_{in} > S_{in}^{H2}$ (or equivalently $S \in ]\lambda_v(D^*), S_H^2[$), $E_1$ returns stable via a second Hopf bifurcation.

## Appendix F. Parameter values used in numerical simulations

All the parameter values used in the numerical simulations are provided in Table 20.
A flocculation model in the chemostat with distinct removal rates and yields 31

<table><caption>Table 20. Parameter values used for system (1.1) when the specific growth rates *f* and *g* are given by (3.1) and *y*<sub>u,v</sub> = 1. The abbreviation *Var* means *Variable*.</caption><thead><tr><th>Parameter<br>Figures</th><th>m<sub>1</sub><br>(h<sup>-1</sup>)</th><th>k<sub>1</sub><br>(g/l)</th><th>m<sub>2</sub><br>(h<sup>-1</sup>)</th><th>k<sub>2</sub><br>(g/l)</th><th>a<br>(l/h/g)</th><th>b<br>(h<sup>-1</sup>)</th><th>α</th><th>β</th><th>m<sub>u</sub><br>(h<sup>-1</sup>)</th><th>m<sub>v</sub><br>(h<sup>-1</sup>)</th></tr></thead><tbody><tr><td>1(a), 13</td><td>4.5</td><td>1</td><td>3</td><td>2.7</td><td>2</td><td>3</td><td>0.8</td><td>0.5</td><td>0.2</td><td>0.25</td></tr><tr><td>1(b-c), 3, 14</td><td></td><td></td><td></td><td></td><td>4</td><td>2</td><td></td><td>0.9</td><td></td><td>0</td></tr><tr><td>2, 4-6, 15-17</td><td></td><td></td><td></td><td></td><td>4</td><td>2</td><td></td><td>1</td><td></td><td>1</td></tr><tr><td>7, 8, 9</td><td></td><td></td><td></td><td></td><td>Var</td><td>Var</td><td></td><td>1</td><td></td><td>1</td></tr><tr><td>10</td><td>5</td><td>2</td><td>5</td><td>3</td><td>Var</td><td>2</td><td>1</td><td>1</td><td>3.25</td><td>1</td></tr><tr><td>11</td><td></td><td></td><td></td><td></td><td>20</td><td>Var</td><td></td><td>1</td><td></td><td>1</td></tr><tr><td>12</td><td></td><td></td><td></td><td></td><td>30</td><td>Var</td><td></td><td>1</td><td></td><td>1</td></tr><tr><td>18-21</td><td>3.5</td><td>2.5</td><td>3</td><td>1.5</td><td>1</td><td>1</td><td>1</td><td>0.75</td><td>0.7</td><td>0.4</td></tr></tbody></table>

References

1. H. L. Smith and P. Waltman, *The Theory of the Chemostat: Dynamics of Microbial Competition*, Cambridge University Press, Cambridge, UK (1995).

2. G. Robledo, F. Grognard and J-L. Gouzé, Global stability for a model of competition in the chemostat with microbial inputs, *Nonlinear Anal.: Real World Appl.*, **13**, 582–598 (2012).

3. R. Fekih-Salem, J. Harmand, C. Lobry, A. Rapaport and T. Sari, Extensions of the chemostat model with flocculation, *J. Math. Anal. Appl.*, **397**, 292–306 (2013).

4. R. Fekih-Salem, A. Rapaport and T. Sari, Emergence of coexistence and limit cycles in the chemostat model with flocculation for a general class of functional responses, *Appl. Math. Modell.*, **40**, 7656–7677 (2016).

5. R. Fekih-Salem and T. Sari, Properties of the chemostat model with aggregated biomass and distinct removal rates, *SIAM J. Appl. Dyn. Syst. (SIADS)*, **18**, 481–509 (2019).

6. R. Fekih-Salem and T. Sari, Operating diagram of a flocculation model in the chemostat, *ARIMA Journal*, **31**, 45–58 (2020).

7. B. Haegeman and A. Rapaport, How flocculation can explain coexistence in the chemostat, *J. Biol. Dyn.*, **2**, 1–13 (2008).

8. N. Abdellatif, R. Fekih-Salem and T. Sari, Competition for a single resource and coexistence of several species in the chemostat, *Math. Biosci. Eng.*, **13**, 631–652 (2016).

9. P. De Leenheer, D. Angeli and E. D. Sontag, Crowding effects promote coexistence in the chemostat, *J. Math. Anal. Appl.*, **319**, 48–60 (2006).

10. R. Fekih-Salem, C. Lobry and T. Sari, A density-dependent model of competition for one resource in the chemostat, *Math. Biosci.*, **286**, 104–122 (2017).

11. L. Imhof and S. Walcher, Exclusion and persistence in deterministic and stochastic chemostat models, *J. Differ. Equations*, **217**, 26–53 (2005).

12. C. Lobry and J. Harmand, A new hypothesis to explain the coexistence of *n* species in the presence of a single resource, *C. R. Biol.*, **329**, 40–46 (2006).

13. C. Lobry, F. Mazenc and A. Rapaport, Persistence in ecological models of competition for a single resource, *C. R. Acad. Sci. Paris Ser. I*, **340**, 199–204 (2005).

14. T. Mtar, R. Fekih-Salem and T. Sari, Interspecific density-dependent model of predator-prey relationship in the chemostat, *Int. J. Biomath.*, **14**, 2050086 (2021).

15. T. Mtar, R. Fekih-Salem and T. Sari, Mortality can produce limit cycles in density-dependent models with a predator-prey relationship, *Discrete and Continuous Dyn. Syst. - B*, **27**, 7445–7467 (2022).

16. B. Bar and T. Sari, The operating diagram for a model of competition in a chemostat with an external lethal inhibitor, *Discrete and Continuous Dyn. Syst. - B*, **25**, 2093–2120 (2020).
17. M. Dellal and B. Bar, Global analysis of a model of competition in the chemostat with internal inhibitor, *Discrete and Continuous Dyn. Syst.* - B, **26**, 1129–1148 (2021).

18. M. Dellal, B. Bar and M. Lakrib, A competition model in the chemostat with allelopathy and substrate inhibition, *Discrete and Continuous Dyn. Syst.* - B, **27**, 2025–2050 (2022).

19. M. Dellal, M. Lakrib and T. Sari, The operating diagram of a model of two competitors in a chemostat with an external inhibitor, *Math. Biosci.*, **302**, 27–45 (2018).

20. N. E. H. Zitouni, M. Dellal and M. Lakrib, Substrate inhibition can produce coexistence and limit cycles in the chemostat model with allelopathy, *J. Math. Biol.*, **87**:7, 35–41 (2023).

21. M. P. Boer, B. W. Kooi and S. A. L. M. Kooijman, Food chain dynamics in the chemostat, *Math. Biosci.*, **150**, 43–62 (1998).

22. G. J. Butler, S. B. Hsu and P. Waltman, Coexistence of competing predators in a chemostat, *J. Math. Biology*, **17**, 133–151 (1983).

23. G. J. Butler and G. S. K. Wolkowicz, Predator-mediated competition in the chemostat, *J. Math. Biol.*, **24**, 167–191 (1986).

24. S.-B. Hsu, C. A. Klausmeier and C.-J. Lin, Analysis of a model of two parallel food chains, *Discrete and Continuous Dyn. Syst.* - B, **12**, 337–359 (2009).

25. G. S. K. Wolkowicz, Successful invasion of a food web in a chemostat, *Math. Biosci.*, **93**, 249–268 (1989).

26. J. W. Costeron, Overview of microbial biofilms, *J Ind Microbiol*, **15**, 137–40 (1995).

27. IWA Task Group on Biofilm Modeling, *Mathematical Modeling of Biofilms*, IWA Publishing, London (2006).

28. D.N. Thomas, S.J. Judd and N. Fawcett, Flocculation modelling: a review, *Water Research*, **33**, 1579–1592 (1999).

29. R. Freter, H. Brickner, J. Fekete, M. M. Vickerman and K. E. Carey, Survival and implantation of Escherichia coli in the intestinal tract, *Infect Immun*, **39**, 686–703 (1983).

30. D. Jones, H. Kojouharov, D. Le and H. L. Smith, The Freter model: A simple model of biofilm formation, *J. Math. Biol.*, **47**, 137–152 (2003).

31. M. Ballyk and H. L. Smith, A model of microbial growth in a plug flow reactor with wall attachment, *Math Biosci.*, **158**, 95–126 (1999).

32. S. Pilyugin and P. Waltman, The Simple Chemostat with Wall Growth, *SIAM J. Appl. Math.*, **59**, 1552–1572 (1999).

33. E. D. Stemmons and H. L. Smith, Competition in a chemostat with wall attachment, *SIAM J. Appl. Math.*, **61**, 567–595 (2000).

34. B. Haegeman, C. Lobry and J. Harmand, Modeling bacteria flocculation as density-dependent growth, *AIChE J*, **53**, 535–539 (2007).

35. A. Rapaport, Properties of the chemostat model with aggregated biomass, *Eur. J. Appl. Math.*, **29**, 972–990 (2018).

36. J. Harmand, C. Lobry, A. Rapaport, and T. Sari, *The Chemostat: Mathematical Theory of Microorganism Cultures*, Chemical Eng. Ser., ISTE-Wiley, New York (2017).

37. S. Shen, G. C. Premier, A. Guwy and R. Dinsdale, Bifurcation and stability analysis of an anaerobic digestion model, *Nonlinear Dynam.*, **48**, 391–408 (2007).

38. M. J. Wade, J. Oakley, S. Harbisher, N. G. Parker and J. Dolfing, MI-Sim: A MATLAB package for the numerical analysis of microbial ecological interactions, *PLoS ONE*, **12**, e0173249 (2017).

39. M. Hanaki, J. Harmand, Z. Mghazli, A. Rapaport, T. Sari and P. Ugalde, Mathematical study of a two-stage anaerobic model when the hydrolysis is the limiting step, *Processes*, **9**, 2050 (2021).
40. Z. Khedim, B. Benyahia, B. Cherki, T. Sari and J. Harmand, Effect of control parameters on biogas production during the anaerobic digestion of protein-rich substrates, *Appl. Math. Model.*, **61**, 351–376 (2018).

41. M. J. Wade, R. W. Pattinson, N. G. Parker and J. Dolfing, Emergent behaviour in a chlorophenol-mineralising three-tiered microbial ‘food web’, *J. Theor. Biol.*, **389**, 171–186 (2016).

42. M. Weedermann, G. S. K. Wolkowicz and J. Sasarag, Optimal biogas production in a model for anaerobic digestion, *Nonlinear Dyn*, **81**, 1097–1112 (2015).

43. A. Xu, J. Dolfing, T. P. Curtis, G. Montague and E. Martin, Maintenance affects the stability of a two-tiered microbial ‘food chain’?, *J. Theor. Biol.*, **276**, 35–41 (2011).

44. A. Dhooge, W. Govaerts, Yu. A. Kuznetsov, H. G. E. Meijer and B. Sautois, New features of the software MatCont for bifurcation analysis of dynamical systems, *Math. Comput. Model. Dyn.*, **14**, 147–175 (2008).

45. MAPLE [Software], *Waterloo Maple Inc.*, Waterloo, Ontario, version 17.0.0.0 (2018)

46. M. Dali Youcef, A. Rapaport and T. Sari, Study of performance criteria of serial configuration of two chemostats, *Math Biosci Eng.*, **17**, 6278–6309 (2020).

47. M. Dali-Youcef, A. Rapaport and T. Sari, Performance study of two serial interconnected chemostats with mortality, *Bull. Math. Biol.* **84**, 110 (2022).

48. Y. Daoud, N. Abdellatif, T. Sari and J. Harmand, Steady state analysis of a syntrophic model: The effect of a new input substrate concentration, *Math. Model. Nat. Phenom.*, **13**, 1–22 (2018).

49. R. Fekih-Salem, Y. Daoud, N. Abdellatif and T. Sari, A mathematical model of anaerobic digestion with syntrophic relationship, substrate inhibition and distinct removal rates, *SIAM J. Appl. Dyn. Syst. (SIADS)*, **20**, 1621–1654 (2021).

50. N. Hamidi, M. Dellal and M. Lakrib, Competition between plasmid-bearing and plasmid-free organisms in a chemostat with distinct removal rates and an external inhibitor, *J. Biol. Syst.* (2024).

51. S. Nouaoura, R. Fekih-Salem, N. Abdellatif and T. Sari, Operating diagrams for a three-tiered microbial food web in the chemostat, *J. Math. Biol.*, **85**, 7445–7467 (2022).

52. T. Sari, Commensalism and syntrophy in the chemostat: a unifying graphical approach, *AIMS Mathematics*, **9**, 18625–18669 (2024).

53. T. Sari and B. Benyahia, The operating diagram for a two-step anaerobic digestion model, *Nonlinear Dyn.*, **105**, 2711–2737 (2021).

54. T. Sari and J. Harmand, A model of a syntrophic relationship between two microbial species in a chemostat including maintenance, *Math. Biosci.*, **275**, 1–9 (2016).

55. T. Sari and M. J. Wade, Generalised approach to modelling a three-tiered microbial food-web, *Math. Biosci.*, **291**, 21–37 (2017).

56. MATCONT [Software], W. Govaerts, Y. A. Kuznetsov and H. G. E. Meijer, *Software version 9.4 (2023)*.

57. Yu A. Kuznetsov, Numerical Normalization Techniques for All Codim 2 Bifurcations of Equilibria in Ode's, *SIAM Journal on Numerical Analysis*, **36**, 1104–1124 (1999).

58. SCILAB [Software], *Enterprises SAS* version 6.1.1 (2021).