# A qualitative study of an anaerobic food-web reveals the importance of
hydrogen for microbial stability

Szymon Sobieszek
Department of Mathematics & Statistics
McMaster University
Hamilton, ON Canada
<EMAIL_ADDRESS>
&Gail S. K. Wolkowicz
Department of Mathematics & Statistics
McMaster University
Hamilton, ON Canada
<EMAIL_ADDRESS>
&Matthew J. Wade
School of Engineering
Newcastle University, Newcastle-upon-Tyne
United Kingdom
<EMAIL_ADDRESS>

###### Abstract

The mathematical analysis of a three-tiered food-web describing anaerobic
chlorophenol mineralisation has suggested the emergence of interesting
dynamical behaviour through its specific ecological interactions, which
include competition, syntrophy and product inhibition. Previous numerical
analyses have revealed the possibility for a Hopf bifurcation occurring
through the interior equilibrium and the role of extraneous substrate inputs
in both mitigating the emergence of periodic solutions and expanding the
desired stable positive steady-state, where full mineralisation occurs. Here
we show that, for a generalised model, the inflow of multiple substrates
results in greater dynamical complexity and prove the occurrence of a
supercritical Hopf bifurcation resulting from variations in these operating
parameters. Further, using numerical estimation, we also show that variations
in the dilution rate can lead to Bogdanov-Takens and Bautin bifurcations.
Finally, we are able to show apply persistence theory for a range of parameter
sets to demonstrate unique persistence in the cases where chlorophenol and
hydrogen are extraneously added to the system, mirroring recent applied
studies highlighting the role of hydrogen in maintaining stable anaerobic
microbial communities.

_K_ eywords Chlorophenol mineralisation; Bifurcation analysis; Hopf
bifurcation; Hydrogen

##  1 Introduction

A mathematical model of the complete anaerobic mineralisation of a generic
monochlorophenol isomer (C6H4ClOH) by a canonical food-web has recently been
described [18]. The system comprises three microbial species (chlorophenol
degrader, phenol degrader, hydrogenotrophic methanogen) whose interactions are
summarised as follows:

  * •

Reductive dehalogenation of chlorophenol in the presence of hydrogen by the
chlorophenol degrader producing phenol [13];

  * •

Phenol is mineralised to acetate and hydrogen via the benzoyl-CoA pathway or a
caproate intermediary [10];

  * •

Anaerobic acetogenic (acetate producing) processes are known to be endergonic
(the reaction results in a net loss of energy to the system). The production
of hydrogen can lead to thermodynamic constraints, or inhibition, if its
partial pressure is high enough. In other words, the reaction becomes
decreasingly exergonic as more hydrogen is produced until it ceases to be
thermodynamically favourable [15, 8]. Hydrogen scavengers such as the
methanogen form a syntrophic partnership with the phenol degrader by
maintaining the hydrogen partial pressure at concentrations low enough for the
mineralisation reaction to proceed;

  * •

Given that the chlorophenol degrader may also act as a syntrophic partner with
the phenol degrader, a competitive interaction between the two hydrogen
utilisers occurs. This positive and negative feedback loops reframes the
ecological network from being a simple food-chain to a more complex food-web
that allows for the possibility of periodicity. This additionally leads to the
possibility of the system reducing to a self-sustaining two-species network in
which the chlorophenol degrader acts as the syntrophic partner to the phenol
degrader.

For reference, the food-web is presented diagrammatically in [18] (Figure 8)
and [14] (Figure 1). The model is a simplified representation of the system at
the population level, ignoring metabolic intermediates and dead-end products
such as methane, which does not contribute to the process dynamics.
Acetoclastic methanogenesis, the conversion of acetate to methane, is also
omitted from the model.

Hydrogen, however, has been shown to play an important role in stability of
anaerobic microbial communities through the effects of inhibition and
competition [1, 4, 17]. Given that external hydrogen addition will maintain
the methanogen population (no washout when the methanogen growth rate is
greater than the combined dilution and decay/maintenance rates) under a wider
operating parameter regime (chlorophenol inflow and dilution rate) [18], a
global analysis of the model can provide deeper insights into the ecological
role of hydrogen through its association with community stability and
criticality of the Hopf bifurcation.

Here, we focus on the mathematical analysis of the model, extending the work
reported in the literature. For example, an analytical approach was taken to
characterise the existence and stability of the system equilibria with and
without inclusion of a microbial decay term using a general representation of
the species growth functions [14]. With no decay, local stability and the
conditions giving rise to asymptotic coexistence of all three species have
been shown analytically, where the possibility of periodic orbits are also not
excluded [6]. However, numerical analysis has suggested the presence of a Hopf
bifurcation emerging through the positive steady-state, with the concentration
of influent chlorophenol as the bifurcating parameter [14].

In this work, we extend the analysis of the model providing a proof of the
existence, uniqueness and stability of six identified equilibria with the
addition of only chlorophenol, and the more general case where all three
substrates as external inputs to the system. The procedure we use allows to
identify the sufficient conditions for the emergence of a Hopf bifurcation in
those inflow concentration parameters. We are also able to prove that the
dimensionless model is uniformly persistent, a new result for the system.

##  2 The model revisited

We first present concisely the original chemostat model using identical
scaling to that given by [18]

$\displaystyle\begin{split}\begin{cases}x_{0}^{\prime}&=-\alpha
x_{0}+\mu_{0}\left(s_{0},s_{2}\right)x_{0}-k_{{}_{A}}x_{0},\\\
x_{1}^{\prime}&=-\alpha
x_{1}+\mu_{1}\left(s_{1},s_{2}\right)x_{1}-k_{{}_{B}}x_{1},\\\
x_{2}^{\prime}&=-\alpha
x_{2}+\mu_{2}\left(s_{2}\right)x_{2}-k_{{}_{C}}x_{2},\\\
s_{0}^{\prime}&=\alpha\left(u_{f}-s_{0}\right)-\mu_{0}\left(s_{0},s_{2}\right)x_{0},\\\
s_{1}^{\prime}&=\alpha(u_{g}-s_{1})+\omega_{0}\mu_{0}\left(s_{0},s_{2}\right)x_{0}-\mu_{1}\left(s_{1},s_{2}\right)x_{1},\\\
s_{2}^{\prime}&=\alpha(u_{h}-s_{2})-\omega_{2}\mu_{0}\left(s_{0},s_{2}\right)x_{0}+\omega_{1}\mu_{1}\left(s_{1},s_{2}\right)x_{1}-\mu_{2}\left(s_{2}\right)x_{2},\end{cases}\end{split}$
(1)

with

$\displaystyle\alpha$ $\displaystyle=\frac{D}{k_{\text{m,ch}}Y_{\text{ch}}},$
(2) $\displaystyle u_{f}$
$\displaystyle=\frac{S_{\text{ch,in}}}{K_{\text{S,ch}}},\quad
u_{g}=\frac{S_{\text{ph,in}}}{K_{\text{S,ph}}},\quad
u_{h}=\frac{S_{\text{H${}_{2}$,in}}}{K_{\text{S,H${}_{2}$}}},$ (3)
$\displaystyle\omega_{0}$
$\displaystyle=\frac{K_{\text{S,ch}}}{K_{\text{S,ph}}}\frac{224}{208}\left(1-Y_{\text{ch}}\right),\quad\omega_{1}=\frac{K_{\text{S,ph}}}{K_{\text{S,H${}_{2}$}}}\frac{32}{224}\left(1-Y_{\text{ph}}\right),\quad\omega_{2}=\frac{16}{208}\frac{K_{\text{S,ch}}}{K_{\text{S,H${}_{2}$}}},$
(4) $\displaystyle\phi_{1}$
$\displaystyle=\frac{k_{\text{m,ph}}Y_{\text{ph}}}{k_{\text{m,ch}}Y_{\text{ch}}},\quad\phi_{2}=\frac{k_{\text{m,H${}_{2}$}}Y_{\text{H${}_{2}$}}}{k_{\text{m,ch}}Y_{\text{ch}}},$
(5) $\displaystyle K_{P}$
$\displaystyle=\frac{K_{\text{S,H${}_{2}$,c}}}{K_{\text{S,H${}_{2}$}}},\quad
K_{I}=\frac{K_{\text{S,H${}_{2}$}}}{K_{I,\text{H${}_{2}$}}},$ (6)
$\displaystyle k_{{}_{A}}$
$\displaystyle=\frac{k_{\text{dec,ch}}}{k_{\text{m,ch}}Y_{\text{ch}}},\quad
k_{{}_{B}}=\frac{k_{\text{dec,ph}}}{k_{\text{m,ch}}Y_{\text{ch}}},\quad
k_{{}_{C}}=\frac{k_{\text{dec,H${}_{2}$}}}{k_{\text{m,ch}}Y_{\text{ch}}}.$ (7)

$\displaystyle\mu_{0}\left(s_{0},s_{2}\right)$
$\displaystyle=\frac{s_{0}}{1+s_{0}}\frac{s_{2}}{K_{P}+s_{2}},\quad\mu_{1}\left(s_{1},s_{2}\right)=\frac{\phi_{1}s_{1}}{1+s_{1}}\frac{1}{1+K_{I}s_{2}},\quad\mu_{2}\left(s_{2}\right)=\frac{\phi_{2}s_{2}}{1+s_{2}},$
(8)

where $\alpha$ is the dilution rate, $u_{f},u_{g},u_{h}$ are the chlorophenol,
phenol and hydrogen inflow concentrations, respectively, and
$k_{A},k_{B},k_{C}$ are the decay (or maintenance) terms. These are scaled to
be dimensionless, as are the other parameters using the scaling provided by
[18]. Briefly, $k_{\mathrm{m,\cdot}}$ are the specific growth rates,
$K_{\mathrm{S,\cdot}}$ are the half-saturation coefficients,
$Y_{\mathrm{\cdot}}$ are the substrate yield coefficients,
$K_{\mathrm{I,H_{2}}}$ is the kinetic inhibition constant of hydrogen on the
phenol degrader, and $k_{\mathrm{dec,\cdot}}$ are the unscaled decay terms.
The numeric values indicate the stoichiometric coefficients given in terms of
units of Chemical Oxygen Demand rather than molarity, as is common for
environmental engineering models. The $\mu_{n}(\cdot),n={0,1,2}$ are the
species growth functions described by double Monod, Monod with product
inhibition, and Monod kinetics, respectively. Subscripts
${}_{\mathrm{ch}},_{\mathrm{ph}},_{\mathrm{H_{2}}}$ relate to the chlorophenol
degrader, phenol degrader, and methanogen, respectively.

For numerical bifurcation analysis given in Section 4.4, we consider the same
parameter values as provided in the original work, as shown in Table 1.

Parameters | Value
---|---
$\omega_{0}$ | 0.1854
$\omega_{1}$ | 1656.69
$\omega_{2}$ | 163.08
$\phi_{1}$ | 1.8875
$\phi_{2}$ | 3.8113
$K_{P}$ | 0.04
$K_{I}$ | 7.1429
Table 1: Parameter regimes for the system (1).

##  3 Reduction of the model

We are able to obtain many theoretical results assuming general forms of the
growth functions provided we assume the death rates of the microbial
populations are insignificant compared to the dilution rate. We thus consider
the following system that is identical to system (1), except that we assume
$k_{i}=0$, $i\in\\{A,B,C\\}$:

$\displaystyle\begin{split}\begin{cases}x_{0}^{\prime}&=-\alpha
x_{0}+\mu_{0}\left(s_{0},s_{2}\right)x_{0},\\\ x_{1}^{\prime}&=-\alpha
x_{1}+\mu_{1}\left(s_{1},s_{2}\right)x_{1},\\\ x_{2}^{\prime}&=-\alpha
x_{2}+\mu_{2}\left(s_{2}\right)x_{2},\\\
s_{0}^{\prime}&=\alpha\left(u_{f}-s_{0}\right)-\mu_{0}\left(s_{0},s_{2}\right)x_{0},\\\
s_{1}^{\prime}&=\alpha(u_{g}-s_{1})+\omega_{0}\mu_{0}\left(s_{0},s_{2}\right)x_{0}-\mu_{1}\left(s_{1},s_{2}\right)x_{1},\\\
s_{2}^{\prime}&=\alpha(u_{h}-s_{2})-\omega_{2}\mu_{0}\left(s_{0},s_{2}\right)x_{0}+\omega_{1}\mu_{1}\left(s_{1},s_{2}\right)x_{1}-\mu_{2}\left(s_{2}\right)x_{2},\end{cases}\end{split}$
(9) $x_{i}\left(0\right)\geq 0,\quad s_{i}\left(0\right)\geq 0,\quad
i\in\left\\{1,2,3\right\\}.$

We assume that $\mu_{0}\left(s_{0},s_{2}\right)$,
$\mu_{1}\left(s_{1},s_{2}\right)$, $\mu_{2}\left(s_{2}\right)$ are
$\mathcal{C}^{1}$ functions that satisfy the following general conditions:

  * •

For all $s_{0}\geq 0$ and $s_{2}\geq 0$, $\mu_{0}\left(0,s_{2}\right)=0$,
$\mu_{0}\left(s_{0},0\right)=0$. As a consequence,
$\partial_{s_{0}}\mu_{0}\left(s_{0},0\right)=0$,
$\partial_{s_{2}}\mu_{0}\left(0,s_{2}\right)=0$. Thus we assume that the
chlorophenol degrader cannot grow in the absence of either chlorophenol or
hydrogen;

  * •

For all $s_{0}>0$ and $s_{2}>0$,
$\partial_{s_{0}}\mu_{0}\left(s_{0},s_{2}\right)>0$,
$\partial_{s_{2}}\mu_{0}\left(s_{0},s_{2}\right)>0$. Thus we assume that the
chlorophenol degrader grows on both chlorophenol and hydrogen;

  * •

For all $s_{2}\geq 0$ and $s_{1}\geq 0$, $\mu_{1}\left(0,s_{2}\right)=0$,
$\partial_{s_{2}}\mu_{1}\left(0,s_{2}\right)=0$. Thus we assume that the
phenol degrader cannot grow in the absence of phenol;

  * •

For all $s_{1}>0$ and $s_{2}>0$,
$\partial_{s_{1}}\mu_{1}\left(s_{1},s_{2}\right)>0$,
$\partial_{s_{2}}\mu_{1}\left(s_{1},s_{2}\right)<0$. Thus we assume that the
supply of phenol results in growth of the phenol degrader, and that hydrogen
inhibits its growth;

  * •

$\mu_{2}\left(0\right)=0$ and $\mu_{2}^{\prime}\left(s_{2}\right)>0$ for all
$s_{2}>0$. Thus we assume that the mathanogen cannot grow without the presence
of hydrogen, and that increasing the supply of hydrogen results in faster
growth of the methanogen.

We use the prototypes $\mu_{0}$, $\mu_{1}$, and $\mu_{2}$ defined in (8),
which satisfy these conditions, when we are able to prove results in general,
and when providing numerical simulations or bifurcation diagrams.

We now prove a lemma that we will use to show global well-posedness of system
(9).

######  Lemma 3.1.

All solutions of system (9) with positive initial conditions remain positive
and bounded for all positive times.

If $x_{i}\left(0\right)=0$, $i\in\left\\{1,2,3\right\\}$, then
$x_{i}\left(t\right)=0$ for all $t\geq 0$.

###### Proof.

Consider any solution $\vec{\varphi}\left(t\right)$ with positive initial
conditions. By existence and uniqueness theory, there cannot be a time
$\bar{t}>0$, such that $x_{i}\left(\bar{t}\right)=0$ for some
$i\in\left\\{1,2,3\right\\}$, since then $x_{i}\left(t\right)\equiv 0$ for all
$t\in\mathbb{R}$, contradicting $x_{i}\left(0\right)>0$. Hence
$x_{i}\left(t\right)>0$ for all $t\geq 0$. Also, if $x_{i}\left(0\right)=0$
for some $i\in\left\\{1,2,3\right\\}$, then there is a solution of system (9)
with $x_{i}\left(t\right)\equiv 0$ for all $t\in\mathbb{R}$. By existence and
uniqueness theory, this is the only solution.

Now, consider $\vec{\varphi}\left(t\right)$, and suppose that there is some
$\bar{t}>0$, such that $s_{0}\left(t\right)>0$ for
$t\in\left[0,\bar{t}\right)$, $s_{0}\left(\bar{t}\right)=0$, and
$s_{1}\left(t\right),s_{2}\left(t\right)\geq 0$ for
$t\in\left[0,\bar{t}\right]$. Then $s_{0}^{\prime}\left(\bar{t}\right)\leq 0$.
However, from system (9), $s_{0}^{\prime}\left(\bar{t}\right)=\alpha u_{f}$.
If $u_{f}>0$, then $s_{0}^{\prime}\left(\bar{t}\right)>0$, a contradiction. If
$u_{f}=0$, then there is a solution of system (9) with
$s_{0}\left(t\right)\equiv 0$, which contradicts uniqueness of solutions. It
follows that $s_{0}\left(t\right)>0$ for all $t\geq 0$.

Next, consider $\vec{\varphi}\left(t\right)$, and suppose that there is some
$\bar{t}>0$, such that $s_{2}\left(t\right)>0$ for
$t\in\left[0,\bar{t}\right)$, $s_{2}\left(\bar{t}\right)=0$, and
$s_{1}\left(t\right)\geq 0$ for $t\in\left[0,\bar{t}\right]$. Then
$s_{2}^{\prime}\left(\bar{t}\right)\leq 0$. However, from system (9),
$s_{2}^{\prime}\left(\bar{t}\right)=\alpha
u_{h}+\omega_{1}\mu_{1}\left(s_{1}\left(\bar{t}\right),0\right)x_{1}$. If
$u_{h}>0$, or $s_{1}\left(\bar{t}\right)>0$, then
$s_{2}^{\prime}\left(\bar{t}\right)>0$, a contradiction. If both $u_{h}=0$,
and $s_{1}\left(\bar{t}\right)=0$, then
$s_{2}^{\prime}\left(\bar{t}\right)=0$, and there is another solution with
$s_{2}\left(t\right)\equiv 0$, which contradicts uniqueness of solutions to
initial value problems. It follows that $s_{2}\left(t\right)>0$ for all $t\geq
0$.

