# Analysis and operating diagram of an interspecific density-dependent model *

TAHANI MTAR † and RADHOUANE FEKIH-SALEM *†‡

**Abstract.** This paper studies a two microbial species model in competition for a single resource in the chemostat including general interspecific density-dependent growth rates with distinct removal rates for each species. We give the necessary and sufficient conditions of existence, uniqueness, and local stability of all steady states. We show that a positive steady state, if it exists, then it is unique and unstable. In this case, the system exhibits a bi-stability where the behavior of the process depends on the initial condition. Our mathematical analysis proves that at most one species can survive which confirms the competitive exclusion principle. We conclude that adding only interspecific competition in the classical chemostat model is not sufficient to show the coexistence of two species even considering mortality in the dynamics of two species. Otherwise, we focus on the study, theoretically and numerically, of the operating diagram which depicts the existence and the stability of each steady state according to the two operating parameters of the process which are the dilution rate and the input concentration of the substrate. Using our mathematical analysis, we construct analytically the operating diagram by plotting the curves that separate their various regions. Our numerical method using MATCONT software validates these theoretical results but it reveals new bifurcations that occur by varying two parameters as Bogdanov-Takens and Zero-Hopf bifurcations. The bifurcation analysis shows that all steady states can appear or disappear only through transcritical bifurcations.

**Key words.** Bi-stability, Bifurcation, Chemostat, MATCONT, Operating Diagram.

**MSC codes.** 34A34, 34D20, 37N25, 92B05

**1. Introduction.** Invented by Monod [33] and almost simultaneously by Novick and Szilard [37] in 1950, the chemostat is a bioreactor fed continuously with a nutrient flow for the cultivation of cells or microorganisms. More precisely, the chemostat represents a device consisting of an enclosure, containing the reaction volume, connected by input to feed the system with resources, and an outlet through which microorganisms are removed. It is the source of several mathematical models describing various biological and ecological interactions between microbial species in the chemostat, for instance: mutualism [45, 46], commensalism [3, 4, 42], syntrophy [8, 13, 14, 40, 48], predator-prey relationship [5, 25], etc. Some of these models have been studied to understand and exploit the anaerobic digestion process used in wastewater and waste treatment. It includes a large number of species that coexist with a very complex relationship (see [14, 32], and the references therein).

The classic chemostat model describes the competition of two or more microbial populations on a single non-reproducing nutrient [44]. It is well known that the *Competitive Exclusion Principle* (CEP) predicts, in a generic case, that only one of them can survive. That is, only the species with the lowest "break-even" concentration survives while all other species will be

<footnote>*Submitted to the editors 2024-01-03.</footnote>

