A model of hepatitis B virus with random interference infection rate

1.   Introduction

A mathematical model of the complete anaerobic mineralisation of a generic monochlorophenol isomer (C$_6$H$_4$ClOH) by a canonical food-web has recently been described ^[1]. The system comprises three microbial species (chlorophenol degrader, phenol degrader, hydrogenotrophic methanogen). Their interactions are summarised as follows:

● Reductive dehalogenation of chlorophenol in the presence of hydrogen by the chlorophenol degrader producing phenol ^[2];

● Phenol is mineralised by a phenol degrader to acetate and hydrogen via the benzoyl-CoA pathway or a caproate intermediary ^[3];

● 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 ^[4,5]. Hydrogen scavengers such as the hydrogenotrophic 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. These positive and negative feedback loops reframe the ecological network from a simple food-chain to a more complex food-web that allows for the possibility of periodicity. This additionally allows for the case where the system reduces to a self-sustaining two-species network in which the chlorophenol degrader acts as the syntrophic partner to the phenol degrader.

A diagrammatic representation of the food-web is presented in 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 do 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 ^[6,7,8]. 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) ^[1], 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 constant using a general representation of the species growth functions ^[9]. With no decay and without considering production of hydrogen by the phenol degrader, (an unusual constraint that appears to make the system entirely theoretical with extraneous addition of hydrogen necessary for complete chlorophenol mineralisation) local stability and the conditions giving rise to asymptotic coexistence of all three species have been shown analytically, with the possibility of periodic orbits not excluded ^[10]. Numerical analysis has suggested the presence of a Hopf bifurcation of the positive steady-state, with the concentration of influent chlorophenol as the bifurcating parameter ^[9].

In this work, we extend the analysis of the model to provide a proof of the existence and uniqueness of system equilibria and their stability for the most general case where up to three external inputs (chlorophenol, phenol, hydrogen) may be included. The procedure we use allows us to identify sufficient conditions for the emergence of a Hopf bifurcation as the inflow concentration parameters vary. We are also able to prove that there is a range of parameters for which the model is uniformly persistent, i.e., all populations avoid extinction, a new result for the system. Two-parameter bifurcation diagrams reveal rich dynamics that include Bogdanov-Takens and Bautin bifurcations and, hence, homoclinic and saddle-node of limit cycle bifurcations.

2.   The model revisited

We first present concisely the original given by ^[1] using the identical scaling:

with

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 ^[1]. 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 $\cdot_\mathrm{ch}$, $\cdot_\mathrm{ph}$, $\cdot_\mathrm{H_2}$ relate to the chlorophenol degrader, phenol degrader, and methanogen, respectively.

For the numerical bifurcation analysis given in Section 4.4, we consider the same parameter values as provided in the original work in ^[1], as shown in Table 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 (2.1), except that we assume $k_{i} = 0$,

$i\in\{A, B, C\}$:

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:

(H1) 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;

(H2) 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;

(H3) For all $s_2\geq 0$, $\mu_1\left(0, s_2\right) = 0$. As a consequence $\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;

(H4) 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;

(H5) $\mu_2\left(0\right) = 0$ and $\mu_2'\left(s_2\right) > 0$ for all $s_2 > 0$. Thus we assume that the methanogen 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 Eq (2.2), which satisfy these conditions, when we are unable to prove results in general, and when providing numerical simulations or bifurcation diagrams.

The proof of the following lemma, needed for well-posedness of system (3.7), is standard, and hence omitted.

Lemma 3.1. All solutions of system (3.1) with positive initial conditions remain positive and bounded for all positive times and if $x_i\left(0\right) = 0$, $i\in\left\{0, 1, 2\right\}$, then $x_i\left(t\right) = 0$ for all $t\geq 0$.

We can further simplify analysis of system (3.1) using the following lemma.

Lemma 3.2. The omega limit set of any solution of system (3.1) is contained in a positively invariant attracting set $\Omega \subset \mathbb{R}^6$ defined as

Proof. By adding the first and the fourth equations of system (3.1), we obtain

hence

which implies that

Similarly, we obtain

and