Finally, consider $\vec{\varphi}\left(t\right)$, and suppose that there is
some $\bar{t}>0$, such that $s_{1}\left(t\right)>0$ for
$t\in\left[0,\bar{t}\right)$, $s_{1}\left(\bar{t}\right)=0$. Then
$s_{1}^{\prime}\left(\bar{t}\right)\leq 0$. However, from system (9),
$s_{1}^{\prime}\left(\bar{t}\right)=\alpha
u_{g}+\omega_{0}\mu_{0}\left(s_{0}\left(\bar{t}\right),s_{2}\left(\bar{t}\right)\right)x_{0}$,
so $s_{1}^{\prime}\left(\bar{t}\right)>0$, a contradiction. It follows that
$s_{1}\left(t\right)>0$ for all $t\geq 0$.

We have thus proved the positivity of solutions and move on to showing the
boundedness of solutions.

By adding the first and the fourth equations of (9), we obtain

$x_{0}^{\prime}+s_{0}^{\prime}=-\alpha\left(x_{0}+s_{0}-u_{f}\right),$

hence

$\left(x_{0}+s_{0}-u_{f}\right)^{\prime}=-\alpha\left(x_{0}+s_{0}-u_{f}\right),$

which implies that

$x_{0}(t)+s_{0}(t)=u_{f}+\left(x_{0}(0)+s_{0}(0)-u_{f}\right)e^{-\alpha t}.$
(10)

Similarly, we obtain

$x_{1}(t)+\omega_{0}s_{0}(t)+s_{1}(t)=\omega_{0}u_{f}+u_{g}+\left(x_{1}(0)+\omega_{0}s_{0}(0)+s_{1}(0)-\omega_{0}u_{f}-u_{g}\right)e^{-\alpha
t},$ (11)

and

$\omega_{2}x_{0}(t)+x_{2}(t)+\omega_{0}\omega_{1}s_{0}(t)+\omega_{1}s_{1}(t)+s_{2}(t)=\omega_{0}\omega_{1}u_{f}+\omega_{1}u_{g}+u_{h}+\\\
\left(\omega_{0}x_{0}(0)+x_{2}(0)+\omega_{0}\omega_{1}s_{0}(0)+\omega_{1}s_{1}(0)+s_{2}(0)-\omega_{0}\omega_{1}u_{f}-\omega_{1}u_{g}-u_{h}\right)e^{-\alpha
t}.$ (12)

Since all terms of the sums in (10), (11) and (12) are positive for all
positive initial conditions, the solutions of (9) are bounded. Also, taking
the limit as $t\to\infty$ in equations (10), (11), and (12) we obtain that

$\displaystyle\lim_{t\to\infty}\left(x_{0}\left(t\right)+s_{0}\left(t\right)\right)$
$\displaystyle=u_{f},$ (13)
$\displaystyle\lim_{t\to\infty}\left(x_{1}\left(t\right)+\omega_{0}s_{0}\left(t\right)+s_{1}\left(t\right)\right)$
$\displaystyle=\omega_{0}u_{f}+u_{g},$ (14)
$\displaystyle\lim_{t\to\infty}\left(\omega_{2}x_{0}\left(t\right)+x_{2}\left(t\right)+\omega_{0}\omega_{1}s_{0}\left(t\right)+\omega_{1}s_{1}\left(t\right)+s_{2}\left(t\right)\right)$
$\displaystyle=\omega_{0}\omega_{1}u_{f}+\omega_{1}u_{g}+u_{h}.$ (15)

Starting with any positive initial conditions, the solutions of system (9)
eventually satisfy

$\displaystyle x_{0}+s_{0}$ $\displaystyle=u_{f},$ (16) $\displaystyle
x_{1}+\omega_{0}s_{0}+s_{1}$ $\displaystyle=\omega_{0}u_{f}+u_{g},$ (17)
$\displaystyle\omega_{2}x_{0}+x_{2}+\omega_{0}\omega_{1}s_{0}+\omega_{1}s_{1}+s_{2}$
$\displaystyle=\omega_{0}\omega_{1}u_{f}+\omega_{1}u_{g}+u_{h}.$ (18)

∎

We call relations (16)-(18) "conservation principles". In other words, system
(9) admits a positively invariant attracting set
$\Omega\subset\mathbb{R}^{6}$, such that

$\displaystyle\begin{split}\Omega=\\{\left(x_{0},x_{1},x_{2},s_{0},s_{1},s_{2}\right)\in\mathbb{R}^{6}:x_{i},s_{i}&\geq
0,i=0,1,2;\\\ x_{0}+s_{0}&=u_{f},\\\
x_{1}+\omega_{0}s_{0}+s_{1}&=\omega_{0}u_{f}+u_{g},\\\
\omega_{2}x_{0}+x_{2}+\omega_{0}\omega_{1}s_{0}+\omega_{1}s_{1}+s_{2}&=\omega_{0}\omega_{1}u_{f}+\omega_{1}u_{g}+u_{h}\\}.\end{split}$
(19)

Using the conservation principles we can compute $s_{0}$, $s_{1}$, and $s_{2}$
as functions of $x_{0}$, $x_{1}$, $x_{2}$

$\displaystyle\begin{split}s_{0}&=-x_{0}+u_{f},\\\
s_{1}&=\omega_{0}x_{0}-x_{1}+u_{g},\\\
s_{2}&=-\omega_{2}x_{0}+\omega_{1}x_{1}-x_{2}+u_{h}.\end{split}$ (20)

Now, we can reduce the analysis of the original system (9) to the analysis of
the following equivalent three-dimensional system on the invariant set
$\Omega$

$\displaystyle\begin{cases}\begin{split}x_{0}^{\prime}&=-\alpha
x_{0}+\mu_{0}\left(-x_{0}+u_{f},-\omega_{2}x_{0}+\omega_{1}x_{1}-x_{2}+u_{h}\right)x_{0},\\\
x_{1}^{\prime}&=-\alpha
x_{1}+\mu_{1}\left(\omega_{0}x_{0}-x_{1}+u_{g},-\omega_{2}x_{0}+\omega_{1}x_{1}-x_{2}+u_{h}\right)x_{1},\\\
x_{2}^{\prime}&=-\alpha
x_{2}+\mu_{2}\left(-\omega_{2}x_{0}+\omega_{1}x_{1}-x_{2}+u_{h}\right)x_{2}.\end{split}\end{cases}$
(21)

From now on, we will study the reduced system (21). We begin by analyzing all
the possible equilibria.

###  3.1 Equilibria of the reduced system and their local stability

The equilibria are found by setting the right hand sides of equations in (21)
equal to zero. Below, we list all the possibilities obtained this way. Since
equations (20) give a one-to-one correspondence of the equilibria of system
(21) with the equilibria of system (9), we also list the corresponding steady
states $\left(x_{0},x_{1},x_{2},s_{0},s_{1},s_{2}\right)$ of the six-
dimensional system in each case.

Types of equilibria of system (21):

  * •

Zero equilibrium ${}_{\left(000\right)}\mathcal{E}=\left(0,0,0\right)$. The
corresponding equilibrium $\prescript{}{\left(000\right)}{\mathbb{E}}$ in the
six-dimensional system:

$\prescript{}{\left(000\right)}{\mathbb{E}}=\left(0,0,0,u_{f},u_{g},u_{h}\right).$
(22)

In this case, all the populations die, hence the only source for the
substrates comes from the inflow rates $u_{f}$, $u_{g}$, and $u_{h}$.

  * •

Boundary equilibria:

    * –

${}_{\left(100\right)}\mathcal{E}=\left(\prescript{}{\left(100\right)}{x_{0}},0,0\right)$,
where $x_{0}=\prescript{}{\left(100\right)}{x_{0}}>0$ is a solution (if it
exists) of

$\mu_{0}(-x_{0}+u_{f},-\omega_{2}x_{0}+u_{h})=\alpha.$ (23)

The corresponding equilibrium $\prescript{}{\left(100\right)}{\mathbb{E}}$ in
the six-dimensional system:

$\prescript{}{\left(100\right)}{\mathbb{E}}=\left(\prescript{}{\left(100\right)}{x_{0}},0,0,-\prescript{}{\left(100\right)}{x_{0}}+u_{f},\omega_{0}\prescript{}{\left(100\right)}{x_{0}}+u_{g},-\omega_{2}\prescript{}{\left(100\right)}{x_{0}}+u_{h}\right).$
(24)

In this case, the only microorganism surviving is the chlorophenol degrader.
It consumes the chlorophenol, hence the value of $s_{0}$ is given as the
balance between this consumption, and the supply inflow $u_{f}$. Since $x_{0}$
produces phenol this value is added to $u_{g}$ in the total phenol amount
$s_{1}$. Since $x_{0}$ consumes hydrogen as well, the value
$\omega_{2}\prescript{}{\left(100\right)}{x_{0}}$ is subtracted from $s_{2}$
as well. This steady state is not desirable because of the phenol build-up in
the system.

    * –

${}_{\left(010\right)}\mathcal{E}=\left(0,\prescript{}{\left(010\right)}{x_{1}},0\right)$,
where $x_{1}=\prescript{}{\left(010\right)}{x_{1}}>0$ is a solution (if it
exists) of

$\mu_{1}(-x_{1}+u_{g},\omega_{1}x_{1}+u_{h})=\alpha.$ (25)

The corresponding equilibrium $\prescript{}{\left(010\right)}{\mathbb{E}}$ in
the six-dimensional system:

$\prescript{}{\left(010\right)}{\mathbb{E}}=\left(0,\prescript{}{\left(010\right)}{x_{1}},0,u_{f},-\prescript{}{\left(010\right)}{x_{1}}+u_{g},\omega_{1}\prescript{}{\left(010\right)}{x_{1}}+u_{h}\right).$
(26)

In this case, only the phenol degrader survives, and hence the value of
$s_{1}$ at the equilibrium is equal to the balance between its consumption and
inflow $u_{g}$. Chlorophenol is not being consumed, hence its total amount
equals the inflow concentration $u_{f}$. Hydrogen is being produced by the
phenol degrader, and also its value is increased by the inflow $u_{h}$.

    * –

${}_{\left(001\right)}\mathcal{E}=\left(0,0,\prescript{}{\left(001\right)}{x_{2}}\right)$,
where $x_{2}=\prescript{}{\left(001\right)}{x_{2}}>0$ is a solution (if it
exists) of

$\mu_{2}(-x_{2}+u_{h})=\alpha.$ (27)

The corresponding equilibrium $\prescript{}{\left(001\right)}{\mathbb{E}}$ in
the six-dimensional system:

$\prescript{}{\left(001\right)}{\mathbb{E}}=\left(0,0,\prescript{}{\left(001\right)}{x_{2}},u_{f},u_{g},-\prescript{}{\left(001\right)}{x_{2}}+u_{h}\right).$
(28)

Here, only the methanogen is present, hence the values of chlotophenol and
phenol are equal the inflow concentrations $u_{f}$ and $u_{g}$, respectively.

    * –

${}_{\left(101\right)}\mathcal{E}=\left(\prescript{}{\left(101\right)}{x_{0}},0,-\omega_{2}\prescript{}{\left(101\right)}{x_{0}}+u_{h}-\mu_{2}^{-1}\left(\alpha\right)\right)$,
where $x_{0}=\prescript{}{\left(101\right)}{x_{0}}>0$ is a solution (if it
exists) of

$\mu_{0}(-x_{0}+u_{f},\mu_{2}^{-1}(\alpha))=\alpha.$ (29)

The corresponding equilibrium $\prescript{}{\left(101\right)}{\mathbb{E}}$ in
the six-dimensional system:

$\prescript{}{\left(101\right)}{\mathbb{E}}=\Big{(}\prescript{}{\left(101\right)}{x_{0}},0,-\omega_{2}\prescript{}{\left(101\right)}{x_{0}}+u_{h}-\mu_{2}^{-1}\left(\alpha\right),-\prescript{}{\left(101\right)}{x_{0}}+u_{f},\\\
\omega_{0}\prescript{}{\left(101\right)}{x_{0}}+u_{g},\mu_{2}^{\left(-1\right)}\left(\alpha\right)\Big{)}.$
(30)

In this case, both chlorophenol degrader and methanogen are present. The lack
of phenol degrader results in phenol build-up. We can also observe competition
for hydrogen between the phenol degrader and methanogen.

    * –

${}_{\left(011\right)}\mathcal{E}=\left(0,\prescript{}{\left(011\right)}{x_{1}},\omega_{1}\prescript{}{\left(011\right)}{x_{1}}+u_{h}-\mu_{2}^{-1}\left(\alpha\right)\right)$,
where $x_{1}=\prescript{}{\left(011\right)}{x_{1}}>0$ is a solution (if it
exists) of

$\mu_{1}\left(-x_{1}+u_{g},\mu_{2}^{-1}\left(\alpha\right)\right)=\alpha.$
(31)

The corresponding equilibrium $\prescript{}{\left(011\right)}{\mathbb{E}}$ in
the six-dimensional system:

$\prescript{}{\left(011\right)}{\mathbb{E}}=\left(0,\prescript{}{\left(011\right)}{x_{1}},\omega_{1}\prescript{}{\left(011\right)}{x_{1}}+u_{h}-\mu_{2}^{-1}\left(\alpha\right),u_{f},-\prescript{}{\left(011\right)}{x_{1}}+u_{g},\mu_{2}^{-1}\left(\alpha\right)\right).$
(32)

This steady state represents a two-tiered food chain, with the phenol degrader
and methanogen present. Hydrogen has an inhibiting effect on the phenol
degrader.

    * –

${}_{\left(110\right)}\mathcal{E}=\left(\prescript{}{\left(110\right)}{x_{0}},\prescript{}{\left(110\right)}{x_{1}},0\right)$,
where $x_{0}=\prescript{}{\left(110\right)}{x_{0}}>0$ and
$x_{1}=\prescript{}{\left(110\right)}{x_{1}}>0$ are solutions of

$\displaystyle\begin{split}\mu_{0}\left(-x_{0}+u_{f},-\omega_{2}x_{0}+\omega_{1}x_{1}+u_{h}\right)&=\alpha,\\\
\mu_{1}\left(\omega_{0}x_{0}-x_{1}+u_{g},-\omega_{2}x_{0}+\omega_{1}x_{1}+u_{h}\right)&=\alpha.\end{split}$
(33)

The corresponding equilibrium $\prescript{}{\left(110\right)}{\mathbb{E}}$ in
the six-dimensional system:

$\prescript{}{\left(110\right)}{\mathbb{E}}=\Big{(}\prescript{}{\left(110\right)}{x_{0}},\prescript{}{\left(110\right)}{x_{1}},0,-\prescript{}{\left(110\right)}{x_{0}}+u_{f},\omega_{0}\prescript{}{\left(110\right)}{x_{0}}-\prescript{}{\left(110\right)}{x_{1}}+u_{g},\\\
-\omega_{2}\prescript{}{\left(110\right)}{x_{0}}+\omega_{1}\prescript{}{\left(110\right)}{x_{1}}+u_{h}\Big{)}.$
(34)

In this case, both the chlorophenol and phenol degraders are present, however
the methanogen is washed out. Thus, full mineralisation to methane is not
possible and, hence, the hydrogen accumulates to some theoretical maximum,
balanced such that the the inhibitory effect on the phenol degrader does not
induce washout, whilst providing enough hydrogen for chlorophenol degrader
activity.

  * •

Positive (interior) equilibrium
${}_{\left(111\right)}\mathcal{E}=\left(\overset{*}{x}_{0},\overset{*}{x}_{1},\overset{*}{x}_{2}\right)$,
where $x_{0}=\overset{*}{x}_{0}>0$, $x_{1}=\overset{*}{x}_{1}>0$, and
$x_{2}=\overset{*}{x}_{2}>0$ are solutions of

$\displaystyle\begin{split}\mu_{0}\left(-x_{0}+u_{f},-\omega_{2}x_{0}+\omega_{1}x_{1}-x_{2}+u_{h}\right)&=\alpha,\\\
\mu_{1}\left(\omega_{0}x_{0}-x_{1}+u_{g},-\omega_{2}x_{0}+\omega_{1}x_{1}-x_{2}+u_{h}\right)&=\alpha,\\\
\mu_{2}\left(-\omega_{2}x_{0}+\omega_{1}x_{1}-x_{2}+u_{h}\right)&=\alpha.\end{split}$
(35)

The corresponding equilibrium $\prescript{}{\left(111\right)}{\mathbb{E}}$ in
the six-dimensional system:

$\prescript{}{\left(111\right)}{\mathbb{E}}=\left(\overset{*}{x}_{0},\overset{*}{x}_{1},\overset{*}{x}_{2},-\overset{*}{x}_{0}+u_{f},\omega_{0}\overset{*}{x}_{0}-\overset{*}{x}_{1}+u_{g},\mu_{2}^{\left(-1\right)}\left(\alpha\right)\right).$
(36)

Here, all species are present, and thus we observe full chlorophenol
mineralisation. For this reason, asymptotic stability of this equilibrium is
the desired operational state.

As it is not clear whether the listed equations have solutions, and if the
solutions are unique, we now derive conditions on the parameters that address
these questions.

####  3.1.1 Existence and uniqueness

Since there are many parameters in system (21), it was not possible to obtain
explicit expressions for some of the equilibria. We did however simplify the
computation by only looking for equilibria in the invariant set $\Omega$. This
assumption is reasonable, since the dynamics of the original system reduces to
the dynamics on the set $\Omega$.

  * •

${}_{\left(000\right)}\mathcal{E}=\left(0,0,0\right)$ equilibrium always
exists.

  * •

