Commensalism and syntrophy in the chemostat: a unifying graphical approach

1.   Introduction

The study of interactions within microbial ecological communities is one of the most important issues in microbial ecology ^{[1, 2, 3, 4, 5]}. An interaction has an impact (neutral, beneficial, or detrimental) on the partner microorganisms involved. Commensalism and syntrophy are among the interactions that benefit the species. A commensalistic association is a relationship in which one microbe derives benefits from the other, while the other microbe is not affected (neither detrimental nor beneficial) by the association. A food chain can be considered a commensal interaction, when the second organism in the food chain lives off the waste products of the first, which in turn is unaffected by the interaction ^{[6, 7]}. In a syntrophic association, both organisms benefit from the presence of the other. A food chain can be seen as a syntrophic interaction when the first organism in the food chain is inhibited by the accumulation of its own waste in the environment. This inhibition is diminished by the second organism, which uses the wastes of the first one as a food source ^{[8, 9, 10, 11, 12, 13, 14]}.

The aim of this work is to present a unified graphical approach of the models of two species having a commensalistic or a syntrophic relationship. Our objective is threefold. First, we show that the graphical concepts introduced by Tilman ^{[15, 16]} to study patterns of competition allow us to study these patterns of commensalism and syntrophy in a unified graphical way. Second, we determine the operating diagram of the model, thus obtaining an extension of all the results on the operating diagrams of the literature. Third, we consider the case including mortality, which has not been addressed in previous studies.

Concerning the first objective, we show that the concepts and graphical methods of Tilman ^{[15, 16]}, complemented by Ballyk and Wolkowicz ^[17], to determine the outcome of competition between two species for two resources are very useful to understand the condition of existence and stability of the equilibria of commensalistic and syntrophic models. Although these systems are not competitive systems, it appears that the concepts introduced in ^{[15, 16, 17]} permit a unifying graphical approach of their study.

Concerning the second objective, we recall that the operating diagram has the operating parameters as its coordinates and the various regions defined within it correspond to qualitatively different asymptotic behaviors. Our model has three operating parameters that are the input concentrations of substrate and the dilution rate $D$. These parameters are control parameters since they are under the control of the experimenter. Apart from these three parameters, which can vary, all other parameters have biological meaning and are fitted using experimental data from ecological and/or biological observations of organisms and substrates. Therefore, the operating diagram is the bifurcation diagram that shows how the system behaves when we vary the control parameters. The importance of the operating diagram for bioreactors was emphasized in ^[18]. Although this diagram was not considered in the classic monograph on the chemostat ^[19], its importance was mentioned by Smith and Waltman, who stated that the operating diagram is probably the most useful tool for discussing the behavior of the model in relation to the parameters [19, p. 252]. The operating diagram is introduced in the recent book on the mathematical theory of the chemostat ^[20]. It is often constructed both in biological literature ^{[14, 21, 22, 23]} and in mathematical literature, in the study of anaerobic digestion ^{[24, 25, 26, 27]}, commensalistic and syntrophic systems ^{[12, 28, 29, 30, 31]}, microbial food-webs ^{[23, 32]}, inhibition ^{[33, 34]}, chemostats in series ^[35], and density-dependent models ^{[36, 37, 38]}.

Concerning the third objective, when there is no mortality, by using the theory of asymptotically autonomous systems, we can reduce the study to that of a model in dimension two so that the study of the stability is easy. In the general case, this reduction cannot be made and the Routh Hurwitz conditions must be used to determine whether the equilibrium is stable or not.

This paper is organized as follows: In Section 2, we present the mathematical model and give the assumptions made on the growth functions. In Section 3, we describe Tilman's graphical method and we show how the existence and stability condition of the equilibria of the model can be easily read on the position of the projections of equilibria into the feasible set. We also construct in this section the operating diagram of the model. In Section 4, we show that the graphical method also applies in the case where the species can be inhibited by their limiting substrate. Some complements on the operating diagrams are given in Section 5. Finally, some conclusions are drawn in Section 6. Global asymptotic stability results are given in Appendix A, and bifurcation diagrams with respect to one of the operating parameters are shown in Appendix B. The technical proofs are reported in Appendix C. In Appendix D, we propose a review of the results of the literature on commensalism and syntrophy.

2.   Mathematical model

We consider the two-species system modeled by

where $X_i$, $i = 1, 2$, represents the concentration of species $i$; $S_j$, $j = 1, 2$, is the concentration of chemical $j$ and $S_j^{in}$ its inlet concentration; $k_i$, $i = 1, 2, 3$, are yield factors, referring to the amount of chemical that is produced or consumed by the growth of a unit amount of the biomass of microbial species; and $D_i$, $i = 1, 2$, represents the disappearance rates of species $i$, and are modeled by

where $D = 1/HRT$ is the dilution rate, $HRT$ being the hydraulic retention time; the non-negative death (or decay) rate parameters $a_1$ and $a_2$ are taken into consideration; and the coefficient $\alpha_i\in[0, 1]$, $i = 1, 2$, represents the biomass proportion that leaves the bioreactor. This coefficient models the decoupling between solids and hydraulic residence time ^{[30, 31]}.

Well-mixed continuous bioreactors (i.e., chemostats) for the culture of multiple species are modeled with the commonly used Monod-type kinetics. For example, in ^[2], the growth rates of the two species, $\mu_i$, are modeled by multiplying substrate limitation terms described as Monod kinetics and inhibition terms as follows:

where $m_i$ is the maximum growth rate, $K_i$ is the half-velocity constant, and $L_{i}$ quantifies the strength of inhibition of substrate $S_j$, $j\neq i$, on species $i$. If $L_{i} = 0$, then there is no inhibition, and so $\mu_i$ depends only on $S_i$.

In Figure 1, we illustrate (2.1) by using the notations and representations of microbial communities proposed by Di and Yang ^[2]. Species $X_2$ is produced by consuming chemical $S_2$, which itself is produced by species $X_1$ through consuming chemical $S_1$, that acts as the substrate of the overall system. System C1, which corresponds to $L_1 = L_2 = 0$ in (2.3), is an example of pure commensalism, where the second population (the commensal population) depends for its growth on the first population (the host), and thus benefits from its interaction, while the host population is not affected by the growth of the commensal population. System C2, which corresponds to $L_1 = 0$ and $L_2 > 0$ in (2.3), is also an example of commensalism since it has a cascade structure and the first population is not affected by the growth of the second population.

On the other hand, systems S1, which corresponds to $L_1 > 0$ and $L_2 = 0$ in (2.3), and S2, which corresponds to $L_1 > 0$ and $L_2 > 0$ in (2.3), are not commensal since the first population is affected by the growth of the second population. For instance, in S1, the first organism is inhibited by the substrate $S_2$ produced by the degradation of $S_1$ by $X_1$. Hence, the extent to which the substrate $S_1$ is degraded by the organism $X_1$ to produce the substrate $S_2$, which is necessary for the growth of $X_2$, depends on the efficiency of the removal of the product $S_2$ by the bacteria $X_2$. Therefore, each population needs the other one for its growth, that is, there is a syntrophic relationship between them.

Remark 1. Case C1 corresponds to the base system depicted on [2, Figure 1]. Cases C2, S1, and S2 correspond, respectively, to cases 1-4, 1-1, and 2-3 shown in [2, Figures 2 and 3]. Notice that ^[2] also considered inhibitions between species, see cases 1-2 and 1-3 in [2, Figure 2], and five other systems obtained by combining them, see [2, Figure 3]. Although very important for microbial ecology, these species inhibitions will not be considered in this paper. The models C1, C2, S1, and S2 have been extensively studied in the literature, see Tables A3 and A4 in Appendix D.

There are many ways of representing inhibition other than (2.3). For example, instead of decreasing the maximum growth rate $m_i$ of the Monod function as in (2.3), we can increase its half-velocity constant (or substrate affinity) $K_i$:

where $L_{i}$ quantifies the strength of inhibition of $S_j$ on species $i$. See also the Kreikenbohm and Bohl function KB1 in Table A7.

In this work, instead of considering particular growth functions such as (2.3), (2.4), or KB1, we consider generalized growth functions characterized by their qualitative behaviors. Assume that $\mu_1, \mu_2: \mathbb{R}_+^2 \rightarrow \mathbb{R}_+$ are of class $\mathcal{C}^1$ such that $\mu_1(0, S_2) = 0$ for $S_2\geq 0$ and $\mu_2(S_1, 0) = 0$ for $S_1\geq 0$. Hence, no growth takes place for species $i = 1, 2$ without substrate $S_i$. This hypothesis is made throughout the paper and will not be repeated. The solutions of (2.1), with prescribed initial conditions, exist, are unique, are bounded, and the non-negative cone is positively invariant. The positive cone is also positively invariant. We make the following assumptions:

${\bf H1.}$ For $S_1 > 0$ and $S_2\geq 0$, we have $\frac{\partial\mu_1}{\partial S_1}(S_1, S_2) > 0$ and $\frac{\partial\mu_1}{\partial S_2}(S_1, S_2)\leq 0$.

${\bf H2.}$ For $S_1\geq 0$ and $S_2 > 0$, we have $\frac{\partial\mu_2}{\partial S_2}(S_1, S_2) > 0$ and $\frac{\partial\mu_2}{\partial S_1}(S_1, S_2)\leq 0$.

Hypotheses H1 and H2 signify that the growth for species $i = 1, 2$ increases with $S_i$, while it is inhibited by the other substrate $S_j$, $j\neq i$ : the first organism is inhibited by its product $S_2$ and the second organism is inhibited by the food $S_1$ of the first organism.

