# On a competition model of two toxin-producing species in the chemostat

Radhouane Fekih-Salem

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

## To cite this version:

Radhouane Fekih-Salem. On a competition model of two toxin-producing species in the chemostat. African Conference on Research in Computer Science and Applied Mathematics - CARI, Nov 2024, Bejaia, Algérie, Algeria. hal-04859900

HAL Id: hal-04859900

https://hal.science/hal-04859900v1

Submitted on 31 Dec 2024

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_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# On a competition model of two toxin-producing species in the chemostat

Radhouane Fekih-Salem¹[0000-0003-1168-4930]

University of Tunis El Manar, National Engineering School of Tunis, LAMSIN, 1002,
Tunis, Tunisia

University of Monastir, Higher Institute of Computer Science of Mahdia, 5147,
Mahdia, Tunisia radhouane.fekih-salem@enit.utm.tn

**Abstract.** Understanding and exploiting the competition and allelopathic effects of toxic substances on marine aquatic ecosystems or biodiversity in the soil or among plants is a major challenge in microbial biology and ecology. In this work, we focus on a model of the chemostat of two toxin-producing microbial species competing for a single nutrient. Each species produces toxin which affects the growth of the other species as well as its growth. Removal rates are distinct and include the specific mortality rate and autotoxicity of each species. Our study gives a complete analysis of the model by analyzing the joined effects of competition, allelopathy, and mortality. Using general monotonic growth rates, the mathematical study establishes the conditions of existence and local stability of all equilibria according to the control parameters. The model can present the multiplicity of the coexistence steady states. Numerical simulations show that the three-dimensional system can exhibit bistability where the asymptotic behavior of the system depends on the initial condition.

**Keywords:** Allelopathy · Chemostat · Coexistence · Stability · Toxicity.

## 1 Introduction

The chemostat is an important laboratory device used for experiments on the controlled growth of microorganisms in microbial ecology and biology, with applications in water treatment, biomass energy recovery, and biotechnology. The mathematical study of the classical chemostat model of several species competing on the same limiting resource shows that only the most competitive species that consumes less substrate to reach its equilibrium survives the competition while all other species are excluded [11,19]. This result is well known as the *Competitive Exclusion Principle* (CEP) which states that only the species which makes optimal use of the resource can survive to the competition. However, the biodiversity observed in nature and microbial ecosystems contradicts the CEP.

To reconcile the mathematical results of the classical chemostat model asserting the CEP and the experimental results and nature showing the biodiversity of microbial species, many recent researches aimed at understanding and
explaining biodiversity in microbial ecosystems by taking into account various mechanisms to promote the coexistence of species. In the literature, we can cite the following coexistence mechanisms: intra- and interspecific competition [1,3], flocculation [7,9,10], density-dependence [8,13,14,16,17], presence of an internal or external inhibitors [4,5], and the references therein. Although these various mechanisms of coexistence have been widely analyzed, the allelopathic effect of the most competitive species producing toxic substances may represent another mechanism of biodiversity between species that has been relatively neglected in the existing literature. In addition, the allelopathic effects of toxic substances on marine aquatic ecosystems or biodiversity in the soil or among plants remain poorly understood despite their environmental importance [20].

Following the modeling of allelopathic competitions in Maynard-Smith [15] based on the formulation of phenomenological Lotka–Volterra type, we consider an original model of two species competing for a single nutrient in the chemostat by producing a toxic substance. The dynamical equations of the mechanistic model are

$$ egin{cases} \dot{S} = D(S_{in} - S) - \mu_1(S) \frac{x_1}{y_1} - \mu_2(S) \frac{x_2}{y_2}, \\ \dot{x}_1 = [\mu_1(S) - \alpha_1 x_1 x_2^2 - D_1]x_1, \\ \dot{x}_2 = [\mu_2(S) - \alpha_2 x_2 x_1^2 - D_2]x_2, \end{cases} \quad (1) $$

where at time $t$, $S(t)$ is the concentration of the substrate; $x_1(t)$ and $x_2(t)$ are the concentrations of the two toxic species; $\mu_1(S)$ and $\mu_2(S)$ are the growth functions of the two species; $D$ and $S_{in}$ are the dilution rate and the input substrate concentration in the chemostat; $y_i$, $i = 1, 2$, is the yield coefficient; $D_i$ is the disappearance rate of the species $i$ that can take the form

$$ D_i = \theta_i D + m_i, \quad \text{where } \theta_i \in [0, 1], \quad \text{and } m_i = \epsilon_i + a_i, $$

with $\theta_i$ is the proportion of the species $i$ leaving the reactor; $m_i$ represents the nonnegative mortality rate of the species $i$ which is the sum of the specific death rate $\epsilon_i$ and the autotoxicity parameters $a_i$. Following Solé et al. [20], the term $\alpha_i x_i x_j^2$, $i, j = 1, 2$, $i \neq j$ is the allelopathic (or phytotoxic) interaction of the species $j$ on the species $i$ with $\alpha_i$ is the phytotoxic coefficient. Thus, the negative effects of allelopathic interactions increase with the concentrations of the two toxic species $x_1$ and $x_2$.

In [2], the author has considered a two-dimensional system of Lotka–Volterra type (without substrate) where two competitive species are also affected by toxic substances. It is shown that the ratio of the toxic substances of the two species has some stabilizing effect on the system. In [20], the simple Lotka–Volterra type model with allelopathic interactions between two marine phytoplankton species was modified such that the allelopathic term depends on the product of the square of the concentration of the target species by the concentration of the toxic species to have a good fit between the model and the experimental data of the toxic alga.