${}_{\left(100\right)}\mathcal{E}=\left(\prescript{}{\left(100\right)}{x_{0}},0,0\right)$.
As mentioned at the beginning of subsection 3.1.1, we are looking for the
equilibria in the feasible set $\Omega$, i.e., where all of the components in
corresponding six-dimensional equilibria are nonnegative. Thus, we want
$x_{0}=\prescript{}{\left(100\right)}{x_{0}}$ to satisfy
$\prescript{}{\left(100\right)}{x_{0}}\in\left(0,u_{f}\right]$, and
$\prescript{}{\left(100\right)}{x_{0}}\leq\frac{u_{h}}{\omega_{2}}$; hence we
consider only
$x_{0}\in\left(0,\min{\left(u_{f},\frac{u_{h}}{\omega_{2}}\right)}\right]$.
For such $x_{0}$ the differentiable mapping
$x_{0}\mapsto\mu_{0}\left(u_{f}-x_{0},-\omega_{2}x_{0}+u_{h}\right)$ is
decreasing, and thus ${}_{\left(100\right)}\mathcal{E}$ exists in $\Omega$ if
and only if $\alpha\in\left[0,\mu_{0}(u_{f},u_{h})\right)$, and when it
exists, it is unique.

  * •

${}_{\left(010\right)}\mathcal{E}=\left(0,\prescript{}{\left(010\right)}{x_{1}},0\right)$.
By a similar argument, we consider $x_{1}\in\left(0,u_{g}\right]$. The
differentiable mapping
$x_{1}\mapsto\mu_{1}(u_{g}-x_{1},\omega_{1}x_{1}+u_{h})$ is decreasing, so
${}_{\left(010\right)}\mathcal{E}$ exists in $\Omega$ if and only if
$\alpha\in\left[0,\mu_{1}(u_{g},u_{h})\right)$, and when it exists, it is
unique.

  * •

${}_{\left(001\right)}\mathcal{E}=\left(0,0,\prescript{}{\left(001\right)}{x_{2}}\right)$.
Once again, we consider only $x_{2}\in(0,u_{h}]$, for which the differentiable
mapping $x_{2}\mapsto\mu_{2}(-x_{2}+u_{h})$ is decreasing, so
${}_{\left(001\right)}\mathcal{E}$ exists in $\Omega$ if and only if
$\alpha\in[0,\mu_{2}(u_{h}))$, and when it exists, it is unique (notice that
if $\mu_{2}(s_{2})=\frac{\phi_{2}s_{2}}{1+s_{2}}$, we have
$\prescript{}{\left(001\right)}{x_{2}}=u_{h}-\frac{\alpha}{\phi_{2}-\alpha}$).

  * •

${}_{\left(101\right)}\mathcal{E}=\left(\prescript{}{\left(101\right)}{x_{0}},0,-\omega_{2}\prescript{}{\left(101\right)}{x_{0}}+u_{h}-\mu_{2}^{-1}\left(\alpha\right)\right)$.
For
$\prescript{}{\left(101\right)}{x_{2}}=-\omega_{2}\prescript{}{\left(101\right)}{x_{0}}+u_{h}-\mu_{2}^{-1}\left(\alpha\right)$,
the restriction $\prescript{}{\left(101\right)}{x_{2}}>0$ gives us the
condition
$\prescript{}{\left(101\right)}{x_{0}}<\frac{u_{h}-\mu_{2}^{-1}(\alpha)}{\omega_{2}}$
(notice that this already implies that
$\prescript{}{\left(101\right)}{x_{0}}\leq\frac{u_{h}}{\omega_{2}}$ is
satisfied), and thus the requirement $\prescript{}{\left(101\right)}{x_{0}}>0$
results in the first condition on $\alpha$, i.e.,
$\alpha<\mu_{2}\left(u_{h}\right)$. We also require that
$\prescript{}{\left(101\right)}{x_{0}}\leq u_{f}$. We therefore only consider
the mapping $x_{0}\mapsto\mu_{0}(-x_{0}+u_{f},\mu_{2}^{-1}(\alpha))$ for
$x_{0}\in\left(0,\min{\left(u_{f},\frac{u_{h}-\mu_{2}^{-1}\left(\alpha\right)}{\omega_{2}}\right)}\right)$.
This differentiable mapping is decreasing, so
${}_{\left(101\right)}\mathcal{E}$ exists in $\Omega$ if and only if

$\alpha\in\left(\mu_{0}\left(u_{f}-\min{\left(u_{f},\frac{u_{h}-\mu_{2}^{-1}\left(\alpha\right)}{\omega_{2}}\right)},\mu_{2}^{-1}\left(\alpha\right)\right),\mu_{0}\left(u_{f},\mu_{2}^{-1}\left(\alpha\right)\right)\right),$
(37)

and $\alpha<\mu_{2}\left(u_{h}\right)$, and when it exists, it is unique.

By solving the equation

$\mu_{0}(-x_{0}+u_{f},\mu_{2}^{-1}(\alpha))=\alpha,$ (38)

we obtain the following explicit formulas for
$\prescript{}{\left(101\right)}{x_{0}}$ and
$\prescript{}{\left(101\right)}{x_{2}}$ with our test prototypes $\mu_{0}$,
$\mu_{1}$, and $\mu_{2}$:

$\displaystyle\prescript{}{\left(101\right)}{x_{0}}$
$\displaystyle=\frac{\alpha\left(1+u_{f}\right)\left(K_{P}+\mu_{2}^{-1}\left(\alpha\right)\right)-u_{f}\mu_{2}^{-1}\left(\alpha\right)}{\alpha\left(K_{P}+\mu_{2}^{-1}\left(\alpha\right)\right)-\mu_{2}^{-1}\left(\alpha\right)},$
(39) $\displaystyle\prescript{}{\left(101\right)}{x_{2}}$
$\displaystyle=\omega_{2}\frac{u_{f}\mu_{2}^{-1}\left(\alpha\right)-\alpha\left(1+u_{f}\right)\left(K_{P}+\mu_{2}^{-1}\left(\alpha\right)\right)}{\alpha\left(K_{P}+\mu_{2}^{-1}\left(\alpha\right)\right)-\mu_{2}^{-1}\left(\alpha\right)}+u_{h}-\mu_{2}^{-1}\left(\alpha\right).$
(40)

  * •

${}_{\left(011\right)}\mathcal{E}=\left(0,\prescript{}{\left(011\right)}{x_{1}},\omega_{1}\prescript{}{\left(011\right)}{x_{1}}+u_{h}-\mu_{2}^{-1}\left(\alpha\right)\right)$.
For
$x_{2}=\omega_{1}\prescript{}{\left(011\right)}{x_{1}}+u_{h}-\mu_{2}^{-1}(\alpha)$
the restriction $\prescript{}{\left(011\right)}{x_{2}}>0$ gives the condition
$x_{1}>\frac{\mu_{2}^{-1}\left(\alpha\right)-u_{h}}{\omega_{1}}$, so the
requirement $\prescript{}{\left(011\right)}{x_{1}}\leq u_{g}$ results in the
first condition on $\alpha$, i.e.,
$\alpha<\mu_{2}\left(\omega_{1}u_{g}+u_{h}\right)$ (notice that this already
implies that $\omega_{1}\prescript{}{\left(011\right)}{x_{1}}+u_{h}\geq 0$ is
satisfied). We therefore only consider the mapping
$x_{1}\mapsto\mu_{1}\left(-x_{1}+u_{g},\mu_{2}^{-1}(\alpha)\right)$ for
$\left(\max\left(0,\frac{\mu_{2}^{-1}\left(\alpha\right)-u_{h}}{\omega_{1}}\right),u_{g}\right]$.
For such $x_{1}$, this differentiable mapping is decreasing, so
${}_{\left(011\right)}\mathcal{E}$ exists in $\Omega$ if and only if
$\alpha\in\left[0,\mu_{1}\left(u_{g}-\max{\left(0,\frac{\mu_{2}^{-1}\left(\alpha\right)-u_{h}}{\omega_{1}}\right)},\mu_{2}^{-1}\left(\alpha\right)\right)\right)$
and $\alpha<\mu_{2}\left(\omega_{1}u_{g}+u_{h}\right)$, and when it exists, it
is unique. By solving the equation

$\mu_{1}\left(-x_{1}+u_{g},\mu_{2}^{-1}\left(\alpha\right)\right)=\alpha,$
(41)

for the prototypes given by (8), we obtain the following explicit formulas for
$\prescript{}{\left(011\right)}{x_{1}}$ and
$\prescript{}{\left(011\right)}{x_{2}}$

$\displaystyle\prescript{}{\left(011\right)}{x_{1}}$
$\displaystyle=\frac{\alpha\left(1+u_{g}\right)\left(1+K_{I}\mu_{2}^{-1}\left(\alpha\right)\right)-u_{g}\phi_{1}}{\alpha\left(1+K_{I}\mu_{2}^{-1}\left(\alpha\right)\right)-\phi_{1}},$
(42) $\displaystyle\prescript{}{\left(011\right)}{x_{2}}$
$\displaystyle=\omega_{1}\frac{\alpha\left(1+u_{g}\right)\left(1+K_{I}\mu_{2}^{-1}\left(\alpha\right)\right)-u_{g}\phi_{1}}{\alpha\left(1+K_{I}\mu_{2}^{-1}\left(\alpha\right)\right)-\phi_{1}}+u_{h}-\mu_{2}^{-1}\left(\alpha\right).$
(43)

  * •

${}_{\left(110\right)}\mathcal{E}=\left(\prescript{}{\left(110\right)}{x_{0}},\prescript{}{\left(110\right)}{x_{1}},0\right)$.
This case is much more complicated since we cannot explicitly compute $x_{0}$
as a function of $x_{1}$, or $x_{1}$ as a function of $x_{0}$ in the same way
as in the previous cases. In this case more than one equilibrium of the form
${}_{\left(110\right)}\mathcal{E}$ can exist. In the case of the growth
functions defined in (8), it was proved in [14] that if $u_{g}=u_{h}=0$, then
there exist at most two equilibria of this form. By using the specific growth
functions (8), the equilibria, given as positive solutions of the following
system of equations

$\begin{split}\frac{-x_{0}+u_{f}}{1-x_{0}+u_{f}}\frac{-\omega_{2}x_{0}+\omega_{1}x_{1}+u_{h}}{K_{P}-\omega_{2}x_{0}+\omega_{1}x_{1}+u_{h}}&=\alpha,\\\
\phi_{1}\frac{\omega_{0}x_{0}-x_{1}+u_{g}}{1+\omega_{0}x_{0}-x_{1}+u_{g}}\frac{1}{1+K_{I}\omega_{0}x_{0}-x_{1}+u_{g}}&=\alpha,\end{split}$
(44)

must also satisfy $x_{0}<u_{f}$ and
$\max\left(0,\frac{\omega_{2}x_{0}-u_{h}}{\omega_{1}}\right)<x_{1}<u_{g}+\omega_{0}x_{0}$.
Notice that the first equation in (44) is linear in $x_{1}$, hence we can
compute it as a function of $x_{0}$, and substitute this expression into the
second equation of (44), obtaining a fourth order polynomial in $x_{0}$. Each
zero of this polynomial, together with the corresponding value of $x_{1}$,
which satisfies the aforementioned conditions, will constitute an equilibrium
$\prescript{}{\left(110\right)}{\mathcal{E}}$ of system (21). Since the
polynomial in $x_{0}$ is of order four, and $x_{1}$ is given as a function of
$x_{0}$, we can have at most four equilibria of the form
$\prescript{}{\left(110\right)}{\mathcal{E}}$.

  * •

$\prescript{}{\left(111\right)}{\mathcal{E}}=\left(\overset{*}{x}_{0},\overset{*}{x}_{1},\overset{*}{x}_{2}\right)$.
For the interior equilibrium, we have to consider two cases, depending on the
sign of $\omega_{2}u_{f}-u_{h}$, since the bounds on the values of $x_{1}$ are
different in each case. We are looking for solutions of system (35), for which
$x_{0}$, $x_{1}$, and $x_{2}$ satisfy

$\displaystyle\begin{split}x_{0}&\in\left(0,u_{f}\right],\\\
x_{1}&\in\left(\max\left(0,\frac{\omega_{2}x_{0}-u_{h}}{\omega_{1}}\right),\omega_{0}x_{0}+u_{g}\right],\\\
x_{2}&\in\left(0,-\omega_{2}x_{0}+\omega_{1}x_{1}+u_{h}\right],\end{split}$
(45)

if $\omega_{2}u_{f}-u_{h}>0$, and

$\displaystyle\begin{split}x_{0}&\in\left(0,u_{f}\right],\\\
x_{1}&\in\left(0,\omega_{0}x_{0}+u_{g}\right],\\\
x_{2}&\in\left(0,-\omega_{2}x_{0}+\omega_{1}x_{1}+u_{h}\right],\end{split}$
(46)

if $\omega_{2}u_{f}-u_{h}<0$. In both cases, if we let
$x_{2}=-\omega_{2}x_{0}+\omega_{1}x_{1}+u_{h}-\mu_{2}^{-1}\left(\alpha\right)$
(which immediately gives us a necessary condition
$\alpha<\sup\limits_{s_{2}\geq 0}\mu_{2}\left(\alpha\right)$), we obtain the
following system for $x_{0}$ and $x_{1}$

$\displaystyle\begin{split}\mu_{0}\left(-x_{0}+u_{f},\mu_{2}^{-1}\left(\alpha\right)\right)&=\alpha,\\\
\mu_{1}\left(\omega_{0}x_{0}-x_{1}+u_{g},\mu_{2}^{-1}\left(\alpha\right)\right)&=\alpha.\end{split}$
(47)

For $x_{0}\in\left(0,u_{f}\right]$ the differentiable mapping
$x_{0}\mapsto\mu_{0}\left(-x_{0}+u_{f},\mu_{2}^{-1}\left(\alpha\right)\right)$
is decreasing, so $x_{0}=\overset{*}{x}_{0}$ exists if and only if
$\alpha\in\left[0,\mu_{0}\left(u_{f},\mu_{2}^{-1}\left(\alpha\right)\right)\right)$
and when this value exists, it is unique. Now consider

$\mu_{1}\left(\omega_{0}\overset{*}{x}_{0}-x_{1}+u_{g},\mu_{2}^{-1}\left(\alpha\right)\right)=\alpha.$
(48)

We have two cases

    * –

$\omega_{2}u_{f}-u_{h}>0$. For
$x_{1}\in\left(\max\left(0,\frac{\omega_{2}\overset{*}{x}_{0}-u_{h}}{\omega_{1}}\right),\omega_{0}\overset{*}{x}_{0}+u_{g}\right]$
the differentiable mapping
$x_{1}\mapsto\mu_{1}\left(\omega_{0}\overset{*}{x}_{0}-x_{1}+u_{g},\mu_{2}^{-1}\left(\alpha\right)\right)$
is decreasing, so $x_{1}=\overset{*}{x}_{1}$ exists if and only if
$\alpha\in\left[0,\mu_{1}\left(\omega_{0}\overset{*}{x}_{0}-\max\left(0,\frac{\omega_{2}\overset{*}{x}_{0}-u_{h}}{\omega_{1}}\right)+u_{g},\mu_{2}^{-1}\left(\alpha\right)\right)\right)$
and when this value exists, it is unique.

    * –

$\omega_{2}u_{f}-u_{h}<0$. Similarly, by considering
$x_{1}\in\left(0,\omega_{0}\overset{*}{x}_{0}+u_{g}\right]$ it follows that
$x_{1}=\overset{*}{x}_{1}$ exists if and only if
$\alpha\in\left[0,\mu_{1}\left(\omega_{0}\overset{*}{x}_{0}+u_{g},\mu_{2}^{-1}\left(\alpha\right)\right)\right)$.

Having $\overset{*}{x}_{0}$ and $\overset{*}{x}_{1}$ defined, we let
$\overset{*}{x}_{2}=-\omega_{2}\overset{*}{x}_{0}+\omega_{1}\overset{*}{x}_{1}+u_{h}-\mu_{2}^{-1}\left(\alpha\right)$.
In order to have $\overset{*}{x}_{2}>0$ we need
$\alpha<\mu_{2}\left(-\omega_{2}\overset{*}{x}_{0}+\omega_{1}\overset{*}{x}_{1}+u_{h}\right)$;
hence $\prescript{}{\left(111\right)}{\mathcal{E}}$ exists in $\Omega$ if and
only if

$\alpha\in\Bigg{[}0,\min\Bigg{(}\mu_{0}\left(u_{f},\mu_{2}^{-1}\left(\alpha\right)\right),\\\
\mu_{1}\left(\omega_{0}\overset{*}{x}_{0}-\max\left(0,\frac{\omega_{2}\overset{*}{x}_{0}-u_{h}}{\omega_{1}}\right)+u_{g},\mu_{2}^{-1}\left(\alpha\right)\right),\\\
\mu_{2}\left(-\omega_{2}\overset{*}{x}_{0}+\omega_{1}\overset{*}{x}_{1}+u_{h}\right)\Bigg{)}\Bigg{)},$
(49)

in the $\omega_{2}u_{f}-u_{h}>0$ case, and

$\alpha\in\Big{[}0,\min\Big{(}\mu_{0}\left(u_{f},\mu_{2}^{-1}\left(\alpha\right)\right),\mu_{1}\left(\omega_{0}\overset{*}{x}_{0}+u_{g},\mu_{2}^{-1}\left(\alpha\right)\right),\mu_{2}\left(-\omega_{2}\overset{*}{x}_{0}+\omega_{1}\overset{*}{x}_{1}+u_{h}\right)\Big{)}\Big{)},$
(50)

in the $\omega_{2}u_{f}-u_{h}<0$ case. Although the conditions on $\alpha$ are
implicit and very complicated, we now know that if
$\prescript{}{\left(111\right)}{\mathcal{E}}$ exists, it is unique. With the
growth functions defined in (8), we can solve the equations (35) explicitly
and obtain the following formulas for the interior equilibrium