Remark 2. Since the partial derivative $\frac{\partial\mu_i}{\partial S_j}(S_1, S_2)$, for $j\neq i$, can be equal to zero, the inhibition of species $i$ by substrate $j\neq i$ is not mandatory. Hence, the assumptions H1 and H2 cover the cases C1, C2, S1, and S2 of Figure 1. For example, C1 corresponds to $\mu_1$ depending only on $S_1$ and increasing for all $S_1 > 0$ and $\mu_2$ depending only on $S_2$ and increasing for all $S_2 > 0$.

To reduce the number of parameters, model (2.1) is further converted to a simpler one where the yield factors $k_i$ are fixed to 1. We scale system (2.1) using the following change of variables and notations:

We obtain the system of differential equations

where the yield factors $k_i$ are fixed to 1 and the growth functions $f_1, f_2: \mathbb{R}_+^2 \rightarrow \mathbb{R}_+$ are given by

Since $\mu_1$ and $\mu_2$ in (2.1) satisfy the conditions of H1 and H2, then the growth functions $f_1$ and $f_2$ in (2.5) have the same qualitative properties.

3.   Results

3.1.   Graphical method

In this section, we explain the graphical method of ^{[15, 16, 17]}.

3.1.1.   The feasible set boundary (FSB)

The "feasible set" $\mathcal{F}$ is the set of points $(s_1, s_2)\in\mathbb{R}_+^2$ where the $(s_1, s_2)$-coordinate of any equilibrium point must be located. We have the following result.

Lemma 1. If $E = (s_1, x_1, s_2, x_2)$ is an equilibrium point of (2.5), then $0 < s_1\leq s_1^{in}$, $0 < s_2\leq s_1^{in}+s_2^{in}-s_1$ and

Proof. The equilibria of the system are the solutions of the following set of equations obtained by setting the right-hand sides of the equations in (2.5) equal to zero

At an equilibrium point we necessarily have $s_1 > 0$ and $s_2 > 0$. Using (3.2)+(3.3) and (3.4)+(3.5)-(3.3), we obtain the equations

By solving these equations, we obtain $x_1$ and $x_2$ as functions of $s_1$ and $s_2$, as given by (3.1). Consequently, the requirement that the components $x_i$ are non-negative imposes that the coordinates $(s_1, s_2)$ of an equilibrium point must satisfy $s_1\leq s_1^{in}$ and $s_1+s_2\leq s_1^{in}+s_2^{in}$. □

Therefore, the feasible set is given by

The boundary of $\mathcal{F}$ consists of two portions of the coordinate axes, together with two curves, namely the feasible set boundary ${\rm FSB}_i$ for population $i$, defined as follows:

These line segments are plotted in green and red, respectively, in Figures 2 and 3.

Remark 3. The formulas (3.1) show that the $x_1$ and $x_2$ components of an equilibrium point are uniquely determined by its $s_1$ and $s_2$ components. This is why we can, without risk of confusion, denote by the same notation an equilibrium point and its projection $(s_1, s_2)$ on the feasible set, as we do in Figures 2 and 3. We will use this abuse of notation in all the figures.

3.1.2.   The zero net growth isocline (ZNGI) and the equilibria

Then we construct the ZNGI for each population and plot them in $(s_1, s_2)$-plane using the same color used for population $i$ as used for its ${\rm FSB}_i$. The ZNGI for population $i$, denoted ${\rm ZNGI}_i$, is the curve of substrate concentrations along which the decline in biomass density is balanced by its growth. Therefore,

The following step is to identify equilibria: Each intersection of a green curve and a red curve in the feasible set corresponds to an equilibrium point, i.e., is the projection on the feasible set of an equilibrium point, as depicted in Figures 2 and 3. More precisely, we have the following result.

Lemma 2. ● The intersection point $\left(s_1^{in}, s_2^{in}\right) = {\rm FSB}_1\cap{\rm FSB}_2$ corresponds to the washout equilibrium $E_0 = \left(s_1^{in}, 0, s_2^{in}, 0\right)$ where both species are extinct.

● Any point $\left(\bar s_1, \bar s_2\right)\in{\rm ZNGI}_1\cap{\rm FSB}_2$ corresponds to a boundary equilibrium $E_1 = \left(\bar s_1, \bar x_1, \bar s_2, 0\right)$, where species 2 is absent and species 1 is present.

● Any point $\left(\tilde{s}_1, \tilde{s}_2\right)\in{\rm FSB}_1\cap{\rm ZNGI}_2$ corresponds to a boundary equilibrium $E_2 = \left(\tilde{s}_1, 0, \tilde{s}_2, \tilde{x}_2\right)$, where species 1 is absent and species 2 is present.

● Any point $\left(s_1^*, s_2^*\right)\in{\rm ZNGI}_1\cap{\rm ZNGI}_2$ lying in the interior $\mathcal{F}^{\mathrm{o}}$ of the feasible set $\mathcal{F}$ corresponds to a coexistence equilibrium $E_c = \left(s_1^*, x_1^*, s_2^*, x_2^*\right)$, where both species are present.

Proof. Recall that an equilibrium point $E = (s_1, x_1, s_2, x_2)$ is a solution of Eqs (3.2)–(3.5). Four cases must be distinguished:

1) Assume that $x_1 = x_2 = 0$, which corresponds to the washout equilibrium point $E_0$. Then, from (3.2) and $x_1 = 0$ it is deduced that $s_1 = s_1^{in}$ and from (3.4) and $x_2 = 0$, it is deduced that $s_2 = s_2^{in}$. Hence $(s_1, s_2) = \left(s_1^{in}, s_2^{in}\right) = {\rm FSB}_1\cap{\rm FSB}_2$.

2) Assume that $x_1 > 0$ and $x_2 = 0$, which corresponds to a boundary equilibrium point $E_1$. Then, from (3.3) and $x_1 > 0$ it is deduced that $f_1(s_1, s_2) = D_1$. Therefore, $(s_1, s_2)\in {\rm ZNGI}_1$. Now using (3.2)+(3.4), we obtain

From $x_2 = 0$ and (3.8), we deduce that $s_1+s_2 = s_1^{ in}+s_2^{ in}$. Now using $x_1 > 0$, from (3.1) we deduce that $s_1 < s_1^{in}$. Therefore, $(s_1, s_2)\in {\rm FSB}_2$. Hence, $(s_1, s_2)\in {\rm ZNGI}_1\cap{\rm FSB}_2$.

3) Assume that $x_1 = 0$ and $x_2 > 0$, which corresponds to a boundary equilibrium point $E_2$. Then, from $x_1 = 0$ and (3.2), we deduce that $s_1 = s_1^{ in}$. Now using $x_2 > 0$, from (3.1) we deduce that $s_2 < s_2^{in}$. Therefore, $(s_1, s_2)\in {\rm FSB}_1$. From $x_2 > 0$ and (3.5), we have $f_2(s_1, s_2) = D_2$. Therefore, $(s_1, s_2)\in{\rm ZNGI}_2$. Hence, $(s_1, s_2)\in{\rm FSB}_1 \cap{\rm ZNGI}_2$.

4) If $x_1 > 0$ and $x_2 > 0$, which corresponds to a coexistence equilibrium point $E_c$, then, from $x_1 > 0$ and (3.3), we have $f_1(s_1, s_2) = D_1$ and from $x_2 > 0$ and (3.5), we have $f_2(s_1, s_2) = D_2$. Therefore, $(s_1, s_2)\in {\rm ZNGI}_1\cap{\rm ZNGI}_2$. The equilibrium point $E_c$ exists if and only if $x_1$ and $x_2$ are positive, i.e., $(s_1, s_2)\in \mathcal{F}^{\mathrm{o}}$. □

3.1.3.   Existence and stability of equilibria

The last step is to determine the conditions of existence and stability of equilibria, which result from their location with respect to the ZNGIs. A stable equilibrium is plotted with a solid circle and an unstable one with a circle, as depicted in Figures 2 and 3. This convention will be used in all figures.

The washout equilibrium $E_0$ is unique and always exists. The number and nature of the boundary and positive equilibria depend on the relative positions of the ZNGIs. The stability of $E_0$, as well as the existence, uniqueness, and stability of the other equilibrium points, require the use of assumptions H1 and H2 on the growth functions $\mu_1$ and $\mu_2$. Since the growth functions $f_1$ and $f_2$ are given by (2.6), they satisfy the following conditions, where we use the notations $f_{ij} = {\partial f_i}/{\partial s_j}$, $i, j = 1, 2$, for the partial derivatives.

In order to describe the ZNGIs, it is convenient to use the concept of break-even concentration.

Definition 1. Let

For $s_2\geq 0$ and $D\in [0, m_1(s_2))$, the break-even concentration is the unique solution $s_1 = \lambda_1(s_2, D)$ of equation $f_1(s_1, s_2) = D$.

Let

For $s_1\geq 0$ and $D\in [0, m_2(s_1))$, the break-even concentration is the unique solution $s_2 = \lambda_2(s_1, D)$ of equation $f_2(s_1, s_2) = D$.

Notice that the existence and uniqueness of $s_1 = \lambda_1(s_2, D)$ follows from (3.9), since $\partial f_1/\partial s_1 > 0$ for all $s_1 > 0$. Using the implicit function theorem, one sees that ${\partial\lambda_1}/{\partial s_2} = -f_{12}/f_{11}\geq 0$. Therefore, $\lambda_1(s_2, D)$ is increasing in $s_2$. Similarly, $s_2 = \lambda_2(s_1, D)$ is well defined and satisfies ${\partial\lambda_2}/{\partial s_1} = -f_{21}/f_{22}\geq 0$. Therefore, it is increasing in $s_1$.

Using the break-even concentrations, the ZNGIs (3.7) are given by

Since ${\rm ZNGI}_1$ is the graph of the increasing function $s_2\mapsto s_1 = \lambda_1(s_2, D_1)$, it divides the positive cone $\mathbb{R}_+^2$ into two connected regions

