1. Introduction {#sec0001}
===============

The mathematical modelling of engineered biological systems has entered a new era in recent years with the expansion and standardisation of existing models aimed at collating disparate components of these processes and provide scientists, engineers and practitioners with the tools to better predict, control and optimise them [@bib0023]. In engineered biological systems, mechanistic modelling reached consensus with the development of the Activated Sludge Models [@bib0009], [@bib0010] for wastewater treatment processes, followed by the Anaerobic Digestion Model No. 1 (ADM1) [@bib0012] a few years later. The development of ADM1 was enabled largely due to the possibilities for better identification and characterisation of functional microbial groups responsible for the chemical transformations within anaerobic digesters. It describes a set of stoichiometric and kinetic functions representing the standard anaerobic processes, remaining the scientific benchmark to the present day. However, there has been a growing awareness that the model should take advantage of improved empirical understanding and extension of biochemical processes included in its structure, to acquire a better trade-off between model realism and complexity [@bib0013].

The full ADM1 model is highly parameterised with a large number of physical, chemical and biological processes described by numerous state variables and algebraic expressions. Whilst suitable for dynamic simulation, more rigorous mathematical analysis of the model is difficult. To the authors knowledge, only numerical investigations are available [@bib0003]. Due to the analytical intractability of the full ADM1, work has been made towards the construction of simpler models that preserve biological meaning whilst reducing the computational effort required to find mathematical solutions of the model equations [@bib0007], [@bib0008].

The most common models used to describe microbial systems are two-tiered models, which take the form of a cascade of two biological reactions where one substrate is consumed by one microorganism to produce a product that serves as the main limiting substrate for a second microorganism. When the second organism has no feedback on the first organism, the system is known as commensalistic [@bib0016], [@bib0020]. The system has a cascade structure and the number of steady-states and their (mathematical) stability as a function of model inputs and parameters may be investigated [@bib0001], [@bib0002], [@bib0019]. When the growth of the first organism is affected by the substrate produced by the second organism the system is known as syntrophic. For instance, if the first organism is inhibited by high concentrations of the product, the extent to which the first substrate is degraded by the first organism depends on the efficiency of the removal of the product by the second organism. The mathematical analysis of such a model is more delicate than for commensalistic models, (see for instance early work by [@bib0004], [@bib0014], [@bib0015], [@bib0027] and the more recent papers [@bib0005], [@bib0011], [@bib0017], [@bib0018], [@bib0022]). An important and interesting extension should be mentioned here: [@bib0025], [@bib0026] analysed an 8-dimensional mathematical model, which includes syntrophy and inhibition, both mechanisms considered by [@bib0002] and [@bib0006].

As an example of this for anaerobic digestion, a previous study investigated the effect of maintenance on the stability of a two-tiered 'food-chain' comprising two species and two substrates [@bib0028]. Maintenance is defined as the energy consumed by an organism that is used for all biological processes other than growth. In [@bib0028] and here, it is analogous to a first-order decay rate constant, or biomass death term. Although the authors were not able to determine the general conditions under which this four dimensional syntrophic consortium was stable, further work has shown that a model with generality can be used to answer the question posed, determining that the two-tiered food-chain is always stable when maintenance is included [@bib0018].

More recently, the model described by [@bib0028] was extended by the addition of a third organism and substrate to create a three-tiered 'food-web' [@bib0024], as shown in [Fig. 1](#fig0001){ref-type="fig"}. In this paper, the existence and stability of the steady-states were determined only numerically. Although the results were important in revealing emergent properties of this extended model, the motivation of this work is to give an analytical study of the model. Moreover, our analysis does not require that growth functions are of the specific form considered and are valid for a large class of growth functions. This is critical as it provides the means by which microbiologists can theoretically test the influence of the growth characteristics of organisms on the properties of the system and the interactions between multiple species.Fig. 1Schematic of the three-tier chlorophenol mineralising food-web indicating the flow and conversion of chemical substrates and products in the system.Fig. 1

Here, we pursue a generalised description and analysis of the model given by [@bib0024]. Chlorophenols are chemicals of importance due to their impact on the environment and to public health, their recalcitrance in food-webs and resistance to aerobic biodegradation via the oxygenase enzyme [@bib0021]. Although we consider the monochlorophenol isomer here, extension to multiple isomeric chlorophenols would be straightforward. It is important to note that, although the particular biological transformation provided is for chlorophenol mineralisation, the structure and methods employed are much more general and apply to any theoretical ecological interactions that may be hypothesised or observed in a microbial community. We therefore stress that this work can provide a good approach for analytically investigating the behaviour of microbial food-webs where numerical parameters are difficult to obtain or uncertain. Ultimately, we demonstrate here the advantage of a generalised approach for mathematical analysis of microbial ecological interactions from both a theoretical perspective and its potential if providing knowledge in applied research where these communities and processes are studied empirically.

The paper is organised as follows. In [Section 2](#sec0002){ref-type="sec"}, we present a description of the model to be investigated, and its reduction in [Section 3](#sec0003){ref-type="sec"}. Model assumptions and notations are provided in [Section 4](#sec0004){ref-type="sec"}. In [Section 5](#sec0005){ref-type="sec"} we demonstrate the existence of the three steady-states and define four interesting cases for specific parameter values that are investigated using the solutions, whilst also indicating the regions of existence of the steady-states for the operating parameter values. We present results on the behaviour of the system whilst varying two control parameters in [Section 6](#sec0006){ref-type="sec"}. In [Section 7](#sec0011){ref-type="sec"} we perform local stability analysis of the steady-states without maintenance and, in [Section 8](#sec0016){ref-type="sec"}, we undertake a comprehensive numerical stability analysis of the four cases for both the model with and without maintenance. We show that our approach leads to the discovery of five operating regions in which one leads to the possibility of instability of the positive steady-state, where all three organisms exist, a fact that has not been reported by [@bib0024]. Indeed, we suggest that a stable limit-cycle can occur at the boundary with this region. Finally, in [Section 9](#sec0021){ref-type="sec"}, we make comment on the role of the kinetic parameters used in the four example cases, in maintaining stability, which points to the importance of the relative aptitude of the two hydrogen consumers in sustaining a viable chlorophenol mineralising community. In the Appendix we describe the numerical method used in [Section 8](#sec0016){ref-type="sec"}, give the assumptions on the growth functions we used and the proofs of the results.

2. The model {#sec0002}
============

The model developed in [@bib0024] has six components, three substrate (chlorophenol, phenol and hydrogen) and three biomass (chlorophenol and phenol degraders, and a hydrogenotrophic methanogen) variables. The substrate and biomass concentrations evolve according to the six-dimensional dynamical of ODEs: $$\begin{aligned}
\frac{dX_{\text{ch}}}{dt} & {= - DX_{\text{ch}} + Y_{\text{ch}}f_{0}\left( S_{\text{ch}},S_{H_{2}} \right)X_{\text{ch}} - k_{\text{dec},\text{ch}}X_{\text{ch}}} \\
\end{aligned}$$$$\begin{aligned}
\frac{dX_{\text{ph}}}{dt} & {= - DX_{\text{ph}} + Y_{\text{ph}}f_{1}\left( S_{\text{ph}},S_{H_{2}} \right)X_{\text{ph}} - k_{\text{dec},\text{ph}}X_{\text{ph}}} \\
\end{aligned}$$$$\begin{aligned}
\frac{dX_{H_{2}}}{dt} & {= - DX_{H_{2}} + Y_{H_{2}}f_{2}\left( S_{H_{2}} \right)X_{H_{2}} - k_{\text{dec},H_{2}}X_{H_{2}}} \\
\end{aligned}$$$$\begin{aligned}
\frac{dS_{\text{ch}}}{dt} & {= D\left( S_{\text{ch},\text{in}} - S_{\text{ch}} \right) - f_{0}\left( S_{\text{ch}},S_{H_{2}} \right)X_{\text{ch}}} \\
\end{aligned}$$$$\begin{array}{ccl}
\frac{dS_{\text{ph}}}{dt} & = & {D\left( S_{\text{ph},\text{in}} - S_{\text{ph}} \right) + \frac{224}{208}\left( 1 - Y_{\text{ch}} \right)f_{0}\left( S_{\text{ch}},S_{H_{2}} \right)X_{\text{ch}}} \\
 & & {- f_{1}\left( S_{\text{ph}},S_{H_{2}} \right)X_{\text{ph}}} \\
\end{array}$$$$\begin{array}{ccl}
\frac{dS_{H_{2}}}{dt} & = & {D\left( S_{H_{2},\text{in}} - S_{H_{2}} \right) + \frac{32}{224}\left( 1 - Y_{\text{ph}} \right)f_{1}\left( S_{\text{ph}},S_{H_{2}} \right)X_{\text{ph}}} \\
 & & {- \frac{16}{208}f_{0}\left( S_{\text{ch}},S_{H_{2}} \right)X_{\text{ch}} - f_{2}\left( S_{H_{2}} \right)X_{H_{2}}} \\
\end{array}$$where *S*~ch~ and *X*~ch~ are the chlorophenol substrate and degrader concentrations, *S*~ph~ and *X*~ph~ for phenol and $S_{H_{2}}$ and $X_{H_{2}}$ for hydrogen; *Y*~ch~, *Y*~ph~ and *Y*~H2~ are the yield coefficients, $224/208\left( 1 - Y_{\text{ch}} \right)$ represents the part of chlorophenol degraded to phenol, and $32/224\left( 1 - Y_{\text{ph}} \right)$ represents the part of phenol that is transformed to hydrogen. Growth functions take Monod form with hydrogen inhibition acting on the phenol degrader and represented in *f*~1~ (see ([Eq. 7](#eq0007){ref-type="disp-formula"})) as a product inhibition term. In [@bib0024], this model is given in dimensionless form that significantly reduces the number of parameters describing the dynamics. In the analysis of the generalised model ([Sections 5](#sec0005){ref-type="sec"} and [7](#sec0011){ref-type="sec"}), we do not assume that the growth functions *f*~0~, *f*~1~ and *f*~2~ have the specific analytical expression shown. We will only assume that the growth functions satisfy properties that are listed in [Section 4](#sec0004){ref-type="sec"}. Therefore, we cannot benefit from the dimensionless rescaling used by [@bib0024], because it uses some kinetics parameters of the specific growth functions ([Eq. 7](#eq0007){ref-type="disp-formula"}), while we work with general *unspecified* growth functions.

In [Section 3](#sec0003){ref-type="sec"} we consider another rescaling that does not use the kinetics parameters. For the numerical analyses given in [Section 8](#sec0016){ref-type="sec"}, we return to an assumed set of Monod growth functions given by ([Eq. 7](#eq0007){ref-type="disp-formula"}), in order to directly compare with the results found by [@bib0024]. $$\begin{array}{rcl}
{f_{0}\left( S_{\text{ch}},S_{H_{2}} \right)} & = & {\frac{k_{m,\text{ch}}S_{\text{ch}}}{K_{S,ch} + S_{\text{ch}}}\frac{S_{H_{2}}}{K_{S,H_{2},c} + S_{H_{2}}}} \\
{f_{1}\left( S_{\text{ph}},S_{H_{2}} \right)} & = & {\frac{k_{m,\text{ph}}S_{\text{ph}}}{K_{S,\text{ph}} + S_{\text{ph}}}\frac{1}{1 + \frac{S_{H_{2}}}{K_{i,H_{2}}}}} \\
{f_{2}\left( S_{H_{2}} \right)} & = & \frac{k_{m,H_{2}}S_{H_{2}}}{K_{S,H_{2}} + S_{H_{2}}} \\
\end{array}$$

Here, apart from the four operating (or control) parameters, which are the inflowing concentrations *S*~ch,\ in~, *S*~ph,\ in~, $S_{H_{2},\text{in}}$ and the dilution rate *D*, that can vary, all others have biological meaning and are fixed depending on the organisms and substrate considered. We use the following notations in ([(1)](#eq0001){ref-type="disp-formula"}, [(2)](#eq0002){ref-type="disp-formula"}, [(3)](#eq0003){ref-type="disp-formula"}, [(4)](#eq0004){ref-type="disp-formula"}, [(5)](#eq0005){ref-type="disp-formula"}, [(6)](#eq0006){ref-type="disp-formula"}): $$\begin{array}{cl}
 & {X_{0} = X_{\text{ch}},\quad X_{1} = X_{\text{ph}},\quad X_{2} = X_{H_{2}}} \\
 & {S_{0} = S_{\text{ch}},\quad S_{1} = S_{\text{ph}},\quad S_{2} = S_{H_{2}}} \\
 & {S_{0}^{\text{in}} = S_{\text{ch},\text{in}},\quad S_{1}^{\text{in}} = S_{\text{ph},\text{in}},\quad S_{2}^{\text{in}} = S_{H_{2},\text{in}}} \\
 & {Y_{0} = Y_{\text{ch}},\quad Y_{1} = Y_{\text{ph}},\quad Y_{2} = Y_{H_{2}}} \\
 & {Y_{3} = \frac{224}{208}\left( 1 - Y_{\text{ch}} \right),\quad Y_{4} = \frac{32}{224}\left( 1 - Y_{\text{ph}} \right),\quad Y_{5} = \frac{16}{208}} \\
 & {a_{0} = k_{\text{dec},\text{ch}},\quad a_{1} = k_{\text{dec},\text{ph}},\quad a_{2} = k_{\text{dec},H_{2}}} \\
\end{array}$$With these notations ([(1)](#eq0001){ref-type="disp-formula"}, [(2)](#eq0002){ref-type="disp-formula"}, [(3)](#eq0003){ref-type="disp-formula"}, [(4)](#eq0004){ref-type="disp-formula"}, [(5)](#eq0005){ref-type="disp-formula"}, [(6)](#eq0006){ref-type="disp-formula"}) can be written as follows: $$\begin{aligned}
\frac{dX_{0}}{dt} & {= - DX_{0} + Y_{0}f_{0}\left( S_{0},S_{2} \right)X_{0} - a_{0}X_{0}} \\
\end{aligned}$$$$\begin{aligned}
\frac{dX_{1}}{dt} & {= - DX_{1} + Y_{1}f_{1}\left( S_{1},S_{2} \right)X_{1} - a_{1}X_{1}} \\
\end{aligned}$$$$\begin{aligned}
\frac{dX_{2}}{dt} & {= - DX_{2} + Y_{2}f_{2}\left( S_{2} \right)X_{2} - a_{2}X_{2}} \\
\end{aligned}$$$$\begin{aligned}
\frac{dS_{0}}{dt} & {= D\left( S_{0}^{\text{in}} - S_{0} \right) - f_{0}\left( S_{0},S_{2} \right)X_{0}} \\
\end{aligned}$$$$\begin{aligned}
\frac{dS_{1}}{dt} & {= D\left( S_{1}^{\text{in}} - S_{1} \right) + Y_{3}f_{0}\left( S_{0},S_{2} \right)X_{0} - f_{1}\left( S_{1},S_{2} \right)X_{1}} \\
\end{aligned}$$$$\begin{aligned}
\frac{dS_{2}}{dt} & {= D\left( S_{2}^{\text{in}} - S_{2} \right) + Y_{4}f_{1}\left( S_{1},S_{2} \right)X_{1} - Y_{5}f_{0}\left( S_{0},S_{2} \right)X_{0} - f_{2}\left( S_{2} \right)X_{2}} \\
\end{aligned}$$Furthermore, we restrict our analysis to the case where we only have one substrate added to the system, such that: $S_{0}^{in} > 0,$$S_{1}^{in} = 0,$ and $S_{2}^{in} = 0$. As shown in [@bib0024], the general case with $S_{1}^{in} > 0,$ and $S_{2}^{in} > 0$ present many important and interesting behaviours and deserves future work.

3. Model reduction {#sec0003}
==================

To ease the mathematical analysis, we can rescale the system ([(8)](#eq0008){ref-type="disp-formula"}, [(9)](#eq0009){ref-type="disp-formula"}, [(10)](#eq0010){ref-type="disp-formula"}, [(11)](#eq0011){ref-type="disp-formula"}, [(12)](#eq0012){ref-type="disp-formula"}, [(13)](#eq0013){ref-type="disp-formula"}) using the following change of variables adapted from [@bib0018]: $$\begin{array}{cl}
 & {x_{0} = \frac{Y_{3}Y_{4}}{Y_{0}}X_{0},\quad x_{1} = \frac{Y_{4}}{Y_{1}}X_{1},\quad x_{2} = \frac{1}{Y_{2}}X_{2}} \\
 & {s_{0} = Y_{3}Y_{4}S_{0},\quad s_{1} = Y_{4}S_{1},\quad s_{2} = S_{2}} \\
\end{array}$$We obtain the following system: $$\begin{aligned}
\frac{dx_{0}}{dt} & {= - Dx_{0} + \mu_{0}\left( s_{0},s_{2} \right)x_{0} - a_{0}x_{0}} \\
\end{aligned}$$$$\begin{aligned}
\frac{dx_{1}}{dt} & {= - Dx_{1} + \mu_{1}\left( s_{1},s_{2} \right)x_{1} - a_{1}x_{1}} \\
\end{aligned}$$$$\begin{aligned}
\frac{dx_{2}}{dt} & {= - Dx_{2} + \mu_{2}\left( s_{2} \right)x_{2} - a_{2}x_{2}} \\
\end{aligned}$$$$\begin{aligned}
\frac{ds_{0}}{dt} & {= D\left( s_{0}^{\text{in}} - s_{0} \right) - \mu_{0}\left( s_{0},s_{2} \right)x_{0}} \\
\end{aligned}$$$$\begin{aligned}
\frac{ds_{1}}{dt} & {= - Ds_{1} + \mu_{0}\left( s_{0},s_{2} \right)x_{0} - \mu_{1}\left( s_{1},s_{2} \right)x_{1}} \\
\end{aligned}$$$$\begin{aligned}
\frac{ds_{2}}{dt} & {= - Ds_{2} + \mu_{1}\left( s_{1},s_{2} \right)x_{1} - \omega\mu_{0}\left( s_{0},s_{2} \right)x_{0} - \mu_{2}\left( s_{2} \right)x_{2}} \\
\end{aligned}$$where the inflowing concentration is: $$s_{0}^{\text{in}} = Y_{3}Y_{4}S_{0}^{\text{in}},$$the growth functions are: $$\begin{array}{ccl}
 & & {\mu_{0}\left( s_{0},s_{2} \right) = Y_{0}f_{0}\left( \frac{s_{0}}{Y_{3}Y_{4}},s_{2} \right)} \\
 & & {\mu_{1}\left( s_{1},s_{2} \right) = Y_{1}f_{1}\left( \frac{s_{1}}{Y_{4}},s_{2} \right)} \\
 & & {\mu_{2}\left( s_{2} \right) = Y_{2}f_{2}\left( s_{2} \right)} \\
\end{array}$$and $$\omega = \frac{Y_{5}}{Y_{3}Y_{4}} = \frac{1}{2\left( 1 - Y_{0} \right)\left( 1 - Y_{1} \right)}$$The benefit of our rescaling is that it permits to fix in ([(14)](#eq0014){ref-type="disp-formula"}, [(15)](#eq0015){ref-type="disp-formula"}, [(16)](#eq0016){ref-type="disp-formula"}, [(17)](#eq0017){ref-type="disp-formula"}, [(18)](#eq0018){ref-type="disp-formula"}, [(19)](#eq0019){ref-type="disp-formula"}) all yield coefficients to one except that denoted by *ω* and defined by ([Eq. 22](#eq0022){ref-type="disp-formula"}), and to discuss the existence and stability with respect to this sole parameter.

Using ([Eq. 21](#eq0021){ref-type="disp-formula"}) and the growth functions ([Eq. 7](#eq0007){ref-type="disp-formula"}), we obtain the model ([(14)](#eq0014){ref-type="disp-formula"}, [(15)](#eq0015){ref-type="disp-formula"}, [(16)](#eq0016){ref-type="disp-formula"}, [(17)](#eq0017){ref-type="disp-formula"}, [(18)](#eq0018){ref-type="disp-formula"}, [(19)](#eq0019){ref-type="disp-formula"}) with the following Monod-type growth functions: $$\begin{array}{ccl}
 & & {\mu_{0}\left( s_{0},s_{2} \right) = \frac{m_{0}s_{0}}{K_{0} + s_{0}}\frac{s_{2}}{L_{0} + s_{2}}} \\
 & & {\mu_{1}\left( s_{1},s_{2} \right) = \frac{m_{1}s_{1}}{K_{1} + s_{1}}\frac{1}{1 + s_{2}/K_{i}}} \\
 & & {\mu_{2}\left( s_{2} \right) = \frac{m_{2}s_{2}}{K_{2} + s_{2}}} \\
\end{array}$$where: $$\begin{array}{cl}
 & {m_{0} = Y_{0}k_{m,\text{ch}},\quad K_{0} = Y_{3}Y_{4}K_{s,\text{ch}},\quad L_{0} = K_{S,H_{2},c}} \\
 & {m_{1} = Y_{1}k_{m,\text{ph}},\quad K_{1} = Y_{4}K_{s,\text{ph}},\quad K_{i} = K_{i,H_{2}}} \\
 & {m_{2} = Y_{2}k_{m,H_{2}},\quad K_{2} = K_{S,H_{2}}} \\
\end{array}$$For the numerical simulations we will use the nominal values in [Table 1](#tbl0001){ref-type="table"} given in [@bib0024].Table 1Nominal parameter values. We use units expressed in Chemical Oxygen Demand (COD) as used by [@bib0012], [@bib0024].Table 1ParametersNominal valuesUnits*k*~*m*,\ ch~29kgCOD~S~/kgCOD~X~/d*K*~*S*,\ ch~0.053kgCOD/m^3^*Y*~ch~0.019kgCOD~X~/kgCOD~S~*k*~*m*,\ ph~26kgCOD~S~/kgCOD~X~/d*K*~*S*,\ ph~0.302kgCOD/m^3^*Y*~ph~0.04kgCOD~X~/kgCOD~S~$k_{m,H_{2}}$35kgCOD~S~/kgCOD~X~/d$K_{S,H_{2}}$2.5$\times 10^{- 5}$kgCOD/m^3^$K_{S,{H_{2},c}}$1.0$\times 10^{- 6}$kgCOD/m^3^$Y_{H_{2}}$0.06kgCOD~X~/kgCOD~S~*k*~dec,\ i~0.02$d^{- 1}$$K_{I,H_{2}}$3.5$\times 10^{- 6}$kgCOD/m^3^

4. Assumptions on the model and notations {#sec0004}
=========================================

Our study does not require that growth functions are of Monod type ([Eq. 23](#eq0023){ref-type="disp-formula"}). Actually, the results are valid for a more general class of growth functions satisfying the following conditions, which concur with those given by ([Eq. 23](#eq0023){ref-type="disp-formula"}): H1For all *s*~0~ \> 0 and *s*~2~ \> 0 then $0 < \mu_{0}\left( s_{0},s_{2} \right) < + \infty$ and $\mu_{0}\left( 0,s_{2} \right) = 0,$$\mu_{0}\left( s_{0},0 \right) = 0$.H2For all *s*~1~ \> 0 and *s*~2~ ≥ 0 then $0 < \mu_{1}\left( s_{1},s_{2} \right) < + \infty$ and $\mu_{1}\left( 0,s_{2} \right) = 0$.H3For all *s*~2~ \> 0 then $0 < \mu_{2}\left( s_{2} \right) < + \infty$ and $\mu_{2}\left( 0 \right) = 0$.H4For all *s*~0~ \> 0 and *s*~2~ \> 0, $$\frac{\partial\mu_{0}}{\partial s_{0}}\left( s_{0},s_{2} \right) > 0,\quad\frac{\partial\mu_{0}}{\partial s_{2}}\left( s_{0},s_{2} \right) > 0.$$H5For all *s*~1~ \> 0 and *s*~2~ \> 0, $$\frac{\partial\mu_{1}}{\partial s_{1}}\left( s_{1},s_{2} \right) > 0,\quad\frac{\partial\mu_{1}}{\partial s_{2}}\left( s_{1},s_{2} \right) < 0.$$H6For all *s*~2~ \> 0, $\frac{d\mu_{2}}{ds_{2}}\left( s_{2} \right) > 0$.H7The function $\left. s_{2}\mapsto\mu_{0}\left( + \infty,s_{2} \right) \right.$ is monotonically increasing and the function $\left. s_{2}\mapsto\mu_{1}\left( + \infty,s_{2} \right) \right.$ is monotonically decreasing.

The proof of the following result is standard and hence omitted. Proposition 1*For non-negative initial conditions, all solutions of the system (Eqs.* [14](#eq0014){ref-type="disp-formula"}*--*[19](#eq0019){ref-type="disp-formula"}*) are bounded and remain non-negative for all t* \> 0*.*

In the following lemma we define the functions *M*~0~(*y, s*~2~), *M*~1~(*y, s*~2~) and *M*~2~(*y*), which are the inverse functions of the functions *s*~0~↦*μ*~0~(*s*~0~, *s*~2~), *s*~1~↦*μ*~1~(*s*~1~, *s*~2~) and *s*~2~↦*μ*~2~(*s*~2~), respectively. Lemma 1*Let s*~2~ ≥ 0 *be fixed. There exists a unique function* $$\left. y \in \left\lbrack 0,\mu_{0}\left( + \infty,s_{2} \right) \right)\mapsto M_{0}\left( y,s_{2} \right) \in \left\lbrack 0, + \infty \right), \right.$$*such that for s*~0~ ≥ 0, *s*~2~ ≥ 0 *and*$y \in \left\lbrack 0,\mu_{0}\left( + \infty,s_{2} \right) \right),$ *we have:* $$\left. s_{0} = M_{0}\left( y,s_{2} \right)\Longleftrightarrow y = \mu_{0}\left( s_{0},s_{2} \right) \right.$$*Let s*~2~ ≥ 0 *be fixed. There exists a unique function* $$\left. y \in \left\lbrack 0,\mu_{1}\left( + \infty,s_{2} \right) \right)\mapsto M_{1}\left( y,s_{2} \right) \in \left\lbrack 0, + \infty \right), \right.$$*such that for s*~1~ ≥ 0, *s*~2~ ≥ 0 *and*$y \in \left\lbrack 0,\mu_{1}\left( + \infty,s_{2} \right) \right),$ *we have:* $$\left. s_{1} = M_{1}\left( y,s_{2} \right)\Longleftrightarrow y = \mu_{1}\left( s_{1},s_{2} \right) \right.$$*There exists a unique function* $$\left. y \in \left\lbrack 0,\mu_{2}\left( + \infty \right) \right)\mapsto M_{2}\left( y \right) \in \left\lbrack 0, + \infty \right), \right.$$*such that, for s*~2~ ≥ 0 *and*$y \in \left\lbrack 0,\mu_{2}\left( + \infty \right) \right)$ *we have:* $$\left. s_{2} = M_{2}\left( y \right)\Longleftrightarrow y = \mu_{2}\left( s_{2} \right) \right.$$

Using the functions *M*~0~, *M*~1~ and *M*~2~ we define now the function *ψ*(*s*~2~, *D*), *F*~1~(*D*), *F*~2~(*D*) and *F*~3~(*D*). Let us observe first that by **H7**, for $D + a_{0} < \mu_{0}\left( + \infty, + \infty \right)$ and $D + a_{1} < \mu_{1}\left( + \infty,0 \right)$ there exist unique values $s_{2}^{0} = s_{2}^{0}\left( D \right)$ and $s_{2}^{1} = s_{2}^{1}\left( D \right),$ see [Fig. 2](#fig0002){ref-type="fig"}(a): $$\mu_{0}\left( + \infty,s_{2}^{0} \right) = D + a_{0},\quad\mu_{1}\left( + \infty,s_{2}^{1} \right) = D + a_{1}$$Assume that *ω* \< 1. Let *D* be fixed such that $s_{2}^{0}\left( D \right)$ and $s_{2}^{1}\left( D \right)$ exist. We consider the function defined on $\left( s_{2}^{0}\left( D \right),s_{2}^{1}\left( D \right) \right)$ by: $$\begin{array}{ccl}
 & & \left. s_{2} \in \left( s_{2}^{0}\left( D \right),s_{2}^{1}\left( D \right) \right)\mapsto\psi\left( s_{2},D \right) \right. \\
 & & {\quad = M_{0}\left( D + a_{0},s_{2} \right) + \frac{M_{1}\left( D + a_{1},s_{2} \right) + s_{2}}{1 - \omega},} \\
\end{array}$$It should be noted that *ψ*(*s*~2~, *D*) \> 0 for $s_{2}^{0}\left( D \right) < s_{2} < s_{2}^{1}\left( D \right)$. The function *ψ* together with the values $s_{2}^{0},$$s_{2}^{1}$ all depend on *D*. However, to avoid cumbersome notation we will use the more precise form *ψ*(*s*~2~, *D*), $s_{2}^{0}\left( D \right)$ and $s_{2}^{1}\left( D \right)$ only if necessary. From ([Eq. 25](#eq0025){ref-type="disp-formula"}), ([Eq. 26](#eq0026){ref-type="disp-formula"}) and ([Eq. 28](#eq0028){ref-type="disp-formula"}) we deduce that: $$M_{0}\left( D + a_{0},s_{2}^{0} \right) = + \infty,\quad M_{1}\left( D + a_{1},s_{2}^{1} \right) = + \infty$$Therefore, we have, see [Fig. 2](#fig0002){ref-type="fig"}(b): $$\lim\limits_{s_{2}\rightarrow s_{2}^{0}}\psi\left( s_{2} \right) = \lim\limits_{s_{2}\rightarrow s_{2}^{1}}\psi\left( s_{2} \right) = + \infty$$Hence, the function *ψ*(*s*~2~), which is positive and tends to $+ \infty$ at the extremities of the interval $\left( s_{2}^{0},s_{2}^{1} \right),$ has a minimum value on this interval. We add the following assumption: H8The function *ψ* has a unique minimum ${\overline{s}}_{2}\left( D \right)$ on the interval $\left( s_{2}^{0}\left( D \right),s_{2}^{1}\left( D \right) \right)$ and $\frac{\partial\psi}{\partial s_{2}}\left( s_{2},D \right)$ is negative on $\left( s_{2}^{0}\left( D \right),{\overline{s}}_{2}\left( D \right) \right)$ and positive on $\left( {\overline{s}}_{2}\left( D \right),s_{2}^{1}\left( D \right) \right),$ respectively.Fig. 2Graphical definitions. (a): $s_{2}^{0}$ and $s_{2}^{1}$. (b) : *ψ*(*s*~2~), ${\overline{s}}_{2},$*F*~1~(*D*) and *F*~2~(*D*). (c): $s_{2}^{\flat}$ and $s_{2}^{\sharp}$.Fig. 2Fig. 3Graphs of $s_{2}^{0}\left( D \right)$ and $s_{2}^{1}\left( D \right)$ (in black) and $M_{2}\left( D + a_{2} \right)$ (in red) and graphical depiction of $I_{1} = \left\lbrack 0,D_{1} \right),$ where *D*~1~ is the solution of $s_{2}^{0}\left( D \right) = s_{2}^{1}\left( D \right),$ and *I*~2~. (a): $I_{2} = \left\lbrack 0,D_{2} \right)$ where *D*~2~ is the solution of $M_{2}\left( D + a_{2} \right) = s_{2}^{1}\left( D \right)$. (b) : $I_{2} = \left\lbrack 0,D_{2} \right)$ where *D*~2~ is the solution of $M_{2}\left( D + a_{2} \right) = s_{2}^{0}\left( D \right)$. (c): *I*~2~ is empty. (d) : $I_{2} = \left( D_{2min},D_{2max} \right)$ where *D*~2*min*~ and *D*~2*max*~ are the solutions of $M_{2}\left( D + a_{2} \right) = s_{2}^{0}\left( D \right)$ and $M_{2}\left( D + a_{2} \right) = s_{2}^{1}\left( D \right),$ respectively. Cases (a)--(d) are obtained with the numerical parameter values listed in [Tables 2](#tbl0002){ref-type="table"} and [3](#tbl0003){ref-type="table"}. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)Fig. 3

The value ${\overline{s}}_{2}$ depends on *D*. However, to avoid cumbersome notation we will use the more precise form ${\overline{s}}_{2}\left( D \right)$ only if necessary.

We consider the function *F*~1~(*D*), which is the infimum of *ψ*(*s*~2~, *D*), with respect to the variable *s*~2~, see [Fig. 2](#fig0002){ref-type="fig"}(b): $$\begin{aligned}
{F_{1}\left( D \right)} & {= \inf\limits_{s_{2} \in {(s_{2}^{0}{(D)},s_{2}^{1}{(D)})}}\psi\left( s_{2},D \right) = \psi\left( {\overline{s}}_{2},D \right)} \\
\end{aligned}$$We also define the functions *F*~2~(*D*) and *F*~3~(*D*): $$\begin{aligned}
{F_{2}\left( D \right)} & {= \psi\left( M_{2}\left( D + a_{2} \right),D \right)} \\
\end{aligned}$$$$\begin{aligned}
{F_{3}\left( D \right)} & {= \frac{\partial\psi}{\partial s_{2}}\left( M_{2}\left( D + a_{2} \right),D \right)} \\
\end{aligned}$$The function *F*~1~(*D*) is defined for *D* ∈ *I*~1~ where: $$I_{1} = \left\{ D \geq 0:s_{2}^{0}\left( D \right) < s_{2}^{1}\left( D \right) \right\}$$The functions *F*~2~(*D*) and *F*~3~(*D*) are defined for *D* ∈ *I*~2~ where: $$I_{2} = \left\{ D \in I_{1}:s_{2}^{0}\left( D \right) < M_{2}\left( D + a_{2} \right) < s_{2}^{1}\left( D \right) \right\}$$For all for *D* ∈ *I*~2~, *F*~1~(*D*) ≤ *F*~2~(*D*). The equality $F_{1}\left( D \right) = F_{2}\left( D \right)$ holds if, and only if, $M_{2}\left( D + a_{2} \right) = {\overline{s}}_{2}\left( D \right)$ that is, $\frac{\partial\psi}{\partial s_{2}}\left( M_{2}\left( D + a_{2} \right) \right) = 0,$ that is if, and only if, $F_{3}\left( D \right) = 0$. We define $$I_{3} = \left\{ D \in I_{2}:F_{3}\left( D \right) < 0 \right\}$$For the Monod-type growth functions ([Eq. 23](#eq0023){ref-type="disp-formula"}), straightforward computations show that the functions *M*~0~, *M*~1~ and *M*~2~ are given respectively by:

For $y \in \left\lbrack 0,\mu_{0}\left( + \infty,s_{2} \right) = \frac{m_{0}s_{2}}{L_{0} + s_{2}} \right),$$$M_{0}\left( y,s_{2} \right) = \frac{K_{0}y}{\frac{m_{0}s_{2}}{L_{0} + s_{2}} - y}$$For $y \in \left\lbrack 0,\mu_{1}\left( + \infty,s_{2} \right) = \frac{m_{1}}{1 + s_{2}/K_{i}} \right),$$$M_{1}\left( y,s_{2} \right) = \frac{K_{1}y}{\frac{m_{1}}{1 + s_{2}/K_{i}} - y}$$For $y \in \left\lbrack 0,\mu_{2}\left( + \infty \right) = m_{2} \right),$$$M_{2}\left( y \right) = \frac{K_{2}y}{m_{2} - y}$$Moreover, we have: $$s_{2}^{0}\left( D \right) = \frac{L_{0}\left( D + a_{0} \right)}{m_{0} - D - a_{0}},\quad s_{2}^{1}\left( D \right) = \frac{K_{i}\left( m_{1} - D - a_{1} \right)}{D + a_{1}}$$$$\psi\left( s_{2},D \right) = \frac{K_{0}\left( D + a_{0} \right)}{m_{0} - D - a_{0}}\frac{L_{0} + s_{2}}{s_{2} - s_{2}^{0}\left( D \right)} + \frac{\frac{K_{1}\left( K_{i} + s_{2} \right)}{s_{2}^{1}\left( D \right) - s_{2}} + s_{2}}{1 - \omega}$$and $$\frac{\partial^{2}\psi}{\partial s_{2}^{2}}\left( s_{2},D \right) = \frac{2K_{0}\left( D + a_{0} \right)}{m_{0} - D - a_{0}}\frac{L_{0} + s_{2}^{0}\left( D \right)}{\left( s_{2} - s_{2}^{0}\left( D \right) \right)^{3}} + \frac{2K_{1}\left( K_{i} + s_{2}^{1}\left( D \right) \right)}{\left( 1 - \omega \right)\left( s_{2}^{1}\left( D \right) - s_{2} \right)^{3}}$$Hence, $\frac{\partial^{2}\psi}{\partial s_{2}^{2}} > 0$ for all $s_{2} \in \left( s_{2}^{0}\left( D \right),s_{2}^{1}\left( D \right) \right),$ so that the function *s*~2~↦*ψ*(*s*~2~, *D*) is convex and, thus, it satisfies assumption **H8**, see [Fig. 2](#fig0002){ref-type="fig"}(b). The minimum ${\overline{s}}_{2}\left( D \right)$ is a solution of an algebraic equation of degree 4 in *s*~2~. Although mathematical software, such as *Maple*, cannot give solutions explicitly with respect to the parameters, ${\overline{s}}_{2}\left( D \right)$ could be obtained analytically since algebraic equations of degree 4 can theoretically be solved by quadratures. We do not try to obtain such an explicit formula. However, if the biological parameters are fixed, the function ${\overline{s}}_{2}\left( D \right)$ and, hence, $F_{1}\left( D \right) = \psi\left( {\overline{s}}_{2}\left( D \right),D \right),$ can be obtained numerically. Since *M*~2~ and *ψ* are given explicitly the functions *F*~2~(*D*) and *F*~3~(*D*) are given explicitly with respect to the biological parameters in ([Eq. 23](#eq0023){ref-type="disp-formula"}). Since $\left. D\mapsto s_{2}^{0}\left( D \right) \right.$ is increasing and $\left. D\mapsto s_{2}^{1}\left( D \right) \right.$ is decreasing, and assuming $s_{2}^{0}\left( 0 \right) < s_{2}^{1}\left( 0 \right),$ the domain of definition *I*~1~ of *F*~1~(*D*) is an interval $I_{1} = \left\lbrack 0,D_{1} \right),$ where *D*~1~ is the solution of $s_{2}^{0}\left( D \right) = s_{2}^{1}\left( D \right),$ see [Fig. 3](#fig0003){ref-type="fig"}. For the domain of definition *I*~2~ of *F*~2~(*D*), several cases can be distinguished. *I*~2~ is an interval $I_{2} = \left\lbrack 0,D_{2} \right),$ where *D*~2~ is the solution of $M_{2}\left( D + a_{2} \right) = s_{2}^{1}\left( D \right),$ see [Fig. 3](#fig0003){ref-type="fig"}(a), or the solution of equation $M_{2}\left( D + a_{2} \right) = s_{2}^{0}\left( D \right),$ see [Fig. 3](#fig0003){ref-type="fig"}(b). *I*~2~ is empty, see [Fig. 3](#fig0003){ref-type="fig"}(c). *I*~2~ is an interval $I_{2} = \left( D_{2min},D_{2max} \right)$ where *D*~2*min*~ and *D*~2*max*~ are the solutions of $M_{2}\left( D + a_{2} \right) = s_{2}^{0}\left( D \right)$ and $M_{2}\left( D + a_{2} \right) = s_{2}^{1}\left( D \right),$ respectively, see [Fig. 3](#fig0003){ref-type="fig"}(d).

5. Existence of steady-states {#sec0005}
=============================

In this paper, as we are restricted only to local stability analysis, any reference to steady-state stability should be considered as local stability. A steady-state of ([(14)](#eq0014){ref-type="disp-formula"}, [(15)](#eq0015){ref-type="disp-formula"}, [(16)](#eq0016){ref-type="disp-formula"}, [(17)](#eq0017){ref-type="disp-formula"}, [(18)](#eq0018){ref-type="disp-formula"}, [(19)](#eq0019){ref-type="disp-formula"}) is obtained by setting the right-hand sides equal to zero: $$\begin{aligned}
{\left\lbrack \mu_{0}\left( s_{0},s_{2} \right) - D - a_{0} \right\rbrack x_{0}} & {= 0} \\
\end{aligned}$$$$\begin{aligned}
{\left\lbrack \mu_{1}\left( s_{1},s_{2} \right) - D - a_{1} \right\rbrack x_{1}} & {= 0} \\
\end{aligned}$$$$\begin{aligned}
{\left\lbrack \mu_{2}\left( s_{2} \right) - D - a_{2} \right\rbrack x_{2}} & {= 0} \\
\end{aligned}$$$$\begin{aligned}
{D\left( s_{0}^{\text{in}} - s_{0} \right) - \mu_{0}\left( s_{0},s_{2} \right)x_{0}} & {= 0} \\
\end{aligned}$$$$\begin{aligned}
{- Ds_{1} + \mu_{0}\left( s_{0},s_{2} \right)x_{0} - \mu_{1}\left( s_{1},s_{2} \right)x_{1}} & {= 0} \\
\end{aligned}$$$$\begin{aligned}
{- Ds_{2} + \mu_{1}\left( s_{1},s_{2} \right)x_{1} - \omega\mu_{0}\left( s_{0},s_{2} \right)x_{0} - \mu_{2}\left( s_{2} \right)x_{2}} & {= 0} \\
\end{aligned}$$A steady-state exists (or is said to be 'meaningful') if, and only if, all its components are non-negative. Lemma 2*The only steady-state of* ([Eqs. 14](#eq0014){ref-type="disp-formula"}*--*[19](#eq0019){ref-type="disp-formula"}*), for which*$x_{0} = 0$ *or*$x_{1} = 0,$ *is the washout steady-state*${\text{SS}1} = \left( x_{0},x_{1},x_{2},s_{0},s_{1},s_{2} \right),$ *where*$x_{0} = 0,$$x_{1} = 0,$$x_{2} = 0,$$s_{0} = s_{0}^{in},$$s_{1} = 0$ *and*$s_{2} = 0$*. This steady-state always exists and is always locally stable.*

Besides the steady-state SS1, the system can have at most two other type of steady-states. •SS2: *x*~0~ \> 0, *x*~1~ \> 0 and $x_{2} = 0,$ where species *x*~2~ is washed out, while species *x*~0~ and and *x*~1~ exist. We show below that, generally, the system can have two steady-states SS2^♭^ and SS2^♯^.•SS3: *x*~0~ \> 0, *x*~1~ \> 0, and *x*~2~ \> 0, where all populations are maintained. We show below that the steady-state SS3 is unique if it exists.

We can state now the necessary and sufficient conditions of existence of SS2 and SS3. Lemma 3*If ω* ≥ 1 *then SS2 does not exist. If ω* \< 1 *then SS2 exists if, and only if,*$s_{0}^{\text{in}} \geq F_{1}\left( D \right)$*. Therefore, a necessary condition for the existence of SS2 is that D* ∈ *I*~1~, *where I*~1~ *is defined by (Eq.* [33](#eq0033){ref-type="disp-formula"}*). If*$s_{0}^{\text{in}} \geq F_{1}\left( D \right)$ *then each solution s*~2~ *of equation* $$\psi\left( s_{2} \right) = s_{0}^{\text{in}},\quad s_{2} \in \left( s_{2}^{0},s_{2}^{1} \right)$$*gives a steady-state*${\text{SS}2} = \left( x_{0},x_{1},x_{2},s_{0},s_{1},s_{2} \right)$ *where* $$\begin{array}{cl}
 & {s_{0} = M_{0}\left( D + a_{0},s_{2} \right),\quad s_{1} = M_{1}\left( D + a_{1},s_{2} \right)} \\
 & {x_{0} = \frac{D}{D + a_{0}}\left( s_{0}^{\text{in}} - s_{0} \right),\quad x_{1} = \frac{D}{D + a_{1}}\left( s_{0}^{\text{in}} - s_{0} - s_{1} \right),\quad x_{2} = 0} \\
\end{array}$$ Lemma 4*If ω* ≥ 1 *then SS3 does not exist. If ω* \< 1 *then SS3 exists if, and only if,*$s_{0}^{\text{in}} > F_{2}\left( D \right)$*. Therefore, a necessary condition of existence of SS3 is that D* ∈ *I*~2~, *where I*~2~ *is defined by (Eq.* [34](#eq0034){ref-type="disp-formula"}*). If*$s_{0}^{\text{in}} > F_{2}\left( D \right)$ *then the steady-state*${\text{SS}3} = \left( x_{0},x_{1},x_{2},s_{0},s_{1},s_{2} \right)$ *is given by* $$\begin{array}{ccl}
 & & {s_{0} = M_{0}\left( D + a_{0},M_{2}\left( D + a_{2} \right) \right),\quad s_{1} = M_{1}\left( D + a_{1},M_{2}\left( D + a_{2} \right) \right),} \\
 & & {s_{2} = M_{2}\left( D + a_{2} \right)} \\
\end{array}$$*called the break-even concentrations, and* $$\begin{array}{ccl}
 & & {x_{0} = \frac{D}{D + a_{0}}\left( s_{0}^{\text{in}} - s_{0} \right),\quad x_{1} = \frac{D}{D + a_{1}}\left( s_{0}^{\text{in}} - s_{0} - s_{1} \right),} \\
 & & {x_{2} = \frac{D}{D + a_{2}}\left( \left( 1 - \omega \right)\left( s_{0}^{\text{in}} - s_{0} \right) - s_{1} - s_{2} \right)} \\
\end{array}$$ Remark 1If $s_{0}^{\text{in}} > F_{1}\left( D \right)$ then ([Eq. 41](#eq0041){ref-type="disp-formula"}) has exactly two solutions denoted by $s_{2}^{\flat}$ and $s_{2}^{\sharp}$ and such that, see [Fig. 2](#fig0002){ref-type="fig"}(c), $$s_{2}^{0} < s_{2}^{\flat} < {\overline{s}}_{2} < s_{2}^{\sharp} < s_{2}^{1}$$If $s_{0}^{\text{in}} = F_{1}\left( D \right)$ then $s_{2}^{0} < s_{2}^{\flat} = {\overline{s}}_{2} = s_{2}^{\sharp} < s_{2}^{1}$.

To these solutions, $s_{2}^{\flat}$ and $s_{2}^{\sharp},$ correspond two steady-states of SS2, which are denoted by SS2^♭^ and SS2^♯^. These steady-states coalesce when $s_{0}^{\text{in}} = F_{1}\left( D \right)$.

Since *F*~1~(*D*) ≤ *F*~2~(*D*), the condition $s_{0}^{\text{in}} > F_{2}\left( D \right)$ for the existence of the positive steady-state SS3 implies that the condition $s_{0}^{\text{in}} > F_{1}\left( D \right)$ for the existence of the two steady-states SS2^♭^ and SS2^♯^ is satisfied. Therefore, if SS3 exists then SS2^♭^ and SS2^♯^ exist and are distinct. If $s_{0}^{\text{in}} = F_{2}\left( D \right)$ then SS3 coalesces with SS2^♭^ if *F*~3~(*D*) \< 0, and with SS2^♯^ if *F*~3~(*D*) \> 0, respectively. Remark 2Using ([Eq. 20](#eq0020){ref-type="disp-formula"}), the conditions $s_{0}^{\text{in}} > F_{1}\left( D \right)$ and $s_{0}^{\text{in}} > F_{2}\left( D \right)$ of existence of the steady-state SS2 and SS3, respectively are equivalent to the conditions $$S_{\text{ch},\text{in}} > \frac{F_{1}\left( D \right)}{Y_{3}Y_{4}}\text{and} S_{\text{ch},\text{in}} > \frac{F_{2}\left( D \right)}{Y_{3}Y_{4}}$$respectively expressed with respect to the inflowing concentration *S*~ch,\ in~. Moreover, for the nominal parameter values given in [Table 1](#tbl0001){ref-type="table"}, we have *ω* ≈ 0.53. Therefore, the steady-states SS2 and SS3 exist if the conditions on the operating parameters stated in [Lemmas 3](#enun0004){ref-type="statement"} and [4](#enun0005){ref-type="statement"}, respectively are satisfied, see [Table 7](#tbl0007){ref-type="table"}.

6. Operating diagram {#sec0006}
====================

The *Operating diagrams* show how the system behaves when we vary the two control parameters *S*~ch,\ in~ and *D* (i.e., substrate inflow rate and medium flow-rate per culture volume) in ([(1)](#eq0001){ref-type="disp-formula"}, [(2)](#eq0002){ref-type="disp-formula"}, [(3)](#eq0003){ref-type="disp-formula"}, [(4)](#eq0004){ref-type="disp-formula"}, [(5)](#eq0005){ref-type="disp-formula"}, [(6)](#eq0006){ref-type="disp-formula"}). This conventional bifurcation plot is used to visualise the existence and stability of steady-states with respect to these operating parameters, as they are the parameters most easily manipulated in a chemostat by which the system under scrutiny can be examined.

According to [Remark 2](#enun0007){ref-type="statement"}, the curve *Γ*~1~ of equation $$S_{\text{ch},\text{in}} = \frac{1}{Y_{3}Y_{4}}F_{1}\left( D \right)$$is the border to which SS2 exists, and the curve *Γ*~2~ of equation $$S_{\text{ch},\text{in}} = \frac{1}{Y_{3}Y_{4}}F_{2}\left( D \right)$$is the border to which SS3 exists, see [Fig. 4](#fig0004){ref-type="fig"}. If we want to plot the operating diagram we must fix the values of the biological parameters. In the remainder of the Section we plot the operating diagrams corresponding to cases (a)--(d) depicted in [Fig. 3](#fig0003){ref-type="fig"}.Fig. 4The curves *Γ*~1~ (black), *Γ*~2~ (red) and *Γ*~3~ (green) for case (a). (i) : regions of steady-state existence, with maintenance. On the right, a magnification for $0 < D < D_{3} = 0.058$ showing the region $\mathcal{J}_{4}$. (ii) : regions of steady-state existence and their stability, without maintenance. On the right, a magnification for $0 < D < D_{3} = 0.078$ showing the regions $\mathcal{J}_{4}$ and $\mathcal{J}_{5}$. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)Fig. 4

6.1. Operating diagram: case (a) {#sec0007}
--------------------------------

This case corresponds to the parameter values used by [@bib0024]. We have seen in [Table 2](#tbl0002){ref-type="table"} that the curves *Γ*~1~ and *Γ*~2~ are defined for *D* \< *D*~1~ and *D* \< *D*~2~, respectively and that they are tangent for $D = D_{3},$ where $D_{1} = 0.432,$$D_{2} = 0.373$ and $D_{3} = 0.058$. Therefore, they separate the operating plane (*S*~ch,\ in~, *D*) into four regions, as shown in [Fig. 4](#fig0004){ref-type="fig"}(i), labelled $\mathcal{J}_{1},$$\mathcal{J}_{2}$ and $\mathcal{J}_{3}$ and $\mathcal{J}_{4}$.Table 2Parameter values of the maintenance terms *a~i~*, $i = 0,1,2,$ for cases (a), (b) and (c) of [Fig. 3](#fig0003){ref-type="fig"}. Unspecified parameter values are as in [Table 1](#tbl0001){ref-type="table"}. The table gives the values of *D*~1~, *D*~2~ and *D*~3~ where $I_{1} = \left\lbrack 0,D_{1} \right),$$I_{2} = \left\lbrack 0,D_{2} \right)$ and $I_{3} = \left\lbrack 0,D_{3} \right)$.Table 2$K_{S,{H_{2},c}}$*a~i~D*~1~*D*~2~*D*~3~(a)$1.0 \times 10^{- 6}$0.020.4320.3730.05800.4520.3930.078(b)$4.0 \times 10^{- 6}$0.020.3290.236$I_{3} = I_{2}$00.3490.256$I_{3} = I_{2}$(c)$7.0 \times 10^{- 6}$0.020.287$I_{2} = \varnothing$00.303$I_{2} = \varnothing$

The results are summarised in [Table 4](#tbl0004){ref-type="table"}, which shows the existence of the steady-states SS1, SS2 and SS3 in the regions of the operating diagram in [Fig. 4](#fig0004){ref-type="fig"}(i).

6.2. Operating diagram: case (b) {#sec0008}
--------------------------------

This case corresponds to the parameter values used by [@bib0024], except that $K_{S,{H_{2},c}}$ is changed from $1.0 \times 10^{- 6}$ to $4.0 \times 10^{- 6}$. We have seen in [Table 2](#tbl0002){ref-type="table"} that the curves *Γ*~1~ and *Γ*~2~ are defined for *D* \< *D*~1~ and *D* \< *D*~2~, respectively and *F*~1~(*D*) \< *F*~2~(*D*) for all *D* \< *D*~2~, where $D_{1} = 0.329$ and $D_{2} = 0.236$. Therefore, they separate the operating plane (*S*~ch,\ in~, *D*) in three regions, as shown in [Fig. 5](#fig0005){ref-type="fig"}(i), labelled $\mathcal{J}_{1},$$\mathcal{J}_{3}$ and $\mathcal{J}_{4}$.Fig. 5The curves *Γ*~1~ (black), *Γ*~2~ (red) and *Γ*~3~ (green) for case (b). (i) : regions of steady-state existence, with maintenance. (ii) : regions of steady-state existence and their stability, without maintenance. On the right, a magnification for 0 \< *D* \< 0.1. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)Fig. 5

The results are summarised in [Table 5](#tbl0005){ref-type="table"}, which shows the existence of the steady-states SS1, SS2 and SS3 in the regions of the operating diagram in [Fig. 5](#fig0005){ref-type="fig"}(i). Note that the region $\mathcal{J}_{2}$ has disappeared.

6.3. Operating diagram: case (c) {#sec0009}
--------------------------------

This case corresponds to the parameter values used by [@bib0024], except that $K_{S,{H_{2},c}}$ is changed from $1.0 \times 10^{- 6}$ to $7.0 \times 10^{- 6}$. We have seen in [Table 2](#tbl0002){ref-type="table"} that the curve *Γ*~1~ is defined for $D < D_{1} = 0.287$ and that *I*~2~ is empty so that SS3 does not exist. Therefore, *Γ*~1~ separates the operating plane (*S*~ch,\ in~, *D*) in two regions, as shown in [Fig. 6](#fig0006){ref-type="fig"}(i), labelled $\mathcal{J}_{1}$ and $\mathcal{J}_{4}$.Fig. 6The curve *Γ*~1~ for case (c). (i) : regions of steady-state existence, with maintenance. (ii) : regions of steady-state existence and their stability, without maintenance. On the right, a magnification for 0 \< *D* \< 0.1.Fig. 6

The results are summarised in [Table 6](#tbl0006){ref-type="table"}, which shows the existence of the steady-states SS1 and SS2 in the regions of the operating diagram in [Fig. 6](#fig0006){ref-type="fig"}(i). Note that the region $\mathcal{J}_{3}$ of existence of SS3 has disappeared.

6.4. Operating diagram: case (d) {#sec0010}
--------------------------------

We end this discussion on the role of kinetic parameters by the presentation of this case, which presents a new behaviour that did not occur in the preceding cases: there exists *D*~2*min*~ such that for *D* \< *D*~2*min*~ the system cannot have a positive steady-state SS3. This case corresponds to the parameter values used by [@bib0024], except that three of them are changed as indicated in [Table 3](#tbl0003){ref-type="table"}. This table shows that the curves *Γ*~1~ and *Γ*~2~ are defined for *D* \< *D*~1~ and *D*~2*min*~ \< *D* \< *D*~2*max*~ and that they are tangent for $D = D_{3},$ where $D_{1} = 0.238,$$D_{2min} = 0.101,$$D_{2max} = 0.198$ and $D_{3} = 0.161$. Therefore, *Γ*~1~ and *Γ*~2~ separate the operating plane (*S*~ch,\ in~, *D*) in four regions, as shown in [Fig. 7](#fig0007){ref-type="fig"}(i), labelled $\mathcal{J}_{1},$$\mathcal{J}_{2},$$\mathcal{J}_{3}$ and $\mathcal{J}_{4}$. The results are summarised in [Table 4](#tbl0004){ref-type="table"}, which shows the existence of the steady-states SS1, SS2 and SS3 in the regions of the operating diagram in [Fig. 7](#fig0007){ref-type="fig"}(i).Fig. 7The curves *Γ*~1~ (black), *Γ*~2~ (red) and *Γ*~3~ (green) for case (d). (i) : regions of steady-state existence, with maintenance. (ii) : regions of steady-state existence and their stability, without maintenance. On the right, a magnification for 0 \< *D* \< 0.1. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)Fig. 7Table 3Parameter values of the maintenance terms *a~i~*, $i = 0,1,2,$ for case (d) of [Fig. 3](#fig0003){ref-type="fig"}: $K_{S,{H_{2},c}} = 1.2 \times 10^{- 5},$$K_{S,H_{2}} = 0.5 \times 10^{- 5}$ and $k_{m,H_{2}} = 5$. Unspecified parameter values are as in [Table 1](#tbl0001){ref-type="table"}. The table gives the values of *D*~1~, *D*~2*min*~, *D*~2*max*~ and *D*~3~ where $I_{1} = \left\lbrack 0,D_{1} \right),$$I_{2} = \left( D_{2min},D_{2max} \right)$ and $I_{3} = \left( D_{2min},D_{3} \right)$.Table 3*a~i~D*~1~*D*~2*min*~*D*~2*max*~*D*~3~(d)0.020.2380.1010.1980.16100.2580.1210.2180.181Table 4Existence of steady-states in the regions of the operating diagrams of [Fig. 4](#fig0004){ref-type="fig"}(i) and [Fig. 7](#fig0007){ref-type="fig"}(i).Table 4RegionSteady-states$\mathcal{J}_{1}$SS1$\mathcal{J}_{2} \cup \mathcal{J}_{4}$SS1, SS2^♭^, SS2^♯^$\mathcal{J}_{3}$SS1, SS2^♭^, SS2^♯^, SS3Table 5Existence of steady-states in the regions of the operating diagram of [Fig 5](#fig0005){ref-type="fig"}(i).Table 5RegionSteady-states$\mathcal{J}_{1}$SS1$\mathcal{J}_{4}$SS1, SS2^♭^, SS2^♯^$\mathcal{J}_{3}$SS1, SS2^♭^, SS2^♯^, SS3Table 6Existence of steady-states in the regions of the operating diagram of [Fig 6](#fig0006){ref-type="fig"}(i).Table 6RegionSteady-states$\mathcal{J}_{1}$SS1$\mathcal{J}_{4}$SS1, SS2^♭^, SS2^♯^Table 7Existence (with or without maintenance) and stability (without maintenance) of steady-states. The conditions for existence and stability are satisfied given that the inequalities are observed.Table 7ExistenceStabilitySS1Always existsAlways stableSS2^♭^$s_{0}^{\text{in}} > F_{1}\left( D \right)$Always unstableSS2^♯^$s_{0}^{\text{in}} > F_{1}\left( D \right)$*F*~3~(*D*) \> 0 and $s_{0}^{in} < F_{2}\left( D \right)$SS3$s_{0}^{\text{in}} > F_{2}\left( D \right)$*F*~3~(*D*) ≥ 0 or*F*~3~(*D*) \< 0 and $F_{4}\left( D,s_{0}^{\text{in}} \right) > 0$

7. Local stability without maintenance {#sec0011}
======================================

We know that SS1 is always stable. However, the analytical study of the stability of SS2 and SS3 is very difficult because the conditions for Routh-Hurwitz in the 6-dimensional case are intractable. For this reason we will consider in this section the question of the stability only for the case without maintenance ($a_{0} = a_{1} = a_{2} = 0$), since the system reduces to a three-dimensional system. The general case will be considered only numerically in [Section 8](#sec0016){ref-type="sec"}. When maintenance is not considered in the model, the steady-states SS1, SS2 and SS3 are given by: 1.${\text{SS}1} = \left( 0,0,0,s_{0}^{\text{in}},0,0 \right)$2.${\text{SS}2} = \left( x_{0},x_{1},0,s_{0},s_{1},s_{2} \right)$ where *s*~2~ a solution of equation $$s_{0}^{\text{in}} = \psi\left( s_{2} \right) = M_{0}\left( D,s_{2} \right) + \frac{M_{1}\left( D,s_{2} \right) + s_{2}}{1 - \omega}$$and $$s_{0} = M_{0}\left( D,s_{2} \right),\;\;\, s_{1} = M_{1}\left( D,s_{2} \right)\quad x_{0} = s_{0}^{\text{in}} - s_{0},\;\;\, x_{1} = s_{0}^{\text{in}} - s_{0} - s_{1}$$3.${\text{SS}3} = \left( x_{0},x_{1},x_{2},s_{0},s_{1},s_{2} \right)$ where $$\begin{aligned}
s_{2} & {= M_{2}\left( D \right),\quad s_{0} = M_{0}\left( D,s_{2} \right),\quad s_{1} = M_{1}\left( D,s_{2} \right)} \\
x_{0} & {= s_{0}^{\text{in}} - s_{0},} \\
x_{1} & {= s_{0}^{\text{in}} - s_{0} - s_{1},_{q}uadx_{2} = \left( 1 - \omega \right)\left( s_{0}^{\text{in}} - s_{0} \right) - s_{1} - s_{2}} \\
\end{aligned}$$ Proposition 2*Let*${\text{SS}2} = \left( x_{0},x_{1},0,s_{0},s_{1},s_{2} \right)$ *be a steady-state. Then SS2 is stable if, and only if, μ*~2~(*s*~2~) \< *D and*$\frac{\partial\psi}{\partial s_{2}} > 0$*.*

Therefore, SS2^♭^ is always unstable and SS2^♯^ is stable if, and only if, *μ*~2~(*s*~2~) \< *D*. This last condition is equivalent to $M_{2}\left( D \right) > s_{2}^{\sharp},$ which implies that *F*~3~(*D*) \> 0. Hence, if SS3 exists then SS2^♯^ is necessarily unstable. Therefore, SS2^♯^ is stable if, and only if, *F*~3~(*D*) \> 0 and SS3 does not exist. Proposition 3*Let*${\text{SS}3} = \left( x_{0},x_{1},x_{2},s_{0},s_{1},s_{2} \right)$ *be a steady-state. If F*~3~(*D*) ≥ 0 *then SS3 is stable as long as it exists. If F*~3~(*D*) \< 0 *then SS3 can be unstable. The instability of SS3 occurs inparticular when s*~2~ *is sufficiently close to*$s_{2}^{\flat},$ *i.e. when SS3 is sufficiently close to* SS2^♭^*.*

The condition *F*~3~(*D*) ≥ 0 is equivalent to $\frac{\partial\psi}{\partial s_{2}}\left( M_{2}\left( D \right) \right) \geq 0,$ namely $s_{2} = M_{2}\left( D \right) \in \left\lbrack {\overline{s}}_{2},s_{2}^{\sharp} \right)$. If $\frac{\partial\psi}{\partial s_{2}} < 0,$ viz. $s_{2} \in \left( s_{2}^{\flat},{\overline{s}}_{2} \right),$ then SS3 can be unstable.

When *D* is such that *F*~3~(*D*) \< 0, the determination of the boundary between the regions of stability and instability of SS3 needs to examine the Routh-Hurwitz condition of stability for SS3. For this purpose we define the following functions. Let ${\text{SS}3} = \left( x_{0},x_{1},x_{2},s_{0},s_{1},s_{2} \right)$ be a steady-state and $$E = \frac{\partial\mu_{0}}{\partial s_{0}},\quad F = \frac{\partial\mu_{0}}{\partial s_{2}},\quad G = \frac{\partial\mu_{1}}{\partial s_{1}},\quad H = - \frac{\partial\mu_{1}}{\partial s_{2}},\quad I = \frac{d\mu_{2}}{ds_{2}}$$evaluated at the steady-state SS3 defined by [(48)](#eq0048){ref-type="disp-formula"}, i.e. for $$s_{2} = M_{2}\left( D \right),\quad s_{0} = M_{0}\left( D,s_{2} \right),\quad s_{1} = M_{1}\left( D,s_{2} \right)$$For *D* ∈ *I*~3~ and $s_{0}^{\text{in}} > F_{2}\left( D \right),$ we define: $$\begin{array}{ccl}
{F_{4}\left( D,s_{0}^{\text{in}} \right)} & = & {\left( EIx_{0}x_{2} + \left\lbrack E\left( G + H \right) - \left( 1 - \omega \right)FG \right\rbrack x_{0}x_{1} \right)f_{2}} \\
 & & {+ \left( Ix_{2} + \left( G + H \right)x_{1} + \omega Fx_{0} \right)GIx_{1}x_{2}} \\
\end{array}$$where $f_{2} = Ix_{2} + \left( G + H \right)x_{1} + \left( E + \omega F \right)x_{0}$. Notice that to compute $F_{4}\left( D,s_{0}^{\text{in}} \right)$ we must replace *x*~0~, *x*~1~, *x*~2~, *s*~0~, *s*~1~ and *s*~2~ by their values at SS3, given by [(48)](#eq0048){ref-type="disp-formula"}. Hence, this function depends on the operating parameters *D* and $s_{0}^{\text{in}}$. For each fixed *D* ∈ *I*~3~, $F_{4}\left( D,s_{0}^{\text{in}} \right)$ is polynomial in $s_{0}^{\text{in}}$ of degree 3 and tends to $+ \infty$ when $s_{0}^{\text{in}}$ tends to $+ \infty$. Therefore, it is necessarily positive for large enough $s_{0}^{\text{in}}$. The values of the operating parameters *D* and $s_{0}^{\text{in}}$ for which $F_{4}\left( D,s_{0}^{\text{in}} \right)$ is positive correspond to the stability of SS3 as shown in the following proposition. Proposition 4*Let*${\text{SS}3} = \left( x_{0},x_{1},x_{2},s_{0},s_{1},s_{2} \right)$ *be a steady-state. If F*~3~(*D*) \< 0 *then SS3 is stable if, and only if,*$F_{4}\left( D,s_{0}^{\text{in}} \right) > 0$*.*

The results on the existence of steady-states (with or without maintenance) of [Lemmas 2](#enun0003){ref-type="statement"}--[4](#enun0005){ref-type="statement"}, and their stability (without maintenance) of [Props 2](#enun0008){ref-type="statement"}--[4](#enun0010){ref-type="statement"}, are summarised in [Table 7](#tbl0007){ref-type="table"}.

7.1. Operating diagram: case (a) {#sec0012}
--------------------------------

This case corresponds to the parameter values used by [@bib0024] but without maintenance. We see from [Table 2](#tbl0002){ref-type="table"} that the curves *Γ*~1~ and *Γ*~2~ of the operating diagram, given by ([Eq. 45](#eq0045){ref-type="disp-formula"}) and ([Eq. 46](#eq0046){ref-type="disp-formula"}), respectively are defined now for *D* \< *D*~1~ and *D* \< *D*~2~, respectively and that they are tangent for $D = D_{3},$ where $D_{1} = 0.452,$$D_{2} = 0.393$ and $D_{3} = 0.078$. Beside these curves, we plot also on the operating diagram of [Fig. 4](#fig0004){ref-type="fig"}(ii), the curve *Γ*~3~ of equation $$F_{4}\left( D,Y_{3}Y_{4}S_{\text{ch},\text{in}} \right) = 0$$According to [Prop. 4](#enun0010){ref-type="statement"}, this curve is defined for $D < D_{3} = 0.078$ and it separates the region of existence of SS3 into two subregions labelled $\mathcal{J}_{3}$ and $\mathcal{J}_{5},$ such that SS3 is stable in $\mathcal{J}_{3}$ and unstable in $\mathcal{J}_{5}$. The other regions $\mathcal{J}_{1},$$\mathcal{J}_{2}$ and $\mathcal{J}_{4}$ are defined as in the previous section. The operating diagram is shown [Fig. 4](#fig0004){ref-type="fig"}(ii). It looks very similar to [Fig. 4](#fig0004){ref-type="fig"}(i), except near the origin, as it is indicated in the magnification for $0 < D < D_{3} = 0.078$. From [Table 7](#tbl0007){ref-type="table"}, we deduce the following result Proposition 5[Table 8](#tbl0008){ref-type="table"} *shows the existence and stability of the steady-states SS1, SS2 and SS3 in the regions of the operating diagram in* [Fig. 4](#fig0004){ref-type="fig"}*(ii).*Table 8Existence and stability of steady-states in the regions of the operating diagrams of [Fig. 4](#fig0004){ref-type="fig"}(ii) and [Fig. 7](#fig0007){ref-type="fig"}(ii)., where S and U indicate stable and unstable steady-states, respectively.Table 8RegionSS1SS2^♭^SS2^♯^SS3$\mathcal{J}_{1}$S$\mathcal{J}_{2}$SUS$\mathcal{J}_{3}$SUUS$\mathcal{J}_{4}$SUU$\mathcal{J}_{5}$SUUU

7.2. Operating diagram: case (b) {#sec0013}
--------------------------------

We see from [Table 2](#tbl0002){ref-type="table"} that the curves *Γ*~1~ and *Γ*~2~ are defined now for $D < D_{1} = 0.349$ and $D < D_{2} = 0.256,$ respectively and that *F*~1~(*D*) \< *F*~2~(*D*) for all *D*. Beside these curves, we plot also on the operating diagram of [Fig. 5](#fig0005){ref-type="fig"}(ii), the curve *Γ*~3~ of equation ([Eq. 50](#eq0050){ref-type="disp-formula"}) which separates the region of existence of SS3 into two subregions labelled $\mathcal{J}_{3}$ and $\mathcal{J}_{5},$ such that SS3 is stable in $\mathcal{J}_{3}$ and unstable in $\mathcal{J}_{5}$. Therefore, the curves *Γ*~1~, *Γ*~2~ and *Γ*~3~ separate the operating plane (*S*~ch,\ in~, *D*) into four regions, as shown in [Fig. 5](#fig0005){ref-type="fig"}(ii), labelled $\mathcal{J}_{1},$$\mathcal{J}_{3},$$\mathcal{J}_{4}$ and $\mathcal{J}_{5}$. Proposition 6[Table 9](#tbl0009){ref-type="table"} *shows the existence and stability of the steady-states SS1, SS2 and SS3 in the regions of the operating diagram in* [Fig. 5](#fig0005){ref-type="fig"}*(ii)*.Table 9Existence and stability of steady-states in the regions of the operating diagram of [Fig. 5](#fig0005){ref-type="fig"}(ii).Table 9RegionSS1SS2^♭^SS2^♯^SS3$\mathcal{J}_{1}$S$\mathcal{J}_{3}$SUUS$\mathcal{J}_{4}$SUU$\mathcal{J}_{5}$SUUU

7.3. Operating diagram: case (c) {#sec0014}
--------------------------------

We see from [Table 2](#tbl0002){ref-type="table"} that *Γ*~1~ is defined for $D < D_{1} = 0.303$ and that *I*~2~ is empty so that SS3 does not exist. Therefore, *Γ*~1~ separates the operating plane (*S*~ch,\ in~, *D*) into two regions, as shown in [Fig. 6](#fig0006){ref-type="fig"}(ii), labelled $\mathcal{J}_{1}$ and $\mathcal{J}_{4}$. Proposition 7[Table 10](#tbl0010){ref-type="table"} *shows the existence and stability of the steady-states SS1, SS2 and SS3 in the regions of the operating diagram in* [Fig. 6](#fig0006){ref-type="fig"}*(ii).*Table 10Existence and stability of steady-states in the regions of the operating diagram of [Fig. 6](#fig0006){ref-type="fig"}(ii).Table 10RegionSS1SS2^♭^SS2^♯^SS3$\mathcal{J}_{1}$S$\mathcal{J}_{4}$SUU

7.4. Operating diagram: case (d) {#sec0015}
--------------------------------

We see in [Table 3](#tbl0003){ref-type="table"} that the curves *Γ*~1~ and *Γ*~2~ are defined for *D* \< *D*~1~ and *D*~2*min*~ \< *D* \< *D*~2*max*~ and that they are tangent for $D = D_{3},$ where $D_{1} = 0.258$ and $D_{2min} = 0.121,$$D_{2max} = 0.218$ and $D_{3} = 0.181$. Beside these curves, we plot also on the operating diagram of [Fig. 7](#fig0007){ref-type="fig"}(ii), the curve *Γ*~3~ defined by ([Eq. 50](#eq0050){ref-type="disp-formula"}), which separates the region of existence of SS3 into two subregions labelled $\mathcal{J}_{3}$ and $\mathcal{J}_{5},$ such that SS3 is stable in $\mathcal{J}_{3}$ and unstable in $\mathcal{J}_{5}$. Therefore, the curves *Γ*~1~, *Γ*~2~ and *Γ*~3~ separate the operating plane (*S*~ch,\ in~, *D*) into five regions, as shown in [Fig. 7](#fig0007){ref-type="fig"}(ii), labelled $\mathcal{J}_{1},$$\mathcal{J}_{2},$$\mathcal{J}_{3},$$\mathcal{J}_{4}$ and $\mathcal{J}_{5}$. Proposition 8[Table 8](#tbl0008){ref-type="table"} *shows the existence and stability of the steady-states SS1, SS2 and SS3 in the regions of the operating diagram in* [Fig. 7](#fig0007){ref-type="fig"}*(ii).*

8. Numerical analysis to illustrate the analytical results {#sec0016}
==========================================================

The aim of this section is to study numerically (the method is explained in [Appendix A](#sec0023){ref-type="sec"}) the existence and stability of the steady-states SS2 and SS3. We obtain numerically the operating diagrams that were described in [Sections 5](#sec0005){ref-type="sec"} and [7](#sec0011){ref-type="sec"}. The results in this section confirm the results on existence of the steady-states obtained in [Section 5](#sec0005){ref-type="sec"} in the cases with or without maintenance and the results of stability obtained in [Section 7](#sec0011){ref-type="sec"} in the case without maintenance. These results permit also to extend the analytical results and elucidate the problem of the local stability of SS2 and SS3, which was left open in [Section 7](#sec0011){ref-type="sec"}.

8.1. Operating diagram: case (a) {#sec0017}
--------------------------------

We endeavoured to find numerically the operating conditions under which SS3 is unstable, previously unreported by [@bib0024]. Given that we have determined analytically in [Proposition 3](#enun0009){ref-type="statement"} that when SS3 is close to SS2^♭^ it becomes unstable, we performed numerical simulations with the parameters defined in [Table 1](#tbl0001){ref-type="table"} over an operating region similar to that shown in [Fig. 2](#fig0002){ref-type="fig"} from [@bib0024] whilst also satisfying our conditions. In [Fig. 8](#fig0008){ref-type="fig"} we show the case when maintenance is excluded. When magnified, we observe more clearly that region $\mathcal{J}_{5}$ does exist for the conditions described above, and also note that the region $\mathcal{J}_{4}$ occurs in a small area between $\mathcal{J}_{1}$ and $\mathcal{J}_{5},$ which corresponds to the results shown in [Fig. 4](#fig0004){ref-type="fig"}(ii), and is in agreement with [Proposition 5](#enun0011){ref-type="statement"}. In [Fig. 9](#fig0009){ref-type="fig"} we confirm that region $\mathcal{J}_{5}$ does exist for the conditions described above, when maintenance is included, but could not be determined analytically, the curve *Γ*~3~ is absent in [Fig. 4](#fig0004){ref-type="fig"}(i). Furthermore, we demonstrate that there is evidence of Hopf bifurcation, which occurs along the boundary of *F*~3~(*D*) for values of *D* \< *D*~3~ by selecting values of *S*~ch,\ in~ (indicated by $\left( \alpha \right) - \left( \delta \right)$ in [Fig. 9](#fig0009){ref-type="fig"}) at a fixed dilution rate of $0.01\mspace{6mu} d^{- 1},$ and running dynamic simulations for 10000 *d*. The three-dimensional phase plots, with the axes representing biomass concentrations, are shown in [Fig. 10](#fig0010){ref-type="fig"}, and show that as *S*~ch,\ in~ approaches $\mathcal{J}_{3}$ from $\mathcal{J}_{5},$ emergent periodic orbits are shown to diminish to a stable limit cycle at the boundary (see  [Appendix B](#sec0024){ref-type="sec"} for details). Subsequently, increasing *S*~ch,\ in~ to $\mathcal{J}_{3}$ results in the orbit reducing to a fixed point equilibrium at SS3.Fig. 8Numerical analysis for the existence and stability of steady-states for case (a), without maintenance. On the right, a magnification for 0 \< *D* \< 0.16.Fig. 8Fig. 9Numerical analysis for the existence and stability of steady-states for case (a), with maintenance. This is a magnification for 0 \< *D* \< 0.1, showing the presence and extent of region $\mathcal{J}_{5}$ undetectable by the analytical method. The coordinates labelled $\left( \alpha \right) - \left( \delta \right)$ are subsequently used to simulate the system dynamics, as shown in the proceeding [Fig. 10](#fig0010){ref-type="fig"}.Fig. 9Fig. 10Three-dimensional phase plane diagrams of the biomass dynamics for $t = 10000\mspace{6mu} d,$ showing initial (green dot) and final (red dot) conditions for a dilution rate, $D = 0.01\mspace{6mu} d^{- 1}$ and chlorophenol input, *S*~ch,\ *in*~ (*kgCOD*/*m*^3^), of (*α*): $S_{\text{ch},in} = 0.01$ - the system converges to SS1, (*β*): $S_{\text{ch},in} = 0.097$ - the system enters a periodic orbit of increasing amplitude, ultimately converging to SS1, (*γ*): $S_{\text{ch},in} = 0.10052$ - the system is close to a stable limit cycle, (*δ*): $S_{\text{ch},in} = 0.16$ - the system undergoes damped oscillations and converges to SS3. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)Fig. 10

8.2. Operating diagram: case (b) {#sec0018}
--------------------------------

Whilst the numerical parameters chosen for this work are taken from the original study [@bib0024], there somewhat arbitrary nature leaves room to explore the impact of the parameters on the existence and stability of the steady-states. Case (b), discussed in [Sections 6.2](#sec0008){ref-type="sec"} and [7.2](#sec0013){ref-type="sec"}, involves a small increase to the half-saturation constant (or inverse of substrate affinity), $K_{S,{H_{2},c}},$ of the chlorophenol degrader on hydrogen. Following the same approach as with the preceding case, we confirm in [Fig. 11](#fig0011){ref-type="fig"}(i) the [Proposition 6](#enun0012){ref-type="statement"} in the scenario without maintenance. Furthermore, the extension of this proposition with maintenance included, corresponding to the existence and stability of all three steady-states given in [Table 9](#tbl0009){ref-type="table"}, is shown in Fig [11](#fig0011){ref-type="fig"}(ii). It shows the region $\mathcal{J}_{5}$ that cannot be obtained analytically (cf. [Fig. 5](#fig0005){ref-type="fig"}(i)). In both cases, region $\mathcal{J}_{2}$ has disappeared, as observed analytically. Additionally, the ideal region $\mathcal{J}_{3},$ where all organisms are present and stable, diminishes.Fig. 11Numerical analysis for the existence and stability of steady-states for case (b). (i) : without maintenance. (ii) : with maintenance.Fig. 11

8.3. Operating diagram: case (c) {#sec0019}
--------------------------------

Here, $K_{S,H_{2},c},$ was further increased and confirm the [Proposition 7](#enun0013){ref-type="statement"}, where the function SS3 never exist and SS2 never stable for the case without maintenance. The extension of this proposition to the case with maintenance, shown in [Table 10](#tbl0010){ref-type="table"}, produce similar results as shown in the comparison of [Fig. 12](#fig0012){ref-type="fig"}(i) and (ii).Fig. 12Numerical analysis for the existence and stability of steady-states for case (c). (i) : without maintenance. (ii) : with maintenance. On the right, a magnification for 0 \< *D* \< 0.1.Fig. 12

8.4. Operating diagram: case (d) {#sec0020}
--------------------------------

With the final investigated scenario, where $k_{m,H_{2}} < k_{m,\text{ch}}$ and $K_{S,H_{2}} < K_{S,{H_{2},c}},$ we observe once again the presence of all operating regions, $\mathcal{J}_{1} - \mathcal{J}_{5},$ without and with maintenance, as shown in [Fig. 13](#fig0013){ref-type="fig"}. It can be seen that regions $\mathcal{J}_{4}$ and $\mathcal{J}_{5}$ increase at low dilution rates across a much larger range of *S*~ch,\ in~ than in the default case (a), and the desired condition (stable SS3) is restricted to a much narrower set of *D*.Fig. 13Numerical analysis for the existence and stability of steady-states for case (d). (i) : without maintenance. (ii) : with maintenance.Fig. 13

As with the previous cases, the numerical analysis for case (d) confirms the [Proposition 8](#enun0014){ref-type="statement"} without maintenance and its extension to the case with maintenance, indicated in [Table 8](#tbl0008){ref-type="table"}.

9. The role of kinetic parameters {#sec0021}
=================================

Finally, we give brief consideration to the characterisation of the four cases discussed in the preceding sections. The main difference between cases (a) or (b) and cases (c) or (d) is that, for small values of *D*, the coexistence steady-state SS3 can exist for cases (a) and (b), but cannot exist for cases (c) or (d). The cases (a) or (b) occur if and only if $s_{2}^{0}\left( 0 \right) < M_{2}\left( 0 \right)$ holds or $s_{2}^{0}\left( 0 \right) = M_{2}\left( 0 \right)$ and $\frac{ds_{2}^{0}}{dD}\left( 0 \right) < \frac{dM_{2}}{dD}\left( 0 \right)$ hold, viz. $$\begin{array}{r}
{\frac{L_{0}a_{0}}{m_{0} - a_{0}} < \frac{K_{2}a_{2}}{m_{2} - a_{2}}\text{or}} \\
\end{array}$$$$\begin{array}{r}
{\frac{L_{0}a_{0}}{m_{0} - a_{0}} = \frac{K_{2}a_{2}}{m_{2} - a_{2}}\text{and}\frac{L_{0}m_{0}}{\left( m_{0} - a_{0} \right)^{2}} < \frac{K_{2}m_{2}}{\left( m_{2} - a_{2} \right)^{2}}} \\
\end{array}$$The cases (c) or (d) occur if and only if $s_{2}^{0}\left( 0 \right) > M_{2}\left( 0 \right)$ holds or $s_{2}^{0}\left( 0 \right) = M_{2}\left( 0 \right)$ and $\frac{ds_{2}^{0}}{dD}\left( 0 \right) > \frac{dM_{2}}{dD}\left( 0 \right)$ hold, viz. $$\begin{array}{r}
{\frac{L_{0}a_{0}}{m_{0} - a_{0}} > \frac{K_{2}a_{2}}{m_{2} - a_{2}}\text{or}} \\
\end{array}$$$$\begin{array}{r}
{\frac{L_{0}a_{0}}{m_{0} - a_{0}} = \frac{K_{2}a_{2}}{m_{2} - a_{2}}\text{and}\frac{L_{0}m_{0}}{\left( m_{0} - a_{0} \right)^{2}} > \frac{K_{2}m_{2}}{\left( m_{2} - a_{2} \right)^{2}}} \\
\end{array}$$

Notice that it is easy to make the difference between case (c) and case (d): the first occurs when $M_{2}\left( D_{1} \right) < s_{2}^{0}\left( D_{1} \right)$ and the second when $M_{2}\left( D_{1} \right) > s_{2}^{0}\left( D_{1} \right)$. Since *D*~1~ is the positive solution of the algebraic quadratic equation $s_{2}^{0}\left( D \right) = s_{2}^{1}\left( D \right),$ it is possible to have an expression for *D*~1~ with respect to the biological parameters. However, this is a complicated expression involving many parameters and the preceding conditions $M_{2}\left( D_{1} \right) < s_{2}^{0}\left( D_{1} \right)$ or $M_{2}\left( D_{1} \right) > s_{2}^{0}\left( D_{1} \right)$ have no biological interpretation. We simply remark here that the function $s_{2}^{0}\left( D \right)$ has a vertical asymptote for $D = m_{0} - a_{0}$ and the function *M*~2~(*D*) has a vertical asymptote for $D = m_{2} - a_{2}$. Therefore, if $m_{0} - a_{0} < m_{2} - a_{2}$ then case (c) occurs, so that a necessary (but not sufficient) condition for case (d) to occur is $m_{0} - a_{0} > m_{2} - a_{2}$. If *m*~2~ is sufficiently small then case (d) can occur.

The observations from the numerical analysis suggest that the role of the chlorophenol degrader as a secondary hydrogen scavenger is critical in maintaining full chlorophenol mineralisation and system stability, particularly at higher dilution rates, as shown by comparing cases (c) and (d) . More significantly, the results coupled with the parameter relationships shown in ([(51)](#eq0051){ref-type="disp-formula"}, [(52)](#eq0052){ref-type="disp-formula"}, [(53)](#eq0053){ref-type="disp-formula"}, [(54)](#eq0054){ref-type="disp-formula"}), highlight the necessary conditions under which the ideal case (SS3 stable) is achieved and, in general, this is a coupling of the two key parameters describing the half-saturation constant and maximum specific growth rates between the two hydrogen competitors.

10. Conclusions {#sec0022}
===============

In this work we have generalised a simplified mechanistic model describing the anaerobic mineralisation of chlorophenol in a two-step food-web. We give conditions for the existence and stability of the steady-states in the case that maintenance is excluded from the system. However, with a decay term present, purely analytical determination of stability was not achievable.

We confirm the findings of previous numerical analysis by [@bib0024] that with chlorophenol as the sole input substrate, three steady-states are possible. However, the analysis goes further and we determine that under certain operating conditions, two of these steady-states (SS2 and SS3) can become stable, whilst SS1 always exists and is always stable. Furthermore, without maintenance we can explicitly determine the stability of the system, and form analytical expressions of the boundaries between the different stability regions.

As the boundary of $\mathcal{J}_{3}$ is not open to analytical determination in the case with maintenance, we determined numerically (substituting the general growth function with the classical Monod-type growth kinetics) the existence and stability of the system over a range of practical operating conditions (dilution rate and chlorophenol input). For comparison and confirmation, we also performed this for the case without maintenance and found the same regions in both cases, with variations only in their shape and extent. For example, whilst the boundary between $\mathcal{J}_{1}$ and $\mathcal{J}_{4}$ terminates at the origin without maintenance, with maintenance it is located at *F*~1~(0)/*Y*~3~*Y*~4~ ≈ 0.0195. More interestingly, the addition of a decay term results in an extension of the SS3 unstable steady-state, reducing the potential for successful chlorophenol demineralisation at relatively low dilution rates and substrate input concentrations. Additionally, we show that at the boundary between $\mathcal{J}_{3}$ and $\mathcal{J}_{5},$ there is numerical evidence of Hopf bifurcation occurring and that a limit cycle in SS3 emerges.

Finally, we gave an example of how the model could be used to probe the system to answer specific questions regarding model parameterisation. Here we have indicated that a switch in dominance between two organisms competing for hydrogen results in the system becoming unstable and a loss in viability. This is perhaps intuitive to microbiologists, but here it has been proven using mathematical analysis and could be used to determine critical limits of the theoretical parameter values in shifting between a stable and unstable system. Whilst parameters are not arbitrary in reality, the potential for microbial engineering or synthetic biology to manipulate the properties of organisms makes this observation all the more pertinent.

Appendix A. Numerical methods {#sec0023}
=============================

We consider sets of operating parameters (*D* and *S*~ch,\ in~) for each of the three steady-states, and using Matlab, the complex polynomials for each steady-state can be solved by substitution of parameter values (see [Table 1](#tbl0001){ref-type="table"}) into the explicit solution. By investigating the signs of the solutions and the eigenvalues, respectively we determine which steady-states are meaningful and stable. By exploring a localised region of suitable operating parameters, we then generate a phase plot showing where each steady-state is stable, bistable or unstable.

Appendix B. Numerical evidence for Hopf Bifurcation {#sec0024}
===================================================

In [Section 8](#sec0016){ref-type="sec"}, we show the operating diagrams with the parameters given in [Table 1](#tbl0001){ref-type="table"}, and determine numerically that as the parameter *S*~ch,\ in~ increases at a fixed dilution rate ($D = 0.01\mspace{6mu} d^{- 1}$), the system bifurcates through several stability domains. We claim that as we cross the boundary between regions $\mathcal{J}_{5}$ and $\mathcal{J}_{3},$ we observe a Hopf bifurcation, and, in $\mathcal{J}_{5},$ close to the boundary with $\mathcal{J}_{3},$ a limit cycle appears. In order to test this numerically, we checked the real parts of the six eigenvalues at each point along the transect shown in Fig [9](#fig0009){ref-type="fig"} (10000 points in total), and plotted their values. [Fig B.14](#fig0014){ref-type="fig"} indicates the conditions for a Hopf bifurcation are satisfied as eigenvalues 2 and 3 both change their sign when passing through the coordinate (0, 0.1034) and the real part of all eigenvalues 1, 4 and 6 remain negative.Fig. B.14Real parts of the eigenvalues determined at $D = 0.01$ and $S_{\text{ch},\text{in}} = \left\lbrack 0.08,0.12 \right\rbrack,$ in the case with maintenance. The red vertical lines indicate the location where the eigenvalue crosses zero. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)Fig. B.14

Appendix C. Proofs {#sec0025}
==================

In this Section we give the proofs of the results.

C1. General case {#sec0026}
----------------

In these proofs, we do not assume that the growth function are of Monod type ([Eq. 23](#eq0023){ref-type="disp-formula"}). We only assume that the growth functions satisfy **H1**--**H8**.

### Proof of [Lemma 1](#enun0002){ref-type="statement"} {#sec0030}

Let *s*~2~ ≥ 0 be fixed. By **H4**, the function $$\left. s_{0} \in \left\lbrack 0, + \infty \right)\mapsto\mu_{0}\left( s_{0},s_{2} \right) \in \left\lbrack 0,\mu_{0}\left( + \infty,s_{2} \right) \right) \right.$$is monotonically increasing. Hence, it has an inverse function denoted by: $$\left. y \in \left\lbrack 0,\mu_{0}\left( + \infty,s_{2} \right) \right)\mapsto M_{0}\left( y,s_{2} \right) \in \left\lbrack 0, + \infty \right), \right.$$such that for all *s*~0~ ≥ 0, *s*~2~ ≥ 0 and $y \in \left\lbrack 0,\mu_{0}\left( + \infty,s_{2} \right) \right)$ ([Eq. 25](#eq0025){ref-type="disp-formula"}) holds.

The proofs for the existence of *M*~1~ and *M*~2~ are similar and are left to the reader.

C2. Existence of steady-states {#sec0027}
------------------------------

### Proof of [Lemma 2](#enun0003){ref-type="statement"} {#sec0031}

Assume first that $x_{0} = 0$. Then, as a consequence of ([Eq. 38](#eq0038){ref-type="disp-formula"}), $s_{0} = s_{0}^{\text{in}}$ and, as a consequence of ([Eq. 39](#eq0039){ref-type="disp-formula"}), we have: $$Ds_{1} + \mu_{1}\left( s_{1},s_{2} \right)x_{1} = 0$$which implies $s_{1} = 0$ and $\mu_{1}\left( s_{1},s_{2} \right)x_{1} = 0$. Therefore, as a consequence of ([Eq. 36](#eq0036){ref-type="disp-formula"}) we have $x_{1} = 0$. Replacing $x_{0} = 0$ and $x_{1} = 0$ in ([Eq. 40](#eq0040){ref-type="disp-formula"}), we have: $$Ds_{2} + \mu_{2}\left( s_{2} \right)x_{2} = 0$$which implies $s_{2} = 0$ and $\mu_{2}\left( s_{2} \right)x_{2} = 0$. Therefore, as a consequence of ([Eq. 37](#eq0037){ref-type="disp-formula"}), $x_{2} = 0$. Hence, the steady-state is SS1.

Assume now that $x_{1} = 0$. Then, as a consequence of ([Eq. 40](#eq0040){ref-type="disp-formula"}), we have: $$Ds_{2} + \omega\mu_{0}\left( s_{0},s_{2} \right)x_{0} + \mu_{2}\left( s_{2} \right)x_{2} = 0$$which implies $s_{2} = 0,$$\mu_{0}\left( s_{0},s_{2} \right)x_{0} = 0$ and $\mu_{2}\left( s_{2} \right)x_{2} = 0$. Therefore, as a consequence of ([Eq. 35](#eq0035){ref-type="disp-formula"}), $x_{0} = 0$. As shown previously this implies that the steady-state is SS1. Evaluated at SS1 the Jacobian matrix of ([(14)](#eq0014){ref-type="disp-formula"}, [(15)](#eq0015){ref-type="disp-formula"}, [(16)](#eq0016){ref-type="disp-formula"}, [(17)](#eq0017){ref-type="disp-formula"}, [(18)](#eq0018){ref-type="disp-formula"}, [(19)](#eq0019){ref-type="disp-formula"}) is: $$\begin{bmatrix}
{- D - a_{0}} & 0 & 0 & 0 & 0 & 0 \\
0 & {- D - a_{1}} & 0 & 0 & 0 & 0 \\
0 & 0 & {- D - a_{2}} & 0 & 0 & 0 \\
0 & 0 & 0 & {- D} & 0 & 0 \\
0 & 0 & 0 & 0 & {- D} & 0 \\
0 & 0 & 0 & 0 & 0 & {- D} \\
\end{bmatrix}$$Thus, SS1 is stable.

### Proof of [Lemma 3](#enun0004){ref-type="statement"} {#sec0032}

Since *x*~0~ \> 0 and *x*~1~ \> 0, then, as a consequence of ([Eqs. 35](#eq0035){ref-type="disp-formula"}, [36](#eq0036){ref-type="disp-formula"}), we have: $$\mu_{0}\left( s_{0},s_{2} \right) = D + a_{0},\quad\mu_{1}\left( s_{1},s_{2} \right) = D + a_{1}$$Hence, using ([Eqs. 25](#eq0025){ref-type="disp-formula"}, [26](#eq0026){ref-type="disp-formula"}), we have: $$s_{0} = M_{0}\left( D + a_{0},s_{2} \right),\quad s_{1} = M_{1}\left( D + a_{1},s_{2} \right)$$Using ([Eqs. 38](#eq0038){ref-type="disp-formula"}, [39](#eq0039){ref-type="disp-formula"}), we have ([Eq. 42](#eq0042){ref-type="disp-formula"}). Using ([Eq. 40](#eq0040){ref-type="disp-formula"}) we obtain: $$- s_{2} + \left( s_{0}^{\text{in}} - s_{0} - s_{1} \right) - \omega\left( s_{0}^{\text{in}} - s_{0} \right) = 0$$If *ω* ≥ 1 this equation has no solution. If *ω* \< 1 this equation is equivalent to: $$s_{0}^{\text{in}} = s_{0} + \frac{s_{1} + s_{2}}{1 - \omega}.$$Using ([Eq. C.1](#eq0055){ref-type="disp-formula"}) we see that *s*~2~ must be a solution of ([Eq. 41](#eq0041){ref-type="disp-formula"}). Since *s*~1~ \> 0 and *s*~2~ \> 0 then, from ([Eq. C.2](#eq0056){ref-type="disp-formula"}) we have necessarily: $$s_{1} + s_{2} = \left( 1 - \omega \right)\left( s_{0}^{\text{in}} - s_{0} \right) > 0$$so that $s_{0}^{\text{in}} - s_{0} > 0$. From ([Eq. 42](#eq0042){ref-type="disp-formula"}) we deduce that *x*~0~ \> 0. Since $s_{0}^{\text{in}} - s_{0} > 0$ and *s*~2~ \> 0 then, from ([Eq. C.2](#eq0056){ref-type="disp-formula"}) we have necessarily: $$\omega\left( s_{0}^{\text{in}} - s_{0} \right) + s_{2} = s_{0}^{\text{in}} - s_{0} - s_{1} > 0$$so that $s_{0}^{\text{in}} - s_{0} - s_{1} > 0$ From ([Eq. 42](#eq0042){ref-type="disp-formula"}) we deduce that *x*~1~ \> 0.

### Proof of [Lemma 4](#enun0005){ref-type="statement"} {#sec0033}

Since *x*~0~ \> 0, *x*~1~ \> 0 and *x*~2~ \> 0, then, as a consequence of ([(35)](#eq0035){ref-type="disp-formula"}, [(36)](#eq0036){ref-type="disp-formula"}, [(37)](#eq0037){ref-type="disp-formula"}), we obtain: $$\mu_{0}\left( s_{0},s_{2} \right) = D + a_{0},\quad\mu_{1}\left( s_{1},s_{2} \right) = D + a_{1},\quad\mu_{2}\left( s_{2} \right) = D + a_{2}$$Hence, using ([(25)](#eq0025){ref-type="disp-formula"}, [(26)](#eq0026){ref-type="disp-formula"}, [(27)](#eq0027){ref-type="disp-formula"}), *s*~0~, *s*~1~ and *s*~2~ are given by ([Eq. 43](#eq0043){ref-type="disp-formula"}). Using ([(38)](#eq0038){ref-type="disp-formula"}, [(39)](#eq0039){ref-type="disp-formula"}, [(40)](#eq0040){ref-type="disp-formula"}) we obtain ([Eq. 44](#eq0044){ref-type="disp-formula"}). For *x*~2~ to be positive it is necessary that *s*~0~, *s*~1~ and *s*~2~ satisfy the condition: $$\left( 1 - \omega \right)\left( s_{0}^{\text{in}} - s_{0} \right) > s_{1} + s_{2},$$If *ω* ≥ 1 this equation has no solution. If *ω* \< 1 this equation is equivalent to the condition: $$s_{0}^{\text{in}} > s_{0} + \frac{s_{1} + s_{2}}{1 - \omega}.$$Using ([Eq. 43](#eq0043){ref-type="disp-formula"}), this condition is the same as $$s_{0}^{\text{in}} > \psi\left( M_{2}\left( D + a_{2} \right) \right) = F_{2}\left( D \right)$$Therefore, from ([Eq. C.3](#eq0057){ref-type="disp-formula"}) we have $s_{0}^{\text{in}} - s_{0} > 0$ and $s_{0}^{\text{in}} - s_{0} - s_{1} > 0,$ so that *x*~0~ \> 0 and *x*~1~ \> 0.

C3. Stability of steady-states {#sec0028}
------------------------------

We use the change of variables: $$z_{0} = s_{0} + x_{0},\quad z_{1} = s_{1} + x_{1} - x_{0},\quad z_{2} = s_{2} + x_{2} + \omega x_{0} - x_{1}$$Therefore, ([Eqs. 14](#eq0014){ref-type="disp-formula"}--[19](#eq0019){ref-type="disp-formula"}), with $a_{0} = a_{1} = a_{2} = 0,$ become: $$\begin{aligned}
\frac{dx_{0}}{dt} & {= - Dx_{0} + \mu_{0}\left( z_{0} - x_{0},z_{2} - \omega x_{0} + x_{1} - x_{2} \right)x_{0}} \\
\end{aligned}$$$$\begin{aligned}
\frac{dx_{1}}{dt} & {= - Dx_{1} + \mu_{1}\left( z_{1} + x_{0} - x_{1},z_{2} - \omega x_{0} + x_{1} - x_{2} \right)x_{1}} \\
\end{aligned}$$$$\begin{aligned}
\frac{dx_{2}}{dt} & {= - Dx_{2} + \mu_{2}\left( z_{2} - \omega x_{0} + x_{1} - x_{2} \right)x_{2}} \\
\end{aligned}$$$$\begin{aligned}
\frac{dz_{0}}{dt} & {= D\left( s_{0}^{\text{in}} - z_{0} \right)} \\
\end{aligned}$$$$\begin{aligned}
\frac{dz_{1}}{dt} & {= - Dz_{1}} \\
\end{aligned}$$$$\begin{aligned}
\frac{dz_{2}}{dt} & {= - Dz_{2}} \\
\end{aligned}$$In the variables (*x*~0~, *x*~1~, *x*~2~, *z*~0~, *z*~1~, *z*~2~) where *z*~0~, *z*~1~ and *z*~2~ are defined by ([Eq. C.4](#eq0058){ref-type="disp-formula"}), the steady-states SS1, SS2 and SS3 are given by: 1.${\text{SS}1} = \left( 0,0,0,s_{0}^{\text{in}},0,0 \right)$2.${\text{SS}2} = \left( x_{0},x_{1},0,s_{0}^{\text{in}},0,0 \right),$ where *x*~0~ and *x*~1~ are defined by ([Eq. 47](#eq0047){ref-type="disp-formula"}).3.${\text{SS}3} = \left( x_{0},x_{1},x_{2},s_{0}^{\text{in}},0,0 \right),$ where *x*~0~, *x*~1~ and *x*~2~ are defined by ([Eq. 48](#eq0048){ref-type="disp-formula"}).

Let $\left( x_{0},x_{1},x_{2},s_{0}^{\text{in}},0,0 \right)$ be a steady-state. The Jacobian matrix of ([Eqs. C.5](#eq0059){ref-type="disp-formula"}--[C.10](#eq0064){ref-type="disp-formula"}) has the block triangular form: $$\mathbf{J} = \begin{bmatrix}
\mathbf{J}_{1} & \mathbf{J}_{2} \\
0 & \mathbf{J}_{3} \\
\end{bmatrix}$$where $$\mathbf{J}_{1} = \left\lbrack \begin{array}{rrr}
{\mu_{0} - D - \left( E + \omega F \right)x_{0}} & {Fx_{0}} & {- Fx_{0}} \\
{\left( G + \omega H \right)x_{1}} & {\mu_{1} - D - \left( G + H \right)x_{1}} & {Hx_{1}} \\
{- \omega Ix_{2}} & {Ix_{2}} & {\mu_{2} - D - Ix_{2}} \\
\end{array} \right\rbrack$$$$\mathbf{J}_{2} = \begin{bmatrix}
{Ex_{0}} & 0 & {Fx_{0}} \\
0 & {Gx_{1}} & {- Hx_{1}} \\
0 & 0 & {Ix_{2}} \\
\end{bmatrix},\quad\mathbf{J}_{3} = \begin{bmatrix}
{- D} & 0 & 0 \\
0 & {- D} & 0 \\
0 & 0 & {- D} \\
\end{bmatrix}$$and $$E = \frac{\partial\mu_{0}}{\partial s_{0}},\quad F = \frac{\partial\mu_{0}}{\partial s_{2}},\quad G = \frac{\partial\mu_{1}}{\partial s_{1}},\quad H = - \frac{\partial\mu_{1}}{\partial s_{2}},\quad I = \frac{d\mu_{2}}{ds_{2}}$$are evaluated at the steady-state.

Since **J** is a block triangular matrix, its eigenvalues are $- D$ (with multiplicity 3) together with the eigenvalues of the 3 × 3 upper-left matrix **J~1~**. Note that we have used the opposite sign for the partial derivative $H = - \partial\mu_{1}/\partial s_{2},$ so that all constants involved in the computation become positive, which will simplify the analysis of the characteristic polynomial of **J~1~**.

### Proof of [Proposition 2](#enun0008){ref-type="statement"} {#sec0034}

Evaluated at SS2, the matrix **J~1~** is: $$\mathbf{J}_{1} = \begin{bmatrix}
{- \left( E + \omega F \right)x_{0}} & {Fx_{0}} & {- Fx_{0}} \\
{\left( G + \omega H \right)x_{1}} & {- \left( G + H \right)x_{1}} & {Hx_{1}} \\
0 & 0 & {\mu_{2} - D} \\
\end{bmatrix}$$Since **J~1~** is a block triangular matrix, its eigenvalues are simply $\mu_{2} - D,$ together with the eigenvalues of the 2 × 2 upper-left matrix. Note that the trace of this 2 × 2 matrix is negative. Hence, its eigenvalues are of negative real part if, and only if, its determinant is positive, that is if, and only if, $$\begin{array}{r}
{E\left( G + H \right) - \left( 1 - \omega \right)FG > 0} \\
\end{array}$$Using $$\frac{\partial M_{0}}{\partial s_{2}} = - \frac{\partial\mu_{0}}{\partial s_{2}}\left\lbrack \frac{\partial\mu_{0}}{\partial s_{0}} \right\rbrack^{- 1} = - F/E$$$$\frac{\partial M_{1}}{\partial s_{2}} = - \frac{\partial\mu_{1}}{\partial s_{2}}\left\lbrack \frac{\partial\mu_{1}}{\partial s_{1}} \right\rbrack^{- 1} = H/G$$we deduce from $$\psi\left( s_{2} \right) = M_{0}\left( D,s_{2} \right) + \frac{M_{1}\left( D,s_{2} \right) + s_{2}}{1 - \omega}$$that $$\frac{\partial\psi}{\partial s_{2}} = \frac{\partial M_{0}}{\partial s_{2}} + \frac{\frac{\partial M_{1}}{\partial s_{2}} + 1}{1 - \omega} = - \frac{F}{E} + \frac{\frac{H}{G} + 1}{1 - \omega}$$Hence, $$\frac{\partial\psi}{\partial s_{2}} = \frac{E\left( G + H \right) - \left( 1 - \omega \right)FG}{\left( 1 - \omega \right)EG}$$Therefore, the condition of stability, ([Eq. C.11](#eq0065){ref-type="disp-formula"}), is equivalent to $\frac{\partial\psi}{\partial s_{2}} > 0$. Hence, we have proven that SS2 is stable if, and only if, *μ*~2~(*s*~2~) \< *D* and $\frac{\partial\psi}{\partial s_{2}} > 0$.

### Proof of [Proposition 3](#enun0009){ref-type="statement"} {#sec0035}

Evaluated at SS3, the matrix **J~1~** is: $$\mathbf{J}_{1} = \begin{bmatrix}
{- \left( E + \omega F \right)x_{0}} & {Fx_{0}} & {- Fx_{0}} \\
{\left( G + \omega H \right)x_{1}} & {- \left( G + H \right)x_{1}} & {Hx_{1}} \\
{- \omega Ix_{2}} & {Ix_{2}} & {- Ix_{2}} \\
\end{bmatrix}$$The characteristic polynomial is given by: $$\begin{array}{r}
{\lambda^{3} + f_{2}\lambda^{2} + f_{1}\lambda + f_{0} = 0} \\
\end{array}$$where $$\begin{aligned}
f_{2} & {= Ix_{2} + \left( G + H \right)x_{1} + \left( E + \omega F \right)x_{0}} \\
\end{aligned}$$$$\begin{aligned}
f_{1} & {= \Delta x_{0}x_{1} + EIx_{0}x_{2} + GIx_{1}x_{2}} \\
\end{aligned}$$$$\begin{aligned}
f_{0} & {= EGIx_{0}x_{1}x_{2}} \\
\end{aligned}$$and $\Delta = E\left( G + H \right) - \left( 1 - \omega \right)FG$.

To satisfy the Routh-Hurwitz criteria, we require *f~i~* \> 0, for $i = 0,1,2$ and $f_{1}f_{2} - f_{0} > 0$. Notice that: $$f_{1}f_{2} - f_{0} = \left( EIx_{0}x_{2} + \Delta x_{0}x_{1} \right)f_{2} + \left( Ix_{2} + \left( G + H \right)x_{1} + \omega Fx_{0} \right)GIx_{1}x_{2}$$We always have *f*~0~ \> 0 and *f*~2~ \> 0.

From ([Eq. C.12](#eq0066){ref-type="disp-formula"}) we deduce that $\Delta = \left( 1 - \omega \right)EG\frac{\partial\psi}{\partial s_{2}}$ . Therefore, if *F*~3~(*D*) ≥ 0, viz. $\frac{\partial\psi}{\partial s_{2}} \geq 0,$ then *Δ* ≥ 0. Hence, *f*~1~ \> 0 and $f_{1}f_{2} - f_{0} > 0,$ so that SS3 is stable.

On the other hand, if $\frac{\partial\psi}{\partial s_{2}} < 0$ and *x*~2~ is very small, which occurs when SS3 is very close to SS2^♭^, then *f*~2~ has the sign of *Δ* since the term with *x*~2~ is negligible compared to the term *Δx*~0~*x*~1~: $$f_{1} = \Delta x_{0}x_{1} + \left( EIx_{0} + GIx_{1} \right)x_{2} < 0$$Thus, SS3 is unstable.

### Proof of [Proposition 4](#enun0010){ref-type="statement"} {#sec0036}

Since we always have *f*~0~ \> 0 and *f*~2~ \> 0, from the previous proof it follows that SS3 is stable if, and only if, $f_{1}f_{2} - f_{0} > 0$. Indeed, this condition implies that we have also *f*~1~ \> *f*~0~/*f*~2~ \> 0. Using ([Eq. 49](#eq0049){ref-type="disp-formula"}) and ([Eq. C.17](#eq0071){ref-type="disp-formula"}), we see that: $$f_{1}f_{2} - f_{0} = F_{4}\left( D,s_{0}^{\text{in}} \right)$$Therefore, the condition $f_{1}f_{2} - f_{0} > 0$ is equivalent to $F_{4}\left( D,s_{0}^{\text{in}} \right) > 0$.

C4. Operating diagrams {#sec0029}
----------------------

### Proof of [Proposition 5](#enun0011){ref-type="statement"} {#sec0037}

We know that SS1 always exists and is stable. We know that SS2^♭^ is unstable if it exists. The curves *Γ*~1~ and *Γ*~2~ of the operating diagram, given by ([Eq. 45](#eq0045){ref-type="disp-formula"}) and ([Eq. 46](#eq0046){ref-type="disp-formula"}), respectively are defined for *D* \< *D*~1~ and *D* \< *D*~2~, respectively and they are tangent for $D = D_{3},$ where $D_{1} = 0.452,$$D_{2} = 0.393$ and $D_{3} = 0.078$. Therefore *F*~3~(*D*) \< 0 for *D* \< *D*~3~ and *F*~3~(*D*) \> 0 for *D* \> *D*~3~. Using [Table 7](#tbl0007){ref-type="table"} and [Remark 2](#enun0007){ref-type="statement"}, we obtain the following results: •$\mathcal{J}_{1}$ is defined by *D* ≥ *D*~1~ or 0 \< *D* \< *D*~1~ and *S*~ch,\ in~ \< *F*~1~(*D*)/*Y*~3~*Y*~4~. Therefore, SS1 is the only existing steady-state in this region.•$\mathcal{J}_{2}$ is defined by *D*~3~ \< *D* \< *D*~1~ and *F*~1~(*D*)/*Y*~3~*Y*~4~ \< *S*~ch,\ in~ \< *F*~2~(*D*)/*Y*~3~*Y*~4~. Therefore, both steady-state SS2 exist and SS2^♯^ is stable since *F*~3~(*D*) \> 0.•$\mathcal{J}_{3}$ is defined by 0 \< *D* \< *D*~2~ and *F*~2~(*D*)/*Y*~3~*Y*~4~ \< *S*~ch,\ in~, and *F*~4~(*D, Y*~3~*Y*~4~*S*~ch,\ in~) \> 0 when 0 \< *D* \< *D*~3~. Therefore, SS3 exists and is stable, both steady-state SS2 exist and SS2^♯^ is unstable since *F*~3~(*D*) \< 0.•$\mathcal{J}_{4}$ is defined by 0 \< *D* \< *D*~3~ and *F*~1~(*D*)/*Y*~3~*Y*~4~ \< *S*~ch,\ in~ \< *F*~2~(*D*)/*Y*~3~*Y*~4~. Therefore, both steady-state SS2 exist and SS2^♯^ is unstable since *F*~3~(*D*) \< 0.•$\mathcal{J}_{5}$ is defined by 0 \< *D* \< *D*~3~, *F*~2~(*D*)/*Y*~3~*Y*~4~ \< *S*~ch,\ in~, and *F*~4~(*D, Y*~3~*Y*~4~*S*~ch,\ in~) \< 0. Therefore, SS3 exists and is unstable and both steady-state SS2 exist and SS2^♯^ is unstable since *F*~3~(*D*) \< 0.

### Proof of [Propositions 6](#enun0012){ref-type="statement"}--[8](#enun0014){ref-type="statement"} {#sec0038}

The result follows from [Table 7](#tbl0007){ref-type="table"} and [Remark 2](#enun0007){ref-type="statement"}. The details are as in the proof of [Proposition 5](#enun0011){ref-type="statement"}.

Matthew J. Wade was funded by the Biotechnology and Biological Sciences Research Council UK (BB/K003240/2 Engineering synthetic microbial communities for biomethane production) and by the Institute of Sustainability, Newcastle University. The authors would like to thank the anonymous reviewers for their valuable and detailed comments.