$\displaystyle\overset{*}{x}_{0}$
$\displaystyle=1+u_{f}+\frac{1}{K_{P}\left(\phi_{2}-\alpha\right)+\alpha-1},$
(51) $\displaystyle\overset{*}{x}_{1}$
$\displaystyle=\omega_{0}\overset{*}{x}_{0}+u_{g}+1+\frac{\phi_{1}}{\alpha\left(1+K_{I}\frac{\alpha}{\phi_{2}-\alpha}\right)-\phi_{1}},$
(52) $\displaystyle\overset{*}{x}_{2}$
$\displaystyle=-\omega_{2}\overset{*}{x}_{0}+\omega_{1}\overset{*}{x}_{1}+u_{h}-\frac{\alpha}{\phi_{2}-\alpha}.$
(53)

####  3.1.2 Local stability results

We now study the local stability of the equilibria by considering the
eigenvalues of the Jacobian evaluated at each equilibrium.

The Jacobian $J$ for system (21) evaluated at $\left(x_{0},x_{1},x_{2}\right)$
has the following form

$J=\begin{bmatrix}\mu_{0}-\alpha+x_{0}\left(-\frac{\partial\mu_{0}}{\partial
s_{0}}-\omega_{2}\frac{\partial\mu_{0}}{\partial
s_{2}}\right)&\omega_{1}x_{0}\frac{\partial\mu_{0}}{\partial
s_{2}}&-x_{0}\frac{\partial\mu_{0}}{\partial s_{2}}\\\
x_{1}\left(\omega_{0}\frac{\partial\mu_{1}}{\partial
s_{1}}-\omega_{2}\frac{\partial\mu_{1}}{\partial
s_{2}}\right)&\mu_{1}-\alpha+x_{1}\left(-\frac{\partial\mu_{1}}{\partial
s_{1}}+\omega_{1}\frac{\partial\mu_{1}}{\partial
s_{2}}\right)&-x_{1}\frac{\partial\mu_{1}}{\partial s_{2}}\\\
-\omega_{2}x_{2}\mu_{2}^{\prime}&\omega_{1}x_{2}\mu_{2}^{\prime}&\mu_{2}-\alpha-
x_{2}\mu_{2}^{\prime}\end{bmatrix}.$ (54)

  * •

For the zero equilibrium $\prescript{}{\left(000\right)}{\mathcal{E}}$, the
corresponding Jacobian $\prescript{}{\left(000\right)}{J}$ has the following
form

$\prescript{}{\left(000\right)}{J}=\begin{bmatrix}\mu_{0}-\alpha&0&0\\\
0&\mu_{1}-\alpha&0\\\ 0&0&\mu_{2}-\alpha\end{bmatrix},$ (55)

and its eigenvalues are $\lambda_{1}=\mu_{0}\left(u_{f},u_{h}\right)-\alpha$,
$\lambda_{2}=\mu_{1}\left(u_{g},u_{h}\right)-\alpha$,
$\lambda_{3}=\mu_{2}\left(u_{h}\right)-\alpha$. This implies that if

    * –

$\alpha>\max\left(\mu_{0}\left(u_{f},u_{h}\right),\mu_{1}\left(u_{g},u_{h}\right),\mu_{2}\left(u_{h}\right)\right)$,
then $\prescript{}{\left(000\right)}{\mathcal{E}}$ is a stable node,

    * –

$\min\left(\mu_{0}\left(u_{f},u_{h}\right),\mu_{1}\left(u_{g},u_{h}\right),\mu_{2}\left(u_{h}\right)\right)<\alpha<\max\left(\mu_{0}\left(u_{f},u_{h}\right),\mu_{1}\left(u_{g},u_{h}\right),\mu_{2}\left(u_{h}\right)\right)$,
then $\prescript{}{\left(000\right)}{\mathcal{E}}$ is a saddle point,

    * –

$\alpha<\min\left(\mu_{0}\left(u_{f},u_{h}\right),\mu_{1}\left(u_{g},u_{h}\right),\mu_{2}\left(u_{h}\right)\right)$,
then $\prescript{}{\left(000\right)}{\mathcal{E}}$ is an unstable node.

  * •

For the boundary equilibrium $\prescript{}{\left(100\right)}{\mathcal{E}}$,
the corresponding Jacobian $\prescript{}{\left(100\right)}{J}$ has the
following form

$\prescript{}{\left(100\right)}{J}=\begin{bmatrix}\prescript{}{\left(100\right)}{x_{0}}\left(-\frac{\partial\mu_{0}}{\partial
s_{0}}-\omega_{2}\frac{\partial\mu_{0}}{\partial
s_{2}}\right)&\omega_{1}\prescript{}{\left(100\right)}{x_{0}}\frac{\partial\mu_{0}}{\partial
s_{2}}&-\prescript{}{\left(100\right)}{x_{0}}\frac{\partial\mu_{0}}{\partial
s_{2}}\\\ 0&\mu_{1}-\alpha&0\\\ 0&0&\mu_{2}-\alpha\end{bmatrix},$ (56)

and its eigenvalues are
$\lambda_{1}=\prescript{}{\left(100\right)}{x_{0}}\left(-\frac{\partial\mu_{0}}{\partial
s_{0}}-\omega_{2}\frac{\partial\mu_{0}}{\partial s_{2}}\right)<0$,
$\lambda_{2}=\mu_{1}\left(\omega_{0}\prescript{}{\left(100\right)}{x_{0}}+u_{g},-\omega_{2}\prescript{}{\left(100\right)}{x_{0}}+u_{h}\right)-\alpha$,
and
$\lambda_{3}=\mu_{2}\left(-\omega_{2}\prescript{}{\left(100\right)}{x_{0}}+u_{h}\right)-\alpha$.
Hence if

    * –

$\alpha>\max\left(\mu_{1}\left(\omega_{0}\prescript{}{\left(100\right)}{x_{0}}+u_{g},-\omega_{2}\prescript{}{\left(100\right)}{x_{0}}+u_{h}\right),\mu_{2}\left(-\omega_{2}\prescript{}{\left(100\right)}{x_{0}}+u_{h}\right)\right)$,
then $\prescript{}{\left(100\right)}{\mathcal{E}}$ is a stable node,

    * –

$\alpha<\max\left(\mu_{1}\left(\omega_{0}\prescript{}{\left(100\right)}{x_{0}}+u_{g},-\omega_{2}\prescript{}{\left(100\right)}{x_{0}}+u_{h}\right),\mu_{2}\left(-\omega_{2}\prescript{}{\left(100\right)}{x_{0}}+u_{h}\right)\right)$,
then $\prescript{}{\left(100\right)}{\mathcal{E}}$ is a saddle point.

  * •

For the boundary equilibrium $\prescript{}{\left(010\right)}{\mathcal{E}}$,
the corresponding Jacobian $\prescript{}{\left(010\right)}{J}$ has the
following form

$\prescript{}{\left(010\right)}{J}=\begin{bmatrix}\mu_{0}-\alpha&0&0\\\
\prescript{}{\left(010\right)}{x_{1}}\left(\omega_{0}\frac{\partial\mu_{1}}{\partial
s_{1}}-\omega_{2}\frac{\partial\mu_{1}}{\partial
s_{2}}\right)&\prescript{}{\left(010\right)}{x_{1}}\left(-\frac{\partial\mu_{1}}{\partial
s_{1}}+\omega_{1}\frac{\partial\mu_{1}}{\partial
s_{2}}\right)&-\prescript{}{\left(010\right)}{x_{1}}\frac{\partial\mu_{1}}{\partial
s_{2}}\\\ 0&0&\mu_{2}-\alpha\end{bmatrix},$ (57)

and its eigenvalues are
$\lambda_{1}=\mu_{0}\left(u_{f},\omega_{1}\prescript{}{\left(010\right)}{x_{1}}+u_{h}\right)-\alpha$,
$\lambda_{2}=\prescript{}{\left(010\right)}{x_{1}}\left(-\frac{\partial\mu_{1}}{\partial
s_{1}}+\omega_{1}\frac{\partial\mu_{1}}{\partial s_{2}}\right)<0$, and
$\lambda_{3}=\mu_{2}\left(\omega_{1}\prescript{}{\left(010\right)}{x_{1}}+u_{h}\right)-\alpha$.
Hence if

    * –

$\alpha>\max\left(\mu_{0}\left(u_{f},\omega_{1}\prescript{}{\left(010\right)}{x_{1}}+u_{h}\right),\mu_{2}\left(\omega_{1}\prescript{}{\left(010\right)}{x_{1}}+u_{h}\right)\right)$,
then $\prescript{}{\left(010\right)}{\mathcal{E}}$ is a stable node,

    * –

$\alpha<\max\left(\mu_{0}\left(u_{f},\omega_{1}\prescript{}{\left(010\right)}{x_{1}}+u_{h}\right),\mu_{2}\left(\omega_{1}\prescript{}{\left(010\right)}{x_{1}}+u_{h}\right)\right)$,
then $\prescript{}{\left(010\right)}{\mathcal{E}}$ is a saddle point.

  * •

For the boundary equilibrium $\prescript{}{\left(001\right)}{\mathcal{E}}$,
the corresponding Jacobian $\prescript{}{\left(001\right)}{J}$ has the
following form

$\prescript{}{\left(001\right)}{J}=\begin{bmatrix}\mu_{0}-\alpha&0&0\\\
0&\mu_{1}-\alpha&0\\\
-\omega_{2}\prescript{}{\left(001\right)}{x_{2}}\mu_{2}^{\prime}&\omega_{1}\prescript{}{\left(001\right)}{x_{2}}\mu_{2}^{\prime}&\mu_{2}-\alpha-\prescript{}{\left(001\right)}{x_{2}}\mu_{2}^{\prime}\end{bmatrix},$
(58)

and its eigenvalues are
$\lambda_{1}=\mu_{0}\left(u_{f},-\prescript{}{\left(001\right)}{x_{2}}+u_{h}\right)-\alpha$,
$\lambda_{2}=\mu_{1}\left(u_{g},-\prescript{}{\left(001\right)}{x_{2}}+u_{h}\right)-\alpha$,
and
$\lambda_{3}=-\prescript{}{\left(001\right)}{x_{2}}\mu_{2}^{\prime}\left(-\prescript{}{\left(001\right)}{x_{2}}+u_{h}\right)<0$.
Hence if

    * –

$\alpha>\max\left(\mu_{0}\left(u_{f},-\prescript{}{\left(001\right)}{x_{2}}+u_{h}\right),\mu_{1}\left(u_{g},-\prescript{}{\left(001\right)}{x_{2}}+u_{h}\right)\right)$,
then $\prescript{}{\left(001\right)}{\mathcal{E}}$ is a stable node,

    * –

$\alpha<\max\left(\mu_{0}\left(u_{f},-\prescript{}{\left(001\right)}{x_{2}}+u_{h}\right),\mu_{1}\left(u_{g},-\prescript{}{\left(001\right)}{x_{2}}+u_{h}\right)\right)$,
then $\prescript{}{\left(001\right)}{\mathcal{E}}$ is a saddle point.

  * •

For the boundary equilibrium $\prescript{}{\left(101\right)}{\mathcal{E}}$,
the corresponding Jacobian $\prescript{}{\left(101\right)}{J}$ has the
following form

$\prescript{}{\left(101\right)}{J}=\begin{bmatrix}\prescript{}{\left(101\right)}{x_{0}}\left(-\frac{\partial\mu_{0}}{\partial
s_{0}}-\omega_{2}\frac{\partial\mu_{0}}{\partial
s_{2}}\right)&\omega_{1}\prescript{}{\left(101\right)}{x_{0}}\frac{\partial\mu_{0}}{\partial
s_{2}}&-\prescript{}{\left(101\right)}{x_{0}}\frac{\partial\mu_{0}}{\partial
s_{2}}\\\ 0&\mu_{1}-\alpha&0\\\
-\omega_{2}\prescript{}{\left(101\right)}{x_{2}}\mu_{2}^{\prime}&\omega_{1}\prescript{}{\left(101\right)}{x_{2}}\mu_{2}^{\prime}&-\prescript{}{\left(101\right)}{x_{2}}\mu_{2}^{\prime}\end{bmatrix},$
(59)

where
$\prescript{}{\left(101\right)}{x_{2}}=-\omega_{2}\prescript{}{\left(101\right)}{x_{0}}+u_{h}-\mu_{2}^{-1}\left(\alpha\right)$.
We immediately obtain one of the eigenvalues
$\lambda_{1}=\mu_{1}\left(\omega_{0}\prescript{}{\left(101\right)}{x_{0}}+u_{g},\mu_{2}^{-1}\left(\alpha\right)\right)-\alpha$.
The other two eigenvalues are given as the solutions of the following
quadratic equation

$\lambda^{2}+a_{1}\lambda+a_{0}=0,$ (60)

where

$\displaystyle a_{1}=$
$\displaystyle\prescript{}{\left(101\right)}{x_{0}}\left(\frac{\partial\mu_{0}}{\partial
s_{0}}+\omega_{2}\frac{\partial\mu_{0}}{\partial
s_{2}}\right)+\prescript{}{\left(101\right)}{x_{2}}\mu_{2}^{\prime}$ (61)
$\displaystyle a_{0}=$
$\displaystyle\prescript{}{\left(101\right)}{x_{0}}\prescript{}{\left(101\right)}{x_{2}}\frac{\partial\mu_{0}}{\partial
s_{0}}\mu_{2}^{\prime}.$ (62)

Since both $a_{1}>0$, and $a_{0}>0$, by the Routh-Hurwitz criterion, all roots
of the equation (60) have negative real parts. Hence if

    * –

$\alpha>\mu_{1}\left(\omega_{0}\prescript{}{\left(101\right)}{x_{0}}+u_{g},\mu_{2}^{-1}\left(\alpha\right)\right)$,
then $\prescript{}{\left(101\right)}{\mathcal{E}}$ is asymptotically stable,

    * –

$\alpha<\mu_{1}\left(\omega_{0}\prescript{}{\left(101\right)}{x_{0}}+u_{g},\mu_{2}^{-1}\left(\alpha\right)\right)$,
then $\prescript{}{\left(101\right)}{\mathcal{E}}$ is a saddle point.

  * •

For the boundary equilibrium $\prescript{}{\left(011\right)}{\mathcal{E}}$,
the corresponding Jacobian $\prescript{}{\left(011\right)}{J}$ has the
following form

$\prescript{}{\left(011\right)}{J}=\begin{bmatrix}\mu_{0}-\alpha&0&0\\\
\prescript{}{\left(011\right)}{x_{1}}\left(\omega_{0}\frac{\partial\mu_{1}}{\partial
s_{1}}-\omega_{2}\frac{\partial\mu_{1}}{\partial
s_{2}}\right)&\prescript{}{\left(011\right)}{x_{1}}\left(-\frac{\partial\mu_{1}}{\partial
s_{1}}+\omega_{1}\frac{\partial\mu_{1}}{\partial
s_{2}}\right)&-\prescript{}{\left(011\right)}{x_{1}}\frac{\partial\mu_{1}}{\partial
s_{2}}\\\
-\omega_{2}\prescript{}{\left(011\right)}{x_{2}}\mu_{2}^{\prime}&\omega_{1}\prescript{}{\left(011\right)}{x_{2}}\mu_{2}^{\prime}&-\prescript{}{\left(011\right)}{x_{2}}\mu_{2}^{\prime}\end{bmatrix},$
(63)

where
$\prescript{}{\left(011\right)}{x_{2}}=\omega_{1}\prescript{}{\left(011\right)}{x_{1}}+u_{h}-\mu_{2}^{-1}(\alpha)$.
We immediately obtain one of the eigenvalues
$\lambda_{1}=\mu_{0}\left(u_{f},\mu_{2}^{-1}\left(\alpha\right)\right)-\alpha$.
The other two eigenvalues are given as the solutions of the following
quadratic equation

$\lambda^{2}+a_{1}\lambda+a_{0}=0,$ (64)

where

$\displaystyle a_{1}=$
$\displaystyle\prescript{}{\left(011\right)}{x_{1}}\left(\frac{\partial\mu_{1}}{\partial
s_{1}}-\omega_{1}\frac{\partial\mu_{1}}{\partial
s_{2}}\right)+\prescript{}{\left(011\right)}{x_{2}}\mu_{2}^{\prime},$ (65)
$\displaystyle a_{0}=$
$\displaystyle\prescript{}{\left(011\right)}{x_{1}}\prescript{}{\left(011\right)}{x_{2}}\frac{\partial\mu_{1}}{\partial
s_{1}}\mu_{2}^{\prime}.$ (66)

Since both $a_{1}>0$, and $a_{0}>0$, by the Routh-Hurwitz criterion, all roots
of the equation (64) have negative real parts. Hence if

    * –

$\alpha>\mu_{0}\left(u_{f},\mu_{2}^{-1}\left(\alpha\right)\right)$, then
$\prescript{}{\left(011\right)}{\mathcal{E}}$ is asymptotically stable,

    * –

$\alpha<\mu_{0}\left(u_{f},\mu_{2}^{-1}\left(\alpha\right)\right)$, then
$\prescript{}{\left(011\right)}{\mathcal{E}}$ is a saddle point.

  * •

For the boundary equilibrium $\prescript{}{\left(110\right)}{\mathcal{E}}$,
the corresponding Jacobian $\prescript{}{\left(110\right)}{J}$ has the
following form

$\prescript{}{\left(110\right)}{J}=\begin{bmatrix}\prescript{}{\left(110\right)}{x_{0}}\left(-\frac{\partial\mu_{0}}{\partial
s_{0}}-\omega_{2}\frac{\partial\mu_{0}}{\partial
s_{2}}\right)&\omega_{1}\prescript{}{\left(110\right)}{x_{0}}\frac{\partial\mu_{0}}{\partial
s_{2}}&-\prescript{}{\left(110\right)}{x_{0}}\frac{\partial\mu_{0}}{\partial
s_{2}}\\\
\prescript{}{\left(110\right)}{x_{1}}\left(\omega_{0}\frac{\partial\mu_{1}}{\partial
s_{1}}-\omega_{2}\frac{\partial\mu_{1}}{\partial
s_{2}}\right)&\prescript{}{\left(110\right)}{x_{1}}\left(-\frac{\partial\mu_{1}}{\partial
s_{1}}+\omega_{1}\frac{\partial\mu_{1}}{\partial
s_{2}}\right)&-\prescript{}{\left(110\right)}{x_{1}}\frac{\partial\mu_{1}}{\partial
s_{2}}\\\ 0&0&\mu_{2}-\alpha\end{bmatrix}.$ (67)