The region ${R_1^-}$ contains the origin, and the region ${R_1^+}$ is possibly empty. Similarly, since ${\rm ZNGI}_2$ is the graph of the increasing function $s_1\mapsto s_2 = \lambda_2(s_1, D_2)$, it divides the positive cone $\mathbb{R}_+^2$ into two connected regions

The region ${R_2^-}$ contains the origin, and the region ${R_2^+}$ is possibly empty.

We say that an equilibrium is stable if it is locally exponentially stable, i.e., the Jacobian matrix has eigenvalues with strictly negative real parts. We can now give the main result.

Proposition 3. Assume that (3.9) and (3.10) are satisfied. The conditions of existence and stability of the equilibria of (2.5) are given in Table 1.

Proof. The proof is given in Appendix C.1. □

In Table 1, $({\rm ZNGI}_1, {\rm ZNGI}_2)$ is the signed angle between ${\rm ZNGI}_1$ and ${\rm ZNGI}_2$ at their intersection point $E_c$. As it is usual, we define the signed angle between two intersecting curves to be the signed angle between the tangents at the point of intersection.

We have thus obtained a complete description of the existence and stability conditions of the equilibria of (2.5), simply by considering the location of their projections on the feasible set $\mathcal{F}$. Existence of the equilibria is determined from the intersections of the ZNGIs and FSBs plotted with different colors. Stability of the equilibria is determined from the location of the projections of the equilibria $E_0$, $E_1$, $E_2$, and $E_c$, in the regions ${R_1^-}$, ${R_1^+}$, ${R_2^+}$ and ${R_2^-}$ of the feasible set $\mathcal{F}$. Let us illustrate these results in the cases shown in Figures 2 and 3.

● On these figures, $E_0$ is unstable since $E_0\notin{R_1^-}\cap{R_2^-}$ and $E_1$ exists since $E_0\in{R_1^+}$. On Figures 2 and 3 (Case S1), $E_1$ is unstable since $E_1\notin{R_2^-}$. On Figure 3 (Case S2), $E_1$ is stable since $E_1\in{R_2^-}$.

● On Figures 2 and 3 (Case S1), $E_2$ exists since $E_0\in{R_2^+}$ and is unstable since $E_2\notin{R_1^-}$. On Figure 3 (Case S2), $E_2$ does not exist since $E_0\notin{R_2^+}$.

● On Figures 2 and 3 (Case S1), $E_c$ is unique and stable since, at $E_c$, $({\rm ZNGI}_1, {\rm ZNGI}_2) < 0$.

● On Figure 3 (Case S2), $E_c^1$ is stable since, at $E_c^1$, $({\rm ZNGI}_1, {\rm ZNGI}_2) < 0$, and $E_c^2$ is unstable since, at $E_c^2$, $({\rm ZNGI}_1, {\rm ZNGI}_2) > 0$.

3.2.   Operating diagram

In this section, we describe the operating diagram of (2.5). Since system (2.5) has three operating parameters, and it is not easy to visualize regions in the three-dimensional operating parameter space, we fix the dilution rate $D$ and we show the operating diagram in the operating plane $\left(s_1^{in}, s_2^{in}\right)$. The effects of $D$ can be shown in a series of operating diagrams. In Section 5, we fix the operating parameter $s_2^{in}$ and we show the operating diagram in the operating plane $\left(s_1^{in}, D\right)$. The effects of $s_2^{in}$ can be shown in a series of operating diagrams.

We fix the growth functions $f_i$ and the parameters $\alpha_i$ and $a_i$, $i = 1, 2$. We begin with the more simple Cases C1, C2, and S1. We see in Figure 2 and Figure 3 (Case S1) that the ZNGIs have a unique intersection point $\left(s_1^*(D), s_2^*(D)\right): = {\rm ZNGI}_1\cap{\rm ZNGI}_2$.

3.2.1.   Commensalism (Cases C1 and C2)

We consider the curves

Even though the $\Gamma_i$ are defined by the same equations as the ZNGIs, see (3.7), it should be noted that the $\Gamma_i$ are curves in the plane of operating parameters $\left(s_1^{in}, s_2^{in}\right)$, while the ${\rm ZNGI}_i$ are curves in the phase plane $(s_1, s_2)$. The curves $\Gamma_1$ and $\Gamma_2$, defined by (3.13), together with the curve $\Gamma_3$ given by

divide the set of operating parameters $\left(s_1^{in}, s_2^{in}\right)$ into five regions denoted $\mathcal{I}_k$, $k = 0, \ldots, 4$, see Figure 4. For the C1 case, we have:

● $\Gamma_1$ is the vertical line with equation $s_1^{in} = \lambda_1(D_1)$,

● $\Gamma_2$ is the horizontal line with equation $s_2^{in} = \lambda_2(D_2)$,

● $\Gamma_3$ is the oblique line with equation $s_1^{in}+s_2^{in} = \lambda_1(D_1)+\lambda_2(D_2)$.

The novelty when we consider the C2 case is that $\Gamma_2$ becomes the curve with equation $s_2^{in} = \lambda_2\left(s_1^{in}, D_2\right)$. We have the following result:

Proposition 4. The condition of existence and stability of the equilibria of the C1 and C2 models in the regions $\mathcal{I}_k$, $k = 0, \ldots, 4$, of Figure 4 (Case C1 or C2) are given in Table 2.

Proof. The proof consists simply of looking at the location of the equilibria in the feasible set and applying Proposition 3. We give the details for the Case C1. Assume that $\left(s_1^{in}, s_2^{in}\right)\in \mathcal{I}_4$. We see in Figure 5 that $E_0$ is unstable since $E_0\notin{R_1^-}\cap{R_2^-}$, and $E_1$ exists since $E_0\in{R_1^+}$ and is unstable since $E_1\notin{R_2^-}$. Moreover, $E_2$ exists since $E_0\in{R_2^+}$ and is unstable since $E_2\notin{R_1^-}$. On the other hand, $E_c$ exists, is unique, and, stable since, at $E_c$, we have $({\rm ZNGI}_1, {\rm ZNGI}_2) < 0$. The proof for global asymptotic stability is given in Appendix A.1. This proves the results depicted in the last row of Table 2. The proofs for the other regions are illustrated in Figure 5. The details of the proof for the C2 case are given in Appendix C.2. □

Although the regions $\mathcal{I}_3$ and $\mathcal{I}_4$ are different in terms of the existence of equilibria, they are colored green because they correspond to the stability of the coexistence equilibrium $E_c$.

3.2.2.   Pure syntrophy (Case S1)

In addition to the curves $\Gamma_1$ and $\Gamma_2$, defined by (3.13), and the curve $\Gamma_3$, defined by (3.14), consider the line segment defined by

The novelty when we consider model S1 is that the curve $\Gamma_1$ becomes the curve with equation $s_1^{in} = \lambda_1\left(s_2^{in}, D_1\right)$, and is therefore distinct from $\Gamma_4$, which is the vertical line with equation $s_1^{in} = s_1^*(D)$. Therefore, in the S1 case, the curves $\Gamma_i$, $i = 1, \ldots, 4$, defined by (3.13)–(3.15), divide the set of operating parameters $\left(s_1^{in}, s_2^{in}\right)$ into six regions denoted $\mathcal{I}_k$, $k = 0, \ldots, 5$, see Figure 4 (Case S1). We have the following result:

Proposition 5. The conditions of existence and stability of the equilibria of the S1 model in the regions $\mathcal{I}_k$, $k = 0, \ldots, 5$, of Figure 4 (Case S1) are given in Table 3.

Proof. The proof is given in Appendix C.3. □

Although the regions $\mathcal{I}_3$, $\mathcal{I}_4$, and $\mathcal{I}_5$ are different in terms of the existence of equilibria, they are colored green because all three correspond to the stability of the coexistence equilibrium $E_c$.

The number of regions in the operating plane $\left(s_1^{in}, s_2^{in}\right)$ depends on $D$. More precisely, let us define

If $0 < D < \min(\delta_1, \delta_2)$, then all regions $\mathcal{I}_k$, $k = 0, \ldots, 4$, appear as shown in Figure 4. If $\delta_1 < \delta_2$ and $\delta_1\leq D < \delta_2$, only $\mathcal{I}_0$ and $\mathcal{I}_2$ appear. If $\delta_2 < \delta_1$ and $\delta_2\leq D < \delta_1$, only $\mathcal{I}_0$ and $\mathcal{I}_1$ appear. Finally, if $D\geq\max(\delta_1, \delta_2)$, then the region $\mathcal{I}_0$ invades the whole plane, see Figure 6.

3.2.3.   The S2 case

To simplify the analysis, we assume that the ZNGIs intersect in two points, as in Figure 3 (Case S2). Cases with more than two intersections can be studied in the same way. More precisely, we make the following assumption.

$\bf H3.$ There exists $\delta_0 > 0$ such that, for all $D\in(0, \delta_0)$, the ZNGIs have two distinct intersection points $\left(s_1^{*1}(D), s_2^{*1}(D)\right)$ and $\left(s_1^{*2}(D), s_2^{*2}(D)\right)$; for $D = \delta_0$ the ZNGIs are tangent, and for $D > \delta_0$ there is no intersection.

In this case, the line segments $\Gamma_3$ and $\Gamma_4$ defined in (3.14) and (3.15), respectively, must be defined for each intersection point of the ZNGIs. More precisely, we consider the curves

The curves $\Gamma_1$ and $\Gamma_2$ defined by (3.13) and $\Gamma_3^j$, $\Gamma_4^j$, $j = 1, 2$, defined by (3.17) divide the set of operating parameters $\left(s_1^{*2}(D), s_2^{*2}(D)\right)$ in nine regions denoted $\mathcal{I}_k$, $k = 0, \ldots, 8$, see Figure 7. According to the value of $D$, some of the regions can be empty, as shown in the typical example studied in Figure 9. These regions are corresponding to different system behaviors, as stated in the following result.