A complete mathematical study of model (1) has been done in [6] in the particular case where $\alpha_2 = 0$, that is, only the second species $x_2$ is toxic for the first one $x_1$. However, in [6], a model of two species competing for a single
resource involving also an intra- and interspecific competition with allelopathic effects in the dynamics of two species has been studied.

In [18], the author has introduced a chemostat model with two microbial species competing for a single resource. The model is a three-dimensional system of ordinary differential equations where a single species that is toxin-producing thereby having an allelopathic effect on the other. However, the growth rates are specified to be of Michaelis-Menten type. In Kengwoung-Keumo [12], an extension of Roy's competition model between non-toxic phytoplankton and toxic phytoplankton was studied by allowing a general class of monotonic growth rates.

This paper aims to provide a complete mathematical and numerical analysis of the model of two toxin-producing microbial species competing for a single nutrient in a chemostat where each species produces toxin which affects the growth of the other species as well as its growth. We consider distinct removal rates which include the specific mortality rate and autotoxicity of each species such that the three-dimensional system cannot be reduced to a two-dimensional one. Our study provides an understanding of the joint effects of competition, allelopathy, and mortality. By allowing a larger class of growth functions, we determine the conditions of existence and local stability of all equilibria according to the control parameters.

The paper is organized as follows. In section 2, we present the general assumption on the growth rates with some preliminary results. In section 3, we analyze the existence and the multiplicity of positive equilibria according to the input concentration of nutrient. In section 4, we study the local stability of all equilibria according to the two operating parameters which are the dilution rate and input concentration of nutrient. In section 5, the numerical simulations show the trajectories in the three-dimensional space. In the last section 6, we conclude with a discussion of the effect of allelopathic interactions between two toxin-producing species on process behavior. The parameter values used in simulations are provided in Appendix A.

## 2 Assumption and preliminary results

In this work, we consider a general assumption on the growth rates of model (1) which are continuously differentiable ($C^1$) functions satisfying the following hypothesis.

(H1) $\mu_i(0) = 0$ and $\mu'_i(S) > 0$, for all $S > 0$ so that $\mu_i(+\infty) < +\infty$ and $\mu'_i(+\infty) = 0$.

That is, each species can only grow in the presence of the nutrient. In addition, the growth rate of each species increases by increasing nutrient concentration.

The following result shows that the model is well defined from a biological point of view by determining the positivity of solutions that are positively bounded where the proof is standard.

**Proposition 1.** Assume that hypothesis (H1) holds for the dynamical system (1). Then, for any nonnegative initial condition, the solution of (1) exists, re-
mains nonnegative and is bounded for all $t \ge 0$. Moreover, the compact 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 (1), with $D_{min} = \min(D, D_1, D_2)$.

### 3 Existence of steady states

By setting the right-hand sides of equations of (1) equal to zero, we can find the equilibria which are the solutions of the set of the following equations

$$ \left\{ \begin{array}{l} 0 = D(S_{in} - S) - \mu_1(S) \frac{x_1}{y_1} - \mu_2(S) \frac{x_2}{y_2}, \\ 0 = [\mu_1(S) - \alpha_1 x_1 x_2^2 - D_1] x_1, \\ 0 = [\mu_2(S) - \alpha_2 x_2 x_1^2 - D_2] x_2. \end{array} \right. \qquad (2) $$

Thus, model (1) has at most four types of equilibria:

- $\mathcal{E}_0 = (S_{in}, 0, 0)$, $(x_1 = x_2 = 0)$: both species are excluded. It always exists.

- $\mathcal{E}_1 = (\tilde{S}_1, \tilde{x}_1, 0)$, $(x_2 = 0, \tilde{x}_1 > 0)$: only the second species is excluded.

- $\mathcal{E}_2 = (\tilde{S}_2, 0, \tilde{x}_2)$, $(x_1 = 0, \tilde{x}_2 > 0)$: only the first species is excluded.

- $\mathcal{E}^* = (S^*, x_1^*, x_2^*)$, $(x_1^* > 0, x_2^* > 0)$: both species coexist.

To establish the existence conditions of equilibria of model (1), we need to define the usual break-even concentrations and the following auxiliary functions.

<table><caption>**Table 1.** Auxiliary functions when $i, j = 1, 2$ and $i \neq j$.</caption><thead><tr><th>Definition</th><th></th></tr></thead><tbody><tr><td>λ<sub>i</sub>(D) is the unique solution of equation μ<sub>i</sub>(S) = D<sub>i</sub>.</td><td></td></tr><tr><td>λ<sub>i</sub>(D) It is defined for all 0 ≤ D &lt; D&#773;<sub>i</sub> := (μ<sub>i</sub>(+∞) - m<sub>i</sub>)/θ<sub>i</sub>.<br>If D ≥ D&#773;<sub>i</sub>, by convention we let λ<sub>i</sub>(D) = +∞.</td><td></td></tr><tr><td>φ<sub>i</sub>(S) It is defined and nonnegative for all D &lt; D&#773;<sub>i</sub> and S ≥ λ<sub>i</sub>(D).</td><td></td></tr><tr><td>X<sub>i</sub>(S) X<sub>i</sub>(S) = [φ<sup>2</sup><sub>j</sub>(S)/φ<sub>i</sub>(S)]<sup>1/3</sup>.<br>It is defined and positive for all D &lt; D&#773;<sub>i</sub> and S &gt; max(λ<sub>1</sub>(D), λ<sub>2</sub>(D)).</td><td></td></tr><tr><td>H(S) H(S) = μ<sub>1</sub>(S)X<sub>1</sub>(S)/y<sub>1</sub> + μ<sub>2</sub>(S)X<sub>2</sub>(S)/y<sub>2</sub>. It is defined and nonnegative,<br>for all D &lt; min(D&#773;<sub>1</sub>, D&#773;<sub>2</sub>) and S ≥ max(λ<sub>1</sub>(D), λ<sub>2</sub>(D)).</td><td></td></tr></tbody></table>

Due to the particularity of the structure of model (1), we will use the concept of steady-state characteristic that was introduced by Lobry et al. [13,14]. The next proposition establishes the necessary and sufficient conditions of existence of all equilibria of model (1).

**Proposition 2.** Let the hypothesis (H1) hold. Then, the four types of equilibria and their existence conditions for the dynamical system (1) are given in Table 2.
<figure><img src="image_3.png" /><figcaption>Fig. 1. Case $\lambda_1 < \lambda_2$: (a) curves of the functions $\mu_i$ and $\varphi_i(S)$, and the corresponding interval $I$; (b) Definition of the functions $X_i(S)$ on their definition domain.</figcaption></figure>

<table><caption>**Table 2.** Components and existence conditions of all equilibria of (1). The functions $H(S)$ and $X_i(S)$, $i = 1, 2$ are defined in Table 1.</caption><thead><tr><th>(S, x<sub>1</sub>, x<sub>2</sub>) components</th><th>Existence conditions</th></tr></thead><tbody><tr><td>E<sub>0</sub> (S<sub>in</sub>, 0, 0)</td><td>Always exists</td></tr><tr><td>E<sub>1</sub> (λ<sub>1</sub>, y<sub>1</sub>D(S<sub>in</sub> - λ<sub>1</sub>)/D<sub>1</sub>), 0)</td><td>S<sub>in</sub> > λ<sub>1</sub></td></tr><tr><td>E<sub>2</sub> (λ<sub>2</sub>, 0, y<sub>2</sub>D(S<sub>in</sub> - λ<sub>2</sub>)/D<sub>2</sub>))</td><td>S<sub>in</sub> > λ<sub>2</sub></td></tr><tr><td>E<sup>*</sup> (S<sup>*</sup>, X<sub>1</sub>(S<sup>*</sup>), X<sub>2</sub>(S<sup>*</sup>)); S<sup>*</sup> is solution of equation D(S<sub>in</sub> - S<sup>*</sup>) = H(S<sup>*</sup>)</td><td>Equation D(S<sub>in</sub> - S<sup>*</sup>) = H(S<sup>*</sup>) has a solution S<sup>*</sup> ∈ I = (max(λ<sub>1</sub>, λ<sub>2</sub>), +∞)</td></tr></tbody></table>