<footnote>**Funding:** This work was supported by the Euro-Mediterranean research network TREASURE (http://www.inra.fr/treasure).</footnote>

<footnote>†University of Tunis El Manar, National Engineering School of Tunis, LAMSIN, Tunisia (tahani.mtar@enit.utm.tn).</footnote>

<footnote>‡University of Monastir, Higher Institute of Computer Science of Mahdia, Tunisia (radhouane.fekih-salem@enit.utm.tn).</footnote>
excluded [44]. Although this theoretical prediction has been corroborated by the experiences of Hansen and Hubell [20], this principle (CEP) is not compatible with the great biodiversity found in nature. Thus, to reconcile the theory and the experimental results, several suggestions and new hypotheses were considered in the classic chemostat model. More precisely, to promote the coexistence of species in microbial ecosystems, many recent researches take into account various mechanisms, such as the inter and intraspecific competition [1, 9], the flocculation of the species or wall attachment [15, 17, 18, 19, 38], the density-dependence of the growth functions [6, 16, 26, 27, 28, 29, 34, 36], etc.

Historically, an extension of the classical chemostat model was considered in a series of papers by Lobry et al. [22, 26, 27, 28, 29]. This model is the following one

$$ (1.1) \qquad \left\{ \begin{aligned} \dot{S} &= D(S_{in} - S) - \sum_{i=1}^{n} \mu_{i}(S, X_{1}, \dots, X_{n})X_{i}, \\ \dot{X}_{i} &= (\mu_{i}(S, X_{1}, \dots, X_{n}) - D_{i})X_{i}, \quad i = 1, \dots, n \end{aligned} \right. $$

where $S$ and $X_i$ denote, respectively, the concentrations of the substrate and the concentration of species $i$ present in the reactor at time $t$; $S_{in}$ and $D$ denote, respectively, the concentration of substrate in the feed bottle and the dilution rate of the chemostat; $\mu_i(S, X_1, \dots, X_n)$ represents the density-dependent growth rate of species $i$ which it depends not only on $S$ but also on the concentrations of species $i = 1, \dots, n$. The function $\mu_i$ is assumed to be increasing in the variable $S$ and decreasing in each variable $X_j$. The removal rate $D_i$ can be interpreted as the sum of the dilution rate and the natural death rate of species $i$. In [27, 28, 29], the authors considered the case where $\mu_i(S, X_i)$ depends only on $S$ and the concentration of species $i$ (that is, only intraspecific interference is present). They proved that the coexistence of all species is possible. A global study of model (1.1) with $D_i = D$ for all species was done in [28]. The case with distinct removal rates ($D_i \neq D$ for all species) was considered in [27, 29]. For particular growth rates, general model (1.1) was studied only numerically in [22, 26]. They demonstrated that the coexistence of species can occur for small enough interspecific interference but it is no longer the case when it is sufficiently large. It was proved how the intraspecific interference is a mechanism of coexistence for competing species in the chemostat.

In this context, El Hajji et al. [11, 12] focused on the effect of interspecific interference on the coexistence of microbial species in the chemostat. They considered a particular case of model (1.1) with $n = 2$. With the same dilution rate, the two-species model is written as follows:

$$ (1.2) \qquad \left\{ \begin{aligned} \dot{S} &= D(S_{in} - S) - \mu_{1}(S, X_{2})X_{1} - \mu_{2}(S, X_{1})X_{2}, \\ \dot{X}_{1} &= (\mu_{1}(S, X_{2}) - D)X_{1}, \\ \dot{X}_{2} &= (\mu_{2}(S, X_{1}) - D)X_{2}, \end{aligned} \right. $$

where the growth function $\mu_i(S, X_j)$ can depends now only on the substrate $S$ and the concentration of the species $X_j$, $i \neq j$. In [11], three types of interspecific interactions between species are considered to explain the coexistence. The first type is the predator-prey relationship which is characterized by the fact that one species (the prey) promotes the growth of the other species (the predator) which in turn inhibits the growth of the prey species. In [11], it has shown that system (1.2) may exhibit the coexistence with a multiplicity of positive steady
states. In this case of interaction, our study in [36] has shown that adding the mortality of species in system (1.2) can destabilize the positive steady state through Hopf bifurcation where the coexistence can be around a stable limit cycle and not only around a positive steady state as the case without mortality. The second type is the obligate mutualistic relationship where each species promotes the growth of the other species. That is, the function $\mu_i$ is increasing in $S$ and $X_j$. Using general monotonic growth rates, the model exhibits the coexistence where the positive steady state is not necessarily unique [11, 12]. Finally, the third type is the mutual inhibitory relationship which is characterized by the fact that each species inhibits the growth of the other species. In [11], it was proved that the coexistence of two species is impossible which confirms the CEP.

Here, we focus attention on the third type of interaction which describes the mutual inhibitory relationship between species. The present paper studies the competition model of two microbial species on a single nutriment (1.2) but with distinct removal rates. By comparing with the existing literature, our study provides an extension of the results in [11] by considering the mortality of each species.

In the first step, we want to understand the effect of decay on the behavior of the density-dependent model (1.2) considered in [11]. As is well known in the literature, the addition of mortality terms of the species can cause changes in the behaviors of the system. For instance in predator-prey models in the chemostat, stable limit cycles [5, 36] or multiple chaotic attractors [24] can be found when the mortality of species not neglected. So, can the simultaneous presence of mortality and interspecific interference induce the stable coexistence of competitors in a chemostat?

In Fekih-Salem et al. [16], a model of competition of two species with inter and intraspecific density-dependent growth functions and distinct removal rates has been considered. However, the particular case with only interspecific density-dependent growth rates has not been studied theoretically or numerically.

Here, we provide a complete mathematical and numerical analysis of the interspecific density-dependent model (1.2) by considering distinct removal rates for each species. By allowing a general class of interspecific density-dependent growth functions, we describe the multiplicity of all steady states and their necessary and sufficient conditions of existence and local stability according to the operating parameters $S_{in}$ and $D$.

In the second step, our principal objective in this work is to provide an important tool which is the operating diagram of the system. It provides a more global vision of the asymptotic behavior of solutions of the system according to the two operating parameters. On the other hand, the two-parameter diagram is also obtained by the numerical continuation method using MATCONT software [31]. This method allows us to detect other types of bifurcations according to two parameters that cannot be detected theoretically. Furthermore, we study the one-parameter bifurcation diagram with respect to the dilution rate $D$ showing the bifurcations of all steady states. Finally, through numerical simulations, we validate our theoretical results. All these results have not been done in the literature.

The paper is structured as follows. First, we present in Section 2 the mathematical model and some general hypotheses about the growth functions. In Section 3, we define some notations and functions used in this paper. Then, we determine the existence, uniqueness, and local stability conditions of all steady states according to the dilution rate and the input
concentration of the substrate. Next, in Section 4, we describe theoretically the operating diagram of the system and then numerically using MATCONT. The study of the one-parameter bifurcation with respect to the dilution rate and the numerical simulations are reported in Section 5. Some conclusions are drawn in Section 6. Appendix A is devoted to the proofs of results presented in Section 3. All the values of parameters used throughout this paper are provided in Appendix B.

**2. Mathematical model.** Taking into account the death rates of two species, the model is described by the following system of ordinary differential equations:

$$ (2.1) \qquad \left\{ \begin{aligned} \dot{S} &= D(S_{in} - S) - \mu_1(S, X_2)X_1 - \mu_2(S, X_1)X_2, \\ \dot{X}_1 &= (Y_1\mu_1(S, X_2) - D_1)X_1, \\ \dot{X}_2 &= (Y_2\mu_2(S, X_1) - D_2)X_2, \end{aligned} \right. $$

where $Y_1 \le 1$ and $Y_2 \le 1$ are yield coefficients. For $i = 1, 2, j = 1, 2, i \ne j$, $\mu_i$ is the density-dependent growth rate of the species $X_i$ that is assumed to be increasing in the variables $S$ and decreasing on $X_j$; $D_i$ is the removal rate of species $X_i$ and can be modeled as in [14] by

$$ D_i = \alpha_i D + a_i, $$

where the coefficient $\alpha_i \in [0, 1]$ is a parameter for decoupling hydraulic retention time and solid retention time; $a_i$ is the non-negative death rate of the species $X_i$.

To ease the mathematical analysis, we first can rescale system (2.1) using the following change of variables $x_i = X_i/Y_i$. Thus, system (2.1) becomes

$$ (2.2) \qquad \left\{ \begin{aligned} \dot{S} &= D(S_{in} - S) - f_1(S, x_2)x_1 - f_2(S, x_1)x_2, \\ \dot{x}_1 &= (f_1(S, x_2) - D_1)x_1, \\ \dot{x}_2 &= (f_2(S, x_1) - D_2)x_2, \end{aligned} \right. $$

where $f_1$ and $f_2$ are defined by

$$ f_1(S, x_2) = Y_1\mu_1(S, Y_2x_2) \quad \text{and} \quad f_2(S, x_1) = Y_2\mu_2(S, Y_1x_1). $$

In what follows, we consider model (2.2). For $i = 1, 2, j = 1, 2, i \ne j$, we make the following general hypotheses on the growth rate $f_i : \mathbb{R}_+^2 \to \mathbb{R}_+$ which is continuously differentiable.

**Hypothesis 2.1.** No growth can occur for species $x_j$ without the substrate $S : f_i(0, x_j) = 0$, for all $x_j \ge 0$.

**Hypothesis 2.2.** The growth rate of species $x_i$ is increasing by the substrate $S$ but is inhibited by the species $x_j : \frac{\partial f_i}{\partial S}(S, x_j) > 0$ and $\frac{\partial f_i}{\partial x_j}(S, x_j) < 0$ for all $S > 0$ and $x_j > 0$.

As model (2.2) represents a biological system, we require all solutions to be nonnegative and bounded for all time when given relevant initial conditions. Thus, solutions of system (2.2) verify the following result. Its proof uses classical tools, see for instance [36].

**Proposition 2.3.** Assume that Hypotheses 2.1 and 2.2 hold. For any nonnegative initial condition, the solution of system (2.2) exists for all $t \ge 0$, remains nonnegative and is bounded. Let $D_{\min} = \min(D, D_1, D_2)$. The set

$$ \Omega = \{(S, x_1, x_2) \in \mathbb{R}_+^3 : S + x_1 + x_2 \le DS_{in}/D_{\min}\} $$
is positively invariant and is a global attractor for the dynamics (2.2).

**3. Existence and stability of steady states.** Our aim in this section is to determine the existence and local stability conditions of all steady states of model (2.2). For convenience, we use the abbreviation LES for Locally Exponentially Stable steady states. In all figures, the straight line $\delta$ is defined by the equation $D_1x_1/D + D_2x_2/D = S_{in}$ and the curves $\gamma_1$ and $\gamma_2$ are in blue and red, respectively. In what follows, we need to define the set

$$ (3.1) \qquad M := \{(x_1, x_2) \in \mathbb{R}_+^2 : D_1x_1/D + D_2x_2/D \leq S_{in}\} $$

and the following quantities:

$$ (3.2) \qquad E = \frac{\partial f_1}{\partial S}, \quad F = \frac{\partial f_2}{\partial S}, \quad G = -\frac{\partial f_1}{\partial x_2}, \quad H = -\frac{\partial f_2}{\partial x_1}. $$

In (3.2), we have used the opposite sign of the partial derivatives $G = -\partial f_1/\partial x_2$ and $H = -\partial f_2/\partial x_1$, such that all constants involved in the computation become positive.

The steady states of (2.2) are the solutions of the set of equations

$$ (3.3) \qquad \begin{cases} 0 = D(S_{in} - S) - f_1(S, x_2)x_1 - f_2(S, x_1)x_2, \\ 0 = (f_1(S, x_2) - D_1)x_1, \\ 0 = (f_2(S, x_1) - D_2)x_2. \end{cases} $$

In the following, we will show that system (2.2) can have at most four steady states labeled below:

* $\mathcal{E}_0 = (S_{in}, 0, 0)$: the washout of two species ($x_1 = x_2 = 0$).

* $\mathcal{E}_1 = (\tilde{S}_1, \tilde{x}_1, 0)$: only the first species $x_1$ survives ($x_1 > 0, x_2 = 0$).

* $\mathcal{E}_2 = (\tilde{S}_2, 0, \tilde{x}_2)$: only the second species $x_2$ survives ($x_1 = 0, x_2 > 0$).

* $\mathcal{E}^* = (S^*, x_1^*, x_2^*)$: the coexistence of two species ($x_1 > 0, x_2 > 0$).

To determine the steady states of (2.2) under Hypotheses 2.1 and 2.2, we need to define some notations and functions in Table 1 and Lemma 3.1, respectively. The proof of Lemma 3.1 is similar to that of [35].

<table><caption>**Table 1**

Break-even concentrations and notations when $i = 1, 2$, $j = 1, 2$, $i \neq j$.</caption><tr><td>$\lambda_i(D)$</td><td>$S = \lambda_i(D)$ is the solution of $f_i(S, 0) = \alpha_i D + a_i$<br>It is defined for $D < (f_i(+\infty, 0) - a_i)/\alpha_i$</td></tr><tr><td>$\tilde{x}_i$</td><td>$\tilde{x}_i$ is the unique solution of $f_i(S_{in} - D_i x_i / D, 0) = D_i$<br>$\tilde{x}_i = D(S_{in} - \lambda_i(D))/D_i$</td></tr><tr><td>$\bar{x}_i$</td><td>$\bar{x}_i$ is the unique solution of $f_j(S_{in} - D_j x_i / D, x_i) = D_j$</td></tr></table>

**Lemma 3.1.** 1. Let $S_{in} > \lambda_1(D)$. The equation $f_1(S_{in} - D_1x_1/D - D_2x_2/D, x_2) = D_1$ defines a smooth decreasing function

$$ F_1 : \begin{cases} [0, \tilde{x}_1] & \to [0, \bar{x}_2] \\ x_1 & \mapsto F_1(x_1) = x_2, \end{cases} $$
such that $F_1(\tilde{x}_1) = 0$, $F_1(0) = \bar{x}_2$ where $\tilde{x}_1$ and $\bar{x}_2$ are defined in Table 1 with

$$ (3.4) \qquad F'_1(x_1) = -\frac{D_1 E}{D_2 E + DG} < 0, \quad \text{for all } x_1 \in [0, \tilde{x}_1]. $$

Furthermore, the graph $\gamma_1$ of $F_1$ lies in the interior of M for all $x_1 \in (0, \tilde{x}_1)$ (see Figure 1).

2. Suppose that $S_{in} > \lambda_2(D)$. The equation $f_2(S_{in} - D_1x_1/D - D_2x_2/D, x_1) = D_2$ defines a smooth decreasing function

$$ F_2 : [0, \tilde{x}_1] \to [0, \tilde{x}_2] \\ x_1 \mapsto F_2(x_1) = x_2, $$

such that $F_2(\bar{x}_1) = 0$, $F_2(0) = \tilde{x}_2$ where $\bar{x}_1$ and $\tilde{x}_2$ are defined in Table 1 with

$$ (3.5) \qquad F'_2(x_1) = -\frac{D_1 F + DH}{D_2 F} < 0, \quad \text{for all } x_1 \in [0, \bar{x}_1]. $$

In addition, the graph $\gamma_2$ of $F_2$ lies in the interior of M for all $x_1 \in (0, \bar{x}_1)$ (see Figure 1).

<figure><img src="image_1.png" /><figcaption>Figure 1. Relative positions of $\tilde{x}_i$ and $\bar{x}_i$, $i = 1, 2$: (a) Case 1 of (3.6) no intersection, (b) Case 2 of (3.6) a unique intersection, (c) Case 3 of (3.6) no intersection. The red [resp. blue] color is for the LES [resp. unstable] steady state.</figcaption></figure>

Note that $\tilde{x}_i$ and $\bar{x}_i$, $i = 1, 2$ defined in Table 1 represent the coordinates of the intersections of the curves $\gamma_1$ and $\gamma_2$ with the coordinates axes. According to their relative positions, there exist three cases that must be distinguished (see Figure 1):

$$ (3.6) \qquad \begin{array}{ll} \text{Case 1:} & \bar{x}_1 < \tilde{x}_1 \quad \text{and} \quad \tilde{x}_2 < \bar{x}_2. \\ \text{Case 2:} & \bar{x}_1 < \tilde{x}_1 \quad \text{and} \quad \bar{x}_2 < \tilde{x}_2. \\ \text{Case 3:} & \bar{x}_1 < \tilde{x}_1 \quad \text{and} \quad \bar{x}_2 < \tilde{x}_2. \end{array} $$

**Lemma 3.2.** When Hypotheses 2.1 and 2.2 hold and $S_{in} > \lambda_i(D)$, for $i = 1, 2$, we have the following equivalences:

$$ (3.7) \qquad \begin{array}{l} \bar{x}_1 < \tilde{x}_1 \Leftrightarrow f_2(\lambda_1(D), D(S_{in} - \lambda_1(D))/D_1) < D_2, \\ \bar{x}_2 < \tilde{x}_2 \Leftrightarrow f_1(\lambda_2(D), D(S_{in} - \lambda_2(D))/D_2) < D_1. \end{array} $$
The components of these steady states and the conditions of their existence according to the operating parameters are given by the following proposition:

**Proposition 3.3.** Assume that *Hypotheses 2.1 and 2.2* hold. The components of all steady states of system (2.2) as well as their necessary and sufficient conditions of existence are given in Table 2. All steady states are unique, when they exist.

<table><caption>**Table 2**

Components of all steady states of model (2.2) with their necessary and sufficient conditions of existence.</caption><thead><tr><th></th><th>Components</th><th>Existence condition</th></tr></thead><tbody><tr><td>ε<sub>0</sub></td><td>S = S<sub>in</sub>, x<sub>1</sub> = x<sub>2</sub> = 0</td><td>Always exists</td></tr><tr><td>ε<sub>1</sub></td><td>S&#x0303;<sub>1</sub> = λ<sub>1</sub>(D), x<sub>1</sub> = x&#x0303;<sub>1</sub>, x<sub>2</sub> = 0</td><td>S<sub>in</sub> &gt; λ<sub>1</sub>(D)</td></tr><tr><td>ε<sub>2</sub></td><td>S&#x0303;<sub>2</sub> = λ<sub>2</sub>(D), x<sub>1</sub> = 0, x<sub>2</sub> = x&#x0303;<sub>2</sub></td><td>S<sub>in</sub> &gt; λ<sub>2</sub>(D)</td></tr><tr><td>ε<sup>*</sup></td><td>S<sup>*</sup> = S<sub>in</sub> - D<sub>1</sub>x<sub>1</sub><sup>*</sup>/D - D<sub>2</sub>x<sub>2</sub><sup>*</sup>/D, (x<sub>1</sub><sup>*</sup>, x<sub>2</sub><sup>*</sup>) is the unique solution of x<sub>2</sub> = F<sub>1</sub>(x<sub>1</sub>) = F<sub>2</sub>(x<sub>1</sub>)</td><td>f<sub>1</sub>(λ<sub>2</sub>(D), D(S<sub>in</sub> - λ<sub>2</sub>(D))/D<sub>2</sub>) &lt; D<sub>1</sub> and f<sub>2</sub>(λ<sub>1</sub>(D), D(S<sub>in</sub> - λ<sub>1</sub>(D))/D<sub>1</sub>) &lt; D<sub>2</sub></td></tr></tbody></table>

As an immediate corollary of Lemma 3.2 and Proposition 3.3, we have the following result establishing the existence and the uniqueness of the coexistence steady state $\mathcal{E}^*$.

**Corollary 3.4.** Assume that *Hypotheses 2.1 and 2.2* hold and $S_{in} > \lambda_i(D)$ for $i = 1, 2$. A coexistence steady state $\mathcal{E}^*$ exists if and only if Case 2 of (3.6) holds. It is unique when it exists.

**Remark 3.5.** If a positive intersection point $(x_1, x_2)$ of $\gamma_1$ and $\gamma_2$ exists, then the tangent of $\gamma_1$ at point $(x_1, x_2)$ is always above the tangent of $\gamma_2$ at point $(x_1, x_2)$.

The stability conditions of all steady states of (2.2) are provided by the following result:

**Proposition 3.6.** Under *Hypotheses 2.1 and 2.2*, the necessary and sufficient conditions for the local asymptotic stability of the steady states $\mathcal{E}_0, \mathcal{E}_1, \mathcal{E}_2$ and $\mathcal{E}^*$ are given in Table 3.

<table><caption>**Table 3**

Necessary and sufficient conditions of local stability of all steady states of (2.2) according to the operating parameters $S_{in}$ and $D$.</caption><thead><tr><th></th><th>Stability condition</th></tr></thead><tbody><tr><td>ε<sub>0</sub></td><td>S<sub>in</sub> &lt; min(λ<sub>1</sub>(D), λ<sub>2</sub>(D))</td></tr><tr><td>ε<sub>1</sub></td><td>f<sub>2</sub>(λ<sub>1</sub>(D), D(S<sub>in</sub> - λ<sub>1</sub>(D))/D<sub>1</sub>) &lt; D<sub>2</sub></td></tr><tr><td>ε<sub>2</sub></td><td>f<sub>1</sub>(λ<sub>2</sub>(D), D(S<sub>in</sub> - λ<sub>2</sub>(D))/D<sub>2</sub>) &lt; D<sub>1</sub></td></tr><tr><td>ε<sup>*</sup></td><td>unstable whenever it exists</td></tr></tbody></table>

From Propositions 3.3 and 3.6, we obtain the following corollary.

**Corollary 3.7.** • $\mathcal{E}_0$ is LES if and only if $\mathcal{E}_1$ and $\mathcal{E}_2$ do not exist.

• A positive steady state $\mathcal{E}^*$ exists if and only if $\mathcal{E}_1$ and $\mathcal{E}_2$ are LES.

**4. Operating diagram.** The operating diagram is a tool to describe how a system behaves when all biological parameters are fixed and the operating (control) parameters are varied, as
they are the most easily manipulated parameters in a chemostat. This operating diagram is very useful to understand the process from the mathematical and biological point of view. It is often constructed in the mathematical literature (see for instance [1, 2, 7, 8, 14, 16, 21, 34, 36, 39, 40, 41]) and the biological literature (see [23, 42, 47, 48]). In the existing literature, this operating diagram have been studied with three different methods. For more details on these various methods, the reader is referred to [36] and the references therein.

The aim of this section is to study theoretically and numerically the operating diagram of (2.2) under variation of the operating parameters which are the concentration of the substrate in the feed bottle $S_{in}$ and the dilution rate $D$. From Tables 2 and 3, we define in Table 4 the set $\Upsilon = \{\Upsilon_1, \Upsilon_2, \Upsilon_1^c, \Upsilon_2^c\}$. This set represents the boundaries of the different regions in $(S_{in}, D)$-plane.

<caption>**Table 4**</caption>

The set $\Upsilon$ of boundaries of the regions in the operating diagram and their corresponding colors used in Figures 2 and 3.

<table><tr><td>$\Upsilon_1 = \{(S_{in}, D) : S_{in} = \lambda_1(D)\}$</td><td>Black</td></tr><tr><td>$\Upsilon_2 = \{(S_{in}, D) : S_{in} = \lambda_2(D)\}$</td><td>Blue</td></tr><tr><td>$\Upsilon_1^c = \{(S_{in}, D) : \tilde{x}_1(S_{in}, D) = \bar{x}_1(S_{in}, D)\}$</td><td>Red</td></tr><tr><td>$\Upsilon_2^c = \{(S_{in}, D) : \tilde{x}_2(S_{in}, D) = \bar{x}_2(S_{in}, D)\}$</td><td>Magenta</td></tr></table>

<figure><img src="image_2.png" /><figcaption>**Figure 2.** MAPLE: (a) operating diagram of (2.2) in the $(S_{in}, D)$-plane. (b) Magnification on the regions $\mathcal{J}_2$ and $\mathcal{J}_4$ when $(S_{in}, D) \in [0, 1] \times [0, 2]$.</figcaption></figure>

In order to construct the operating diagram, we consider the following specific growth rates satisfying Hypotheses 2.1 and 2.2:

$$ (4.1) \qquad f_1(S, x_2) = \frac{m_1 S}{K_1 + S + \beta_1 x_2}, \quad f_2(S, x_1) = \frac{m_2 S}{K_2 + S + \beta_2 x_1}, $$

where $m_1$ and $m_2$ denote the maximum growth rates; $K_1$ and $K_2$ denote the Michaelis-Menten constants; $\beta_1$ [resp. $\beta_2$] represents the inhibition factor due to $x_2$ [resp. $x_1$] for the growth of the species $x_1$ [resp. $x_2$]. Except for the operating parameters $S_{in}$ and $D$, all these parameters cannot be easily manipulated by the biologist. Thus, we fix it as provided in Table 11. Note
that the construction of the operating diagram is similar for any other specific growth rate satisfying Hypotheses 2.1 and 2.2.

From Tables 2 and 3 establishing the existence and local stability conditions of all steady states, we can state the following result determining theoretically the operating diagram where the various functions and the corresponding curves are defined in Tables 1 and 4.

**Proposition 4.1.** For the specific growth rates $f_1$ and $f_2$ defined in (4.1) and the biological parameter values provided in Table 11, the existence and the local stability of all steady states of system (2.2) in the six regions $\mathcal{J}_k$, $k = 0, \dots, 5$ of the operating diagram shown in Figure 2 are described in Table 5.

<caption>**Table 5**</caption>

Definitions of the regions of the operating diagram in Figures 2 and 3 when $i = 1, 2$. The letter S [resp. U)] means stable [resp. unstable] steady state. The absence of a letter means that the corresponding steady state does not exist.

<table><thead><tr><th>Region</th><th>Condition 1</th><th>Condition 2</th><th>Color</th><th>ε<sub>0</sub></th><th>ε<sub>1</sub></th><th>ε<sub>2</sub></th><th>ε*</th></tr></thead><tbody><tr><td>J<sub>0</sub></td><td>S<sub>in</sub> > min(λ<sub>1</sub>(D), λ<sub>2</sub>(D))</td><td>x&#772;<sub>1</sub>(S<sub>in</sub>, D) < x&#772;<sub>1</sub>(S<sub>in</sub>, D)</td><td>White</td><td>S</td><td></td><td></td><td></td></tr><tr><td>J<sub>1</sub></td><td>λ<sub>1</sub>(D) < S<sub>in</sub> < λ<sub>2</sub>(D)</td><td>x&#772;<sub>2</sub>(S<sub>in</sub>, D) < x&#772;<sub>2</sub>(S<sub>in</sub>, D)</td><td>Green</td><td>U</td><td>S</td><td></td><td></td></tr><tr><td>J<sub>2</sub></td><td>λ<sub>2</sub>(D) < S<sub>in</sub> < λ<sub>1</sub>(D)</td><td>x&#772;<sub>1</sub>(S<sub>in</sub>, D) < x&#772;<sub>1</sub>(S<sub>in</sub>, D)</td><td>Pink</td><td>U</td><td></td><td>S</td><td></td></tr><tr><td>J<sub>3</sub></td><td>S<sub>in</sub> > max(λ<sub>1</sub>(D), λ<sub>2</sub>(D))</td><td>x&#772;<sub>i</sub>(S<sub>in</sub>, D) < x&#772;<sub>i</sub>(S<sub>in</sub>, D)</td><td>Yellow</td><td>U</td><td>S</td><td>S</td><td>U</td></tr><tr><td>J<sub>4</sub></td><td>S<sub>in</sub> > max(λ<sub>1</sub>(D), λ<sub>2</sub>(D))</td><td>x&#772;<sub>2</sub>(S<sub>in</sub>, D) < x&#772;<sub>2</sub>(S<sub>in</sub>, D),<br>x&#772;<sub>1</sub>(S<sub>in</sub>, D) > x&#772;<sub>1</sub>(S<sub>in</sub>, D)</td><td>Pink</td><td>U</td><td>U</td><td>S</td><td></td></tr><tr><td>J<sub>5</sub></td><td>S<sub>in</sub> > max(λ<sub>1</sub>(D), λ<sub>2</sub>(D))</td><td>x&#772;<sub>2</sub>(S<sub>in</sub>, D) > x&#772;<sub>2</sub>(S<sub>in</sub>, D),<br>x&#772;<sub>1</sub>(S<sub>in</sub>, D) < x&#772;<sub>1</sub>(S<sub>in</sub>, D)</td><td>Green</td><td>U</td><td>S</td><td>U</td><td></td></tr></tbody></table>

As shown in Figure 2, the $(S_{in}, D)$-plane is divided into six regions. The region $J_0$ corresponds to the washout of two species. The regions $J_1$ and $J_5$ [resp. $J_2$ and $J_4$] correspond to the competitive exclusion of the second [resp. first] species. The region $J_3$ corresponds to the bi-stability with convergence to one of the two boundary steady states.

In the following, we used MATCONT [31] to numerically analyze the two-parameter diagram of system (2.2) for the set of the biological parameter values of Figure 2, which is provided in Table 11. Indeed, MATCONT [31] is a MATLAB numerical continuation package for the bifurcation study of continuous and discrete parameterized systems of ODEs. Thanks to the numerous features and functions, it allows us to compute curves of steady states and limit cycles as well as their local asymptotic behaviors and bifurcations as Branch Points (BP) or transcritical bifurcations, Limit Points (LP) or saddle-node (or fold) bifurcations, Hopf points (H), Limit Point of Cycles (LPC), Cusp bifurcations (CP), Bogdanov-Takens bifurcations, etc. All these curves are determined under variation of one or more system parameters by numerical continuation. For more details, the reader is referred to [10] and the references therein.

Figure 3 illustrates the operating diagram obtained numerically using MATCONT which is identical to one determined theoretically in Figure 2. However, MATCONT can detect other two-parameter bifurcations as Zero-Hopf Point (ZH) and Bogdanov-Takens Point (BT). These types of bifurcation are not detected in the theoretical study of the operating diagram. In Figure 3, the intersection point of all curves in the set $\mathcal{Y}$ corresponds to a two-parameter
<figure><img src="image_3.png" /><figcaption>**Figure 3.** MATCONT: (a) Operating diagram of (2.2). (b) Magnification of the two-parameter bifurcations at points ZH and BT when $(S_{in}, D) \in [0, 1] \times [0, 2]$.</figcaption></figure>

bifurcation of type Zero-Hopf Point (ZH). Moreover, MATCONT was able to detect a second Zero-Hopf Point (ZH) and a new two-parameter bifurcation of type Bogdanov-Takens (BT) with normal form coefficient $(a, b) = (2.56 \ 10^{-6}, 0.0909)$. Table 6 summarizes the critical operating parameters and the state at BT and ZH bifurcations.

<caption>**Table 6**</caption>

Operating parameters and state for BT and ZH at the bifurcation points in Figure 3. The abbreviation BT [resp. ZH] means a Bogdanov-Takens [resp. Zero-Hopf] bifurcation.

<table><thead><tr><th>Bifurcation</th><th>Parameter (S<sub>in</sub>, D)</th><th>State (S, x<sub>1</sub>, x<sub>2</sub>)</th></tr></thead><tbody><tr><td>ZH</td><td>(0.418289, 0.361063)</td><td>(0.418289, 0, 0)</td></tr><tr><td>BT</td><td>(0.2, 2.9 10<sup>-6</sup>)</td><td>(0.2, 0, 0)</td></tr><tr><td>ZH</td><td>(0.387348, 0.10468)</td><td>(0.387348, 0, 0)</td></tr></tbody></table>

**5. Bifurcation analysis and numerical simulations.** In this section, we will analyze the one-parameter bifurcation diagram of system (2.2). This diagram shows the various types of bifurcations and the asymptotic behavior of the process by varying the dilution rate $D$ as the bifurcating parameter where the input concentration $S_{in}$ is fixed. In a similar way, we can study the one-parameter bifurcation diagram where $S_{in}$ is the bifurcation parameter. Note that we have chosen $S_{in} = 1$ to maximize the passage number through the various regions of the operating diagram in Figures 2 and 3. The following result establishes the nature of all bifurcations that occur by passing from one region to another in the operating diagram in Figures 2 and 3.

**Proposition 5.1.** For the specific growth rates $f_1$ and $f_2$ (4.1) and the set of the biological parameter values provided in Table 11, the bifurcation analysis of all steady states of system (2.2) by crossing the various curves in the set $\Upsilon$ defined in Table 4, is summarized in Table 7.

From Table 7 and as we shall see later, the passage from one region to another through the curves of the set $\Upsilon$ always corresponds to a transcritical bifurcation where steady states coalesce and change stability. Figure 4 represents the one-parameter bifurcation diagram of
<caption>**Table 7**</caption>

Nature of bifurcations of steady states of (2.2) by crossing to the surfaces of Υ. The letter TB means a transcritical bifurcation.

<table><thead><tr><th>Transition</th><th>Υ</th><th>Bifurcation</th><th>Steady state</th></tr></thead><tbody><tr><td>J<sub>0</sub> to J<sub>1</sub> or J<sub>2</sub> to J<sub>4</sub></td><td>Υ<sub>1</sub></td><td>TB</td><td>E<sub>0</sub> = E<sub>1</sub></td></tr><tr><td>J<sub>0</sub> to J<sub>2</sub> or J<sub>1</sub> to J<sub>5</sub></td><td>Υ<sub>2</sub></td><td>TB</td><td>E<sub>0</sub> = E<sub>2</sub></td></tr><tr><td>J<sub>3</sub> to J<sub>4</sub></td><td>Υ<sub>1</sub><sup>c</sup></td><td>TB</td><td>E<sub>1</sub> = E<sub>*</sub></td></tr><tr><td>J<sub>5</sub> to J<sub>3</sub></td><td>Υ<sub>2</sub><sup>c</sup></td><td>TB</td><td>E<sub>2</sub> = E<sub>*</sub></td></tr></tbody></table>

system (2.2) where the ω-limit set is projected to the S coordinate, showing all the bifurcations that occur when we vary D. Similarly, we can obtain the one-parameter bifurcation diagram with $x_1$ or $x_2$ on the y-axis. This diagram corresponds to a vertical line of equation $S_{in} = 1$ in the operating diagram in Figures 2 and 3. A magnification of the one-parameter bifurcation diagram is illustrated in Figure 4(b) showing the transcritical bifurcations occur at $\sigma_3$ and $\sigma_4$ which describe the appearance and disappearance of the unique positive steady state $\mathcal{E}^*$.

<figure><img src="image_4.png" /><figcaption>**Figure 4.** (a) One-parameter bifurcation diagram of (2.2) in the variable S with D as the bifurcation parameter and $S_{in} = 1$. (b) Magnification when $D \in [\sigma_3, \sigma_4]$. Blue dashed curves correspond to unstable steady states and red solid curves to LES steady states. The green solid diamonds represent the transcritical bifurcations.</figcaption></figure>

First, by decreasing D, the one-parameter bifurcation diagram in Figure 4 illustrates the transcritical bifurcation occurring at $D = \sigma_1 = 4$ between $\mathcal{E}_0$ and $\mathcal{E}_1$. Indeed, the washout $\mathcal{E}_0$ loses its stability while the steady state of exclusion of the second species $\mathcal{E}_1$ emerges LES. Decreasing D further, the steady state of exclusion of the first species $\mathcal{E}_2$ appears unstable at $D = \sigma_2 \approx 1.0666$ via a second transcritical bifurcation with the steady state $\mathcal{E}_0$. Next, a third transcritical bifurcation occurs at $D = \sigma_3 \approx 0.6447$ that gives birth to an unstable positive steady state $\mathcal{E}^*$ where $\mathcal{E}_2$ changes its stability and becomes LES. Finally, decreasing D further, once again a fourth transcritical bifurcation occurs at $D = \sigma_4 \approx 0.352$ between $\mathcal{E}^*$ and $\mathcal{E}_1$ where the steady state $\mathcal{E}^*$ disappears while $\mathcal{E}_1$ remains unstable. In the following proposition, we summarize the study of the one-parameter bifurcation diagram in D corresponding to the operating diagram in Figures 2 and 3 when the input concentration $S_{in}$ is fixed.
**Proposition 5.2.** Let $S_{in} = 1$. Using the biological parameter values provided in Table 11 and the specific growth rates (4.1), the existence and stability of steady states of (2.2) according to the control parameter D are given in Table 9 where the bifurcation values $\sigma_i$, $i = 1, \dots, 4$ and the corresponding nature of the bifurcations are defined in Table 8.

<table><caption>**Table 8**

Critical values $\sigma_i$, $i = 1, \dots, 4$ of D and the corresponding nature of bifurcations when $S_{in} = 1$.</caption><thead><tr><th>Definition</th><th>Value</th><th>Bifurcation</th></tr></thead><tbody><tr><td>$\sigma_1 = f_1(S_{in}, 0)$</td><td>4</td><td>TB</td></tr><tr><td>$\sigma_2 = f_2(S_{in}, 0)$</td><td>1.0666</td><td>TB</td></tr><tr><td>$\sigma_3$ is the solution of equation $\tilde{x}_2(S_{in}, D) = \bar{x}_2(S_{in}, D)$</td><td>0.6447</td><td>TB</td></tr><tr><td>$\sigma_4$ is the solution of equation $\tilde{x}_1(S_{in}, D) = \bar{x}_1(S_{in}, D)$</td><td>0.352</td><td>TB</td></tr></tbody></table>

<table><caption>**Table 9**

Existence and stability of all steady states according to D where $\sigma_i$, $i = 1, \dots, 4$ are defined in Table 8.</caption><thead><tr><th>Interval of D</th><th>$\mathcal{E}_0$</th><th>$\mathcal{E}_1$</th><th>$\mathcal{E}_2$</th><th>$\mathcal{E}^*$</th></tr></thead><tbody><tr><td>$(0, \sigma_4)$</td><td>U</td><td>U</td><td>S</td><td></td></tr><tr><td>$(\sigma_4, \sigma_3)$</td><td>U</td><td>S</td><td>S</td><td>U</td></tr><tr><td>$(\sigma_3, \sigma_2)$</td><td>U</td><td>S</td><td>U</td><td></td></tr><tr><td>$(\sigma_2, \sigma_1)$</td><td>U</td><td>S</td><td></td><td></td></tr><tr><td>$(\sigma_1, +\infty)$</td><td>S</td><td></td><td></td><td></td></tr></tbody></table>

In what follows, we illustrate numerically different cases of stability of steady states using SCILAB [43]. Using the same set of parameter values in Figures 2 and 3, the numerical simulations in Figure 5 illustrate the trajectories in the three-dimensional phase space $(S, x_1, x_2)$. Figure 5(a) shows that the solutions of system (2.2) converge to $\mathcal{E}_1$ for several positive initial conditions when $(S_{in}, D) = (1, 0.7) \in \mathcal{J}_5$. In this case, there exist three steady states given by

$$ \mathcal{E}_0 = (1, 0, 0), \quad \mathcal{E}_1 \approx (0.4607, 0.4015, 0), \quad \mathcal{E}_2 \approx (0.6666, 0, 0.4242), $$

where $\mathcal{E}_1$ is LES while $\mathcal{E}_0$ and $\mathcal{E}_2$ are unstable. The numerical simulations were done for several positive initial conditions permitting us to conjecture the convergence to $\mathcal{E}_1$ for any positive initial condition. Next, we choose $(S_{in}, D) = (1, 0.2) \in \mathcal{J}_4$. System (2.2) admits the following three steady states

$$ \mathcal{E}_0 = (1, 0, 0), \quad \mathcal{E}_1 \approx (0.3987, 0.1431, 0), \quad \mathcal{E}_2 \approx (0.3157, 0, 0.4561), $$

where $\mathcal{E}_2$ is LES while $\mathcal{E}_0$ and $\mathcal{E}_1$ are unstable. Figure 5(b) illustrates the convergence towards $\mathcal{E}_2$ for several positive initial conditions. Finally, when $(S_{in}, D) = (1, 0.5)$ is chosen in the bistability region $\mathcal{J}_3$, the washout steady state is given by $\mathcal{E}_0 = (1, 0, 0)$ and the other steady states are given by

$$ \mathcal{E}_1 \approx (0.4354, 0.3136, 0), \quad \mathcal{E}_2 \approx (0.5142, 0, 0.5396), \quad \mathcal{E}^* = (0.5181, 0.1491, 0.2371), $$
<figure><img src="image_5.png" /><figcaption>**Figure 5.** SCILAB: the three-dimensional space $(S, x_1, x_2)$ of (2.2) according to regions in the operating diagram of Figures 2 and 3 when $S_{in} = 1$. (a) Convergence to $\mathcal{E}_1$ for $(S_{in}, D) \in \mathcal{J}_5$. (b) Convergence to $\mathcal{E}_2$ for $(S_{in}, D) \in \mathcal{J}_4$. (c) Bistability of $\mathcal{E}_1$ and $\mathcal{E}_2$ when $(S_{in}, D) \in \mathcal{J}_3$.</figcaption></figure>

where $\mathcal{E}_0$ and $\mathcal{E}^*$ are unstable while $\mathcal{E}_1$ and $\mathcal{E}_2$ are stable. In fact, the system exhibits bistability with two basins of attraction, which are separated by the stable manifold of saddle point $\mathcal{E}^*$. One basin attracts the solutions to the steady state $\mathcal{E}_1$ and the other to the steady state $\mathcal{E}_2$, as shown in Figure 5(c). The three cases of (3.6) and their correspondences in the three cases considered in Figure 5 (with the stable steady states and the regions of the operating diagram) are summarized in Table 10.

<table><caption>**Table 10**
The three cases of (3.6) and their correspondences illustrated in Figure 5.</caption><thead><tr><th>Case</th><th>Region</th><th>Stable steady state</th><th>Figure</th></tr></thead><tbody><tr><td>Case 1</td><td>J<sub>5</sub></td><td>E<sub>1</sub></td><td>Figure 5(a)</td></tr><tr><td>Case 2</td><td>J<sub>3</sub></td><td>E<sub>1</sub> and E<sub>2</sub></td><td>Figure 5(c)</td></tr><tr><td>Case 3</td><td>J<sub>4</sub></td><td>E<sub>2</sub></td><td>Figure 5(b)</td></tr></tbody></table>

**6. Conclusion.** In this work, we have extended the intra-specific density-dependent model of [11] where two species compete for a single resource by allowing distinct removal rates. Our
study considers the effect of mortality terms, which were neglected in the theoretical analysis of [11]. The three-dimensional system cannot be reduced to a two-dimensional one.

On the other hand, this model was studied in Fekih-Salem et al. [16] in the case of intra- and interspecific density-dependent growth rates. However, the study of the particular case with only interspecific interference between the two populations of microorganisms was not investigated in Fekih-Salem et al. [16]. Our study has allowed us to understand whether intraspecific interference is necessary to guarantee persistence between microbial species.

Under general interspecific density-dependent growth functions, we gave an exhaustive study of the behavior of system (2.2). More precisely, the necessary and sufficient conditions of existence and stability of all steady states are determined according to the control parameters $S_{in}$ and $D$. We have shown that system (2.2) admits at most four types of steady states: the washout steady state $\mathcal{E}_0$ which always exists, two boundary steady states $\mathcal{E}_1$ and $\mathcal{E}_2$ in which one of species degenerates and the other persists and a coexistence steady state of two species $\mathcal{E}^*$. We have proved that whenever $\mathcal{E}^*$ exists, it is unique and unstable. In fact, the existence and the uniqueness of the positive steady state $\mathcal{E}^*$ depends on the relative position of the values $\tilde{x}_i$ and $\bar{x}_i$, $i = 1, 2$ (as shown in Figure 1). In addition, a positive steady state exists if and only if both boundary steady states are LES. In this case, the model exhibits the bi-stability.

Our result on the instability of the positive steady state in this paper can be deduced from those in [16]. More precisely, [16, condition (28)] holds in the particular case without intraspecific interference and [16, Theorem 1] implies in this case that a positive steady state, if it exists, is always unique and unstable.

Using our mathematical analysis of model (2.2), we have determined, first, theoretically the operating diagram to provide a complete description of the asymptotic behavior of the process with respect to the control parameters $S_{in}$ and $D$. Our description of the operating diagram in the $(S_{in}, D)$-plane using the specific growth functions (4.1) shows that this diagram is divided into six regions: $\mathcal{J}_0$ corresponds to the washout of two species; the regions $\mathcal{J}_1$ and $\mathcal{J}_5$ [resp. $\mathcal{J}_2$ and $\mathcal{J}_4$] correspond to the competitive exclusion of the second [resp. first] species; the bi-stability region $\mathcal{J}_3$ of the boundary steady states with the existence of a unique unstable positive steady state.

On the other hand, we have established numerically the operating diagram of (2.2) using the software MatCont [31]. We found the same operating diagram obtained theoretically from the existence and stability conditions. However, this numerical method was able to detect new types of two-parameter bifurcations: BT, which corresponds to a Bogdanov-Takens bifurcation, and ZH, which corresponds to a Zero-Hopf bifurcation (see Figure 3).

When only the dilution rate $D$ varies, the study of the one-parameter bifurcation diagram has illustrated the nature of bifurcations of all steady states by passing through the various regions of the operating diagram in Figures 2 and 3. It was shown that all steady states can appear or disappear only through transcritical bifurcations. The numerical simulations have shown in three-dimensional phase space $(S, x_1, x_2)$ that the system can exhibit bi-stability with convergence either to the exclusion of the first species or to the exclusion of the second species.

Finally, our mathematical analysis proves that the outcome of competition satisfies the CEP which predicts that only one species can exist. Thus, our study allows us to say that
the addition of mortality terms of the species in the interspecific density-dependent model with a mutual inhibitory relationship, could not lead to the coexistence around a stable positive steady state or stable limit cycle. Comparing our results with those in [16], we can conclude that interspecific interference is not sufficient to guarantee coexistence between species, contrary to intraspecific interference.

## Appendix A. Proofs.

**Proof of Lemma 3.2.** Under Hypothesis 2.2, it is easy to see that the function $x_1 \mapsto f_2(S_{in} - D_1x_1/D, x_1) - D_2$ is decreasing. From the definition of $\bar{x}_1$ in Table 1, we have

$$ \bar{x}_1 < \tilde{x}_1 \Leftrightarrow f_2(S_{in} - D_1\tilde{x}_1/D, \tilde{x}_1) - D_2 < f_2(S_{in} - D_1\bar{x}_1/D, \bar{x}_1) - D_2 = 0. $$

From the definition of $\tilde{x}_1$ in Table 1, we obtain the first equivalence in (3.7). Similar arguments apply to the second equivalence in (3.7). ■

**Proof of Proposition 3.3.** • When $x_1 = x_2 = 0$, from the first equation of (3.3), one has $S = S_{in}$. So, $\mathcal{E}_0 = (S_{in}, 0, 0)$, it always exists.

• For the boundary steady state $\mathcal{E}_i$, $x_i > 0$ and $x_j = 0$, $i = 1, 2$, $j = 1, 2$, $j \neq i$. From system (3.3), the components $\tilde{S}_i$ and $\tilde{x}_i$ are the solutions of

$$ (A.1) \qquad D(S_{in} - \tilde{S}_i) = D_i \tilde{x}_i $$

$$ (A.2) \qquad f_i(\tilde{S}_i, 0) = D_i. $$

Since the function $S \mapsto f_i(S, 0)$ is increasing (under Hypothesis 2.2), it follows from equation (A.2) and Table 1 that

$$ \tilde{S}_i = \lambda_i(D). $$

Replacing $\tilde{S}_i$ by this expression in (A.1), we can see that

$$ \tilde{x}_i = D(S_{in} - \lambda_i(D))/D_i, $$

which is positive if and only if $S_{in} > \lambda_i(D)$. If $\mathcal{E}_i$ exists, then it is unique.

• For the positive steady state $\mathcal{E}^*$, $x_1 > 0$ and $x_2 > 0$. From system (3.3), the components $S^*$, $x_1^*$ and $x_2^*$ are the solutions of the set of equations

$$ (A.3) \qquad D(S_{in} - S) = D_1x_1 + D_2x_2, $$

$$ (A.4) \qquad f_1(S, x_2) = D_1, $$

$$ (A.5) \qquad f_2(S, x_1) = D_2. $$

Hence, from (A.3), we obtain

$$ (A.6) \qquad S^* = S_{in} - D_1x_1^*/D - D_2x_2^*/D. $$

Replacing expression (A.6) in (A.4)-(A.5), one has $(x_1 = x_1^*, x_2 = x_2^*)$ is a solution of

$$ (A.7) \qquad \left\{ \begin{array}{l} f_1(S_{in} - D_1x_1/D - D_2x_2/D, x_2) = D_1, \\ f_2(S_{in} - D_1x_1/D - D_2x_2/D, x_1) = D_2. \end{array} \right. $$
From (A.6), one see that $S^*$ is positive if and only if $D_1x_1^*/D + D_2x_2^*/D < S_{in}$. Thus, (A.7) has a positive solution in the interior of the set $M$ defined by (3.1), that is, the curves $\gamma_1$ of $F_1$ and $\gamma_2$ of $F_2$, respectively, has a positive intersection in the interior of $M$ (from Lemma 3.1). More precisely, $(x_1 = x_1^*, x_2 = x_2^*)$ is a solution of equation

$$x_2 = F_1(x_1) \quad \text{and} \quad x_2 = F_2(x_1).$$

Let $F(x_1)$ be the function defined by

$$(A.8) \qquad F(x_1) := F_1(x_1) - F_2(x_1).$$

Using the derivatives (3.4) of the function $F_1$ and (3.5) of $F_2$, respectively, a straightforward calculation shows that

$$F_1'(x_1) - F_2'(x_1) = \frac{D(DGH + D_1FG + D_2EH)}{D_2F(D_2E + DG)} > 0, \quad \text{for all} \quad x_1 < \min(\tilde{x}_1, \bar{x}_1).$$

Therefore, the function $F$ is increasing from $F(0) = \bar{x}_2 - \tilde{x}_2$ to $F(\bar{x}_1) = F_1(\bar{x}_1) - F_1(\tilde{x}_1)$ since $F_2(\bar{x}_1) = 0 = F_1(\tilde{x}_1)$. Hence, $F(\bar{x}_1) > 0$ if and only if $\bar{x}_1 < \tilde{x}_1$. Consequently, there exists an $x_1 = x_1^* \in (0, \bar{x}_1)$ satisfying (A.8) if and only if

$$\bar{x}_1 < \tilde{x}_1 \quad \text{and} \quad \bar{x}_2 < \tilde{x}_2.$$

If $\tilde{x}_1$ exists, then it is unique. This completes the proof.

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

**Proof of Proposition 3.6.** Here, we study the stability conditions of all steady states of (2.2). The Jacobian matrix of (2.2) at $(S, x_1, x_2)$ takes the form

$$J = \begin{pmatrix} -D - x_1E - x_2F & -f_1(S, x_2) + x_2H & x_1G - f_2(S, x_1) \\ x_1E & f_1(S, x_2) - D_1 & -x_1G \\ x_2F & -x_2H & f_2(S, x_1) - D_2 \end{pmatrix}$$

where $E, F, G$ and $H$ are given by (3.2).

* At $\mathcal{E}_0 = (S_{in}, 0, 0)$, the characteristic polynomial is

$$P(\lambda) = (-D - \lambda)(f_1(S_{in}, 0) - D_1 - \lambda)(f_2(S_{in}, 0) - D_2 - \lambda).$$

Therefore, $\mathcal{E}_0$ is LES if and only if $f_i(S_{in}, 0) < D_i$, for $i = 1, 2$. Using the definition of $\lambda_i(D)$ in Table 1, we conclude that $\mathcal{E}_0$ is LES if and only if $S_{in} < \lambda_i(D)$, for $i = 1, 2$.

* At $\mathcal{E}_1 = (\lambda_1(D), D(S_{in} - \lambda_1(D))/D_1, 0)$, the characteristic polynomial is

$$P(\lambda) = (f_2(\lambda_1(D), D(S_{in} - \lambda_1(D))/D_1) - D_2 - \lambda)(\lambda^2 + \beta_1\lambda + \beta_2),$$

where $\beta_1 = D(1 + E(S_{in} - \lambda_1(D))/D_1)$ and $\beta_2 = DE(S_{in} - \lambda_1(D))$. Since $\beta_1 > 0$ and $\beta_2 > 0$, the real parts of the roots of the quadratic factor are negative. Therefore, $\mathcal{E}_1$ is LES if and only if $f_2(\lambda_1(D), D(S_{in} - \lambda_1(D))/D_1) < D_2$.
* Similarly, at $\mathcal{E}_2 = (\lambda_2(D), 0, D(S_{in} - \lambda_2(D))/D_2)$, the characteristic polynomial is

$$P(\lambda) = (f_1(\lambda_2(D), D(S_{in} - \lambda_2(D))/D_2) - D_1 - \lambda)(\lambda^2 + \delta_1\lambda + \delta_2),$$

where $\delta_1 = D(1 + F(S_{in} - \lambda_2(D))/D_2)$ and $\delta_2 = DF(S_{in} - \lambda_2(D))$. Since $\delta_1 > 0$ and $\delta_2 > 0$, the real parts of the roots of the quadratic factor are negative. Therefore, $\mathcal{E}_2$ is LES if and only if $f_1(\lambda_2(D), D(S_{in} - \lambda_2(D))/D_2) < D_1$.

* The characteristic polynomial of the Jacobian matrix $J$ evaluated at $\mathcal{E}^*$ is given by:

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

where

$$ \begin{aligned} c_1 &= D + Ex_1^* + Fx_2^*, & c_2 &= D_1Ex_1^* + D_2Fx_2^* - (GH + FG + EH)x_1^*x_2^*, \\ c_3 &= -(DGH + D_1FG + D_2EH)x_1^*x_2^*. \end{aligned} $$

Hence, the coefficient $c_3 < 0$. According to the Routh-Hurwitz criterion, $\mathcal{E}^*$ is always unstable, whenever it exists. ■

**Appendix B. Parameters used in numerical simulations.** All the values of the parameters used in the numerical simulations are provided in Table 11. For the numerical simulations, we have used MAPLE [30] to plot Figures 1, 2, and 4, MATCONT [31] for Figure 3 and SCILAB [43] for Figure 5.

<table><caption>**Table 11**

Parameter values used for model (2.2) when the growth rates $f_1$ and $f_2$ are given by (4.1).</caption><thead><tr><th>Parameter</th><th>m<sub>1</sub></th><th>K<sub>1</sub></th><th>β<sub>1</sub></th><th>m<sub>2</sub></th><th>K<sub>2</sub></th><th>β<sub>2</sub></th><th>α<sub>1</sub></th><th>α<sub>2</sub></th><th>a<sub>1</sub></th><th>a<sub>2</sub></th><th>S<sub>in</sub></th><th>D</th></tr></thead><tbody><tr><td>Figs. 1(a), 5(a)</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td>1</td><td>0.7</td></tr><tr><td>Figs. 1(b), 5(c)</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td>1</td><td>0.5</td></tr><tr><td>Figs. 1(c), 5(b)</td><td>4</td><td>1.5</td><td>1.2</td><td>2.2</td><td>2</td><td>0.1</td><td>0.2</td><td>0.5</td><td>0.8</td><td>0.2</td><td>1</td><td>0.2</td></tr><tr><td>Figs. 2, 3</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td>var</td><td>var</td></tr><tr><td>Fig. 4</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td>1</td><td>var</td></tr></tbody></table>

REFERENCES

[1] 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).

[2] 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).

[3] B. Benyahia, T. Sari, B. Cherki and J. Harmand, Bifurcation and stability analysis of a two step model for monitoring anaerobic digestion processes, J. Proc. Control., **22**, 1008–1019 (2012).

[4] O. Bernard, Z. Hadj-Sadok, D. Dochain, A. Genovesi and J-P. Steyer, Dynamical model development and parameter identification for an anaerobic wastewater treatment process, Biotechnol. Bioeng., **75**, 424–438 (2001).

[5] M. P. Boer, B. W. Kooi and S. A. L. M. Kooijman, Food chain dynamics in the chemostat, Math. Biosci., **150**, 43–62 (1998).
[6] F. Borsali and K. Yadi, Contribution to the study of the effect of the interspecificity on a two nutrients competition model, *Int. J. Biomath.*, **8**, 1550008, 17 pp. (2015).

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

[8] 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).

[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] A. Dhooge, W. Govaerts, Y. 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).

[11] M. El Hajji, How can inter-specific interferences explain coexistence or confirm the competitive exclusion principle in a chemostat?, *Int. J. Biomath.*, **11**, 1850111, 20 pp. (2018).

[12] M. El Hajji, J. Harmand, H. Chaker and C. Lobry, Association between competition and obligate mutualism in a chemostat, *J. Biol. Dyn.*, **3**, 635–647 (2009).

[13] M. El Hajji, F. Mazenc and J. Harmand, A mathematical study of a syntrophic relationship of a model of anaerobic digestion process, *Mathematical Biosciences and Engineering*, **7**, 641–656 (2010).

[14] 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).