Proposition 6. The conditions of existence and stability of the equilibria of (2.5) in the regions $\mathcal{I}_k$, $k = 0, \ldots, 8$, shown in Figure 7 are given in the table shown in this figure.

Proof. The proof consists simply in looking at the location of the equilibria in the feasible set and applying Proposition 3. Assume that $\left(s_1^{*2}(D), s_2^{*2}(D)\right)\in \mathcal{I}_8$. We see in Figure 8 that $E_0$ is unstable since $E_0\notin{R_1^-}\cap{R_2^-}$, and $E_1$ does not exist since $E_0\notin{R_1^+}$. Moreover, $E_2$ exists since $E_0\in{R_2^+}$ and is stable since $E_2\in{R_1^-}$. On the other hand, $E_c^1$ is stable since, at $E_c^1$, we have $({\rm ZNGI}_1, {\rm ZNGI}_2) < 0$ and $E_c^2$ is unstable since, at $E_c^2$, we have $({\rm ZNGI}_1, {\rm ZNGI}_2) > 0$. This proves the results depicted in the last row of the table in Figure 7. Similarly, if $\left(s_1^{*2}(D), s_2^{*2}(D)\right)\in \mathcal{I}_7$, we see in Figure 8 that $E_0$ is stable since $E_0\in{R_1^-}\cap{R_2^-}$, $E_1$ does not exist since $E_0\notin{R_1^+}$, and $E_2$ does not exist since $E_0\notin{R_2^+}$. On the other hand, $E_c^1$ is stable since, at $E_c^1$, we have $({\rm ZNGI}_1, {\rm ZNGI}_2) < 0$, and $E_c^2$ is unstable since, at $E_c^2$, we have $({\rm ZNGI}_1, {\rm ZNGI}_2) > 0$. This proves the results described in the penultimate row of the table in Figure 7. The proofs for the other regions are shown in Figure 8. □

Although the regions $\mathcal{I}_3$, $\mathcal{I}_4$, and $\mathcal{I}_5$ are different in terms of the existence of equilibria, they are colored green because all three correspond to the stability of the coexistence equilibrium $E_c^1$.

We now illustrate how the operating diagram behaves when $D$ is varied.

To do this, we consider the growth functions defined by (2.3), and the values of the biological parameters used are listed in Table 4. The operating diagrams corresponding to various values of $D$ are shown in Figure 9. This figure shows that, as $D$ is increased, the green regions $\mathcal{I}_3$, $\mathcal{I}_4$, and $\mathcal{I}_5$, corresponding to the stability of $E_c^1$, shrink until they disappear completely at $D = \delta_0$. At this value of $D$, a saddle-node bifurcation occurs, in which $E_c^1$ and $E_c^2$ collide and annihilate each other. For $D > \delta_0$, the red region $\mathcal{I}_0$, corresponding to the stability of $E_0$, increases in size until it invades the entire operating parameter plane.

4.   Self-inhibition

4.1.   Classification of models with self inhibition

In Figure 1, we considered only the case without self inhibition of a species $i$ by its limiting substrate $S_i$, as shown by the examples (2.3) and (2.4), or the general assumptions H1 and H2. In this section, our aim is to show that our graphical method applies without added difficulty to the case where self-inhibition is possible. For example, to represent self-inhibition, we can replace the Monod functions in (2.3) with Haldane functions and obtain

or in (2.4) and obtain

See also the Kreikenbohm and Bohl function KB2 in Table A7.

In Figure 10, we illustrate the commensalistic systems with one self inhibition, the syntrophic systems with one self inhibition, and the systems with two self inhibitions. Some of the twelve systems shown in Figure 10 were considered in the literature, see Tables A5 and A6 in Appendix D.

We will not try to propose a general definition of inhibition, as such a definition would probably not cover all the cases of interest for applications. We will not make a full description of every case in Figure 10 as we have done when there is no self-inhibition. Our aim is above all to show that our qualitative graphical method is applicable in each of these cases, and we will illustrate this in Cases ${\rm C1}_2$ and ${\rm S1}_2$ which have been particularly studied in the literature, see Tables A5 and A6.

4.2.   Commensalism with inhibition of the second species by its limiting substrate

We consider Case ${\rm C1}_2$ in Figure 10. The model takes the form

We assume that $f_1$ satisfies the following condition:

We assume that $f_2$ satisfies the following condition:

We now define the break-even concentrations $\lambda_1(D)$, $\lambda_2^1(D)$, and $\lambda_2^2(D)$ of (4.1).

Definition 2. Let

For $D\in [0, m_1)$, the break-even concentration is the unique solution $s_1 = \lambda_1(D)$ of equation $f_1(s_1) = D$.

Let

For $D\in [0, m_2)$, the break-even concentration are the solutions $\lambda_2^1(D)$ and $\lambda_2^2(D)$ of equation $f_2(s_2) = D$, such that $0 < \lambda_2^1(D) < s_2^m < \lambda_2^2(D)\leq +\infty$.

Note that for $D$ close to $m_2$, we necessarily have two solutions, i.e.,

However, for $D$ small enough, there may be only one solution, $\lambda_2^1(D)$, and the second solution $\lambda_2^2(D)$ does not exist. The ZNGIs of (4.1) are given by

and

Note that, in this case, the ${R_2^{-}}$ region has two connected components: ${R_2^{-}} = {R_{21}^{-}}\cup {R_{22}^{-}}$, where

Moreover, we have

What is new compared to the C1 case (see Figure 2, Case C1) is that, since ${\rm ZNGI}_2$ has two components, we have a multiple crossing of ${\rm ZNGI}_2$ with ${\rm ZNGI}_1$ and ${\rm ZNGI}_2$ with ${\rm FSB}_1$. The system can have up to six equilibria, see Figure 11. The existence and stability conditions of these equilibria are summarized in Table 5.

To construct the operating diagram, we consider

● the vertical line $\Gamma_1$ of equation $s_1^{in} = \lambda_1(D_1)$,

● the horizontal lines $\Gamma_2^j$ of equations $s_2^{in} = \lambda_2^j(D_2)$, $j = 1, 2$,

● the oblique lines $\Gamma_3^j$ of equations $s_1^{in}+s_2^{in} = \lambda_1(D_1)+\lambda_2^j(D_2)$, $j = 1, 2$.

These lines divide the set of operating parameters $\left(s_1^{in}, s_2^{in}\right)$ in nine regions denoted $\mathcal{I}_k$, $k = 0, \ldots, 8$, as depicted in Figure 12 (Case ${\rm C1}_2$).

Proposition 7. The conditions of existence and stability of the equilibria of (4.1) in the regions $\mathcal{I}_k$ of Figure 12 (Case ${\rm C1}_2$) are given in Table 6.

Proof. The proof is given in Figure 13. Assume that $\left(s_1^{in}, s_2^{in}\right)\in \mathcal{I}_8$. We see in Figure 13 that $E_0$ is unstable since $E_0\notin{R_1^-}\cap{R_2^-}$, $E_1$ exists since $E_0\in{R_1^+}$ and is stable since $E_1\in{R_2^-}$, $E_2^1$ and $E_2^2$ exist since $E_0\in{R_{22}^{-}}$, and $E_2^1$ is unstable since $E_2^1\notin{R_1^-}$. Moreover, $E_c^1$ and $E_c^2$ exist, $E_c^1$ is stable, and $E_c^2$ is unstable. The proofs for all other regions are similar. The proof for global asymptotic stability is given in Appendix A.1. □

For a more detailed study of model ${\rm C1}_2$ and information on how the operating diagram changes when the biological parameters are changed, as well as operating diagrams in the operating plane $\left(s_1^{in}, D\right)$, where $s_2^{in}$ is kept constant, we refer the reader to ^{[31, 39]}.

4.3.   Syntrophy with inhibition of the second species by its limiting substrate

We consider the Case ${\rm S1}_2$ in Figure 10. The model takes the form

We assume that $f_1(s_1, s_2)$ and $f_2(s_2)$ satisfy (3.9) and (4.3), respectively.

The break-even concentrations $\lambda_1(s_2, D)$, $\lambda_2^1(D)$, and $\lambda_2^2(D)$ of (4.5) are defined by Definitions 1 and 2, respectively. The ${\rm ZNGI}_1$ of (4.1) is given by ${\rm ZNGI}_1 = \left\{(s_1, s_2)\in\mathbb{R}_+^2:s_1 = \lambda_1(s_2, D_1)\right\}$, while ${\rm ZNGI}_2$ is given by (4.4). As for the commensalistic case, the system can have up to six equilibria, see Figure 11, whose conditions of existence and stability are also given by Table 5.

What is new compared to the ${\rm C1}_2$ case is that the vertical line $\Gamma_1$ of the commensalistic model becomes now the curve with equation $s_1^{in} = \lambda_1\left(s_2^{in}, D_1\right)$. The horizontal lines $\Gamma_2^j$ are not changed, and the oblique lines $\Gamma_3^j$ now have the equations $s_1^{in}+s_2^{in} = \lambda_1\left(\lambda_2^j(D_2), D_1\right)+\lambda_2^j(D_2)$, $j = 1, 2$. We must also consider the vertical lines $\Gamma_4^j$ with equations $s_1 = \lambda_1\left(\lambda_2^j(D_2), D_1\right)$, $j = 1, 2$. All these lines divide the set of operating parameters $\left(s_1^{in}, s_2^{in}\right)$ into twelve regions denoted $\mathcal{I}_k$, $k = 0, \ldots, 11$, depicted in Figure 12 (Case ${\rm S1}_2$).

Proposition 8. The existence and stability of the equilibria of (4.5) in the regions $\mathcal{I}_k$, $k = 0, \ldots, 11$ of Figure 12 (Case ${\rm S1}_2$) are given in Table 7.