By taking the limit as $t \to \infty$ in Eqs (3.3)–(3.5) we obtain

which yields the desired result.

The three equalities in the definition of $\Omega$ given in system (3.2) are called "conservation principles". Using them, we can compute $s_0$, $s_1$, and $s_2$ as functions of $x_0$, $x_1$, $x_2$:

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

From now on, we will study the reduced system (3.7). We begin by analysing 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 Eq (3.7) equal to zero. Below, we list all the possibilities obtained this way. Since Eq (3.6) give a one-to-one correspondence of the equilibria of system (3.7) with the equilibria of system (3.1), 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. Each equilibrium is preceded by a subscript, indicating which population of microorganisms survives. That is, subscript $(100)$ implies that $x_0 > 0$, $x_1 = x_2 = 0$, $(110)$ that $x_0$, $x_1 > 0$, $x_2 = 0$, and so on.

Types of equilibria of system (3.7):

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

In this case, all the populations die out, hence the only source for the substrates comes from the inflows $u_f$, $u_g$ and $u_h$.

● Boundary equilibria:

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

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

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$ to obtain the total phenol amount $s_1$. Since $x_0$ consumes hydrogen as well, the value $\omega_{2_{\left(100\right)}}{x_0}$ is subtracted from $s_2$. This steady state is not desirable, because of the phenol build-up in the system.

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

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

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, _{\left(001\right)}{x_2}\right)$, where $x_2 = _{\left(001\right)}{x_2} > 0$ is a solution (if it exists) of

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

Here, only the methanogen is present, hence the values of chlorophenol and phenol are equal to the the inflow concentrations $u_f$ and $u_g$, respectively.

- $_{\left(101\right)}\mathcal{E} = \left(_{\left(101\right)}{x_0}, 0, -\omega_{2_{\left(101\right)}}{x_0}+u_h-\mu_2^{-1}\left(\alpha\right)\right)$,

where $x_0 = _{\left(101\right)}{x_0} > 0$ is a solution (if it exists) of

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

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

- $_{\left(011\right)}\mathcal{E} = \left(0, _{\left(011\right)}{x_1}, \omega_{1 _{\left(011\right)}}{x_1} + u_h -\mu_2^{-1}\left(\alpha\right)\right)$,

where $x_1 = _{\left(011\right)}{x_1} > 0$ is a solution (if it exists) of

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

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(_{\left(110\right)}{x_0}, _{\left(110\right)}{x_1}, 0\right)$, where $x_0 = _{\left(110\right)}{x_0} > 0$ and $x_1 = _{\left(110\right)}{x_1} > 0$ are solutions of

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

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

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

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.

3.2.   Existence, uniqueness, and local stability

The issues regarding the conditions on existence, uniqueness, and local stability of equilibria of the system have been fully answered in two recent preprints ^[11,12]. The system considered in these works is more general, as it allows non-zero death rates. It is related to our system (3.1) by a linear change of variables, that is, by replacing $x_0$ by $\omega_0 \omega_1 x_0$, $s_0$ by $\omega_0 \omega_1 s_0$, $x_1$ by $\omega_1 s_1$, $s_1$ by $\omega_1 s_1$, and replacing $\omega_0$, $\omega_1$, and $\omega_2$ by a single parameter $\omega = \frac{\omega_2}{\omega_0\omega_1}$. It follows that all results described in ^[11,12] apply to system (3.1) under this rescaling. Specifically, it has been shown that all equilibria, except $_{\left(110\right)}{\mathbb{E}}$, are unique, when they exist. Also, system (3.1) can have at most two equilibria of the form $_{\left(110\right)}{\mathbb{E}}$. Note that, following ^[1], different steady-state notation is adopted in ^[11,12], i.e., $_{\left(000\right)}{\mathbb{E}} = \mathrm{SS1}$, $_{\left(001\right)}{\mathbb{E}} = \mathrm{SS2}$, $_{\left(100\right)}{\mathbb{E}} = \mathrm{SS3}$, $_{\left(110\right)}{\mathbb{E}} = \mathrm{SS4}$, $_{\left(101\right)}{\mathbb{E}} = \mathrm{SS5}$, $_{\left(111\right)}{\mathbb{E}} = \mathrm{SS6}$, $_{\left(010\right)}{\mathbb{E}} = \mathrm{SS7}$, $_{\left(011\right)}{\mathbb{E}} = \mathrm{SS8}$. For more details, we refer readers to the aforementioned preprints, as this paper is focused more on possible bifurcations and persistence analysis of system (3.7).