We immediately obtain one eigenvalue $\lambda_{1}=\mu_{2}-\alpha$. The other
two eigenvalues are solutions of the following quadratic equation

$\lambda^{2}+a_{1}\lambda+a_{0}=0,$ (68)

where

$\displaystyle a_{1}$
$\displaystyle=\prescript{}{\left(110\right)}{x_{0}}\left(\frac{\partial\mu_{0}}{\partial
s_{0}}+\omega_{2}\frac{\partial\mu_{0}}{\partial
s_{2}}\right)+\prescript{}{\left(110\right)}{x_{1}}\left(\frac{\partial\mu_{1}}{\partial
s_{1}}-\omega_{1}\frac{\partial\mu_{1}}{\partial s_{2}}\right),$ (69)
$\displaystyle a_{0}$
$\displaystyle=\prescript{}{\left(110\right)}{x_{0}}\prescript{}{\left(110\right)}{x_{1}}\left(\frac{\partial\mu_{0}}{\partial
s_{0}}+\omega_{2}\frac{\partial\mu_{0}}{\partial
s_{2}}\right)\left(\frac{\partial\mu_{1}}{\partial
s_{1}}-\omega_{1}\frac{\partial\mu_{1}}{\partial s_{2}}\right)$ (70)
$\displaystyle\phantom{{}=}-\omega_{1}\prescript{}{\left(110\right)}{x_{0}}\prescript{}{\left(110\right)}{x_{1}}\frac{\partial\mu_{0}}{\partial
s_{2}}\left(\omega_{0}\frac{\partial\mu_{1}}{\partial
s_{1}}-\omega_{2}\frac{\partial\mu_{1}}{\partial s_{2}}\right).$

We have $a_{1}>0$, and

$a_{0}>0\iff\frac{\partial\mu_{0}}{\partial
s_{0}}\frac{\partial\mu_{1}}{\partial
s_{1}}-\omega_{1}\frac{\partial\mu_{0}}{\partial
s_{0}}\frac{\partial\mu_{1}}{\partial
s_{2}}+\left(\omega_{2}-\omega_{0}\omega_{1}\right)\frac{\partial\mu_{0}}{\partial
s_{2}}\frac{\partial\mu_{1}}{\partial s_{1}}>0.$ (71)

Hence if

    * –

$\alpha>\mu_{2}$ and $\frac{\partial\mu_{0}}{\partial
s_{0}}\frac{\partial\mu_{1}}{\partial
s_{1}}-\omega_{1}\frac{\partial\mu_{0}}{\partial
s_{0}}\frac{\partial\mu_{1}}{\partial
s_{2}}+\left(\omega_{2}-\omega_{0}\omega_{1}\right)\frac{\partial\mu_{0}}{\partial
s_{2}}\frac{\partial\mu_{1}}{\partial s_{1}}>0$ (where all the functions are
evaluated at the steady state), then
$\prescript{}{\left(110\right)}{\mathcal{E}}$ is asymptotically stable,

    * –

$\alpha<\mu_{2}$ or $\frac{\partial\mu_{0}}{\partial
s_{0}}\frac{\partial\mu_{1}}{\partial
s_{1}}-\omega_{1}\frac{\partial\mu_{0}}{\partial
s_{0}}\frac{\partial\mu_{1}}{\partial
s_{2}}+\left(\omega_{2}-\omega_{0}\omega_{1}\right)\frac{\partial\mu_{0}}{\partial
s_{2}}\frac{\partial\mu_{1}}{\partial s_{1}}<0$, then
$\prescript{}{\left(110\right)}{\mathcal{E}}$ is unstable.

  * •

For the interior equilibrium $\prescript{}{\left(111\right)}{\mathcal{E}}$,
the corresponding Jacobian $\prescript{}{\left(111\right)}{J}$ has the
following form

$\prescript{}{\left(111\right)}{J}=\begin{bmatrix}\overset{*}{x}_{0}\left(-\frac{\partial\mu_{0}}{\partial
s_{0}}-\omega_{2}\frac{\partial\mu_{0}}{\partial
s_{2}}\right)&\omega_{1}\overset{*}{x}_{0}\frac{\partial\mu_{0}}{\partial
s_{2}}&-\overset{*}{x}_{0}\frac{\partial\mu_{0}}{\partial s_{2}}\\\
\overset{*}{x}_{1}\left(\omega_{0}\frac{\partial\mu_{1}}{\partial
s_{1}}-\omega_{2}\frac{\partial\mu_{1}}{\partial
s_{2}}\right)&\overset{*}{x}_{1}\left(-\frac{\partial\mu_{1}}{\partial
s_{1}}+\omega_{1}\frac{\partial\mu_{1}}{\partial
s_{2}}\right)&-\overset{*}{x}_{1}\frac{\partial\mu_{1}}{\partial s_{2}}\\\
-\omega_{2}\overset{*}{x}_{2}\mu_{2}^{\prime}&\omega_{1}\overset{*}{x}_{2}\mu_{2}^{\prime}&-\overset{*}{x}_{2}\mu_{2}^{\prime}\end{bmatrix},$
(72)

and its eigenvalues are solutions of the following cubic equation

$\lambda^{3}+a_{2}\lambda^{2}+a_{1}\lambda+a_{0}=0,$ (73)

where

$\displaystyle a_{2}$
$\displaystyle=-\overset{*}{x}_{0}\left(-\frac{\partial\mu_{0}}{\partial
s_{0}}-\omega_{2}\frac{\partial\mu_{0}}{\partial
s_{2}}\right)-\overset{*}{x}_{1}\left(-\frac{\partial\mu_{1}}{\partial
s_{1}}+\omega_{1}\frac{\partial\mu_{1}}{\partial
s_{2}}\right)+\overset{*}{x}_{2}\mu_{2}^{\prime},$ (74) $\displaystyle a_{1}$
$\displaystyle=\overset{*}{x}_{1}\frac{\partial\mu_{1}}{\partial
s_{1}}\left(\overset{*}{x}_{0}\frac{\partial\mu_{0}}{\partial
s_{0}}-\left(\omega_{0}\omega_{1}-\omega_{2}\right)\overset{*}{x}_{0}\frac{\partial\mu_{0}}{\partial
s_{2}}+\overset{*}{x}_{2}\mu_{2}^{\prime}\right)+\overset{*}{x}_{0}\frac{\partial\mu_{0}}{\partial
s_{0}}\left(-\omega_{1}\overset{*}{x}_{1}\frac{\partial\mu_{1}}{\partial
s_{2}}+\overset{*}{x}_{2}\mu_{2}^{\prime}\right),$ (75) $\displaystyle a_{0}$
$\displaystyle=\overset{*}{x}_{0}\overset{*}{x}_{1}\overset{*}{x}_{2}\frac{\partial\mu_{0}}{\partial
s_{0}}\frac{\partial\mu_{1}}{\partial s_{1}}\mu_{2}^{\prime}.$ (76)

By the Routh-Hurwitz criterion, all eigenvalues have negative real parts if
and only if $a_{2}>0$, $a_{0}>0$, and $a_{2}a_{1}>a_{0}$. We have $a_{2}>0$,
and $a_{0}>0$ always, and

$\begin{split}a_{2}a_{1}>a_{0}\iff\left[-\overset{*}{x}_{0}\left(-\frac{\partial\mu_{0}}{\partial
s_{0}}-\omega_{2}\frac{\partial\mu_{0}}{\partial
s_{2}}\right)-\overset{*}{x}_{1}\left(-\frac{\partial\mu_{1}}{\partial
s_{1}}+\omega_{1}\frac{\partial\mu_{1}}{\partial
s_{2}}\right)+\overset{*}{x}_{2}\mu_{2}^{\prime}\right]\cdot\\\
\left[\overset{*}{x}_{1}\frac{\partial\mu_{1}}{\partial
s_{1}}\left(\overset{*}{x}_{0}\frac{\partial\mu_{0}}{\partial
s_{0}}-\left(\omega_{0}\omega_{1}-\omega_{2}\right)\overset{*}{x}_{0}\frac{\partial\mu_{0}}{\partial
s_{2}}+\overset{*}{x}_{2}\mu_{2}^{\prime}\right)+\overset{*}{x}_{0}\frac{\partial\mu_{0}}{\partial
s_{0}}\left(-\omega_{1}\overset{*}{x}_{1}\frac{\partial\mu_{1}}{\partial
s_{2}}+\overset{*}{x}_{2}\mu_{2}^{\prime}\right)\right]\\\
-\overset{*}{x}_{0}\overset{*}{x}_{1}\overset{*}{x}_{2}\frac{\partial\mu_{0}}{\partial
s_{0}}\frac{\partial\mu_{1}}{\partial s_{1}}\mu_{2}^{\prime}>0.\end{split}$
(77)

Thus $\prescript{}{\left(111\right)}{\mathcal{E}}$ is asymptotically stable if
$a_{2}a_{1}>a_{0}$, and is unstable if $a_{2}a_{1}<a_{0}$.

We gather the results concerning the existence and local stability of the
equilibria of (21) in table 2. All the functions in the "Local stability"
column are evaluated at the corresponding steady states, and all the symbols
are given by equations (78)-(84). The conditions for existence of the
$\prescript{}{\left(110\right)}{\mathcal{E}}$ equilibrium are discussed in
subsection 3.1.1.

Equilibrium |  Existence |  Local stability
---|---|---
${}_{\left(000\right)}\mathcal{E}$ |  Always |  $\alpha>\max(\mu_{0},\mu_{1},\mu_{2})$
|  |
${}_{\left(100\right)}\mathcal{E}$ |  $\alpha<\mu_{0}\left(u_{f},u_{h}\right)$ |  $\alpha>\max\left(\mu_{1},\mu_{2}\right)$
|  |
${}_{\left(010\right)}\mathcal{E}$ |  $\alpha<\mu_{1}\left(u_{g},u_{h}\right)$ |  $\alpha>\max\left(\mu_{0},\mu_{2}\right)$
|  |
$\prescript{}{\left(001\right)}{\mathcal{E}}$ |  $\alpha<\mu_{2}\left(u_{h}\right)$ |  $\alpha>\max\left(\mu_{0},\mu_{1}\right)$
|  |
$\prescript{}{\left(101\right)}{\mathcal{E}}$ |  $\prescript{}{\left(101\right)}{\star}$ |  $\alpha>\mu_{1}$
|  |
$\prescript{}{\left(011\right)}{\mathcal{E}}$ |  $\prescript{}{\left(011\right)}{\star}$ |  $\alpha>\mu_{0}$
|  |
$\prescript{}{\left(110\right)}{\mathcal{E}}$ |  - |  $\alpha>\mu_{2}$ and  $\prescript{}{\left(110\right)}{a_{0}}>0$
|  |
$\prescript{}{\left(111\right)}{\mathcal{E}}$ |  $\prescript{}{\left(111\right)}{\star}$ |  $\prescript{}{\left(111\right)}{a_{2}}\prescript{}{\left(111\right)}{a_{1}}-\prescript{}{\left(111\right)}{a_{0}}>0$,
|  |
Table 2: Equilibria of system (21) together with their local stability.

$\displaystyle\prescript{}{\left(101\right)}{\star}$
$\displaystyle=\alpha\in\left(\mu_{0}\left(u_{f}-\min{\left(u_{f},\frac{u_{h}-\mu_{2}^{-1}\left(\alpha\right)}{\omega_{2}}\right)},\mu_{2}^{-1}\left(\alpha\right)\right),\mu_{0}\left(u_{f},\mu_{2}^{-1}\left(\alpha\right)\right)\right),$
(78) $\displaystyle\prescript{}{\left(011\right)}{\star}$
$\displaystyle=\alpha\in\left[0,\mu_{1}\left(u_{g}-\max{\left(0,\frac{\mu_{2}^{-1}\left(\alpha\right)-u_{h}}{\omega_{1}}\right)},\mu_{2}^{-1}\left(\alpha\right)\right)\right),$
(79) $\displaystyle\phantom{{}=}\text{ and
}\alpha<\mu_{2}\left(\omega_{1}u_{g}+u_{h}\right),$
$\displaystyle\prescript{}{\left(110\right)}{a_{0}}$
$\displaystyle=\frac{\partial\mu_{0}}{\partial
s_{0}}\frac{\partial\mu_{1}}{\partial
s_{1}}-\omega_{1}\frac{\partial\mu_{0}}{\partial
s_{0}}\frac{\partial\mu_{1}}{\partial
s_{2}}+\left(\omega_{2}-\omega_{0}\omega_{1}\right)\frac{\partial\mu_{0}}{\partial
s_{2}}\frac{\partial\mu_{1}}{\partial s_{1}},$ (80)
$\displaystyle\prescript{}{\left(111\right)}{\star}$
$\displaystyle=\alpha\in\Bigg{[}0,\min\Bigg{(}\mu_{0}\left(u_{f},\mu_{2}^{-1}\left(\alpha\right)\right),\mu_{1}\Bigg{(}\omega_{0}\overset{*}{x}_{0}-\max\left(0,\frac{\omega_{2}\overset{*}{x}_{0}-u_{h}}{\omega_{1}}\right)$
(81)
$\displaystyle\phantom{{}=}+u_{g},\mu_{2}^{-1}\left(\alpha\right)\Bigg{)},\mu_{2}\left(-\omega_{2}\overset{*}{x}_{0}+\omega_{1}\overset{*}{x}_{1}+u_{h}\right)\Bigg{)}\Bigg{)},$
$\displaystyle\prescript{}{\left(111\right)}{a_{2}}$
$\displaystyle=-\overset{*}{x}_{0}\left(-\frac{\partial\mu_{0}}{\partial
s_{0}}-\omega_{2}\frac{\partial\mu_{0}}{\partial
s_{2}}\right)-\overset{*}{x}_{1}\left(-\frac{\partial\mu_{1}}{\partial
s_{1}}+\omega_{1}\frac{\partial\mu_{1}}{\partial
s_{2}}\right)+\overset{*}{x}_{2}\mu_{2}^{\prime},$ (82)
$\displaystyle\prescript{}{\left(111\right)}{a_{1}}$
$\displaystyle=\overset{*}{x}_{1}\frac{\partial\mu_{1}}{\partial
s_{1}}\left(\overset{*}{x}_{0}\frac{\partial\mu_{0}}{\partial
s_{0}}-\left(\omega_{0}\omega_{1}-\omega_{2}\right)\overset{*}{x}_{0}\frac{\partial\mu_{0}}{\partial
s_{2}}+\overset{*}{x}_{2}\mu_{2}^{\prime}\right)$ (83)
$\displaystyle\phantom{{}=}+\overset{*}{x}_{0}\frac{\partial\mu_{0}}{\partial
s_{0}}\left(-\omega_{1}\overset{*}{x}_{1}\frac{\partial\mu_{1}}{\partial
s_{2}}+\overset{*}{x}_{2}\mu_{2}^{\prime}\right),$
$\displaystyle\prescript{}{\left(111\right)}{a_{0}}$
$\displaystyle=\overset{*}{x}_{0}\overset{*}{x}_{1}\overset{*}{x}_{2}\frac{\partial\mu_{0}}{\partial
s_{0}}\frac{\partial\mu_{1}}{\partial s_{1}}\mu_{2}^{\prime}.$ (84)

##  4 Analysis of the full system

###  4.1 Periodic orbits on the faces

We begin by ruling out possibility of having a periodic orbit in one of the
invariant faces of $\Omega$. This can be done by using general forms of
prototypes $\mu_{0}$, $\mu_{1}$ and $\mu_{2}$.

First, consider system (21) on the part of $\Omega$ with $x_{0}=0$, i.e.,

$\displaystyle\begin{cases}\begin{split}x_{1}^{\prime}&=\left(-\alpha+\mu_{1}\left(-x_{1}+u_{g},\omega_{1}x_{1}-x_{2}+u_{h}\right)\right)x_{1},\\\
x_{2}^{\prime}&=\left(-\alpha+\mu_{2}\left(\omega_{1}x_{1}-x_{2}+u_{h}\right)\right)x_{2}.\end{split}\end{cases}$
(85)

The domain for system (85) is given by the following set $\Omega_{12}$:

$\Omega_{12}=\left\\{\left(x_{1},x_{2}\right)\in\mathbb{R}^{2}:0\leq x_{1}\leq
u_{g},0\leq x_{2}\leq u_{h}+\omega_{1}x_{1}\right\\}.$ (86)

Notice that no periodic orbit can intersect the axes $x_{1}=0$ or $x_{2}=0$
since they are invariant. Now, let us define an auxiliary function

$\varphi_{0}\left(x_{1},x_{2}\right)=\frac{1}{x_{1}x_{2}},\quad\left(x_{1},x_{2}\right)\in\Omega_{12}\setminus\left(\left\\{x_{1}=0\right\\}\cup\left\\{x_{2}=0\right\\}\right).$
(87)

Then

$\nabla\cdot\left(\varphi_{0}\left(-\alpha+\mu_{1}\right)x_{1},\varphi_{0}\left(-\alpha+\mu_{2}\right)x_{2}\right)=\frac{-\partial_{s_{1}}\mu_{1}+\omega_{1}\partial_{s_{2}}\mu_{1}}{x_{2}}-\frac{\mu^{\prime}_{2}}{x_{1}}<0$
(88)

for all $\left(x_{1},x_{2}\right)$ in the domain of $\varphi_{0}$. Thus, by
the Dulac’s Criterion [9], there are no periodic orbits in the $x_{1}x_{2}$
face.