Proof. The proof for the regions $\mathcal{I}_k$, $k = 0, \ldots, 8$ is the same as the proof of Proposition 7 and is given in Appendix C.1. The proof for the new regions $\mathcal{I}_k$, $k = 9, 10$, and 11, which do not exist in the commensalistic case, are given in Figure 14. Assume that $\left(s_1^{in}, s_2^{in}\right)\in \mathcal{I}_{11}$. We see in Figure 14 that $E_0$ is stable since $E_0\in{R_1^-}\cap{R_2^-}$, and $E_1$ does exist since $E_0\notin{R_1^+}$. Moreover, $E_2^1$ and $E_2^2$ exist since $E_0\in{R_{22}^{-}}$, and $E_2^1$ is unstable since $E_2^1\notin{R_1^-}$. On the other hand, $E_c^1$ and $E_c^2$ exist, $E_c^1$ is stable, and $E_c^2$ is unstable. The proofs for the other regions is similar. □

The most important novelty on the asymptotic behavior of ${\rm S1}_2$ compared to ${\rm C1}_2$ is the appearance of the white regions $\mathcal{I}_{10}$ and $\mathcal{I}_{11}$ of bistability of the washout equilibrium $E_0$ and the coexistence equilibrium $E_c^1$. For a more detailed study of model ${\rm S1}_2$ and information on how the operating diagram changes when the biological parameters are changed, as well as operating diagrams in the operating plane $\left(s_1^{in}, D\right)$, where $s_2^{in}$ is kept constant, we refer the reader to ^[29].

5.   Operating diagrams where $s_2^{in}$ is kept constant

In the previous section, we determined the operating diagram when the operating parameter $D$ is kept constant and the operating parameters $s_1^{in}$ and $s_2^{in}$ vary. For the results to be useful in practice, it is necessary to determine the operating diagram when $D$ is one of the operating parameters that can vary, as this operating parameter is the most commonly used control parameter in laboratories. Since $s_1^{in}$ is the inlet concentration of the substrate of the overall system, a useful representation is to keep constant the inlet concentration $s_2^{in}$ of the substrate $s_2$ produced by the first reaction and to describe the operating diagram in the operating plane $\left(s_1^{in}, D\right)$. The effects of $s_2^{in}$ are shown in a series of operating diagrams.

The construction of this diagram can be easily deduced from the construction of the operating diagram in the operating plane $\left(s_1^{in}, s_2^{in}\right)$, and $D$ is fixed. The study carried out in the plane $\left(s_1^{in}, s_2^{in}\right)$ showed the existence of regions such that the system behaves in a certain way when the operating parameters are chosen in these regions. It is then sufficient to see how the boundaries of these regions, i.e., the $\Gamma_i$ and $\Gamma_i^j$ curves, are written in the operating plane $\left(s_1^{in}, D\right)$ and which regions of this plane they delimit. In the following sections we illustrate this construction in models C1, C2, and S1, as well as in a typical example of model S2.

5.1.   Commensalism and pure syntrophy

Figure 15 shows the regions $\mathcal{I}_k$, previously identified in Figure 4, in the operating plane $(s_1^{in}, D)$, where $s_2^{in}$ is kept constant. The number of regions depends on $s_2^{in}$. Let us describe the operating diagram in the case $\delta_1 > \delta_2$, where $\delta_1$ and $\delta_2$ are defined by (3.16). The case $\delta_1\leq \delta_2$ is similar and is left to the reader.

Note that for C1, $\Gamma_1$ is the curve with equation $D = \left(f_1\left(s_1^{in}\right)-a_1\right)/\alpha_1$ and $\Gamma_2$ is the horizontal line with equation $D = \left(f_2\left(s_2^{in}\right)-a_2\right)/\alpha_2$, which is not empty if and only if $s_2^{in} > \lambda_2(a_2)$. On the other hand, $\Gamma_3$ is the curve with equation $s_1^{in}+s_2^{in} = \lambda_1(D_1)+\lambda_2(D_2)$. The novelty when we consider model C2 is that $\Gamma_2$ becomes the curve of equation $f_2\left(s_1^{in}, s_2^{in}\right) = D_2$. This curve is not empty if and only if $f_2(0, s_2^{in}) > a_2$, which is equivalent to $s_2^{in} > \lambda_2(0, a_2)$. The novelty when we consider model S1 is that the curve $\Gamma_1$ becomes the curve with equation $s_1^{in} = \lambda_1\left(s_2^{in}, D_1\right)$, and is therefore distinct from $\Gamma_4$, which is the curve with equation $s_1^{in} = \lambda_1(\lambda_2(D_2), D_1)$. For the study of the intersection point of $\Gamma_i$ curves of the model C2, we need to define the following function: $D\in(0, \delta_2)\mapsto\xi(D): = \lambda_2(\lambda_1(D_1), D_2).$ Using (3.9) and (3.10), we see that $\xi'(D) > 0.$ Hence, $\xi$ is an increasing function from $\xi(0) = \lambda_2(\lambda_1(a_1), a_2)$ to $\xi(\delta_2) = +\infty.$ It then has an inverse function $\xi^{-1}$. The curves $\Gamma_i$ intersect at $P\left(s_2^{in}\right)$ defined by

Therefore, for C1 and $s_2^{in} > \lambda_2(a_2)$, the five regions $\mathcal{I}_k$, $k = 0, \ldots, 4$, appear as shown in Figure 15. Similarly for C2 and $s_2^{in} > \lambda_2(\lambda_1(a_1), a_2)$, the five regions $\mathcal{I}_k$ appear, while for S1 and $s_2^{in} > \lambda_2(a_2)$, the sixth region $\mathcal{I}_5$ appears, see Figure 15 (Case S1).

If $0\leq s_2^{in}\leq \lambda_2(a_2)$, for models C1 and S1, only the $\mathcal{I}_0$, $\mathcal{I}_1$, and $\mathcal{I}_3$ regions appear, as depicted in Figure 15 (Cases C1 and S1). For C2, two cases must be distinguished: If $0\leq s_2^{in}\leq \lambda_2(0, a_2)$, only the $\mathcal{I}_0$, $\mathcal{I}_1$, and $\mathcal{I}_3$ regions appear, as shown in Figure 15 (Case C2). If $\lambda_2(0, a_2) < s_2^{in}\leq \lambda_2(\lambda_1(a_1), a_2)$, a small $\mathcal{I}_2$ region also appears close to the origin, as shown in Figure 16.

5.2.   An illustrative example in the S2 case

We consider the model S2 of Table 1. The curves $\Gamma_1$, $\Gamma_2$, $\Gamma_3^i$, and $\Gamma_4^i$, $i = 1, 2$, which are the boundaries of the various regions, have been derived analytically. Indeed, since $s_2^{in}$ is kept fixed, and using (2.2), we have the following equations for these curves:

● $\Gamma_1$ is the curve with equation $D = \frac{f_1\left(s_1^{in}, s_2^{in}\right)-a_1}{\alpha_1}$.

● $\Gamma_2$ is the curve with equation $D = \frac{f_2\left(s_1^{in}, s_2^{in}\right)-a_2}{\alpha_2}$.

● $\Gamma_3^i$ and $\Gamma_4^i$, $i = 1, 2$, as defined by (3.17), are also curves in the operating plane $\left(s_1^{in}, D\right)$.

In addition to these curves we must also consider the horizontal line defined by

where $\delta_0$ is the value of $D$ for which the ZNGIs are tangent, see Assumption H3 and Figure 9. These curves divide the operating plane $(s_1^{in}, D)$ into nine regions corresponding to the nine regions depicted in Figures 7 and 9. Therefore, the operating diagrams can be drawn by plotting all these curves.

It is not possible to give a general qualitative description of the operating diagram as we did in Section 3.2.1 for models C1, C2, and S1. In Figures 17 and 18, we present the specific example given in Table 4. Increasing $s_2^{in}$ from $s_2^{in} = 0$ to $s_2^{in} = 7$, we draw the curves $\Gamma_j$, $j = 0, 1, 2$, $\Gamma_3^i$, and $\Gamma_4^i$, $i = 1, 2$, and color the $\mathcal{I}_k$ regions they delimit with the colors already used in Figure 7.

We briefly explain how these figures are obtained. For $s_2^{in} = 0$, the curves $\Gamma_2$, $\Gamma_4^1$, and $\Gamma_4^2$ are empty and only regions $\mathcal{I}_0$, $\mathcal{I}_1$, $\mathcal{I}_3$, and $\mathcal{I}_6$ appear, as shown in the panel $s_2^{in} = 0$ in Figure 17. When $s_2^{in} = 1$, the curves $\Gamma_2$ and $\Gamma_4^1$ are not empty and the regions $\mathcal{I}_2$, $\mathcal{I}_4$, and $\mathcal{I}_5$ also appear. See the panel $s_2^{in} = 1$ and its enlargement in Figure 17. If $s_2^{in}$ is increased to $s_2^*(\delta_0)$, corresponding to the tangency of the ZNGIs, then no new region appears. If $s_2^{in}$ is increased further, then the curve $\Gamma_4^2$ appears, while the curve $\Gamma_3^1$ disappears, defining two new regions $\mathcal{I}_7$ and $\mathcal{I}_8$. Therefore, the nine regions of the operating diagram in Figure 7 appear, as shown in the panel $s_2^{in} = 4$, and its enlargement, in Figure 17 and the panel $s_2^{in} = 5$, in Figure 18. If $s_2^{in}$ is increased to $s_2^{in} = 6.315$, then the nine regions continue to appear. The value $s_2^{in} = 6.315$ corresponds to the tangency of the curves $\Gamma_0$ and $\Gamma_1$, and is given by the equation $f_1(+\infty, s_2^{in}) = \alpha_1\delta_0+a_1$. If $s_2^{in}$ is increased further, then the region $\mathcal{I}_1$ disappears, see the panel $s_2^{in} = 7$ in Figure 18.