<figure><img src="image_4.png" /><figcaption>Fig. 2. Case $\lambda_1 < \lambda_2$: Existence and local stability of positive equilibria of model (1). (a) $S_{in} = 10 < S^c \approx 14.928$, (b) $S_{in} = 20 > S^c$. The red [resp. blue] color corresponds to LES [resp. unstable] equilibrium point.</figcaption></figure>
*Proof.* To obtain the equilibria of system (1), we set the right-hand sides of the equations in (1) equal to zero :

$$0 = D(S_{in} - S) - \mu_{1}(S)x_{1}/y_{1} - \mu_{2}(S)x_{2}/y_{2}, \quad (3)$$

$$0 = [\mu_{1}(S) - \alpha_{1}x_{1}x_{2}^{2} - D_{1}]x_{1}, \quad (4)$$

$$0 = [\mu_{2}(S) - \alpha_{2}x_{1}^{2}x_{2} - D_{2}]x_{2}. \quad (5)$$

- For $\mathcal{E}_0$, we have $x_1 = x_2 = 0$. Then, the first equation (3) implies that in $S = S_{in}$. Hence, we obtain $\mathcal{E}_0$ which always exists.

- For $\mathcal{E}_1$, we have $x_2 = 0$ and $x_1 > 0$. Hence, the two first equations (3) and (4) imply that $S = \lambda_1$ and $x_1 = Dy_1(S_{in} - \lambda_1)/D_1$. Thus, $\mathcal{E}_1$ exists if and only if $S_{in} > \lambda_1$.

- For $\mathcal{E}_2$, we have $x_1 = 0$ and $x_2 > 0$. Hence, equations (3) and (5) imply that $S = \lambda_2$ and $x_2 = Dy_2(S_{in} - \lambda_2)/D_2$. Thus, $\mathcal{E}_2$ exists if and only if $S_{in} > \lambda_2$.

- For $\mathcal{E}^*$, we have $x_1^* > 0$ and $x_2^* > 0$. Thus, a positive steady state $(S^*, x_1^*, x_2^*)$ must be a solution of the following system

$$0 = D(S_{in} - S^*) - \mu_{1}(S^*)x_{1}^*/y_{1} - \mu_{2}(S^*)x_{2}^*/y_{2}, \quad (6)$$

$$0 = \mu_{1}(S^{*}) - \alpha_{1}x_{1}^{*}x_{2}^{*2} - D_{1}, \quad (7)$$

$$0 = \mu_{2}(S^{*}) - \alpha_{2}x_{1}^{*2}x_{2}^{*} - D_{2}. \quad (8)$$

From equation (7) and the definition of $\varphi_i(S)$ in Table 1, we have

$$x_{1}^{*} = \frac{\varphi_{1}(S^{*})}{x_{2}^{*2}}, \quad \text{for all } S^{*} \in [\lambda_{1}, +\infty). \quad (9)$$

Similarly, from equation (8), we have

$$x_{2}^{*} = \frac{\varphi_{2}(S^{*})}{x_{1}^{*2}}, \quad \text{for all } S^{*} \in [\lambda_{2}, +\infty). \quad (10)$$