4.   Analysis of the full system

4.1.   Periodic orbits on the faces

We begin by ruling out the 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$.

Lemma 4.1. There are no periodic orbits on the invariant faces of the set $\Omega$ of system (3.7).

Proof. First, consider system (3.7) on the part of $\Omega$ with $x_0 = 0$, i.e..

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

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

Then

for all $\left(x_1, x_2\right)$ in the domain of $\varphi_0$, by Conditions (H4) and (H5). Thus, by the Dulac's Criterion ^[13], there are no periodic orbits in the $x_1x_2$ face.

Now, consider system (3.7) on the part of $\Omega$ with $x_1 = 0$, i.e..

defined on $\Omega_{02}$ given by

For auxiliary function

we have, using Conditions (H2) and (H5),

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_0x_2$ face.

Finally, consider system (3.7) on the part of $\Omega$ with $x_2 = 0$, i.e..

defined on $\Omega_{01}$ given by

Analogously to the previous cases, we define

and note that by Conditions (H2) and (H4)

Since the axes $x_0 = 0$, $x_1 = 0$ are invariant, by Dulac's Criterion, there are no periodic orbits in the $x_0x_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 ^[9]. The occurrence of a stable periodic orbit in system (3.7) represents a situation in which the concentrations of all three populations of microorganisms oscillate indefinitely, and as a consequence, the substrate concentrations fluctuate.

The characteristic polynomial of the Jacobian corresponding to the interior equilibrium is given by

where

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 by Conditions (H2), (H4) and (H5)), but $a_1$ might possibly change sign. Notice first that since the polynomial is of 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

In that case, we also have

Hence, the eigenvalues are $\lambda_1 = -a_2$ and $\lambda_{2, 3} = \pm\sqrt{a_1}i$. In the following theorem, we prove that under certain nondegeneracy conditions, system (3.7) exhibits a Hopf bifurcation when the eigenvalues $\lambda_{1, 2}$ cross the imaginary axis. We explore the possibility of a Hopf bifurcation in either of the three parameters $u_f$, $u_g$, or $u_h$. Additionally, the effect of $\alpha$ on possible bifurcations in the system is explored numerically, by constructing bifurcation diagrams in Subsection 4.4. To derive our results, we use specific forms of the growth functions (2.2) and the values of the parameters given in Table 1.

Theorem 4.1. Consider system (3.7) with the prototypes given by Eq (2.2) and with the values of parameters from Table 1 fixed. Define $f:\mathbb{R}_+^4\to \mathbb{R}$ as $f(u_f, u_g, u_h, \alpha) = a_2a_1-a_0$ where $a_0$, $a_1$, $a_2$ are given in Eqs (4.2)–(4.4). Furthermore, 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 fixed $u_g = u_g^\ast$, $u_h = u_h^\ast$, and $\alpha = \alpha^\ast$, denote the resulting third-order polynomial $f(u_f, u_g^\ast, u_h^\ast, \alpha^\ast)$ as

for fixed $u_f = u_f^\ast$, $u_h = u_h^\ast$, and $\alpha = \alpha^\ast$, denote the resulting third-order polynomial $f(u_f^\ast, u_g, u_h^\ast, \alpha^\ast)$ as

and for fixed $u_f = u_f^\ast$, $u_g = u_g^\ast$, and $\alpha = \alpha^\ast$, denote the second-order polynomial $f(u_f^\ast, u_g^\ast, u_h, \alpha^\ast)$ as

There exists a $\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$ and in addition:

I.

then there is a Hopf bifurcation in $u_f$ when $u_f = u_f^\ast$, or

II.

then there is a Hopf bifurcation in $u_g$ when $u_g = u_g^\ast$, or

III.