6.   Discussion and conclusions

In this work we have studied the general model (2.5) of commensalism or syntrophy. This model contains a large number of models in the existing literature. The case where $f_1$ does not depend on $s_2$, or $f_2$ does not depend on $s_1$, were studied in the literature when the removal rates $D_i$ are not equal to the dilution rate $D$, see ^{[12, 28, 29, 39, 40]}. However, the case where $f_1$ depends also on $s_2$, and $f_2$ depends also on $s_1$ were studied in the literature only when $D_1 = D_2 = D$ so that the system can be reduced to a planar system, see ^[11]. Our mathematical analysis of the model has revealed several possible behaviors: Proposition 3 provides a complete theoretical description of asymptotic behavior of the system. In order that the results can be useful in practice, one must have a description of the system with respect to the operating parameters. The study of bifurcations according to the operating parameters $D$, $s_1^{in}$, and $s_2^{in}$ is the most meaningful one for the laboratory model since the experimenter can easily vary these parameters.

Our results contain all the results from the literature as special cases, and more importantly present them in a unified way. The advantage of this unified presentation is that it also shows the similarities between the models and emphasizes the new behaviors that can emerge when new assumptions are introduced. For example:

● Extending the model from Case C1 (i.e., $f_i(s_i)$ depends only on $s_i$, $i = 1, 2$) to Case C2 (i.e., $f_2(s_1, s_2)$ is allowed to depend on $s_1$ as well) introduces no new system behaviors: the same number of equilibria, the same asymptotic behaviors. The only modification is that the regions of the operating diagram are slightly modified, see Figure 4 (Cases C1 and C2) and Table 2.

● However, extending the model from Case C1 to Case S1 (i.e., $f_1(s_1, s_2)$ is allowed to depend on $s_2$ as well) introduces new system behaviors: although the system retains the same number of equilibria, a new region appears in the operating diagram, see Figure 4 (Case S1) and Table 3.

● On the other hand, extending the model from the C2 case to the S2 case (i.e., $f_1(s_1, s_2)$ is allowed to depend on $s_2$ as well and $f_2(s_1, s_2)$ is allowed to depend on $s_1$ as well) introduces the most novelties in the system's behavior, the most important being the possibility of multiple coexistence equilibria, as well as a variety of asymptotic system behavior, including bistability, see Figure 7.

● When self-inhibition is allowed in the model, our graphical method makes it easy to compare the system's asymptotic behaviors and to highlight the contribution of syntrophy (Case ${\rm S1}_2$) versus commensalism (Case ${\rm C1}_2$). Figure 12 and Table 7 show how syntrophy brings out the possibility of bistability between the washout equilibrium $E_0$ and the coexistence equilibrium $E_c^1$ (see the white regions). This bistability does not occur in the commensal (Case ${\rm C1}_2$) model, as shown in Figure 12 and Table 6.

In the case without self-inhibition, the bistability phenomenon necessarily requires that the function $f_1$ depends on $s_2$ and the function $f_2$ depends on $s_2$. Indeed, the bistability of a positive equilibrium and a boundary equilibrium, where one of the species (or both) is extinct, is possible only for the S2 model and not for the S1, C1, or C2 models. If one of the species is unaffected by the food of its partner, the system is unable to undergo bistability.

The growth functions KB1 and KB2 of Kreikenbohm and Bohl ^{[9, 41]} do not satisfy our assumptions since they are not $\mathcal{C}^1$, but only piecewise $\mathcal{C}^1$, see Table A7. It should be interesting to extend our graphical method to those cases. Moreover, the KB1 or KB2 functions are identically 0 until a threshold value and then become positive. Such functions often appear in biological literature since the growth of a species requires that the limiting substrate exceeds a certain threshold. The extension of our method in these cases, will be considered in a future work.

As we pointed out in Remark 1, Di and Yang ^[2] also considered models of the form (2.1) with inhibition by the other species, in which, for example, $\mu_1(S_1, S_2)$ is replaced by $\mu_1(S_1, x_2)$ and $\mu_2(S_1, S_2)$ is replaced by $\mu_2(x_1, S_2)$, where $\mu_i$ is increasing in $S_i$ and decreasing in $x_j$, for $i, j = 1, 2$ ($j\neq i$). Such so-called density-dependent growth functions have been considered in the context of competition between species for a single limiting substrate, see [20, Chapter 4] and ^{[36, 37, 38, 42]}. When considering the combined effects of several interacting substrates and species density dependence, Tilman's graphical method developed in this article may prove useful in clarifying and unifying the many studies on the subject.

Use of AI tools declaration

The author declares he have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgments

The author would like to thank Radhouane Fekih-Salem and Claude Lobry for their remarks which improved the presentation considerably. The author thanks the TREASURE Euro-Mediterranean research network for partial financial support. This work was presented as part of the CIMPA summer school Digital Green: mathematical biology and theoretical ecology. The author would like to thank the three anonymous reviewers for their careful reading and constructive comments.

Conflict of interest

The author declares that he have no competing interest.

Appendix

A.   Global results

Local stability analysis only implies that the solutions starting near an asymptotically stable equilibrium converge toward this equilibrium. Hence, one cannot make assertions about the eventual outcome that are global in the sense that they are independent of the initial conditions. However, in some cases, one can reduce the four-dimensional system (2.5) to a two-dimensional system, whose global study is in general more easy. Thanks to Thieme's theory ^{[43, 44]}, we can deduce the global asymptotic stability of the initial four-dimensional system from the global asymptotic stability of the reduced two-dimensional one. For details and complements on how to use Thieme's theory, see [20, Appendix A] or [19, Appendix F]. This reduction is possible for the commensalistic models in the general case when the removal rates are not equal to the dilution rate and also in the syntrophic models when they are (i.e., $D_1 = D_2 = D$).

A.1.   Commensalism

We consider the commensalistic model ${\rm C2}_{12}$ in Figure 10. System (2.5) becomes

where $f_1$ is not assumed to be necessarily increasing in $s_1$, and similarly $f_2$ is not assumed to be necessarily increasing in $s_2$. This system contains all commensalistic models C1, C2, ${\rm C1}_1$, ${\rm C1}_2$, ${\rm C2}_1$, ${\rm C2}_2$, and ${\rm C1}_{12}$ as particular cases. For example, ${\rm C1}_2$ is obtained when $f_1$ satisfies (4.2) and $f_2$ depends only on $s_2$ and satisfies (4.3). The important thing is that system (A.1) has a cascade structure. Indeed, if $(s_1(t), x_1(t), s_2(t), x_2(t))$ is a solution of (A.1), then $(s_1(t), x_1(t))$ is a solution of the two-dimensional system

and $(s_2(t), x_2(t))$ is a solution of the non autonomous two-dimensional system

System (A.2) is a classical chemostat system. Every solution (except for a set of initial conditions of measure 0) converges toward an equilibrium $\left(s_1^*, x_1^*\right)$, possibly the washout equilibrium $\left(s_1^{in}, 0\right)$. Therefore, system (A.3) is an asymptotically autonomous system whose limiting system is

System (A.4) is also a classical chemostat system and its solutions (except for a set of initial conditions of measure 0) converge toward an equilibrium $\left(s_2^*, x_2^*\right)$. Using Thieme's theory we can conclude on the global asymptotic stability of the system (A.1). This proves the global asymptotic behavior shown in Tables 2 and 6. For details and complements on how these kinds of arguments are conducted, we refer the reader to ^{[31, 39]} for the ${\rm C1}_{2}$ case in Figure 10.

A.2.   Syntrophy with the same removal rates

In this section, we consider the case where $D_1 = D_2 = D$. System (2.5) becomes

We consider the change of variables $z_1 = s_1+x_1$ and $z_2 = s_2 -x_1+x_2$. System (A.5) becomes

Since $z_1(t)$ and $z_2(t)$ exponentially converge toward $s_1^{in}$ and $s_2^{in}$, respectively, (A.6) is an asymptotically autonomous system whose limiting system is

where $\phi_1$ and $\phi_2$ are defined by

Using Thieme's theory, the asymptotic behaviour of the solutions of the reduced model (A.7) is informative for the complete system (A.6). For details and complements on how these kinds of arguments are conducted, we refer the reader to ^[11] for the ${\rm S2}$ model when $D_1 = D_2 = D$.

B.   Bifurcation diagrams

Along $\Gamma_0$, defined by (5.1), i.e., when $D = \delta_0$, a saddle-node bifurcation, in which $E_c^1$ and $E_c^2$ collide and annihilate each other, can occur, see Figures 9, 17, and 18. The transcritical bifurcations occurring along $\Gamma_1$, $\Gamma_2$, $\Gamma_3^j$, and $\Gamma_4^j$, $j = 1, 2$, are summarized in Table A1. The result is a straightforward consequence of Proposition 3.

We illustrate bifurcations by constructing in Figure A1 the one-parameter bifurcation diagram, showing the equilibria as a function of the parameter $D$ when $s_1^{in}$ and $s_2^{in}$ are fixed. We use the growth functions given in Table 4 and take $s_1^{in} = 10$, while the value for $s_2^{in}$ is given in the figure. The behavior of the system is easily deduced from the operating diagrams shown in Figures 17 and 18. Indeed, it corresponds to the behavior on the vertical line $s_1^{in} = 10$ of these diagrams. For example, if we take $s_2^{in} = 4$, we see in Figure 17 (panel $s_2^{in} = 4$) that there exist four bifurcation values $d_0 < d_1 < d_2 < d_3$ such that:

● If $D > d_3$, then $\left(s_1^{in} = 10, D\right)\in \mathcal{I}_0$. Hence, $E_0$ is the only equilibrium and is stable.

● If $d_3 > D > d_2$, then $\left(s_1^{in} = 10, D\right)\in \mathcal{I}_1$. Hence, $E_0$ and $E_1$ are the only equilibria, $E_1$ is stable, and $E_0$ unstable.