Now, consider system (21) on the part of $\Omega$ with $x_{1}=0$, i.e.,

$\displaystyle\begin{cases}\begin{split}x_{0}^{\prime}&=\left(-\alpha+\mu_{0}\left(-x_{0}+u_{f},-\omega_{2}x_{0}-x_{2}+u_{h}\right)\right)x_{0},\\\
x_{2}^{\prime}&=\left(-\alpha+\mu_{2}\left(-\omega_{2}x_{0}-x_{2}+u_{h}\right)\right)x_{2},\end{split}\end{cases}$
(89)

defined on $\Omega_{02}$ given by

$\Omega_{02}=\left\\{\left(x_{0},x_{2}\right)\in\mathbb{R}^{2}:0\leq
x_{0}\leq\min\left\\{u_{f},\frac{u_{h}-x_{2}}{\omega_{2}}\right\\},0\leq
x_{2}\leq-\omega_{2}x_{0}+u_{h}\right\\}.$ (90)

For auxiliary function

$\varphi_{1}\left(x_{0},x_{2}\right)=\frac{1}{x_{0}x_{2}},\quad\left(x_{0},x_{2}\right)\in\Omega_{02}\setminus\left(\left\\{x_{0}=0\right\\}\cup\left\\{x_{2}=0\right\\}\right)$
(91)

we have

$\nabla\cdot\left(\varphi_{1}\left(-\alpha+\mu_{1}\right)x_{0},\varphi_{1}\left(-\alpha+\mu_{2}\right)x_{2}\right)=\frac{-\partial_{s_{0}}\mu_{0}-\omega_{2}\partial_{s_{2}}\mu_{0}}{x_{2}}-\frac{\mu^{\prime}_{2}}{x_{0}}<0.$
(92)

This, together with the fact that the axes $x_{0}=0$, $x_{2}=0$ are invariant,
shows that there are no periodic orbits in the $x_{0}x_{2}$ face.

Finally, consider system (21) on the part of $\Omega$ with $x_{2}=0$, i.e.,

$\displaystyle\begin{cases}\begin{split}x_{0}^{\prime}&=\left(-\alpha+\mu_{0}\left(-x_{0}+u_{f},-\omega_{2}x_{0}+\omega_{1}x_{1}+u_{h}\right)\right)x_{0},\\\
x_{1}^{\prime}&=\left(-\alpha+\mu_{1}\left(\omega_{0}x_{0}-x_{1}+u_{g},-\omega_{2}x_{0}+\omega_{1}x_{1}+u_{h}\right)\right)x_{1}.\end{split}\end{cases}$
(93)

defined on $\Omega_{01}$ given by

$\Omega_{01}=\left\\{\left(x_{0},x_{1}\right)\in\mathbb{R}^{2}:0\leq x_{0}\leq
u_{f},\max\left\\{0,\frac{\omega_{2}x_{0}-u_{h}}{\omega_{1}}\right\\}\leq
x_{1}\leq\omega_{0}x_{0}+u_{g}\right\\}.$ (94)

Analogously to the previous cases, we define

$\varphi_{2}\left(x_{0},x_{1}\right)=\frac{1}{x_{0}x_{1}},\quad\left(x_{0},x_{1}\right)\in\Omega_{01}\setminus\left(\left\\{x_{0}=0\right\\}\cup\left\\{x_{1}=0\right\\}\right)$
(95)

and compute

$\nabla\cdot\left(\varphi\left(-\alpha+\mu_{0}\right)x_{0},\varphi\left(-\alpha+\mu_{1}\right)x_{1}\right)=\frac{-\partial_{s_{0}}\mu_{0}-\omega_{2}\partial_{s_{2}}\mu_{0}}{x_{1}}+\frac{-\partial_{s_{1}}\mu_{1}+\omega_{1}\partial_{s_{2}}\mu_{1}}{x_{0}}<0.$
(96)

Since the axes $x_{0}=0$, $x_{1}=0$ are invariant, by Dulac’s Criterion, there
are no periodic orbits in the $x_{0}x_{1}$ face.

###  4.2 Hopf Bifurcation

In this work we are especially interested in developing a more systematic
approach to studying the Hopf bifurcation of the interior equilibrium. That
Hopf bifurcation that occurs in this model was previously observed numerically
[14]. The occurrence of a stable periodic orbit in system (21) represents a
situation in which all three populations of microorganisms oscillate
indefinitely, and as a consequence, the substrate concentrations fluctuate.
The characteristic polynomial of the Jacobian
$\prescript{}{\left(111\right)}{J}$ corresponding to the interior equilibrium
is given by

$\lambda^{3}+a_{2}\lambda^{2}+a_{1}\lambda+a_{0}=0,$ (97)

where

$\displaystyle a_{2}$
$\displaystyle=-\overset{*}{x}_{0}\left(-\frac{\partial\mu_{0}}{\partial
s_{0}}-\omega_{2}\frac{\partial\mu_{0}}{\partial
s_{2}}\right)-\overset{*}{x}_{1}\left(-\frac{\partial\mu_{1}}{\partial
s_{1}}+\omega_{1}\frac{\partial\mu_{1}}{\partial
s_{2}}\right)+\overset{*}{x}_{2}\mu_{2}^{\prime},$ (98) $\displaystyle a_{1}$
$\displaystyle=\overset{*}{x}_{1}\frac{\partial\mu_{1}}{\partial
s_{1}}\left(\overset{*}{x}_{0}\frac{\partial\mu_{0}}{\partial
s_{0}}-\left(\omega_{0}\omega_{1}-\omega_{2}\right)\overset{*}{x}_{0}\frac{\partial\mu_{0}}{\partial
s_{2}}+\overset{*}{x}_{2}\mu_{2}^{\prime}\right)+\overset{*}{x}_{0}\frac{\partial\mu_{0}}{\partial
s_{0}}\left(-\omega_{1}\overset{*}{x}_{1}\frac{\partial\mu_{1}}{\partial
s_{2}}+\overset{*}{x}_{2}\mu_{2}^{\prime}\right),$ (99) $\displaystyle a_{0}$
$\displaystyle=\overset{*}{x}_{0}\overset{*}{x}_{1}\overset{*}{x}_{2}\frac{\partial\mu_{0}}{\partial
s_{0}}\frac{\partial\mu_{1}}{\partial s_{1}}\mu_{2}^{\prime},$ (100)

and the coefficients $a_{2}$, $a_{1}$, and $a_{0}$ depend on the parameters
$u_{f}$, $u_{g}$, $u_{h}$, and $\alpha$. The coefficients $a_{2}$ and $a_{0}$
are sign-definite (they are both positive), and $a_{1}$ might possibly change
sign. Let us first notice that since the polynomial has order three, a real
eigenvalue always exists. By the Routh-Hurwitz criterion, the above polynomial
has a pair of purely imaginary eigenvalues if and only if

$a_{2}a_{1}=a_{0}\quad\text{(which implies $a_{1}>0$).}$ (101)

In that case we also have

$\lambda^{3}+a_{2}\lambda^{2}+a_{1}\lambda+a_{0}=\lambda^{3}+a_{2}\lambda^{2}+a_{1}\lambda+a_{2}a_{1}=\left(\lambda+a_{2}\right)\left(\lambda^{2}+a_{1}\right).$
(102)

Hence the eigenvalues are $\lambda_{1}=-a_{2}$ and
$\lambda_{2,3}=\pm\sqrt{a_{1}}i$. Since eigenvalues are continuous functions
of the parameters, we can see that if there is some
$\left(u_{f},u_{g},u_{h},\alpha\right)=\left(u_{f}^{\ast},u_{g}^{\ast},u_{h}^{\ast},\alpha^{\ast}\right)$
such that $a_{2}a_{1}=a_{0}$, then if we denote by $\lambda_{1}$ the always
present real eigenvalue, there is some $\delta>0$ such that for
$\left\|\left(u_{f},u_{g},u_{h},\alpha\right)-\left(u_{f}^{\ast},u_{g}^{\ast},u_{h}^{\ast},\alpha^{\ast}\right)\right\|<\delta$
we always have $\lambda_{1}<0$. By lemma $5.1$, section $5.2$ in [9], this
implies the existence of a parameter-dependent, smooth, attracting, two
dimensional, center manifold $W^{c}_{\left(u_{f},u_{g},u_{h},\alpha\right)}$.
In the following analysis, we consider only parameters that are in the
$\delta$-neighborhood of $\left(u_{f},u_{g},u_{h},\alpha\right)$ in order to
ensure that the real eigenvalue $\lambda_{1}$ is negative.

By the Routh-Hurwitz criterion (since we just shown that when
$\prescript{}{\left(111\right)}{\mathcal{E}}$ exists that $a_{0}$ and $a_{2}$
are always positive), a Hopf bifurcation occurs when the expression
$a_{2}a_{1}-a_{0}$ changes sign as a parameter varies. This ensures that the
real part of a pair of complex eigenvalues with nonzero imaginary part passes
through $0$ and hence changes sign. This is related to the transversality
condition: the derivative of the real part of the eigenvalue with respect to
the bifurcation parameter evaluated at the critical value when the real parts
are zero is non-zero. We check this condition for specific forms of the
functions $\mu_{0}$, $\mu_{1}$, and $\mu_{2}$. For the prototypes proposed in
(8), and with the values of parameters from Table 1 fixed, the function

$f\left(u_{f},u_{g},u_{h},\alpha\right)=a_{2}a_{1}-a_{0}$ (103)

is an algebraic function in $u_{f}$, $u_{g}$, $u_{h}$, and $\alpha$. Hence
fixing all the parameters except one makes the function $f$ a polynomial, the
order of which depends on the choice of the free parameter. Specifically if we
choose:

  * •

$\alpha$ \- free parameter $\implies$ $f$ has order $27$,

  * •

$u_{f}$ \- free parameter $\implies$ $f$ has order $3$,

  * •

$u_{g}$ \- free parameter $\implies$ $f$ has order $3$,

  * •

$u_{h}$ \- free parameter $\implies$ $f$ has order $2$.

We construct bifurcation diagrams to explore the possibility of Hopf
bifurcations as $\alpha$ varies in subsection 4.4, and begin our theoretical
analysis by choosing $u_{f}$ as the free parameter. We have

$f_{u_{f}}\left(u_{f}\right)=a_{2}a_{1}-a_{0}=b_{3}u_{f}^{3}+b_{2}u_{f}^{2}+b_{1}u_{f}+b_{0},$
(104)

and we assume that there is a value $u_{f}=u_{f}^{\ast}$ such that
$f_{u_{f}}(u_{f}^{\ast})=0$. We want to find conditions on the coefficients of
$f_{u_{f}}$ that guarantee that the derivative of $a_{2}a_{1}-a_{0}$ with
respect to $u_{f}$ is not equal to zero when $a_{2}a_{1}-a_{0}=0$, i.e., that
$u_{f}^{\ast}$ is not a local extremum of $f_{u_{f}}$. We have

$f^{\prime}_{u_{f}}\left(u_{f}\right)=3b_{3}u_{f}^{2}+2b_{2}u_{f}+b_{1}.$
(105)

The necessary condition for $u_{f}^{\ast}$ to be a local extremum for
$f_{u_{f}}$ is $f^{\prime}_{u_{f}}(u_{f}^{\ast})=0$, that is

$u_{f}^{\ast}=\frac{-b_{2}\pm\sqrt{b_{2}^{2}-3b_{1}b_{3}}}{3b_{3}}.$ (106)

We can derive sufficient conditions for $u_{f}^{\ast}$ to be an extremum (for
example by computing the second derivative of $f_{u_{f}}$), but the condition
(106) is already very restrictive and will be sufficient for our work. We have
thus obtained a sufficient condition for a Hopf bifurcation.

If we choose $u_{g}$ as the free parameter, we have

$f_{u_{g}}\left(u_{h}\right)=c_{3}u_{g}^{3}+c_{2}u_{g}^{2}+c_{1}u_{g}+c_{0},$
(107)

and with the assumption that $f_{u_{g}}\left(u_{g}^{\ast}\right)=0$, by a
similar analysis as in the previous case, we obtain an analogous sufficient
condition for a Hopf bifurcation in $u_{g}$.

Finally, if we choose $u_{h}$ as the free parameter, we have

$f_{u_{h}}\left(u_{h}\right)=d_{2}u_{h}^{2}+d_{1}u_{h}+d_{0}.$ (108)

Once again, we assume that $f_{u_{h}}\left(u_{h}^{\ast}\right)=0$, that is

$u_{h}^{\ast}=\frac{-d_{1}\pm\sqrt{d_{1}^{2}-4d_{2}d_{0}}}{2d_{2}}.$ (109)

Here, $u_{h}^{\ast}$ is the local extremum if and only if the discriminant of
equation (108) is zero, i.e., if $d_{1}^{2}-4d_{2}d_{0}=0$.

We summarize our results in the following theorem.

######  Theorem 1.

Consider system (21) with the prototypes given by (8) and with the values of
parameters from Table 1 fixed. Assume that there exists a point
$\left(u_{f},u_{g},u_{h},\alpha\right)=\left(u_{f}^{\ast},u_{g}^{\ast},u_{h}^{\ast},\alpha^{\ast}\right)$
such that
$f\left(u_{f}^{\ast},u_{g}^{\ast},u_{h}^{\ast},\alpha^{\ast}\right)=0$ for $f$
defined in (103). Then
$f_{u_{f}}\left(u_{f}\right)=f\left(u_{f},u_{g}^{\ast},u_{h}^{\ast},\alpha^{\ast}\right)$,
$f_{u_{g}}\left(u_{g}\right)=f\left(u_{f}^{\ast},u_{g},u_{h}^{\ast},\alpha^{\ast}\right)$,
and
$f\left(u_{h}\right)=f\left(u_{f}^{\ast},u_{g}^{\ast},u_{h},\alpha^{\ast}\right)$
are given by the equations (104), (107), and (108), respectively. Also, there
exists $\delta>0$ such that if
$\left\|\left(u_{f},u_{g},u_{h},\alpha\right)-\left(u_{f}^{\ast},u_{g}^{\ast},u_{h}^{\ast},\alpha^{\ast}\right)\right\|<\delta$,
then

  1. I.

if

$u_{f}^{\ast}\neq\frac{-b_{2}\pm\sqrt{b_{2}^{2}-3b_{1}b_{3}}}{3b_{3}},$ (110)

then there is a Hopf bifurcation in $u_{f}$ at $u_{f}=u_{f}^{\ast}$,

  2. II.

if

$u_{g}^{\ast}\neq\frac{-c_{2}\pm\sqrt{c_{2}^{2}-3c_{1}c_{3}}}{3c_{3}},$ (111)

then there is a Hopf bifurcation in $u_{g}$ at $u_{g}=u_{g}^{\ast}$,

  3. III.

if

$d_{1}^{2}-4d_{2}d_{0}\neq 0,$ (112)

then there is a Hopf bifurcation in $u_{h}$ at $u_{h}=u_{h}^{\ast}$.

We now illustrate the theoretical results with the numerical simulations. To
approximate values of the equilibria we used Maple software [11], rounding all
the values to $6$ significant digits. For the choice of parameters given in
Table 1 and

$\alpha=0.01,\quad u_{f}=0.5,\quad u_{g}=0.0006,$ (113)

it follows that the only zero of (108) occurs for $u_{h}=0.102520$. In Figure
1 we plot the phase space for $u_{h}=0.05$ (just before the Hopf bifurcation)
using the ode15s solver from [12]. For this set of parameters, we have the
following approximate values of the equilibria

$\begin{split}\prescript{}{\left(000\right)}{\mathcal{E}}&=\left(0,0,0\right)\quad\text{(unstable)},\\\
\prescript{}{\left(100\right)}{\mathcal{E}}&=\left(0.000299015,0,0\right)\quad\text{(stable)},\\\
\prescript{}{\left(001\right)}{\mathcal{E}}&=\left(0,0,0.0473693\right)\quad\text{(unstable)},\\\
\prescript{}{\left(110\right)}{\mathcal{E}}^{\left(1\right)}&=\left(0.0545964,0.00534486,0\right)\quad\text{(unstable)},\\\
\prescript{}{\left(110\right)}{\mathcal{E}}^{\left(2\right)}&=\left(0.489465,0.0487183,0\right)\quad\text{(unstable)},\\\
\prescript{}{\left(111\right)}{\mathcal{E}}&=\left(0.306611,0.0520205,36.2277\right)\quad\text{(unstable)},\end{split}$
(114)

Figure 1: Phase space of system (21) with $\alpha=0.01$, $u_{f}=0.5$,
$u_{g}=0.0006$, and $u_{h}=0.05$.

We can see that there exists a stable periodic orbit in the system, but
depending on the initial conditions, the solution might also converge to the
boundary equilibrium $\prescript{}{\left(100\right)}{\mathcal{E}}$. Thus in
this case we observe bistability.

We now repeat the simulations for $u_{h}=0.3$ (after the predicted Hopf
bifurcation), presented in Figure 2. With all the other parameters set to the
same values as in the previous case, we have the following equilibria

$\begin{split}\prescript{}{\left(000\right)}{\mathcal{E}}&=\left(0,0,0\right)\quad\text{(unstable)},\\\
\prescript{}{\left(100\right)}{\mathcal{E}}&=\left(0.00183202,0,0\right)\quad\text{(stable)},\\\
\prescript{}{\left(001\right)}{\mathcal{E}}&=\left(0,0,0.297369\right)\quad\text{(unstable)},\\\
\prescript{}{\left(110\right)}{\mathcal{E}}^{\left(1\right)}&=\left(0.0528596,0.00502299,0\right)\quad\text{(unstable)},\\\
\prescript{}{\left(110\right)}{\mathcal{E}}^{\left(2\right)}&=\left(0.489467,0.0485697,0\right)\quad\text{(unstable)},\\\
\prescript{}{\left(111\right)}{\mathcal{E}}&=\left(0.306611,0.0520205,36.4777\right)\quad\text{(stable)}.\end{split}$
(115)