By replacing $x_2^*$ by its expression (10) in (9), we obtain

$$x_{1}^{*} = \frac{\varphi_{1}(S^{*})}{\varphi_{2}^{2}(S^{*})}x_{1}^{4}, \quad \text{for all } S^{*} \in I := (\max(\lambda_{1}, \lambda_{2}), +\infty).$$

Consequently, we deduce that

$$x_{1}^{*} = \left[ \frac{\varphi_{2}^{2}(S^{*})}{\varphi_{1}(S^{*})} \right]^{\frac{1}{3}} =: X_{1}(S^{*}), \quad \text{for all } S^{*} \in I. \quad (11)$$

Using expression (11) of $x_1^*$ in formula (10) of $x_2^*$, we obtain

$$x_{2}^{*} = \left[ \frac{\varphi_{1}^{2}(S^{*})}{\varphi_{2}(S^{*})} \right]^{\frac{1}{3}} =: X_{2}(S^{*}), \quad \text{for all } S^{*} \in I. \quad (12)$$

Replacing $x_i^*$, $i = 1, 2$ with their expressions $X_i(S^*)$, it follows that $S^*$ is a solution of the equation $D(S_{in} - S^*) = H(S^*)$.
In the following, we study the functions $X_i(S)$, $i = 1, 2$ and the function $H(S)$ defined in Table 1 to determine the multiplicity of positive equilibria.

**Proposition 3.** Let $\eta_i$ be the positive constant defined by

$$ \eta_i = \left[ \frac{\alpha_j(\mu_i(+\infty) - D_i)^2}{\alpha_i^2(\mu_j(+\infty) - D_j)} \right]^{\frac{1}{2}}, \quad i, j = 1, 2 \quad i \neq j. \quad (13) $$

We have

$$ \lim_{S \to +\infty} X_i(S) = \eta_j. $$

1. When $\lambda_1 < \lambda_2$, $X_i(S)$, $i = 1, 2$, are defined and positive for all $S \in (\lambda_2, +\infty)$ so that $X_1(\lambda_2) = 0$ and we have (see Fig. 1)

$$ \lim_{S \to \lambda_2} X_2(S) = +\infty. $$

2. When $\lambda_1 > \lambda_2$, $X_i(S)$, $i = 1, 2$, are defined and positive for all $S \in (\lambda_1, +\infty)$ so that $X_2(\lambda_1) = 0$ and we have

$$ \lim_{S \to \lambda_1} X_1(S) = +\infty. $$

*Proof.* From the definitions of the functions $\varphi_i(S)$ and $X_i(S)$ in Table 1, it is easy to see that $X_i(+\infty) = \eta_j$.

1. When $\lambda_1 < \lambda_2$, we have $X_1(\lambda_2) = 0$ and $X_2(\lambda_2) = +\infty$ since $\varphi_2(\lambda_2) = 0$.

2. Conversely, when $\lambda_1 > \lambda_2$, we have $X_2(\lambda_1) = 0$ and $X_1(\lambda_1) = +\infty$ since $\varphi_1(\lambda_1) = 0$.

Next, we study the function $H$ on its definition domain $I = (\max(\lambda_1, \lambda_2), +\infty)$. According to the definition of the function $H$ in Table 1, we obtain

$$ H'(S) = \frac{1}{y_1}(\mu'_1(S)X_1(S) + \mu_1(S)X'_1(S)) + \frac{1}{y_2}(\mu'_2(S)X_2(S) + \mu_2(S)X'_2(S)). \quad (14) $$

Replacing $X_i(S)$ with its definition in Table 1, a direct calculation yields that

$$ H'(S) = \frac{\mu'_1 F_1 + \mu'_2 F_2}{(\varphi_1 \varphi_2)^{\frac{2}{3}}}, \quad (15) $$

where