● At $D = d_3$, a transcritical bifurcation occurs in which $E_0$ and $E_1$ collide and exchange stability.

● If $d_2 > D > d_1$, then $\left(s_1^{in} = 10, D\right)\in \mathcal{I}_6$. Hence, $E_0$, $E_1$, $E_c^1$, and $E_c^2$ are the equilibria, $E_1$ and $E_c^1$ are stable, and $E_0$ and $E_c^2$ are unstable.

● At $D = d_2$, a saddle-node bifurcation occurs in which $E_c^1$ and $E_c^2$ collide and annihilate.

● If $d_1 > D > d_0$, then $\left(s_1^{in} = 10, D\right)\in \mathcal{I}_3$. Hence, $E_0$, $E_1$, and $E_c^1$ are the only equilibria, $E_c^1$ is stable, and $E_0$ and $E_1$ are unstable.

● At $D = d_1$, a transcritical bifurcation occurs in which $E_1$ and $E_c^2$ collide and $E_1$ becomes unstable.

● If $d_0 > D > 0$, then $\left(s_1^{in} = 10, D\right)\in \mathcal{I}_4$. Hence, $E_0$, $E_1$, $E_2$, and $E_c^1$ are the equilibria, $E_c^1$ is stable, and $E_0$, $E_1$, and $E_2$ are unstable.

● At $D = d_0$, a transcritical bifurcation occurs in which $E_2$ and $E_0$ collide and $E_2$ disappears.

C.   Proofs

C.1.   Proof of Proposition 3

Recall that an equilibrium is said to be stable if it is locally exponentially stable, i.e., the Jacobian matrix has eigenvalues with strictly negative real parts. We begin by proving the following result.

Proposition 9. System (2.5) can have four types of equilibria:

1) The washout equilibrium $E_0 = (s_1^{in}, 0, s_2^{in}, 0)$, which always exists. It is stable if and only if

2) A boundary equilibrium $E_1 = (\bar s_1, \bar x_1, \bar s_2, 0)$, where $\bar{s}_1$ is a solution of the equation

and

It is unique if it exists. It exists if and only if

It is stable if and only if

3) A boundary equilibrium $E_2 = (\tilde{s}_1, 0, \tilde{s}_2, \tilde{x}_2)$, where

It exists if and only if

It is stable if and only if

4) Coexistence equilibria $E_c = \left(s_1^*, x_1^*, s_2^*, x_2^*\right)$, where $(s_1^*, s_2^*)$ is a solution of the system of equations

and

It exists if and only if $(s_1^*, s_2^*)\in\mathcal{F}^{\mathrm{o}}$. It is stable if and only if

Proof. We begin by the conditions for the existence of equilibria. Using Lemma 2, for a boundary equilibrium $E_1 = (\bar s_1, \bar x_1, \bar s_2, 0)$, we have $({\bar s}_1, {\bar s}_2) = {\rm ZNGI}_1\cap{\rm FSB}_2$. This condition is equivalent to

From the second formula in (C.12), we have $\bar s_2 = s_1^{ in}+s_2^{ in}-\bar s_1$, which is the first formula in (C.3). Replacing $\bar s_2$ in the first formula of (C.12), we have

Therefore, $\bar s_1$ is a solution of Eq (C.2). The $x_1$ component is then given by (3.1), which proves the second formula in (C.3). Eq (C.2) is equivalent to $\psi_1(s_1) = D_1$, where

We have

Using (3.9), we have $\psi_1'(s_1) > 0$ for all $s_1 > 0$. Therefore, Eq (C.2) has at most one solution. Hence, if it exists, $E_1$ is unique. $E_1$ exists if and only if equation $\psi_1(s_1) = D_1$ has a solution in the interval $\left(0, s_1^{ in}\right)$. Since $\psi_1(0) = 0$ and $\psi_1\left(s_1^{ in}\right) = f_1\left(s_1^{ in}, s_2^{ in}\right)$, the solution exists if and only if $f_1\left(s_1^{ in}, s_2^{ in}\right) > D_1$, which proves (C.4).

Using Lemma 2, for a boundary equilibrium $E_2 = (\tilde{s}_1, 0, \tilde{s}_2, \tilde{x}_2)$ we have $({\tilde s}_1, {\tilde s}_2) = {\rm FSB}_1\cap{\rm ZNGI}_2$. This condition is equivalent to

The first formula in (C.13) is the first formula in (C.6). Replacing $\tilde s_1$ in the second formula of (C.13), we have $f_2(s_1^{in}, \tilde s_2) = D_2$. Therefore, using Definition 1 we have $\tilde{s}_2 = \lambda_2\left(s_1^{in}, D_2\right)$, which proves the second formula in (C.6). The $x_2$ component is then given by (3.1), which proves the third formula in (C.6). $E_2$ exists if and only if the equation $f_2\left(s_1^{in}, s_2\right) = D_2$ has a solution in the interval $\left(0, s_2^{ in}\right)$. Since $f_2\left(s_1^{in}, 0\right) = 0$, the solution exists if and only if $f_2\left(s_1^{ in}, s_2^{ in}\right) > D_2$, which proves (C.7).

Using Lemma 2, for $E_c = (s_1^*, x_1^*, s_2^*, x_2^*)$ we have $(s_1^*, s_2^*)\in{\rm ZNGI}_1\cap{\rm ZNGI}_2\cap \mathcal{F}^{\mathrm{o}}$. This condition is equivalent to

Therefore, $(s_1^*, s_2^*)$ is a solution of (C.9), lying in the interior $\mathcal{F}^{\mathrm{o}}$ of the feasible set. The $x_1$ and $x_2$ components are then given by (3.1), which proves (C.10).

This ends the proof of the existence conditions in the proposition. The local stability of an equilibrium point of (2.5) depends on the sign of the real parts of the eigenvalues of the corresponding Jacobian matrix of system (2.5). The Jacobian matrix is the matrix of the partial derivatives of the right-hand side of system (2.5), with respect to the state variables, evaluated at the given equilibrium point $(x_1, x_2, s_1, s_2)$:

where $f_i$ and $f_{ij} = \frac{\partial f_i}{\partial s_j}$ are evaluated at the components of the equilibrium point.

At $E_0$, $x_1 = 0$ and $x_1 = 0$. Hence, the Jacobian matrix (C.14) becomes

Its eigenvalues are $f_1\left(s_1^{in}, s_2^{in}\right)-D_1$, $f_2\left(s_1^{in}, s_2^{in}\right)-D_2$ and $-D$. Therefore, $E_0$ is stable if and only if $f_1\left(s_1^{in}, s_2^{in}\right) < D_1$ and $f_2\left(s_1^{in}, s_2^{in}\right) < D_2$ which proves (C.1).

At $E_1$, $x_2 = 0$ and $x_1 > 0$ so that $f_1 = D_1$. Evaluated at $E_1$, the Jacobian matrix (C.14) becomes

Its characteristic polynomial is $P_1(\lambda) = (\lambda+D)\left(\lambda-f_2\left(\overline{s}_1, \overline{s}_2\right)+D_2\right) \left(\lambda^2+c_1\lambda+c_2\right),$ where $c_1 = D+(f_{11} - f_{12})\left(\overline{s}_1, \overline{s}_2\right) x_1$ and $c_2 = D_1(f_{11}-f_{12})\left(\overline{s}_1, \overline{s}_2\right)x_1$. The eigenvalues of $J_1$ are $-D$ and $f_2\left(\overline{s}_1, \overline{s}_2\right)-D_2$, together with the roots of the quadratic polynomial in $P_1(\lambda)$. Since $f_{11} > 0$ and $f_{12}\leq 0$, one has $c_1 > 0$ and $c_2 > 0$. Hence, the roots of the quadratic polynomial have negative real parts. Therefore, $E_1$ is stable if and only if $f_2\left(\overline{s}_1, \overline{s}_2\right) < D_2$, which proves (C.5).

At $E_2$, $x_1 = 0$ and $x_2 > 0$ so that $f_2 = D_2$. Evaluated at $E_2$, the Jacobian matrix (C.14) becomes

Its characteristic polynomial is $P_2(\lambda) = (\lambda+D)\left(\lambda-f_1\left(\tilde{s}_1, \tilde{s}_2\right)+D_1\right) \left(\lambda^2+c_1\lambda+c_2\right),$ where $c_1 = D+ f_{22}\left(\tilde{s}_1, \tilde{s}_2\right) x_2$ and $c_2 = D_2 f_{22}\left(\tilde{s}_1, \tilde{s}_2\right)x_2$. The eigenvalues of $J_2$ are $-D$ and $f_1\left(\tilde{s}_1, \tilde{s}_2\right)-D_1$, together with the roots of the quadratic polynomial in $P_2(\lambda)$. Since $f_{22} > 0$, one has $c_1 > 0$ and $c_2 > 0$. Hence, the roots of the quadratic polynomial have negative real parts. Therefore, $E_2$ is stable if and only if $f_1\left(\tilde{s}_1, \tilde{s}_2 \right) < D_1$, which proves (C.8).

At $E_c$, $x_1 > 0$ and $x_2 > 0$ so that $f_1 = D_1$ and $f_2 = D_2$. Evaluated at $E_c$, the Jacobian matrix (C.14) becomes

where $f_i$ and $f_{ij}$ are evaluated at $\left(s_1^*, s_2^*\right)$. Its characteristic polynomial is $P_c(\lambda) = \lambda^4+c_1\lambda^3+c_2\lambda^2+ c_3\lambda+c_4,$ where

The eigenvalues of $J_c$ have negative real parts if and only if the Routh-Hurwitz conditions

are satisfied. We use the following notations:

Using (3.9) and (3.10), we have $A > 0$, $C\geq 0$, and $B+C = f_{11}-f_{12} > 0$. The coefficients $c_i$ can be written as follows:

Note that $c_1 > 0$. The condition $c_4 > 0$ is equivalent to $B > 0$. If $B > 0$, then we have $c_3 > 0$. Straightforward computations show that $r_1$ can be written $r_1 = pq+r$, where

and

Hence, $r_1 > 0$ if $B > 0$. Therefore, the conditions (C.15) are satisfied if and only if $B > 0$, which is equivalent to $\left(f_{11}f_{22}-f_{12}f_{21}\right)\left(s_1^*, s_2^*\right) > 0$. This proves (C.11). □

The conditions of existence and stability of equilibria are summarized in Table A2. Let us now prove Proposition 3.

Proof of Proposition 3. The conditions $f_1\left(s_1^{in}, s_2^{in}\right) < D_1$ and $f_1\left(s_1^{in}, s_2^{in}\right) < D_2$ of stability of $E_0$ are equivalent to $E_0\in{R_1^-}\cap{R_2^-}$. The condition $f_1\left(s_1^{in}, s_2^{in}\right) > D_1$ of existence of $E_1$ is equivalent to $E_0\in{R_1^+}$. Its condition $f_2(\bar s_1, \bar s_2) < D_2$ of stability is equivalent to $E_1\in{R_2^-}$. The condition $f_2\left(s_1^{in}, s_2^{in}\right) > D_2$ of existence of $E_2$ is equivalent to $E_0\in{R_2^+}$. Its condition $f_1(\tilde s_1, \tilde s_2) < D_1$ of stability is equivalent to $E_2\in{R_1^-}$. The condition of existence of $E_c$ is $\left(s_1^*, s_2^*\right)\in\mathcal{F}^{\mathrm{o}}$. Using ${\partial\lambda_1}/{\partial s_2} = -f_{12}/f_{11}$ and ${\partial\lambda_2}/{\partial s_1} = -f_{21}/f_{22}$, the condition $(f_{11}f_{22}-f_{12}f_{21})\left(s_1^*, s_2^*\right) > 0$ of stability of $E_c$ is equivalent to

This condition means that the signed angle between between the tangent of ${\rm ZNGI}_1$ and the tangent of ${\rm ZNGI}_2$ at the point of intersection $\left(s_1^*, s_2^*\right)$ is negative. □

C.2.   Proof of Proposition 4 (Case C2)

The proof is given in Figure A2. Assume that $\left(s_1^{in}, s_2^{in}\right)\in \mathcal{I}_4$. We see in Figure A2 that $E_0$ is unstable since $E_0\notin{R_1^-}\cap{R_2^-}$, and $E_1$ exists since $E_0\in{R_1^+}$ and is unstable since $E_1\notin{R_2^-}$. Moreover, $E_2$ exists since $E_0\in{R_2^+}$ and is unstable since $E_2\notin{R_1^-}$. On the other hand, $E_c$ exists, is unique, and is stable since, at $E_c$, $({\rm ZNGI}_1, {\rm ZNGI}_2) < 0$. The proof for global asymptotic stability is given in Section 6. This proves the results depicted in the last row of Table 2. The proofs for the other regions are illustrated in Figure A2.

C.3.   Proof of Proposition 5

The proof is given in Figure A3. Assume that $\left(s_1^{in}, s_2^{in}\right)\in \mathcal{I}_5$. We see in Figure A3 that $E_0$ is unstable since $E_0\notin{R_1^-}\cap{R_2^-}$, $E_1$ does not exist since $E_0\notin{R_1^+}$, and $E_2$ exists since $E_0\in{R_2^+}$, and is unstable since $E_2\notin{R_1^-}$. Moreover, $E_c$ exists, is unique and is stable since, at $E_c$, $({\rm ZNGI}_1, {\rm ZNGI}_2) < 0$. This proves the results depicted in the last row of Table 3. The proofs for the other regions are illustrated in Figure A3. When $D_1 = D_2 = D$, the proof for global asymptotic stability is given in Appendix A.2.

C.1.   Proof of Proposition 8 ($\mathcal{I}_k$ regions, $k=0, \ldots, 8$)

The proof is given in Figure A4. Assume that $(s_1^{in}, s_2^{in})\in \mathcal{I}_8$. We see in Figure A4 that $E_0$ is unstable since $E_0\notin{R_1^-}\cap{R_2^-}$, $E_1$ exists since $E_0\in{R_1^+}$ and is stable since $E_1\in{R_2^-}$, $E_2^1$ and $E_2^2$ exist since $E_0\in{R_{22}^{-}}$, $E_2^1$ is unstable since $E_2^1\notin{R_1^-}$, and $E_c^1$ and $E_c^2$ exist. The proofs for all other regions are similar.

D.   Review of models of commensalism and syntrophy

In this section, we review the main results of the existing literature concerning the systems described in Figures 1 and 10 and briefly summarize their contributions.

Table A3 presents the main studies of the commensalistic models C1 and C2, while Table A4 presents the studies for the syntrophic models S1 and S2. The studies for the model including self-inhibition are presented in Tables A5 and A6. The growth functions used in these tables are defined in Table A7. For details and complements, the reader can consult the review papers ^{[12, 29, 45]}.

To our knowledge, Reilly ^[6] was the first to propose a mathematical study of the pure commensalistic model C1, with Monod growth functions and without decay terms of the species. He also considered more complicated commensalistic systems when feedback inhibition and feedforward activation occur to explain oscillations observed in experimental data.

The work of Stephanopoulos ^[7] is an important contribution to commensalism. He considered general growth functions and investigated the case of nonmonotone growth functions. In the case of equal removal rates, he reduced the system to a planar system and gave a complete qualitative description of the Cases C1, ${\rm C1}_1$, ${\rm C1}_2$, and ${\rm C1}_{12}$, see Figures 3 and 4 in ^[7].

The classical two-step (acidogenesis-methanisation) dynamical representation of anaerobic digestion processes, see [46, Eq (1.85)], has the structure of the pure commensalistic model C1 or its extension ${\rm C1}_2$, obtained by adding an inhibition by the second substrate. It was used by Simeonov and Stoyanov ^[40] for the control of the process who considered C1, with Monod growth functions and including decay terms.

An important contribution to the modeling of anaerobic digestion as the commensalistic system ${\rm C1}_2$ is the model of Bernard et al. ^[47] with a Monod function for $\mu_1$ and a Haldane function for $\mu_2$. Their model, sometimes referred to as AMOCO or AM2, included the term $\alpha$ which represents the fraction of the biomass in the liquid phase. Simeonov and Diop ^[48] studied this model for $\alpha = 1$. Sbarciog et al. ^[49] and Weedermann ^[50] studied the model with general growth function and $\alpha = 1$, while the interesting case where $0 < \alpha\leq 1$ was studied by Benyahia et al. ^[39]. Bayen and Gajardo ^[51] studied the steady state optimization of the biogas production of AM2, while the operating diagram of the model is described in ^[31]. For more details and complements on anaerobic digestion, we refer the reader to the recent review ^[52].

System S1 with a syntrophic relationship between the species was considered by Wilkinson ^[13] with an inhibited Monod growth function for $\mu_1$ and a Monod growth function for $\mu_2$, while Kreikenbohm and Bohl ^[9] considered another form of feedback of the second substrate on the first species, see Table A7. Burchard ^[8] extended the results of ^{[9, 13]} to a large class of general growth functions. He highlighted conditions under which there is persistence or extinction. El Hajji et al. ^[10] also obtained results for general growth functions. The studies in ^{[8, 9, 10, 13]} have shown that, if it exists, the positive coexistence equilibrium is unique and stable, and there is no other stable equilibrium.

Another model of syntrophic relationship in anaerobic digestion was considered by El Hajji ^[53].

All these studies did not include the decay terms of the species. Xu et al. ^[14] considered the effects of the decay terms and studied the model S1 with an inhibited Monod growth function for $\mu_1$ and a Monod growth function for $\mu_2$, and $S_2^{in} = 0$. These authors leaved unanswered the question of the stability of the positive steady state as long as it exists. Sari and Harmand ^[12] considered S1 with general growth functions and proved that the positive steady state is stable whenever it exists. Daoud et al. ^[28] extended the results of ^[12] to the case $S_2^{in} > 0$, showing the appearance of a new boundary equilibrium where the species $X_1$ is absent, while the species $X_2$ is present.

System ${\rm S1}_2$, which includes an inhibition of the second species by the second substrate, was considered by Harvey et al. ^[54] in the case without decay terms, where $\mu_1(S_1, S_2)$ is the product of an increasing function of $S_1$ and a decreasing function of $S_2$, which is a particular case of M-I-functions, while $\mu_2(S_2)$ is an H-function, see [54, Figure 2]. The more general case of M-I-function, also including decay terms of the species, was considered by Fekih-Salem et al. ^{[29, 55]}. In this case, the stability analysis is much more delicate since the system cannot be reduced to a planar system.

System S2, with two inhibitions, was considered by Sari et al. ^[11], who showed that, in contrast with the case of S1 with only an inhibition of the second substrate on the first species, a multiplicity of positive equilibria can occur.

An example of system ${\rm S2}_1$ with three different inhibitions was considered by Kreikenbohm and Bohl ^[41]. The mathematical analysis of this model shows the occurrence of bistability as in the case of the system S2, without the additional inhibition of the first species by its limiting substrate.

Other models for which $\mu_1$ and $\mu_2$ depend both on $(S_1, S_2)$ and, in addition, $X_1$ and $X_2$ are in competition on substrate $S_1$, exhibiting the multiplicity of positive equilibrium points, can be found in ^[56]. An important and interesting extension should be mentioned here: ^[26] proposed an 8-dimensional mathematical model, which includes syntrophy and inhibition, and both mechanisms considered by ^[47] and by ^[10]. The effects of decay terms are considered by ^[27].