Figure 2: Phase space of system (21) with $\alpha=0.01$, $u_{f}=0.5$,
$u_{g}=0.0006$, and $u_{h}=0.3$.

We can see that the Hopf bifurcation occurs between $u_{h}=0.05$ and
$u_{h}=0.3$. The stable periodic orbit is no longer present, and the interior
equilibrium is stable. Once again, the boundary equilibrium
$\prescript{}{\left(100\right)}{\mathcal{E}}$ is stable, and thus we observe
bistability in the system.

It is worth to notice, that increasing $u_{h}$ had a stabilizing effect on the
system (actually, this type of behaviour applies also to $u_{f}$ and $u_{g}$).
This result is especially important in the context of the modeled phenomenon,
since the most desirable situation happens when the production of methane is
not fluctuating. Variable rates of gas production can result in decreased
productivity of the biogas plant.

###  4.3 Persistence

The notion of persistence is particularly important in modeling biological
phenomena. Roughly speaking, we say that a system is persistent if all the
species with positive initial populations survive. The formal definition is as
follows.

######  Definition 4.1.

The system

$\begin{split}x_{i}^{\prime}&=x_{i}f_{i}\left(x_{1},x_{2},\ldots,x_{n}\right),\\\
x_{i}\left(0\right)&=x_{i0}\geq 0,\quad i=1,2,\ldots,n,\end{split}$ (116)

is said to be weakly persistent if

$\limsup\limits_{t\rightarrow\infty}x_{i}\left(t\right)>0,\quad
i=1,2,\ldots,n$ (117)

for every trajectory with positive initial conditions, and is said to be
persistent if

$\liminf\limits_{t\rightarrow\infty}x_{i}\left(t\right)>0,\quad
i=1,2,\ldots,n$ (118)

for every trajectory with positive initial conditions. This system is said to
be uniformly persistent if there exists a positive number $\epsilon$ such that

$\liminf\limits_{t\rightarrow\infty}x_{i}\left(t\right)\geq\epsilon,\quad
i=1,2,\ldots,n$ (119)

for every trajectory with positive initial conditions.

To prove that system (21) is persistent, we will use the Butler-McGehee lemma
[16] repeatedly.

######  Lemma 4.1.

Suppose that $x^{\ast}$ is a hyperbolic equilibrium point of the system

$\begin{split}x^{\prime}&=f(x),\\\ x\left(0\right)&=x_{0},\end{split}$ (120)

with $x\in\mathbb{R}^{n}$ and $f:\mathbb{R}^{n}\to\mathbb{R}^{n}$, where $f$
is continuously differentiable. Suppose also that $x^{\ast}$ is in
$\omega\left(x_{0}\right)$, the omega limit set of
$\gamma^{+}\left(x_{0}\right)$ (the positive semi-orbit through $x_{0}$), but
is not the entire omega limit set. Then $\omega\left(x_{0}\right)$ has
nontrivial (i.e., different from $x^{\ast}$) intersection with the stable and
unstable manifolds of $x^{\ast}$.

As we already noticed in section 4.2, there are values of the parameters where
one of the boundary equilibria and the interior equilibrium are both
asymptotically stable, and hence system (21) is not persistent, even though an
interior equilibrium point exists. We will thus focus on the cases for which
no boundary equilibrium point of system (21) is stable.

######  Theorem 2.

Let system (21) have the following equilibria configuration (as represented
schematically in Figure 3):

Equilibrium |  Number of eigenvalues with positive real part |  Number of eigenvalues with negative real part
---|---|---
$\prescript{}{\left(000\right)}{\mathcal{E}}$ |  2 |  1
$\prescript{}{\left(100\right)}{\mathcal{E}}$ |  1 |  2
$\prescript{}{\left(001\right)}{\mathcal{E}}$ |  1 |  2
$\prescript{}{\left(110\right)}{\mathcal{E}}$ |  1 |  2

Then system (21) is persistent.

Figure 3: Schematic representation of the equilibria configuration occurring
in the hypothesis of theorem 2. Black arrows represent stable and unstable
manifolds of each of the equilibrium (marked by the orange dots). In the
example for the parameters we select in model (21) to illustrate this theorem
(see (121)), there is an asymptotically stable interior equilibrium (as
shown). However, this is not necessary for the proof of theorem 2.

###### Proof.

Since planes $x_{0}x_{1}$, $x_{0}x_{2}$ and $x_{1}x_{2}$ are invariant, we
know where the stable and unstable manifolds of the boundary equilibria lie.
This is represented in a schematic way in Figure 3. Keeping this picture in
mind should make the following argument much more transparent. Assume that a
solution
$\vec{x}\left(t\right)=\left(x_{0}\left(t\right),x_{1}\left(t\right),x_{2}\left(t\right)\right)$
with an initial condition
$\vec{x}^{\left(0\right)}=\left(x_{0}^{\left(0\right)},x_{1}^{\left(0\right)},x_{2}^{\left(0\right)}\right)$,
where $x_{i}^{\left(0\right)}>0$, $i=1,2,3$, is given. First, suppose that
$\prescript{}{\left(000\right)}{\mathcal{E}}$ belongs to
$\omega\left(\gamma^{+}\left(\vec{x}^{\left(0\right)}\right)\right)$, the
omega limit set of $\gamma^{+}\left(\vec{x}^{\left(0\right)}\right)$. Since
$\prescript{}{\left(000\right)}{\mathcal{E}}$ is a saddle point with one-
dimensional stable manifold restricted to the $x_{1}$-axis, it is not the
entire omega limit set
$\omega\left(\gamma^{+}\left(\vec{x}^{\left(0\right)}\right)\right)$. Hence,
by lemma 4.1, there is a point
$x^{\ast}\neq\prescript{}{\left(000\right)}{\mathcal{E}}$ in both
$\omega\left(\gamma^{+}\left(\vec{x}^{\left(0\right)}\right)\right)$ and
$W^{s}\left(\prescript{}{\left(000\right)}{\mathcal{E}}\right)$, the stable
manifold of $\prescript{}{\left(000\right)}{\mathcal{E}}$. The entire orbit
through any point in an omega limit set is also in the omega limit set. The
stable manifold of $\prescript{}{\left(000\right)}{\mathcal{E}}$ is the
$x_{1}$-axis, and the $x_{1}$-axis is unbounded. We have already proven in
section 3 that all orbits of system (21) are bounded, and hence the omega
limit set of any orbit of (21) is bounded. This contradicts the existence of
such an $x^{\ast}$ and thus
$\prescript{}{\left(000\right)}{\mathcal{E}}\notin\omega\left(\gamma^{+}\left(\vec{x}^{\left(0\right)}\right)\right)$.

Now, suppose that
$\prescript{}{\left(001\right)}{\mathcal{E}}\in\omega\left(\gamma^{+}\left(\vec{x}^{\left(0\right)}\right)\right)$.
Since $\prescript{}{\left(001\right)}{\mathcal{E}}$ is a saddle point with
two-dimensional stable manifold restricted to the $x_{0}x_{1}$-plane,
$\left\\{\prescript{}{\left(001\right)}{\mathcal{E}}\right\\}$ is not the
entire omega limit set
$\omega\left(\gamma^{+}\left(\vec{x}^{\left(0\right)}\right)\right)$. Thus,
using lemma 4.1, there is a point
$x^{\ast}\in\omega\left(\gamma^{+}\left(\vec{x}^{\left(0\right)}\right)\right)\cap
W^{s}\left(\prescript{}{\left(001\right)}{\mathcal{E}}\right)\setminus\left\\{\prescript{}{\left(001\right)}{\mathcal{E}}\right\\}$.
Since the stable manifold
$W^{s}\left(\prescript{}{\left(001\right)}{\mathcal{E}}\right)$ lies entirely
in the $x_{1}x_{2}$-plane, and the entire orbit through $x^{\ast}$ is in
$\omega\left(\gamma^{+}\left(\vec{x}^{\left(0\right)}\right)\right)$, by the
analysis in subsection 4.1, this orbit becomes unbounded in backward time.
This contradiction shows that
$\prescript{}{\left(001\right)}{\mathcal{E}}\notin\omega\left(\gamma^{+}\left(\vec{x}^{\left(0\right)}\right)\right)$.

Now, suppose that
$\prescript{}{\left(100\right)}{\mathcal{E}}\in\omega\left(\gamma^{+}\left(\vec{x}^{\left(0\right)}\right)\right)$.
Similarly as in the previous cases, this implies that there exists a point
$x^{\ast}\in\omega\left(\gamma^{+}\left(\vec{x}^{\left(0\right)}\right)\right)\cap
W^{s}\left(\prescript{}{\left(100\right)}{\mathcal{E}}\right)\setminus\left\\{\prescript{}{\left(100\right)}{\mathcal{E}}\right\\}$.
This time the stable manifold
$W^{s}\left(\prescript{}{\left(100\right)}{\mathcal{E}}\right)$ is two-
dimensional and lies entirely in the $x_{0}x_{2}$-plane. By the analysis in
subsection 4.1, the entire orbit through $x^{\ast}$ (which belongs to
$\omega\left(\gamma^{+}\left(\vec{x}^{\left(0\right)}\right)\right)$) becomes
unbounded in backward time or its closure contains
$\prescript{}{\left(001\right)}{\mathcal{E}}$. This contradiction proves that
$\prescript{}{\left(100\right)}{\mathcal{E}}\notin\omega\left(\gamma^{+}\left(\vec{x}^{\left(0\right)}\right)\right)$.

Now, suppose that
$\prescript{}{\left(110\right)}{\mathcal{E}}\in\omega\left(\gamma^{+}\left(\vec{x}^{\left(0\right)}\right)\right)$.
Again, $\left\\{\prescript{}{\left(110\right)}{\mathcal{E}}\right\\}$ is not
the entire omega limit set
$\omega\left(\gamma^{+}\left(\vec{x}^{\left(0\right)}\right)\right)$, so there
exists a point
$x^{\ast}\in\omega\left(\gamma^{+}\left(\vec{x}^{\left(0\right)}\right)\right)\cap
W^{s}\left(\prescript{}{\left(110\right)}{\mathcal{E}}\right)\setminus\left\\{\prescript{}{\left(110\right)}{\mathcal{E}}\right\\}$.
This point lies in the $x_{0}x_{1}$-plane, since
$W^{s}\left(\prescript{}{\left(110\right)}{\mathcal{E}}\right)$ is two-
dimensional and is entirely contained in this plane. As in the previous cases,
the entire orbit through $x^{\ast}$ is in
$\omega\left(\gamma^{+}\left(\vec{x}^{\left(0\right)}\right)\right)$. Since
there are no periodic orbits in the $x_{0}x_{1}$ face, and since
$\left\\{\prescript{}{\left(100\right)}{\mathcal{E}}\right\\}\notin\omega\left(\gamma^{+}\left(\vec{x}^{\left(0\right)}\right)\right)$,
the orbit becomes unbounded in backward time. This contradiction proves that
$\prescript{}{\left(110\right)}{\mathcal{E}}\notin\omega\left(\gamma^{+}\left(\vec{x}^{\left(0\right)}\right)\right)$.

Finally, consider any
$\widehat{x}=\left(\widehat{x}_{0},\widehat{x}_{1},\widehat{x}_{2}\right)$,
such that $\widehat{x}_{i}=0$ for at least one $i=1,2,3$, and suppose that
$\widehat{x}\in\omega\left(\gamma^{+}\left(\vec{x}^{\left(0\right)}\right)\right)$.
Then, the entire orbit through $\widehat{x}$ is in
$\omega\left(\gamma^{+}\left(\vec{x}^{\left(0\right)}\right)\right)$. But
since this orbit lies entirely in either $x_{0}x_{1}$, $x_{1}x_{2}$, or
$x_{0}x_{2}$ face, it converges to one of the boundary equilibria. This
implies that this boundary equilibrium is in
$\omega\left(\gamma^{+}\left(\vec{x}^{\left(0\right)}\right)\right)$, and this
possibility has been eliminated in the previous part of the proof.

We have therefore proven that

$\liminf\limits_{t\rightarrow\infty}x_{i}\left(t\right)>0,\quad i=1,2,3,$

i.e., that system (21) is persistent. ∎

An example satisfying the assumptions of theorem 2 occurs for

$\alpha=0.0002,\quad u_{f}=0.6,\quad u_{g}=0,\quad u_{h}=0.1.$ (121)

Persistence can also be observed with the addition of phenol, i.e., with
$u_{g}>0$.

######  Theorem 3.

Let system (21) have the following equilibria configuration (as represented
schematically in Figure 4):

Equilibrium |  Number of eigenvalues with positive real part |  Number of eigenvalues with negative real part
---|---|---
$\prescript{}{\left(000\right)}{\mathcal{E}}$ |  2 |  1
$\prescript{}{\left(100\right)}{\mathcal{E}}$ |  1 |  2
$\prescript{}{\left(001\right)}{\mathcal{E}}$ |  1 |  2
$\prescript{}{\left(011\right)}{\mathcal{E}}$ |  1 |  2
$\prescript{}{\left(110\right)}{\mathcal{E}}$ |  1 |  2

Then system (21) is persistent.

Figure 4: Schematic representation of the equilibria configuration occurring
in the hypothesis of theorem 3. Black arrows represent stable and unstable
manifolds of each of the equilibrium (marked by the orange dots). In the
example for the parameters we select in model (21) to illustrate this theorem
(see (122)), there is an asymptotically stable interior equilibrium (as
shown). However, this is not necessary for the proof of theorem 3.

###### Proof.

The idea behind this proof is very similar to the method presented in the
proof of theorem 2. Let
$\vec{x}\left(t\right)=\left(x_{0}\left(t\right),x_{1}\left(t\right),x_{2}\left(t\right)\right)$
be a solution of (21) with an initial condition
$\vec{x}^{\left(0\right)}=\left(x_{0}^{\left(0\right)},x_{1}^{\left(0\right)},x_{2}^{\left(0\right)}\right)$,
where $x_{i}^{\left(0\right)}>0$, $i=1,2,3$. Since the stable and unstable
manifolds of the all of the equilibria, except
$\prescript{}{\left(001\right)}{\mathcal{E}}$ and
$\prescript{}{\left(011\right)}{\mathcal{E}}$, have the same configuration as
in the hypothesis of theorem 2, the argument eliminating them from the omega
limit set of $\gamma^{+}\left(\vec{x}^{\left(0\right)}\right)$ is exactly the
same and we only need to focus on
$\prescript{}{\left(001\right)}{\mathcal{E}}$ and
$\prescript{}{\left(011\right)}{\mathcal{E}}$ equilibria.

Suppose that
$\prescript{}{\left(001\right)}{\mathcal{E}}\in\omega\left(\gamma^{+}\left(\vec{x}^{\left(0\right)}\right)\right)$.
Since $\prescript{}{\left(001\right)}{\mathcal{E}}$ is a saddle point with
one-dimensional stable manifold restricted to the $x_{2}$-axis, we have
$\omega\left(\gamma^{+}\left(\vec{x}^{\left(0\right)}\right)\right)\setminus\left\\{\prescript{}{\left(001\right)}{\mathcal{E}}\right\\}\neq\emptyset$.
Hence, by lemma 4.1, there is a point
$x^{\ast}\in\omega\left(\gamma^{+}\left(\vec{x}^{\left(0\right)}\right)\right)\cap
W^{s}\left(\prescript{}{\left(001\right)}{\mathcal{E}}\right)\setminus\left\\{\prescript{}{\left(001\right)}{\mathcal{E}}\right\\}$.
The entire orbit through $x^{\ast}$, which also belongs to
$\omega\left(\gamma^{+}\left(\vec{x}^{\left(0\right)}\right)\right)$, either
becomes unbounded in backward time, or converges to the
$\prescript{}{\left(000\right)}{\mathcal{E}}$ equilibrium. Since all orbits of
system (21) are bounded, and
$\prescript{}{\left(000\right)}{\mathcal{E}}\notin\omega\left(\gamma^{+}\left(\vec{x}^{\left(0\right)}\right)\right)$,
we obtain a contradiction. Hence
$\prescript{}{\left(001\right)}{\mathcal{E}}\notin\omega\left(\gamma^{+}\left(\vec{x}^{\left(0\right)}\right)\right)$.

Now, suppose that
$\prescript{}{\left(011\right)}{\mathcal{E}}\in\omega\left(\gamma^{+}\left(\vec{x}^{\left(0\right)}\right)\right)$.
Since $\prescript{}{\left(011\right)}{\mathcal{E}}$ is a saddle point with
two-dimensional stable manifold restricted to the $x_{1}x_{2}$-plane (it is
repelling into the interior), we have
$\omega\left(\gamma^{+}\left(\vec{x}^{\left(0\right)}\right)\right)\setminus\left\\{\prescript{}{\left(011\right)}{\mathcal{E}}\right\\}\neq\emptyset$.
By using lemma 4.1, there exists a point
$x^{\ast}\in\omega\left(\gamma^{+}\left(\vec{x}^{\left(0\right)}\right)\right)\cap
W^{s}\left(\prescript{}{\left(011\right)}{\mathcal{E}}\right)\setminus\left\\{\prescript{}{\left(011\right)}{\mathcal{E}}\right\\}$.
The entire orbit through $x^{\ast}$, which also belongs to
$\omega\left(\gamma^{+}\left(\vec{x}^{\left(0\right)}\right)\right)$, either
becomes unbounded in backward time, or converges to
$\prescript{}{\left(000\right)}{\mathcal{E}}$, or
$\prescript{}{\left(001\right)}{\mathcal{E}}$ (we have previously shown in
subsection 4.1 that there are no periodic orbits in the $x_{1}x_{2}$ face).
Since we have already proven that
$\prescript{}{\left(000\right)}{\mathcal{E}}\notin\omega\left(\gamma^{+}\left(\vec{x}^{\left(0\right)}\right)\right)$,
and
$\prescript{}{\left(001\right)}{\mathcal{E}}\notin\omega\left(\gamma^{+}\left(\vec{x}^{\left(0\right)}\right)\right)$,
we obtain a contradiction, which proves that
$\prescript{}{\left(011\right)}{\mathcal{E}}\notin\omega\left(\gamma^{+}\left(\vec{x}^{\left(0\right)}\right)\right)$.