then there is a Hopf bifurcation in $u_h$ when $u_h = u_h^\ast$.

Proof. 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 $f\left(u_f^\ast, u_g^\ast, u_h^\ast, \alpha^\ast\right) = 0$ (i.e., $a_2a_1 = a_0$), then if we denote by $\lambda_1$ the real eigenvalue (that always exists), 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 ^[13], 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$-neighbourhood 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 when $_{\left(111\right)}{\mathcal{E}}$ exists, $a_0$ and $a_2$ are always positive), a Hopf bifurcation occurs when the function $f$ 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 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$, given in Eq (2.2). We begin our analysis by choosing $u_f$ as the free parameter. We have

for some constants $b_i$, $i\in\{0, 1, 2, 3\}$, and by the hypothesis, for $u_f = u_f^\ast$, we have $f\left(u_f^\ast, u_g^\ast, u_h^\ast, \alpha^\ast\right) = 0$. We want to find conditions on the coefficients of $f$ that guarantee that its derivative with respect to $u_f$ is not equal to zero when $u_f = u_f^\ast$. We have

Thus $\partial_{u_f}f\left(u_f, u_g^\ast, u_h^\ast, \alpha^\ast\right) \neq 0$ if and only if

We have thus obtained a sufficient condition for a Hopf bifurcation.

If we choose $u_g$ as the free parameter, we have