[15] 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).

[16] 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).

[17] 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).

[18] 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).

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

[20] S.R. Hansen and S.P. Hubbell, Single-nutrient microbial competition: Qualitative agreement between experimental and theoretically forecast outcomes, *Science*, **286**, 1491–1493 (1980).

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

[22] J. Harmand, A. Rapaport, D. Dochain and C. Lobry, Microbial ecology and bioprocess control: Opportunities and challenges, *Journal of Process Control*, **18**, 865–875 (2008).

[23] 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).

[24] B. W. Kooi and M. P. Boer, Chaotic behaviour of a predator-prey system in the chemostat, *Dyn. Contin. Discrete Impuls. Syst. Ser. B Appl. Algorithms.*, **10**, 259–272 (2003).

[25] B. Li and Y. Kuang, Simple food chain in a chemostat with distinct removal rates, *J. Math. Anal. Appl.*, **242**, 75–92 (2000).

[26] 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).

[27] C. Lobry and F. Mazenc, Effect on persistence of intra-specific competition in competition models, *Electron. J. Diff. Equ.*, **125**, 1–10 (2007).

[28] 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).

[29] C. Lobry, A. Rapaport and F. Mazenc, Sur un modèle densité-dépendant de compétition pour une ressource, C. R. Biol., **329**, 63-70 (2006).

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

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