$$ \left\{ \begin{array}{l} F_1(S) = \frac{2}{3} (\varphi_1 \varphi_2)^{\frac{1}{3}} \frac{\mu_2}{y_2 \alpha_1} + \varphi_1^{-\frac{2}{3}} \varphi_2^{\frac{4}{3}} (\frac{2}{3} \mu_1 - D_1)/(\alpha_1 y_1), \\ F_2(S) = \frac{2}{3} (\varphi_1 \varphi_2)^{\frac{1}{3}} \frac{\mu_1}{y_1 \alpha_2} + \varphi_2^{-\frac{2}{3}} \varphi_1^{\frac{4}{3}} (\frac{2}{3} \mu_2 - D_2)/(\alpha_2 y_2). \end{array} \right. \quad (16) $$

Note that the sign of $H'(S)$ can change and consequently the function $H(S)$ can be nonmonotonic (see Fig. 2). According to the position of $\lambda_i$, Figure 2 shows that the function $H$ is nonmonotonic on its domain of definition $I$. Furthermore, the positive equilibrium point is given by the intersection of the curve of $H$ with the line $\Delta$ of equation $\delta(S) = D(S_{in} - S)$.

The following result establishes the multiplicity of positive equilibria according to the value of $S_{in}$.
**Proposition 4.** If Hypothesis (H1) holds, then there exists a least one solution $S_c^*$ on I of equation $H'(S) = -D$. For each solution $S_c^*$, there exists a critical values of $S_{in}$ denoted $S_{in}^c$ so that the curve of $H$ is tangent to the line $\Delta^c$ (see Fig. 2) :

$$S_{in}^{c} = S_{c}^{*} + \frac{H(S_{c}^{*})}{D}. \qquad (17)$$

Assume that $S_c^*$ is unique. Then,

- If $S_{in} < S_{in}^c$, system (1) has no positive equilibrium.

- If $S_{in} > S_{in}^c$, there exist two positive equilibria denoted $E_1^*$ and $E_2^*$.

When $S_c^*$ is not unique, then generically the system can have an odd number of positive equilibria.

*Proof.* Let $\bar{\lambda} := \max(\lambda_1, \lambda_2)$. Let $\phi$ be the function defined by

$$\phi(S) := H(S) - \delta(S), \quad \text{for all } S \in I := (\bar{\lambda}, +\infty).$$

Hence,

$$\phi'(S) := H'(S) - \delta'(S) = H'(S) + D, \quad \text{for all } S \in I.$$

From Hypothesis (H1) and the expression of $H'$ in (14), we have

$$\lim_{S \to \bar{\lambda}} H'(S) = -\infty, \quad \lim_{S \to +\infty} H'(S) = 0.$$

In addition, the function $\phi'(S)$ is continuous on $I = (\bar{\lambda}, +\infty)$ so that $\phi'(\bar{\lambda}) = -\infty$ and $\phi'(+\infty) = D$. Thus, there exists a least one solution $S_c^*$ on I of equation $\phi'(S) = 0$ or equivalently $H'(S) = -D$.

The curve of $H$ is tangent to the line $\Delta^c$ when there exists $S_{in}^c$ so that $H(S_c^*) = \delta^c(S_c^*) = D(S_{in}^c - S_c^*)$. Thus, we obtain the expression of $S_{in}^c$ in (17).

In other hand, we have $\phi(S)$ is continuous on $(\bar{\lambda}, S_{in}]$ so that $\phi(\bar{\lambda}) = +\infty$ and $\phi(S_{in}) = H(S_{in}) > 0$. Assume that $S_c^*$ is unique so that

$$\phi(S_c^*) = H(S_c^*) - \delta(S_c^*) = \delta^c(S_c^*) - \delta(S_c^*).$$

If $S_{in} < S_{in}^c$, then $\phi(S_c^*) = D(S_{in}^c - S_c^*) - D(S_{in} - S_c^*) > 0$ so that the equation $\phi(S) = 0$, or equivalently the equation $H(S) = \delta(S)$ has no solution on I. Conversely, if $S_{in} > S_{in}^c$, then $\phi(S_c^*) < 0$ so that the equation $H(S) = \delta(S)$ has two solutions on I. Finally, if $S_c^*$ is not unique, then generically the equation $H(S) = \delta(S)$ can have an odd number of solutions.

# 4 Stability of steady states

In this section, we determine the conditions of local stability of all equilibria of model (1). Note that all reference to the stability of a given equilibrium should be considered as local exponential stability, i.e. all the real parts of the eigenvalues of the Jacobian matrix are negative. We show the next result.
<table><caption>**Table 3.** Local stability conditions of all equilibria of (1). The functions $H(S)$ and $X_i(S)$, $i = 1, 2$ are defined in Table 1, while $c_4$ is defined in Table 4.</caption><thead><tr><th colspan="2">Local stability conditions</th></tr></thead><tbody><tr><td>$\mathcal{E}_0$</td><td>$S_{in} < \min(\lambda_1, \lambda_2)$</td></tr><tr><td>$\mathcal{E}_1$</td><td>$\lambda_1 < \lambda_2$</td></tr><tr><td>$\mathcal{E}_2$</td><td>$\lambda_2 < \lambda_1$</td></tr><tr><td>$\mathcal{E}^*$</td><td>$H'(S^*) < -D$ and condition $c_4 > 0$ holds</td></tr></tbody></table>

<caption>**Table 4.** Routh-Hurwitz coefficient $c_4$ for $\mathcal{E}^*$ where all function are evaluated at $(S^*, x_1^*, x_2^*)$ components.</caption>

$$
\begin{align*}
c_4 &= a_1\mu_1''(S) + a_2\mu_2''(S) + a_3\mu_1'(S)\mu_2'(S) + a_4\mu_1'(S) + a_5\mu_2'(S) + a_6, \\
a_1 &= (\alpha_1\varphi_1 + \alpha_2\varphi_2 + y_1\mu_1)x_1^2/y_1^2, \\
a_2 &= (\alpha_1\varphi_1 + \alpha_2\varphi_2 + y_2\mu_2)x_2^2/y_2^2, \\
a_3 &= (2\alpha_1\varphi_1 + 2\alpha_2\varphi_2 + y_1\mu_1 + y_2\mu_2)x_1x_2/(y_1y_2) \\
a_4 &= (\alpha_1^2\varphi_1^2 + \alpha_2^2\varphi_2^2 + y_1\alpha_1\mu_1\varphi_1 + 2y_1\alpha_2\mu_2\varphi_1 + 2D\alpha_1\varphi_1 + 2D\alpha_2\varphi_2 + 2\alpha_1\alpha_2\varphi_1\varphi_2 \\
&\quad + y_1D\mu_1)x_1/y_1 \\
a_5 &= (\alpha_1^2\varphi_1^2 + \alpha_2^2\varphi_2^2 + y_2\alpha_2\mu_2\varphi_2 + 2y_2\alpha_1\mu_1\varphi_2 + 2D\alpha_1\varphi_1 + 2D\alpha_2\varphi_2 + 2\alpha_1\alpha_2\varphi_1\varphi_2 \\
&\quad + y_2D\mu_2)x_2/y_2 \\
a_6 &= (\alpha_1\varphi_1 + \alpha_2\varphi_2)(D^2 + D(\alpha_1\varphi_1 + \alpha_2\varphi_2) - 3\alpha_1\alpha_2\varphi_1\varphi_2).
\end{align*}
$$

**Proposition 5.** Assume that hypothesis (H1) holds. Then, the local stability conditions of the four types of equilibria of system (1) are given in Table 3.

*Proof.* To investigate the asymptotic stability of all steady states of model (1), we shall use the abbreviation LES for Locally Exponentially Stable. The Jacobian matrix of model (1) at a steady state $(S, x_1, x_2)$ is given by the $3 \times 3$ matrix:

$$J = \begin{bmatrix} -D - \mu_1'(S) \frac{x_1}{y_1} - \mu_2'(S) \frac{x_2}{y_2} & -\frac{\mu_1(S)}{y_1} & -\frac{\mu_2(S)}{y_2} \\ \mu_1'(S)x_1 & \mu_1(S) - 2\alpha_1 x_1 x_2^2 - D_1 & -2\alpha_1 x_1^2 x_2 \\ \mu_2'(S)x_2 & -2\alpha_2 x_1 x_2^2 & \mu_2(S) - 2\alpha_2 x_1^2 x_2 - D_2 \end{bmatrix}.$$

– For $\mathcal{E}_0 = (S_{in}, 0, 0)$, the Jacobian matrix is given by

$$J_{E_0} = \begin{bmatrix} -D & -\mu_1(S_{in})/y_1 & -\mu_2(S_{in})/y_2 \\ 0 & \mu_1(S_{in}) - D_1 & 0 \\ 0 & 0 & \mu_2(S_{in}) - D_2 \end{bmatrix}.$$

The three eigenvalues are on the elements on the diagonal of $J_{E_0}: -D$, $\mu_1(S_{in}) - D_1$ and $\mu_2(S_{in}) - D_2$. They are strictly negative if and only if $S_{in} < \lambda_i$, pour $i = 1, 2$.

– For $\mathcal{E}_1 = (\lambda_1, \bar{x}_1, 0)$, the existence condition is $S_{in} > \lambda_1$ and the Jacobian matrix is given by

$$J_{E_1} = \begin{bmatrix} -D - \mu_1'(\lambda_1)\bar{x}_1/y_1 - D_1/y_1 & -\mu_2'(\lambda_1)/y_2 & 0 \\ \mu_1'(\lambda_1)\bar{x}_1 & 0 & 0 \\ 0 & 0 & \mu_2(\lambda_1) - D_2 \end{bmatrix}.$$
Thus, the characteristic polynomial is given by

$$P(\lambda) = (\mu_2(\lambda_1) - D_2 - \lambda)(\lambda^2 + c_1\lambda + c_2).$$

where

$$c_1 = D + \mu'_1(\lambda_1)\bar{x}_1/y_1 \quad \text{and} \quad c_2 = D_1\mu'_1(\lambda_1)\bar{x}_1/y_1.$$

Since $c_1 > 0$ and $c_2 > 0$, the real parts of the roots of the quadratic factor are negative. Therefore, $\mathcal{E}_1$ is LES if and only if $\mu_2(\lambda_1) < D_2$, that is, $\lambda_1 < \lambda_2$.

- For $\mathcal{E}_2 = (\lambda_2, 0, \bar{x}_2)$, the existence condition is $S_{in} > \lambda_2$ and the Jacobian matrix is given by

$$J_{E_2} = \begin{bmatrix} -D - \mu'_2(\lambda_2)\bar{x}_2/y_2 & -\mu_1(\lambda_2)/y_1 & -D_2/y_2 \\ 0 & \mu_1(\lambda_2) - D_1 & 0 \\ \mu'_2(\lambda_2)\bar{x}_2 & 0 & 0 \end{bmatrix}.$$

Thus, the characteristic polynomial is given by

$$P(\lambda) = (\mu_1(\lambda_2) - D_1 - \lambda)(\lambda^2 + c_1\lambda + c_2).$$

where

$$c_1 = D + \mu'_2(\lambda_2)\bar{x}_2/y_2 \quad \text{and} \quad c_2 = D_2\mu'_2(\lambda_2)\bar{x}_2/y_2.$$

Since $c_1 > 0$ and $c_2 > 0$, the real parts of the roots of the quadratic factor are negative. Therefore, $\mathcal{E}_2$ is LES if and only if $\mu_1(\lambda_2) < D_1$, that is, $\lambda_2 < \lambda_1$.

- For $E^* = (S^*, x_1^*, x_2^*)$ the Jacobian matrix is given by

$$J_{E^*} = \begin{bmatrix} -m_{11} & -m_{12} & -m_{13} \\ m_{21} & -m_{22} & -m_{23} \\ m_{31} & -m_{32} & -m_{33} \end{bmatrix}$$

where

$$\left\{ \begin{array}{l} m_{11} = D + \mu'_1(S^*) \frac{x_1^*}{y_1} + \mu'_2(S^*) \frac{x_2^*}{y_2}, \quad m_{12} = \frac{\mu_1(S^*)}{y_1}, \quad m_{13} = \frac{\mu_2(S^*)}{y_2}, \\ m_{21} = \mu'_1(S^*) x_1^*, \quad m_{22} = -\mu_1(S^*) + 2\alpha_1 x_1^*(x_2^*)^2 + D_1, \quad m_{23} = 2\alpha_1 x_1^{*2} x_2^*, \\ m_{31} = \mu'_2(S^*) x_2^*, \quad m_{32} = 2\alpha_2 x_1^{*2} x_2^*, \quad m_{33} = -\mu_2(S^*) + 2\alpha_2 x_1^{*2} x_2^* + D_2. \end{array} \right. \qquad (18)$$

From (7), we have

$$m_{22} = \alpha_1 x_1^*(x_2^*)^2,$$

which is positive. From (8), we have

$$m_{33} = \alpha_2 x_1^{*2} x_2^*,$$

which is positive. Thus, all $m_{ij}$ are positive for all $i, j = 1, 2, 3$. The characteristic polynomial is given by

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

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

According to the Routh–Hurwitz criterion, $E_1$ is LES if and only if

$$c_1 > 0, \ c_3 > 0 \quad \text{and} \quad c_4 = c_1c_2 - c_3 > 0. \tag{20}$$

Since all $m_{ij}$ are positive, then $c_1 > 0$. Next, we will show that the sign of $c_3$ is given by the position of the curve of function $H$ with respect to the line $\Delta$ of equation $\delta(S) = D(S_{in} - S)$. More precisely, we will establish the link between the sign of $c_3$ or equivalently the determinant of the Jacobian matrix $J_{E^*}$ at $E^* = (S^*, x_1^*, x_2^*)$ and $H'(S^*) + D$. In fact, we will show that

$$c_3 = -\det(J_{E^*}) = (-3\alpha_1\alpha_2\varphi_1(S^*)\varphi_2(S^*)) (D + H'(S^*)) \tag{21}$$

Indeed, replacing $x_1^*$ and $x_2^*$ by (11) and (12), respectively, we have

$$
\begin{align*}
c_3 = -3\alpha_1\alpha_2\varphi_1\varphi_2 \left[ D + \frac{1}{y_1}\mu'_1\left(\frac{\varphi_2^2}{\varphi_1}\right)^{\frac{1}{3}} + \frac{1}{y_2}\mu'_2\left(\frac{\varphi_1^2}{\varphi_2}\right)^{\frac{1}{3}} \right] \\
&\quad + \mu'_1\left[\frac{\varphi_2^2}{\varphi_1}\right]^{\frac{1}{3}} \left[ \frac{\alpha_2}{y_1}\mu_1\varphi_2 - \frac{2\alpha_2}{y_2}\mu_2\varphi_1 \right] + \mu'_2\left[\frac{\varphi_1^2}{\varphi_2}\right]^{\frac{1}{3}} \left[ \frac{\alpha_1}{y_2}\mu_2\varphi_1 - \frac{2\alpha_1}{y_1}\mu_1\varphi_2 \right],
\end{align*}
$$

A straightforward calculation yields that

$$c_3 = (-3\alpha_1\alpha_2\varphi_1\varphi_2) \left[ D + \frac{\mu'_1 F_1 + \mu'_2 F_2}{(\varphi_1\varphi_2)^{\frac{2}{3}}} \right],$$

where $F_1(S^*)$ and $F_2(S^*)$ are defined by (16). From the expression of $H'$ in (15), we obtain the formula of $c_3$ in (21). Since $\varphi_i(S^*)$ is positive for all $S^* \in I$, then $c_3$ is positive if and only if $H'(S^*) < -D$. Finally, for $c_4$, one can view after simplification that

$$
\begin{align*}
c_4 &= m_{11}(m_{12}m_{21} + m_{13}m_{31} + m_{11}m_{22} + m_{11}m_{33}) + m_{22}(m_{12}m_{21} \\
&\quad + m_{11}m_{22} + m_{11}m_{33}) + m_{33}(m_{13}m_{31} + m_{11}m_{22} + m_{11}m_{33}) \\
&\quad + m_{21}m_{32}m_{13} + m_{31}m_{12}m_{23} + (m_{22} + m_{33})(m_{22}m_{33} - m_{23}m_{32}).
\end{align*}
$$

Replacing the $m_{ij}$, by their expressions in (18), a direct calculation shows that the expression of $c_4$ can be written by (4). We conclude that if $H'(S^*) < -D$ and condition $c_4 > 0$ hold, then all the conditions of the Routh–Hurwitz criterion are satisfied so that the positive equilibrium is LES when it exists.

**Remark 1.** 1. If $H'(S^*) > -D$, then $c_3 < 0$, that is, a necessary condition of the Routh–Hurwitz criterion does not hold so that the positive equilibrium $\mathcal{E}^*$ is unstable.
2. Since $\varphi_i(S^*)$, $i = 1, 2$ and $a_j$, $j = 1, \dots, 5$ are positive for all $S^* \in I$, then from the expression of $c_4$ in (4) we see that $c_4$ is positive if

$$ (D^2 + D(\alpha_1\varphi_1 + \alpha_2\varphi_2) - 3\alpha_1\alpha_2\varphi_1\varphi_2) > 0. $$

However, we do not have at our disposal a necessary and sufficient condition of the stability of the positive equilibrium $\mathcal{E}^*$. The complete characterization of local stability of $\mathcal{E}^*$ remains an open problem.

## 5 Numerical simulations

Next, to illustrate the numerical simulations for model (1), we choose the following specific growth rates of Monod-type satisfying Hypothesis (H1):

$$ \mu_1(S) = \frac{\mu_1^m S}{k_1 + S}, \quad \mu_2(S) = \frac{\mu_2^m S}{k_2 + S}, \qquad (22) $$

where $\mu_1^m$ and $\mu_2^m$ are the maximum growth rates; $k_1$ and $k_2$ are the Michaelis-Menten constants. Fix the set of the parameter as the values provided in Table 5. Fig. 3(a) shows the three-dimensional space $(S, x_1, x_2)$ when $S_{in} = 10 < S_{in}^c \approx 14.928$ so that there is any coexistence equilibrium point. In this case, the only equilibrium points that exist are $\mathcal{E}_0$, $\mathcal{E}_1$ and $\mathcal{E}_2$. Since $\lambda_1 \approx 0.857 < \lambda_2 \approx 2.316 < S_{in} = 10$, the exclusion equilibrium of the second species $\mathcal{E}_1$ is LES while the washout equilibrium $\mathcal{E}_0$ and the exclusion equilibrium of the first species $\mathcal{E}_2$ are unstable. For different positive initial conditions, the numerical simulations show the convergence towards $\mathcal{E}_1$. Note that the three equilibria are given by

$$ \mathcal{E}_0 = (10, 0, 0), \quad \mathcal{E}_1 \approx (0.857, 1.524, 0) \quad \text{and} \quad \mathcal{E}_2 \approx (2.316, 0, 2.794). $$

<figure><img src="image_5.png" /><figcaption>Fig. 3. The three-dimensional space $(S, x_1, x_2)$ of system (1) (a) $S_{in} = 10 < S_{in}^c \approx 14.928$, all the solutions convergence to $\mathcal{E}_1$. (b) $S_{in} = 20 > S_{in}^c$, the system exhibits bistability with either convergence to $\mathcal{E}_1$ or $\mathcal{E}_1^*$.</figcaption></figure>
Fig. 3(b) shows that system (1) exhibits bistability between the coexistence equilibrium point $\mathcal{E}_1^*$ and the one of exclusion of the second species $\mathcal{E}_1$, in the case $S_{in} = 20 > S_{in}^c$. In this case, there exist five steady states whose boundary equilibria are

$$ \mathcal{E}_0 = (20, 0, 0), \quad \mathcal{E}_1 \approx (0.857, 3.19, 0), \quad \mathcal{E}_2 \approx (2.316, 0, 6.431), $$

while the coexistence equilibria are

$$ \mathcal{E}_1^* \approx (2.359, 0.051, 6.142) \text{ and } \mathcal{E}_2^* \approx (4.123, 0.454, 2.565). $$

# 6 Conclusion

In this study, we propose an original model of two toxin-producing microbial species competing for a single nutrient in a chemostat. The toxin produced by each species affects the growth of the other species as well as their growth. We allow a large class of monotonic growth functions with distinct removal rates including the specific mortality rates and autotoxicity of each species. The main contribution of this study is to bring out the mutual effects on the coexistence of the competition and allelopathy of two microbial species, which are not studied in the existing literature.

First, we establish the necessary and sufficient conditions of existence of all equilibria of model (1) according to the operating parameters $D$ and $S_{in}$. Then, we show that the model can present the multiplicity of the coexistence equilibria depending on the input concentration of nutrient $S_{in}$. Indeed, model (1) can have up to four types of equilibria: the washout equilibrium denoted by $\mathcal{E}_0$, two types of boundary equilibria, denoted by $\mathcal{E}_1$ and $\mathcal{E}_2$, where the first or the second competitor is extinct, respectively, and a positive equilibrium point, denoted by $\mathcal{E}^*$, where the two microbial populations coexist. When they exist, all equilibria are unique, except the coexistence equilibrium $\mathcal{E}^*$.

Then, we determine the necessary and sufficient conditions of the washout and the boundary equilibria $\mathcal{E}_0$, $\mathcal{E}_1$ and $\mathcal{E}_2$ according to the operating parameters. For the coexistence equilibrium $\mathcal{E}^*$, we provide only a sufficient condition of stability. The complete characterization of local stability of $\mathcal{E}^*$ remains an open problem that deserves further attention and will be the object of future work. As shown in the numerical simulations, model (1) can also exhibit bistability between the coexistence equilibrium and that of exclusion of the second or the first species.

# A Parameter values for system (1)

In Table 5, we provide the biological parameter values of the specific growth rates $\mu_1$ and $\mu_2$ of Monod-type given by (22).

Disclosure of Interests. The author has no competing interests to declare that are relevant to the content of this article.
<table><caption>**Table 5.** Parameter values for system (1).</caption><thead><tr><th>PARAMETER</th><th>μ<sub>1</sub><sup>m</sup></th><th>k<sub>1</sub></th><th>μ<sub>2</sub><sup>m</sup></th><th>k<sub>2</sub></th><th>α<sub>1</sub></th><th>α<sub>2</sub></th><th>y<sub>1</sub></th><th>y<sub>2</sub></th><th>θ<sub>1</sub></th><th>θ<sub>2</sub></th><th>m<sub>1</sub></th><th>m<sub>2</sub></th><th>D</th></tr></thead><tbody><tr><td>FIGS. 1-2-3</td><td>4</td><td>2</td><td>3</td><td>4</td><td>0.5</td><td>0.8</td><td>0.2</td><td>0.4</td><td>0.8</td><td>0.9</td><td>0.4</td><td>0.2</td><td>1</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. J. Chattopadhyay, Effect of toxic substances on a two-species competitive system, *Ecol. Modell.*, **84**, 287–289 (1996).

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

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

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

6. R. Fekih-Salem, Analysis of an intra-and interspecific interference model with allelopathic competition., *hal-04099563v1*, (2023).

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

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

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

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

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

12. J-J. Kengwoung-Keumo, Competition between a nonallelopathic phytoplankton and an allelopathic phytoplankton species under predation, *Math. Biosci. Eng.*, **13**, 787–812 (2016).

13. C. Lobry and 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. C. Lobry and 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).

15. J. Maynard-Smith, *Models in Ecology*, Cambridge University Press, London (1974).

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

18. S. Roy, The coevolution of two phytoplankton species on a single resource: Allopathy as a pseudo-mixotrophy, *Theor. Popul. Biol.*, **75**, 68–75 (2009).

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

20. J. Solé, E. García-Ladona, P. Ruardij and M. Estrada, Modelling allelopathy among marine algae, *Ecol. Modell.*, **373–384**, 7045–7073 (2005).