Finally, consider any
$\widehat{x}=\left(\widehat{x}_{0},\widehat{x}_{1},\widehat{x}_{2}\right)$,
such that $\widehat{x}_{i}=0$ for at least one $i=1,2,3$, and suppose that
$\widehat{x}\in\omega\left(\gamma^{+}\left(\vec{x}^{\left(0\right)}\right)\right)$.
Then, the entire orbit through $\widehat{x}$ is in
$\omega\left(\gamma^{+}\left(\vec{x}^{\left(0\right)}\right)\right)$. But
since this orbit lies entirely in either $x_{0}x_{1}$, $x_{1}x_{2}$, or
$x_{0}x_{2}$ face, it converges to one of the boundary equilibria. This
implies that this boundary equilibrium is in
$\omega\left(\gamma^{+}\left(\vec{x}^{\left(0\right)}\right)\right)$, and this
possibility has been eliminated in the previous part of the proof. ∎

An example satisfying the assumptions of theorem 3 occurs for

$\alpha=0.0002,\quad u_{f}=0.6,\quad u_{g}=0.00015,\quad u_{h}=0.1.$ (122)

######  Remark 1.

Interestingly enough, in many cases of models describing biological phenomena,
persistence already implies uniform persistence. The rigorous results were
obtained in [2]. In our context, the key theorem from [2] states that if
$\mathcal{F}$ is a dynamical system for which $\mathbb{R}^{n}_{+}$ and
$\partial\mathbb{R}^{n}_{+}$ are invariant, then $\mathcal{F}$ is uniformly
persistent provided that

  1. 1.

$\mathcal{F}$ is dissipative (meaning that $\forall x\in\mathbb{R}^{n}_{+}$
$\omega\left(x\right)\neq\emptyset$ and
$\bigcup_{x\in\mathbb{R}^{n}_{+}}\omega\left(x\right)$ has compact closure),

  2. 2.

$\mathcal{F}$ is weakly persistent,

  3. 3.

$\partial\mathcal{F}$ (the restriction of $\mathcal{F}$ to the boundary
$\partial\mathbb{R}^{n}_{+}$) is "isolated",

  4. 4.

$\partial\mathcal{F}$ is "acyclic".

These results can be easily modified so that we consider the flow
$\mathcal{F}$ on $\Omega$ defined in (19). Although $\partial\Omega$ is not
invariant, the theorem from [2], as explained in [3], can be modified so that
it applies in the case when $\partial\Omega$ is the union of two sets
$\Omega_{1}$ and $\Omega_{2}$, for which $\mathcal{F}$ is invariant on
$\Omega_{1}$ and $\Omega_{2}$ is repelling into the interior of $\Omega$,
provided that conditions $3.$ and $4.$ are satisfied for the restriction of
$\mathcal{F}$ to $\Omega_{1}$. In our case, the positively invariant set
$\Omega$, on which we analyze system (21) is bounded, hence condition $1.$ is
satisfied. Condition $2.$ holds by theorem 2 (persistence implies weak
persistence). In our context condition 3. is satisfied, because all the
boundary equilibria are hyperbolic, and hence each one is the maximal
invariant set in a neighbourhood of itself. Also, their union forms a covering
of the omega limit sets of $\Omega_{1}$. Condition 4. is satisfied because the
boundary equilibria are not cyclically linked, i.e., there is no cyclic chain
of heteroclinic orbits joining them. Thus, we have shown not only persistence,
but also uniform persistence of system (21) in the case of theorem 2 and
theorem 3.

We have thus proven the following theorem:

######  Theorem 4.

Under the hypotheses of theorems 2 and 3, system (21) is uniformly persistent.

We finish this subsection by extending the uniform persistence to the original
six dimensional system (1). Notice that if
$\vec{X_{0}}=\left(x_{0}(0),x_{1}(0),x_{2}(0),s_{0}(0),s_{1}(0),s_{2}(0)\right)$
with $x_{i}(0)\geq 0$, $s_{i}(0)\geq 0$, $i=1,2,3$, then we necessarily must
have $\omega\left(\vec{X_{0}}\right)\in\Omega$. Otherwise, there would exist a
point $\vec{p}\in\mathbb{R}_{+}\setminus\Omega$ and a sequence of times
$\left(t_{n}\right)$ with $t_{n}\rightarrow\infty$ for which the corresponding
solution converges to $\vec{p}$. This would mean that $\Omega$ is not globally
attracting, which was proven in section 3. Also, if
$\omega\left(\vec{X_{0}}\right)$ has a point lying in a face with one of the
$x_{i}$, $i=1,2,3$ equal zero, then the entire orbit through that point would
also be in $\omega\left(\vec{X_{0}}\right)$. Thus, if the assumptions of
either theorem 2 or theorem 3 hold, the omega limit set
$\omega\left(\vec{X_{0}}\right)$ is entirely contained in the interior of
$\Omega$. We have thus proven the following theorem:

######  Theorem 5.

Under the hypotheses of theorems 2 and 3, the six dimensional system (1) is
uniformly persistent.

###  4.4 Bifurcation diagrams

As previously stated in section 4.2, we now study numerically effects on the
qualitative behaviour of system (21) when considering $\alpha$ as the
bifurcation parameter. Throughout this section, we assume that parameters
$\omega_{0}$, $\omega_{1}$, $\omega_{2}$, $\phi_{1}$, $\phi_{2}$, $K_{P}$, and
$K_{I}$ are fixed at the values given in Table 1.

We now fix the following parameters

$u_{f}=2,\quad u_{g}=0,\quad u_{h}=0,$ (123)

and plot a one-parameter bifurcation diagram in $\alpha$, with $x_{0}$ on the
$y$-axis. All simulations were performed using [19].

Figure 5: One-parameter bifurcation diagram of system (21) with $\alpha$ as
the bifurcation parameter and $u_{f}=2$, $u_{g}=0$, $u_{h}=0$. Figure 6:
Close-up of on the one-parameter bifurcation diagram represented in Figure 5.

We can see that as $\alpha$ decreases, there is a saddle-node bifurcation,
resulting in two equilibria
$\prescript{}{\left(110\right)}{\mathcal{E}}^{\left(1\right)}$ and
$\prescript{}{\left(110\right)}{\mathcal{E}}^{\left(2\right)}$ appearing (both
unstable). Next, there is a transcritical bifurcation with the
$\prescript{}{\left(110\right)}{\mathcal{E}}^{\left(1\right)}$ equilibrium,
which results in the positive equilibrium coming into the interior of the
admissible region $\Omega$. After that, a saddle-node of limit cycles
bifurcation occurs, which gives birth to a stable and unstable periodic
orbits. The $\prescript{}{\left(111\right)}{\mathcal{E}}$ equilibrium
(unstable), undergoes a Hopf bifurcation, and as a consequence it becomes
asymptotically stable, and the stable periodic orbits disappears. Since these
bifurcations occur for a narrow range of $\alpha$, a close-up is presented in
Figure 6. Stable periodic orbit represents a case in which all three
populations oscillate indefinitely, and hence the production of methane
fluctuates. As already mentioned in section 4.2, this situation is not a
desirable one, because it might result it decreased productivity of the biogas
plant. The unstable periodic orbits acts as a separatrix, giving the border of
the basin of attraction of two asymptotically stable equilibria in the case of
bistability.

Since by the conservation principles (20), $s_{0}=u_{f}-x_{0}$, the
bifurcation diagram in $\alpha$ with $s_{0}$ on the $y$-axis is similar to the
one presented in Figure 5. The amount of chlorophenol in the system is
inversely proportional to the concentration of the phenol degrader. As the
dilution rate $\alpha$ decreases, concentration of the chlorophenol degrader
in the interior equilibrium increases, and as a consequence, the amount of
chlorophenol decreases. It thus suggests, that operating on lower dilution
rates results in the most desirable dynamics, i.e., an asymptotically stable
interior equilibrium and fast chlorophenol removal.

To extend the previous analysis, we now fix the following parameters

$u_{g}=0,\quad u_{h}=0.1,$ (124)

and plot a two-parameter bifurcation diagram of system (21), choosing $\alpha$
and $u_{f}$ as the bifurcation parameters. Each region of the diagram is
labeled and the corresponding dynamics are represented schematically in
figures around it. Black dashed curve corresponds to saddle-node of equilibria
bifurcation (LP), black solid curve represents saddle-node of limit cycles
bifurcation (SNLC), black dotted curve denotes Hopf bifurcation (HB), and grey
solid curves represent transcritical bifurcations (BP). We also denote the
predicted heteroclinic bifurction by a grey dashed curve which lies very close
to the Hopf curve.

Figure 7: Two-parameter bifurcation diagram of system (21) with $u_{g}=0$ and
$u_{h}=0.1$. Figure 8: Close-up on region I of Figure 7. The SNLC curve
intersects the HB curve at Bautin bifurcation. This results in the change of
criticality of the Hopf bifurcation from supercritical (on the left) to
subcritical (on the right). Figure 9: Close-up on region II of Figure 7.

We can see that varying two parameters at the same time can lead to a much
more complicated dynamics than in the case of one-parameter bifurcations.
There is a generalized Hopf bifurcation, at the point at which the Hopf curve
intersects the saddle-node of limit cycles curve. This is the point where the
criticality of the Hopf bifurcation changes from supercritical to subcritical,
looking from left to right. The unstable periodic orbit disappears through a
heteroclinic bifurcation. There are two heteroclinic orbits that form a cycle
that joins the two equilibria in the $x_{0}x_{1}$ face, then passes into the
interior, and then goes back to the boundary in the $x_{0}x_{1}$ face. The
point at which the Hopf, homoclinic, and saddle-node of limits cycles curves
intersect, represents the Bogdanov-Takens bifurcation.

From the biological viewpoint, the most interesting dynamics is occurs in
regions $5$ and $iii$. There, the interior equilibrium is asymptotically
stable. In the case of region $5$ we also observe bistability with the
$\prescript{}{\left(100\right)}{\mathcal{E}}$ equilibrium. In region $iii$,
there is uniform persistence, and thus an interior compact attractor is
present. As was previously anticipated by the analysis of the one-parameter
bifurcation diagram, operating at low dilution rates is the most desirable
approach. If $\alpha$ is small enough, it is possible to remain in region
$iii$, even for high inflow rate $u_{f}$.

##  5 Conclusions

In this work we have generalised the approach presented in [14] by including
multiple substrate inflow into the chemostat, while maintaining generality (in
most cases) with respect to the exact form of the growth functions. We
observed that allowing the inflow of multiple substrates resulted in much more
complex dynamics of the system. For example, eight steady states are possible.
Previously, the theoretical results were limited to existence and uniqueness
of up to three equilibria (when chlorophenol was the only input substrate),
and to numerical evidence suggesting that the model should be subjected to a
more detailed analysis. We also observed that external addition of substrates
can result in bistability - two equilibria can simultaneously be
asymptotically stable. As well, there can be an orbitally asymptotically
stable periodic orbit with all of the populations surviving and an
asymptotically stable equilibrium with only chlorophenol degrader population
surviving.

We have also confirmed the findings of the previous analysis in [14], where
numerical evidence of the occurrence of a supercritical Hopf bifurcation was
given. Theoretical conditions for the existence of a Hopf bifurcation were
provided in the case of specific forms of the growth functions. Varying any
one of the three parameters: chlorophenol, phenol, and hydrogen inflow rates,
was shown to result in a Hopf bifurcation. Theoretical results for varying the
dilution rate as the bifurcation parameter has been left for future work.
However, we have observed numerically, that varying this parameter can result
in a Hopf bifurcation and a saddle-node of limit cycles bifurcation. Our
numerical investigations also showed that increasing the inflow rate of the
substrates has a stabilizing effect on the entire system. From a biological
engineering point of view, i.e. a bioreactor treating a monochlorophenol rich
waste stream, instability would typically be undesirable in terms of process
performance. Therefore, identification of control strategies to avoid periodic
behaviour is an important output of this work.

Another result, particularly important for engineering applications, concerns
the persistence of the system for a range of parameter sets. Knowing when the
microbial populations survive is again crucial from a process control
perspective, and it is one of the main theoretical results of this work. We
have proven that in two configurations of equilibria (in both cases all the
boundary equilibria are saddle points) we observe not only persistence, but
also uniform persistence, a much stronger result. These situations occur when
there is an inflow of all three substrates, but also when phenol addition is
not considered (i.e., when $u_{g}=0$).

Although we now know much more about the dynamics of the system, it is not
fully understood. This follows from the numerical results provided by the two-
parameter bifurcation diagrams. The analyses reveals that varying the dilution
rate and the chlorophenol inflow simultaneously, can lead to a Bogdanov-
Takens, or Bautin (generalized Hopf) bifurcations. Also, for the cases of
bistability, where both a boundary and the interior equilibrium are
asymptotically stable, it is of great importance to empiricists to have an
estimation of the basins of attraction of these equilibria. This result is
usually difficult to obtain theoretically, however numerical estimations are
possible. Another factor that is of interest would be the inclusion of
stochasticity in the model. In practice, even if the interior equilibrium is
globally asymptotically stable, one of the microorganisms may become extinct.
This might occur when a population is very small, and the stochastic noise
effects result in the solution curve reaching one of the invariant faces of
the admissible region.

There has been resistance to the idea that simplified models, of the type
described here, are too remote from the systems they represent to be of worth
to practitioners. Without experimental results to compare against model
predictions, this case becomes stronger. However, we can look to emerging
disciplines such as synthetic biology to help bridge the theoretical and the
applied [7]. Recent studies have shown that synthetically derived anaerobic
communities are able to confirm model predictions and provide insight into the
ecology and dynamics of microbial communities that are relevant in practice
[5]. We believe this work provides a basis by which experimental studies
describing a chlorophenol mineralising food-web could be undertaken.

##  6 Acknowledgments

M.J.W. acknowledges the support from the European Union’s Horizon 2020
research and innovation programme under the Marie Skłodowska-Curie Grant
Agreement No. 702408 (DRAMATIC).

## References

  * [1] I. Bassani, P. G. Kougias, L. Treu, and I. Angelidaki, Biogas upgrading via hydrogenotrophic methanogenesis in two-stage continuous stirred tank reactors at mesophilic and thermophilic conditions, Environ. Sci. Technol., 49 (2015), pp. 12585–12593.
  * [2] G. J. Butler, H. I. Freedman, and P. Waltman, Uniformly persistent systems, Proceedings of the American Mathematical Society, 96 (1986), pp. 425–430.
  * [3] G. J. Butler and G. S. K. Wolkowicz, Predator — mediated coexistence in a chemostat: Coexistence and competition reversal, Mathematical Modelling, 8 (1987), pp. 781 – 785.
  * [4] J. Chen, M. J. Wade, J. Dolfing, and O. S. Soyer, Increasing sulfate levels show a differential impact on synthetic communities comprising different methanogens and a sulfate reducer, J. Royal Soc. Interface, 16 (2019), p. 20190129.
  * [5] H. Delattre, J. Chen, M. Wade, and O. Soyer, Thermodynamic modelling of synthetic communities predicts minimum free energy requirements for sulfate reduction and methanogenesis, bioRxiv, (2019), p. 857276.
  * [6] M. El Hajji, N. Chorfi, and M. Jleli, Mathematical modelling and analysis for a three-tiered microbial food web in a chemostat, Electron. J. Differ. Eq., 255 (2017), pp. 1–13.
  * [7] M. El Karoui, M. Hoyos-Flight, and L. Fletcher, Future Trends in Synthetic Biology — A Report, Front. Bioeng. Biotechnol., 7 (2019), p. 175.
  * [8] T. Großkopf and O. Soyer, Microbial diversity arising from thermodynamic constraints, ISME J., 10 (2016), pp. 2725–2733.
  * [9] Y. Kuznestov, Elements of Applied Bifurcation Theory, vol. 112, Springer-Verlag, New York, 2004.
  * [10] L. Levén, K. Nyberg, and A. Schn urer, Conversion of phenols during anaerobic digestion of organic solid waste – a review of important microorganisms and impact of temperature, J. Env. Manage., 95 (2012), pp. 99–103.
  * [11] Maple, version 18.02, Waterloo Maple Inc., Waterloo, Ontario, 2018.
  * [12] MATLAB, version 9.5.0.944444 (R2018b), The MathWorks Inc., Natick, Massachusetts, 2018.
  * [13] C. Mazur, W. Jones, and C. Tebes-Stevens, H2 consumption during the microbial reductive dehalogenation of chlorinated phenols and tetrachloroethene, Biodegradation, 14 (2003), pp. 285–295.
  * [14] T. Sari and M. J. Wade, Generalised approach to modelling a three-tiered microbial food-web, Mathematical Biosciences, 291 (2017), pp. 21–37.
  * [15] B. Schink, Energetics of syntrophic cooperation in methanogenic degradation, Microbiol. Mol. Biol. Rev., 61 (1997), pp. 262–280.
  * [16] H. L. Smith and P. Waltman, The Theory of the Chemostat: Dynamics of Microbial Competition, Cambridge Studies in Mathematical Biology, Cambridge University Press, 1995.
  * [17] N. W. Smith, P. R. Shorten, E. H. Altermann, N. C. Roy, and W. C. McNabb, Hydrogen cross-feeders of the human gastrointestinal tract, Gut Microbes, 10 (2019), pp. 270–288.
  * [18] M. J. Wade, R. W. Pattinson, N. G. Parker, and J. Dolfing, Emergent behaviour in a chlorophenol-mineralising three-tiered microbial ’food web’, Journal of Theoretical Biology, 389 (2016), pp. 171–186.
  * [19] XPPAUT, version 8.0, Dr. Bard Ermentrout, Dept of Mathematics, University of Pittsburgh, Pittsburgh PA, 2016.