for some constants $c_i$, $i\in\{0, 1, 2, 3\}$, and with the assumption that $f\left(u_f^\ast, u_g^\ast, u_h^\ast, \alpha^\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

for some constants $d_i$, $i\in\{0, 1, 2\}$. Once again, we assume that $f\left(u_f^\ast, u_g^\ast, u_h^\ast, \alpha^\ast\right) = 0$, that is

Here, $\partial_{u_h}f\left(u_f^\ast, u_g^\ast, u_h^\ast, \alpha^\ast\right)\neq 0$ if and only if $d_1^2-4d_2d_0\neq 0$.

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

it follows that the only zero of Eq (4.5) occurs for $u_h = 0.102520$. In Figure 2, we plot phase space for $u_h = 0.05$ (just before the Hopf bifurcation) using the ode15s solver from Matlab ^[15]. For this set of parameters, we have the following approximate values of the equilibria

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 $_{\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 3. With all the other parameters set to the same values as in the previous case, we have the following equilibria

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 $_{\left(100\right)}{\mathcal{E}}$ is stable, and thus we observe bistability in the system.

It is worth noticing that increasing $u_h$ had a stabilising 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 phenomenon modelled, 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 modelling 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

is said to be weakly persistent if

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

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

for every trajectory with positive initial conditions.

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

Lemma 4.2. Assume that $x^\ast$ is a hyperbolic equilibrium point of the system

with $x\in\mathbb{R}^n$ and $f:\mathbb{R}^n\to\mathbb{R}^n$, where $f$ is continuously differentiable. Assume 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 (3.7) 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 (3.7) is stable.

Theorem 4.2. Let system (3.7) have the following equilibria configuration (as represented schematically in Figure 4):

Then system (3.7) is persistent.

Proof. Since planes $x_0x_1$, $x_0x_2$, and $x_1x_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 4. 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 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 = 0, 1, 2$, is given. We show that no point of the omega limit set has a zero coordinate using proof by contradiction.

First, assume that $_{\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 $_{\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.2, there is a point $x^\ast\neq_{\left(000\right)}{\mathcal{E}}$ in both $\omega\left(\gamma^+\left(\vec x^{\left(0\right)}\right)\right)$ and $W^s\left(_{\left(000\right)}{\mathcal{E}}\right)$, the stable manifold of $_{\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 $_{\left(000\right)}{\mathcal{E}}$ is the $x_1$-axis, and the $x_1$-axis is unbounded. By Lemma 3.1, all orbits of system (3.7) are bounded, and hence the omega limit set of any orbit of system (3.7) is bounded. This contradicts the existence of such an $x^\ast$. Thus $_{\left(000\right)}{\mathcal{E}}\notin\omega\left(\gamma^+\left(\vec x^{\left(0\right)}\right)\right)$.

Next, assume that $_{\left(001\right)}{\mathcal{E}}\in\omega\left(\gamma^+\left(\vec x^{\left(0\right)}\right)\right)$. Since $_{\left(001\right)}{\mathcal{E}}$ is a saddle point with two-dimensional stable manifold restricted to the $x_0x_1$-plane, $\left\{_{\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, by Lemma 4.2, there is a point $x^\ast\in \omega\left(\gamma^+\left(\vec x^{\left(0\right)}\right)\right) \cap W^s\left(_{\left(001\right)}{\mathcal{E}}\right)\setminus\left\{_{\left(001\right)}{\mathcal{E}}\right\}$. Since the stable manifold $W^s\left(_{\left(001\right)}{\mathcal{E}}\right)$ lies entirely in the $x_1x_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 $_{\left(001\right)}{\mathcal{E}}\notin\omega\left(\gamma^+\left(\vec x^{\left(0\right)}\right)\right)$.

Now, assume that $_{\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(_{\left(100\right)}{\mathcal{E}}\right)\setminus\left\{_{\left(100\right)}{\mathcal{E}}\right\}$. This time the stable manifold $W^s\left(_{\left(100\right)}{\mathcal{E}}\right)$ is two-dimensional and lies entirely in the $x_0x_2$-plane. By the analysis in Subsection 4.1, the entire orbit through $x^\ast \;\Big ($which belongs to $\omega\left(\gamma^+\left(\vec x^{\left(0\right)}\right)\right) \Big)$ becomes unbounded in backward time or its closure contains $_{\left(001\right)}{\mathcal{E}}$. This contradiction proves that $_{\left(100\right)}{\mathcal{E}}\notin\omega\left(\gamma^+\left(\vec x^{\left(0\right)}\right)\right)$.

Now, assume that $_{\left(110\right)}{\mathcal{E}}\in\omega\left(\gamma^+\left(\vec x^{\left(0\right)}\right)\right)$. Again, $\left\{_{\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(_{\left(110\right)}{\mathcal{E}}\right)\setminus\left\{_{\left(110\right)}{\mathcal{E}}\right\}$. This point lies in the $x_0x_1$-plane, since $W^s\left(_{\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_0x_1$ face, and since $\left\{_{\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 $_{\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 = 0, 1, 2$, and assume 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_0x_1$, $x_1x_2$, or $x_0x_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

i.e., that system (3.7) is persistent.

An example satisfying the assumptions of Theorem 4.2 occurs for

Persistence can also be observed with the addition of phenol, i.e., with $u_g > 0$.

Theorem 4.3. Let system (3.7) have the following equilibria configuration (as represented schematically in Figure 5):

Then system (3.7) is persistent.

Proof. The idea behind this proof is very similar to the method presented in the proof of Theorem 4.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 system (3.7) with 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 = 0, 1, 2$. Since the stable and unstable manifolds of the all of the equilibria, except $_{\left(001\right)}{\mathcal{E}}$ and $_{\left(011\right)}{\mathcal{E}}$, have the same configuration as in the hypothesis of Theorem 4.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 $_{\left(001\right)}{\mathcal{E}}$ and $_{\left(011\right)}{\mathcal{E}}$.

Assume that $_{\left(001\right)}{\mathcal{E}}\in\omega\left(\gamma^+\left(\vec x^{\left(0\right)}\right)\right)$. Since $_{\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\{_{\left(001\right)}{\mathcal{E}}\right\}\neq\emptyset$. Hence, by Lemma 4.2, there is a point $x^\ast\in\omega\left(\gamma^+\left(\vec x^{\left(0\right)}\right)\right)\cap W^s\left(_{\left(001\right)}{\mathcal{E}}\right)\setminus\left\{_{\left(001\right)}{\mathcal{E}}\right\}$. The entire orbit through $x^\ast$ also belongs to $\omega\left(\gamma^+\left(\vec x^{\left(0\right)}\right)\right)$ and either becomes unbounded in backward time or converges to $_{\left(000\right)}{\mathcal{E}}$. Since all orbits of system (3.7) are bounded and $_{\left(000\right)}{\mathcal{E}}\notin\omega\left(\gamma^+\left(\vec x^{\left(0\right)}\right)\right)$, we obtain a contradiction. Hence, $_{\left(001\right)}{\mathcal{E}}\notin\omega\left(\gamma^+\left(\vec x^{\left(0\right)}\right)\right)$.

Next, assume that $_{\left(011\right)}{\mathcal{E}}\in\omega\left(\gamma^+\left(\vec x^{\left(0\right)}\right)\right)$. Since $_{\left(011\right)}{\mathcal{E}}$ is a saddle point with two-dimensional stable manifold restricted to the $x_1x_2$-plane (it is repelling into the interior), we have $\omega\left(\gamma^+\left(\vec x^{\left(0\right)}\right)\right)\setminus\left\{_{\left(011\right)}{\mathcal{E}}\right\}\neq\emptyset$. By Lemma 4.2, there exists a point $x^\ast\in\omega\left(\gamma^+\left(\vec x^{\left(0\right)}\right)\right)\cap W^s\left(_{\left(011\right)}{\mathcal{E}}\right) \setminus\left\{_{\left(011\right)}{\mathcal{E}}\right\}$. The entire orbit through $x^\ast$ belongs to $\omega\left(\gamma^+\left(\vec x^{\left(0\right)}\right)\right)$ and either becomes unbounded in backward time or converges to $_{\left(000\right)}{\mathcal{E}}$ or $_{\left(001\right)}{\mathcal{E}}$. We have previously shown in Subsection 4.1 that there are no periodic orbits in the $x_1x_2$ face. Since we have already proven that $_{\left(000\right)}{\mathcal{E}}\notin\omega\left(\gamma^+\left(\vec x^{\left(0\right)}\right)\right)$ and $_{\left(001\right)}{\mathcal{E}}\notin\omega\left(\gamma^+\left(\vec x^{\left(0\right)}\right)\right)$, we obtain a contradiction. This proves that $_{\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 = 0, 1, 2$, and assume 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 the $x_0x_1$, $x_1x_2$, or $x_0x_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)$. But this possibility has been eliminated in the previous part of the proof.

An example satisfying the assumptions of Theorem 4.3 occurs for

Remark 1. Interestingly enough, in many cases of models describing biological phenomena, persistence already implies uniform persistence. The rigorous results were obtained in ^[17]. In our context, the key theorem from ^[17] 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. $\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. $\mathcal{F}$ is weakly persistent,

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

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

These results can be easily modified ^[18] so that we consider the flow $\mathcal{F}$ on $\Omega$ and $\partial\Omega$, which are positively invariant. The following theorem is proven using results of ^[17,18], by checking conditions 1–4 on $\Omega$ and $\partial \Omega$.

Theorem 4.4. Under the hypotheses of Theorem 4.2 or Theorem 4.3, system (3.7) is uniformly persistent.

Proof. In our case, the positively invariant set $\Omega$, on which we analyse system (3.7) is bounded. Hence, condition $1$ is satisfied. Condition $2$ holds by Theorem 4.2 (persistence implies weak persistence). In our context condition 3 is satisfied, because the omega limit set on $\partial \Omega$ consists of a finite number of hyperbolic equilibria. 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 (3.7) in the cases when the hypotheses of Theorem 4.2 or Theorem 4.3 are satisfied.

We finish this subsection by extending the uniform persistence to the original six dimensional system (2.1).

Theorem 4.5. Under the hypotheses of Theorem 4.2 or Theorem 4.3, the six dimensional system (2.1) is uniformly persistent.

Proof. 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) > 0$, $s_i(0) > 0$, $i = 0, 1, 2$, then we necessarily must have $\omega\left({\vec X_0}\right)\in\Omega$. Otherwise, there would exist a point $\vec p\in\mathbb{R}^6_+\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 Lemma 3.2. Moreover, there is only a finite number of equilibria on $\partial\Omega$, with all of their stable manifolds lying entirely in $\partial\Omega$. Thus, if the assumptions of either Theorem 4.2 or Theorem 4.3 hold, by Theorem 4.4 and its proof, the corresponding three-dimensional system is uniformly persistent and the boundary $\partial\Omega$ is acyclic. Hence, we can apply the results of Thieme ^[19] to conclude that the six-dimensional system (3.1) is also uniformly persistent.

4.4.   Bifurcation diagrams

As previously stated in Section 4.2, we now study effects on the qualitative behaviour of system (3.7) numerically. We consider $\alpha$ as a 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

and plot a one-parameter bifurcation diagram in $\alpha$, with $x_0$ on the $y$-axis. All simulations were performed using ^[20].

In Figure 6, we can see that as $\alpha$ decreases, there is a saddle-node bifurcation, resulting in two equilibria $_{\left(110\right)}{\mathcal{E}}^{\left(1\right)}$ and $_{\left(110\right)}{\mathcal{E}}^{\left(2\right)}$ appearing (both unstable). Next, there is a transcritical bifurcation with the $_{\left(110\right)}{\mathcal{E}}^{\left(1\right)}$ equilibrium. This results in the positive equilibrium moving into the interior of the admissible region, $\Omega$. After that, a saddle-node of limit cycles bifurcation occurs that gives birth to stable and unstable periodic orbits. The $_{\left(111\right)}{\mathcal{E}}$ equilibrium (unstable), undergoes a Hopf bifurcation, and as a consequence it becomes asymptotically stable, and the stable periodic orbit disappear. Since these bifurcations occur for a narrow range of $\alpha$, a close-up is presented in Figure 7. Existence of a 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 by the biogas plant. The unstable periodic orbit 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 (3.6), $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 6. The amount of chlorophenol in the system is inversely proportional to the concentration of the phenol degrader. As the dilution rate $\alpha$ decreases, the 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:

and plot a two-parameter bifurcation diagram of system (3.7), choosing $\alpha$ and $u_f$ as the bifurcation parameters. Each region of the diagram is labeled and the corresponding dynamics are represented schematically in the surrounding figures in phase space. Black dashed curves correspond to saddle-node of equilibria bifurcations (LP), black solid curves represent saddle-node of limit cycles bifurcations (SNLC), black dotted curves depict Hopf bifurcations (HB), and grey solid curves represent transcritical bifurcations (BP). We also depict the predicted heteroclinic bifurcation by a grey dashed curve which lies very close to the BP curve.

We can see in Figures 8–10, that varying two parameters at the same time can demonstrate much more complicated dynamics than in the case of one-parameter bifurcation diagrams. There is a generalised 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, moving 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_0x_1$ face, then passes into the interior, and then goes back to the boundary in the $x_0x_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 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 $_{\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 ^[9] by including multiple substrate inflow in 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. We also observed that external addition of substrates (phenol and hydrogen) 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 the chlorophenol degrader population surviving.

We have also confirmed the findings of the previous analysis in ^[9], 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 stabilising effect on the entire system. From a biological engineering point of view, in 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 aspect 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 and a Bautin (generalised 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.

Anecdotally, there has been some resistance to using reduced order models, of the type described here. Some researchers consider them too removed from the systems they represent to be of practical use. Without experimental results to compare against model predictions, they remain mostly theoretical constructs limited by the underlying assumptions and simplifications imposed. However, we can look to emerging disciplines such as synthetic biology to help bridge the theoretical and the applied ^[21]. 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 ^[22]. For example, we have shown that, in a region of unstable dynamics for a given set of operational parameters and initial conditions, addition of extraneous hydrogen results in an asymptotically stable equilibrium by ensuring survival of the two hydrogen consumers. Hydrogen is a key substrate in methanogenic systems, but process stabilisation by its addition has only recently been suggested ^[6,7]. We believe this work provides a basis from which experimental studies describing microbial food-webs with complex, high-order interactions can be conducted with a posteriori knowledge of the anticipated behaviour.

Acknowledgments

Matthew J. Wade 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). Gail S.K. Wolkowicz acknowledges support from the Natural Sciences and Engineering Council of Canada (NSERC) Discovery Grant with accelerator supplement.

Conflict of interest

The Authors declare no conflict of interest in this work.