[32] T. Meadows, M. Weedermann, and G. S. K. Wolkowicz, Global analysis of a simplified model of anaerobic digestion and a new result for the chemostat, SIAM J. Appl. Math., **79**, 668-689 (2019).

[33] J. Monod, La technique de culture continue. Théorie et applications, Ann. Inst. Pasteur, **79**, 390-410 (1950).

[34] 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).

[35] T. Mtar, R. Fekih-Salem and T. Sari, Mathematical analysis of an interspecific density-dependent model with mutual inhibition and distinct removal rates, CARI'2022, African Conference on Research in Computer Science and Applied Mathematics. (2022).

[36] 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).

[37] A. Novick and L. Szilard, Description of the chemostat, Science, **112**, 715-716 (1950).

[38] S. S. Pilyugin and P. Waltman, The simple chemostat with wall growth, SIAM J. Appl. Math, **59**, 1552-1572 (1999).

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

[40] 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).

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

[42] M. Sbarciog, M. Loccfufer and E. Noldus, Determination of appropriate operating strategies for anaerobic digestion systems, Biochem. Eng. J., **51**, 180-188 (2010).

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

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

[45] J. H. Vandermeer and D. H. Bouchere, Varieties of mutualistic interaction in population models, J. Theor. Biol., **74**, 549-558 (1978).

[46] S. Vet, L. Gelens and D. Gonze, Mutualistic cross-feeding in microbial systems generates bistability via an Allee effect, Sci Rep., **10**, 7763 (2020).

[47] 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).

[48] 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).