Generalised approach to modelling a three-tiered microbial food-web

Ecological modelling of increasingly more complex microbial populations is necessary to reflect the highly functional and diverse behaviour inherent to many systems found in reality. Anaerobic digestion is one such process that has benefitted from the application of mathematical analysis not only for characterising the biological dynamics, but also to investigate emergent behaviour not apparent by simulation alone. Nevertheless, the standard modelling approach has been to describe biological systems using sets of differential equations whose kinetics are generally described by some empirically derived function of growth. The drawbacks of this are two-fold; the growth functions are derived from empirical studies that may not be representative of the system to be modelled and whose parameters may not have a mechanistic meaning, and mathematical analysis is restricted by a conformity to an assumption of the dynamics. Here, we attempt to address these challenges by investigating a generalised form of a three-tier chlorophenol mineralising food-web previously only analysed numerically. We examine the existence and stability of the identified steady-states and find that, without a decay term, the system may be characterised analytically. However, it is necessary to perform numerical analysis for the case when maintenance is included, but in both cases we verify the discovery of two important phenomena; i) the washout steady-state is always stable, and ii) the two other steady-states can be unstable according to the initial conditions and operating parameters.


Introduction
The mathematical modelling of engineered biological systems has entered a new era in recent years with the expansion and standardisation of existing models aimed at collating disparate components of these processes and provide scientists, engineers and practitioners with the tools to better predict, control and optimise them. These forms of mechanistic models emerged initially with the Activated Sludge Models [1,2] for wastewater treatment processes, followed by the Anaerobic Digestion Model No. 1 (ADM1) [3] a few years later. The development of ADM1 was enabled largely due to the possibilities for better identification and characterisation of functional groups responsible for the discrete degradation steps operating in series within anaerobic digesters. It describes a set of fairly complex stoichiometric and kinetic functions representing the standard anaerobic process, remaining the scientific benchmark to the present day, despite an general understanding of its limitations in describing all necessary biochemical transformations. Indeed, there has been a growing argument that the model should take advantage of improved empirical understanding and extension of biochemical processes included in its structure, to acquire a better trade-off between model realism and complexity [4].
It has previously been shown that the simplification or reduction in model complexity can preserve biolog-ical meaning whilst reducing the computational effort required to find mathematical solutions of the model equations [5,6,7]. Whilst simpler models are approximations of real systems, it can be beneficial to consider a reduced model to better understand biological phenomena of sub-processes without the need to consider extraneous system parameters and variables, which tend to make mathematical analysis intractable and cumbersome. Nevertheless, even with gross simplification of a biological system based on a set of ordinary differential equations (ODEs) of relatively low dimensionality, analytical techniques are unable to provide general solutions for the system and numerical methods must suffice.
As an example of this, for anaerobic digestion, a previous study investigated the effect of maintenance on the stability of a two-tiered 'food-chain' comprising two species and two substrates [8]. Although the authors were not able to determine the general conditions under which this four dimensional syntrophic consortium was stable, further work has shown that a model with generality can be used to answer the question posed, determining that the two-tiered food-chain is always stable when maintenance is included [9].
In a more recent example, the model described by [8] was extended by the addition of a third organism and substrate to create a three-tiered 'food-web' [10]. In this model, the stability of some steady-states could be determined analytically, but due to the complexity of the Jacobian matrix for certain steady-states, local solutions were necessary using numerical analysis, when considering the full system behaviour. Although the results were important in revealing emergent properties of this extended model, the motivation of this work is to determine whether the approach carried out in [9], can be applied to the three-tiered model from [10], to provide some general properties of that system. The paper is organised as follows. In Section 2, we present a description of the model to be investigated, before providing an alternative reduction of its structure than that given by [10], in Section 3. With Section 4 we demonstrate the existence of the three steadystates and define four interesting cases for specific parameter values that are investigated using the analytical solutions, whilst also indicating the regions of existence of the steady-states for the operating parameter values (dilution rate and substrate input concentration). In Section 5 we perform local stability analysis of the steady-states without maintenance, and in Section 6, perform a comprehensive numerical stability analysis of the four cases for both the model with and without a decay constant. We show that our approach leads to the discovery of five operating regions, in which one leads to the possibility of instability of the positive steady state, where all three organisms exist, a fact that has not be reported by [10]. Indeed, we prove that a stable limit-cycle can occur in this region. Finally, in Section 7, we make comment on the role of the kinetic parameters used in the four example cases, in maintaining stability, which points to the importance of the relative aptitude of the two hydrogen consumers in sustaining a viable chlorophenol mineralising community. In the Appendix we describe the numerical method used in Section 6 and we give the proofs of the results.

The model
The model developed in [10] has six components, three substrate (chlorophenol, phenol and hydrogen) and three biomass (chlorophenol, phenol and hydrogen degraders) variables. The substrate and biomass concentrations evolve according to the six-dimensional dynamical of ODEs where S ch and X ch are the chlorophenol substrate and biomass concentrations, S ph and X ph those for phenol and S H2 and X H2 those for hydrogen; Y ch , Y ph and Y H2 are the yield coefficients, 224/208 (1 − Y ch ) represents the part of chlorophenol degraded to phenol, and 32/224 (1 − Y ph ) represents the part of phenol that is transformed to hydrogen. Growth functions take Monod form with hydrogen inhibition acting on the phenol degrader and represented in f 1 (see Eq. 7) as a product inhibition term.
Here, apart from the four operating (or control) parameters, which are the inflowing concentrations S ch,in , S ph,in , S H2,in and the dilution rate D, that can vary, all others have biological meaning and are fixed depending on the organisms and substrate considered. We use the following simplified notations in (Eqs. 1-6) With these notations Eqs. 1-6 can be written as follows In [10], this model is reduced to a dimensionless form that significantly reduces the number of parameters describing the dynamics. In this paper we do not assume that the growth functions f 0 , f 1 and f 2 have the specific analytical expression (Eq. 7). We will only assume that the growth functions satisfy properties that are listed in Appendix C. Therefore, we cannot benefit from the dimensionless rescaling used by [10], because this rescaling uses some kinetics parameters of the specific growth functions (Eq. 7), while we work with general unspecified growth functions. In Section 3 we consider another rescaling that does not use the kinetics parameters. Furthermore, we restrict our analysis to the case where we only have one substrate addition to the system, such that: S in 0 > 0, S in 1 = 0, and S in 2 = 0.

Model reduction
To ease the mathematical analysis, we can rescale the system (Eqs. 8-13) using the following change of variables adapted from [9]: We obtain the following system where the inflowing concentration is the growth functions are and The benefit of our rescaling is that it permits to fix in Eqs. 14-19 all yield coefficients to one except that denoted by ω and defined by (Eq. 22), and to discuss the existence and stability with respect to this sole parameter.
Using Eq. 21 and the growth functions (Eq. 7), we obtain the model (Eqs. 14-19) with the following Monod-type growth functions For the numerical simulations we will use the nominal values in Table 1 given in [10]. 3.5×10 −6 kgCOD/m 3 Table 1 Nominal parameter values.

Existence of steady-states
A steady-state of Eqs. 14-19 is obtained by setting the right-hand sides equal to zero: A steady-state exists (or is said to be 'meaningful') if, and only if, all its components are non-negative.
From the previous Lemma we deduce that besides the steady-state SS1, the system can have at most two other steady-states.

2
> 0 for all s 2 ∈ (s 0 2 (D), s 1 2 (D)), so that the function s 2 → ψ(s 2 , D) is convex and, thus, it has a unique minimum s 2 (D), see Fig. 1 Let ω < 1. We define the function as shown in Fig. 1 (b). The minimum s 2 (D) is a solution of an algebraic equation of degree 4 in s 2 . Although mathematical software, such as Maple, cannot give its solutions explicitly with respect to the parameters, s 2 (D) could be obtained analytically since algebraic equations of degree 4 can theoretically be solved by quadratures. We do not try to obtain such an explicit formula. However, if the biological parameters are fixed, the function s 2 (D) and, hence, . We denote by the set on which F 1 (D) is defined.
Let ω < 1. We define the functions Since M 2 and ψ are given explicitly by Eq. 33 and Eq. 34, respectively, the functions F 2 (D) and F 3 (D) are given explicitly with respect to the biological parameters in Eq. 23. The functions F 2 (D) and , that is to say, for D such that We denote by the subset of I 1 on which F 2 (D) and F 3 (D) are defined. For all for D ∈ I 2 , F 1 (D) ≤ F 2 (D). The equality F 1 (D) = F 2 (D) holds if, and only if, M 2 (D + a 2 ) = s 2 (D) that is, dψ ds2 (M 2 (D + a 2 )) = 0. Therefore, holds if, and only if, F 3 (D) = 0. We define Since D → s 0 2 (D) is increasing and D → s 1 2 (D) is decreasing, and assuming s 0 For the domain of definition I 2 of F 2 (D), several cases can be distinguished. I 2 is an interval  Table 2 and 3.
We can state now the necessary and sufficient conditions of existence of SS2 and SS3.

Lemma 2
If ω ≥ 1 then SS2 does not exist. If ω < 1 then SS2 exists if, and only if, s in 0 ≥ F 1 (D). Therefore, a necessary condition for the existence of SS2 is that D ∈ I 1 , where I 1 is defined by Eq. 37. If s in 0 ≥ F 1 (D) then each solution s 2 of equation If ω ≥ 1 then SS3 does not exist. If ω < 1 then SS3 exists if, and only if, s in 0 > F 2 (D). Therefore, a necessary condition of existence of SS3 is that D ∈ I 2 , where I 2 is defined by Eq. 40. If s in 0 > F 2 (D) then the steady-state SS3 = (x 0 , x 1 , x 2 , s 0 , s 1 , s 2 ) is given by and then Eq. 41 has exactly two solutions denoted by s 2 and s 2 and such that, see Fig.  1 To these solutions, s 2 and s 2 , correspond two steady-states of SS2, which are denoted by SS2 and SS2 . These steady-states coalesce when s in 0 = F 0 (D). Since F 1 (D) ≤ F 2 (D), the condition s in 0 > F 2 (D) for the existence of the positive steady-state SS3 implies that the condition s in 0 > F 2 (D) for the existence of the two steady-states SS2 and SS2 is satisfied. Therefore, if SS3 exists then SS2 and SS2 exist and are distinct. If Remark 2 Using Eq. 20, the conditions s in 0 > F 1 (D) and s in 0 > F 2 (D) of existence of the steady-state SS2 and SS3 respectively are equivalent to the conditions respectively, expressed with respect to the inflowing concentration S ch,in .
Our aim now is to describe the operating diagram : The operating diagram shows how the system behaves when we vary the two control parameters S ch,in and D in Eqs. 1-6. According to Remark 2, the curve Γ 1 of equation is the border to which SS2 exists, and the curve Γ 2 of equation is the border to which SS3 exists, see Fig. 3. If we want to plot the operating diagram we must fix the values of the biological parameters. In the remainder of the Section we plot the operating diagrams corresponding to cases (a)-(d) depicted in Fig. 2.

Operating diagram: case (a)
This case corresponds to the parameter values used by [10]. We have seen in Table 2 that the curves Γ 1 and Γ 2 are defined for D < D 1 and D < D 2 , respectively and that they are tangent for D = D 3 , where D 1 = 0.432, D 2 = 0.373 and D 3 = 0.058. Therefore, they separate the operating plane (S ch,in , D) into four regions, as shown in Fig. 3(i), labelled J 1 , J 2 and J 3 and J 4 . The results are summarised in Table 4, which shows the existence of the steady-states SS1, SS2 and SS3 in the regions of the operating diagram in Fig. 3(i).

Region
Steady states Table 4 Existence of steady-states in the regions of the operating diagrams of Fig. 3(i) and Fig. 6(i).

Operating diagram: case (b)
This case corresponds to the parameter values used by [10], except that K S,H2,c is changed from 1.0 × 10 −6 to 4.0 × 10 −6 . We have seen in Table 2 that the curves Γ 1 and Γ 2 are defined for D < D 1 and D < D 2 , respectively and F 1 (D) < F 2 (D) for all D < D 2 , where D 1 = 0.329 and D 2 = 0.236. Therefore, they separate the operating plane (S ch,in , D) in three regions, as shown in Fig. 4(i), labelled J 1 , J 3 and J 4 .
The results are summarised in Table 5, which shows the existence of the steady-states SS1, SS2 and SS3 in the regions of the operating diagram in Fig. 4(i). Note that the region J 2 has disappeared.
The curve Γ 1 for case (c). (i) : regions of steady-state existence, with maintenance. (ii) : regions of steady-state existence. without maintenance and their stability. On the right, a magnification for 0 < D < 0.1.

Operating diagram: case (c)
This case corresponds to the parameter values used by [10], except that K S,H2,c is changed from 1.0 × 10 −6 to 7.0 × 10 −6 . We have seen in Table 2 that the curve Γ 1 is defined for D < D 1 = 0.287 and that I 2 is empty so that SS3 does not exist. Therefore, Γ 1 separates the operating plane (S ch,in , D) in two regions, as shown in Fig. 5(i), labelled J 1 and J 4 . The results are summarised in Table 6, which shows the existence of the steady-states SS1 and SS2 in the regions of the operating diagram in Fig. 5(i). Note that the region J 3 of existence of SS3 has disappeared.

Operating diagram: case (d)
We end this discussion on the role of kinetic parameters by the presentation of this case, which presents a new behaviour that did not occur in the preceding cases: there exists D 2min such that for D < D 2min the system cannot have a positive steady-state SS3. This case corresponds to the parameter values used by [10], except that three of them are changed as indicated in Table 3. This table shows that the curves Γ 1 and Γ 2 are defined for D < D 1 and D 2min < D < D 2max and  Table 4, which shows the existence of the steady-states SS1, SS2 and SS3 in the regions of the operating diagram in Fig. 6(i).

Stability of steady-states
We know that SS1 is always stable. The analytical study of the stability of SS2 and SS3 is very difficult because the conditions for Routh-Hurwitz in the 6-dimensional case are intractable. For this reason we will consider in Section 5 the question of the stability only in the case without maintenance, since the system reduces to a 3dimensional. The general case will be considered only numerically in Section 6.

Local stability without maintenance
When maintenance is not considered in the model, the steady-states SS1, SS2 and SS3 are given by The instability of SS3 occurs in particular when s 2 is sufficiently close to s 2 , that is to say SS3 is sufficiently close to SS2 .
When D is such that F 3 (D) < 0, the determination of the boundary between the regions of stability and instability of SS3 needs to examine the Routh-Hurwitz condition of stability for SS3. For this purpose we define the following functions. Let SS3 = (x 0 , x 1 , x 2 , s 0 , s 1 , s 2 ) be a steady-state. Let For D ∈ I 3 and s in 0 > F 2 (D), we define where f 2 = Ix 2 +(G+H)x 1 +(E+ωF )x 0 . Notice that to compute F 4 D, s in 0 , we must replace x 0 , x 1 , x 2 , s 0 , s 1 and s 2 by their values at SS3, given by (48). Hence, this function depends on the operating parameters D and Existence Stability SS1 Always exists Always stable SS2 Always unstable SS2   Table 7.

Operating diagram: case (a)
This case corresponds to the parameter values used by [10] but without maintenance. We see from Table  2 that the curves Γ 1 and Γ 2 of the operating diagram, given by Eq. 45 and Eq. 46, respectively, are defined now for D < D 1 = 0.452 and D < D 2 = 0.393, respectively and that they are tangent for D = D 3 = 0.078. Beside these curves, we plot also on the operating diagram of Fig. 3(ii), the curve Γ 3 of equation According to Prop. 3, this curve is defined for D < D 3 = 0.078 and it separates the region of existence of SS3 into two subregions labelled J 3 and J 5 , such that SS3 is stable in J 3 and unstable in J 5 . The other regions J 1 , J 2 and J 4 are defined as in the previous section. The operating diagram is shown Fig. 3(ii). It looks very similar to Fig. 3(i), except near the origin, as it is indicated in the magnification for 0 < D < D 3 = 0.078. From Table 7, we deduce the following result Proposition 4 Table 8 shows the existence and stability of the steady-states SS1, SS2 and SS3 in the regions of the operating diagram in Fig. 3(ii). Table 8 Existence and stability of steady-states in the regions of the operating diagrams of Fig. 3(ii) and Fig. 6(ii). Table 9 Existence and stability of steady-states in the regions of the operating diagram of Fig. 4(ii). Table 10 Existence and stability of steady-states in the regions of the operating diagram of Fig. 5(ii).

Operating diagram: case (b)
We see from Table 2 that the curves Γ 1 and Γ 2 are defined now for D < D 1 = 0.349 and D < D 2 = 0.256, respectively and that F 1 (D) < F 2 (D) for all D. Beside these curves, we plot also on the operating diagram of Fig. 4(ii), the curve Γ 3 of equation (Eq. 50) which separates the region of existence of SS3 into two subregions labelled J 3 and J 5 , such that SS3 is stable in J 3 and unstable in J 5 . Therefore, the curves Γ 1 , Γ 2 and Γ 3 separate the operating plane (S ch,in , D) into four regions, as shown in Fig. 4(ii), labelled J 1 , J 3 , J 4 and J 5 .
Proposition 5 Table 9 shows the existence and stability of the steady-states SS1, SS2 and SS3 in the regions of the operating diagram in Fig. 4(ii)

Operating diagram: case (c)
We see from Table 2 that Γ 1 is defined for D < D 1 = 0.303 and that I 2 is empty so that SS3 does not exist. Therefore, Γ 1 separates the operating plane (S ch,in , D) into two regions, as shown in Fig. 5(ii), labelled J 1 and J 4 . Table 10 shows the existence and stability of the steady-states SS1, SS2 and SS3 in the regions of the operating diagram in Fig. 5(ii).

Operating diagram: case (d)
We see in Table 3 that the curves Γ 1 and Γ 2 are defined for D < D 1 and D 2min < D < D 2max and that they are tangent for D = D 3 , where D 1 = 0.258 and D 2min = 0.121, D 2max = 0.218 and D 3 = 0.181. Beside these curves, we plot also on the operating diagram of Fig. 6(ii), the curve Γ 3 defined by Eq. 50, which separates the region of existence of SS3 into two subregions labelled J 3 and J 5 , such that SS3 is stable in J 3 and unstable in J 5 . Therefore, the curves Γ 1 , Γ 2 and Γ 3 separate the operating plane (S ch,in , D) into five regions, as shown in Fig. 6(ii), labelled J 1 , J 2 , J 3 , J 4 and J 5 .
Proposition 7 Table 8 shows the existence and stability of the steady-states SS1, SS2 and SS3 in the regions of the operating diagram in Fig. 6(ii).

Numerical analysis to confirm and extend the analytical results
The aim of this section is to study numerically (the method is explained in Appendix A) the existence and stability of the steady-states SS2 and SS3. We obtain numerically the operating diagrams that were described in Sections 4 and 5. The results in this section confirm the results on existence of the steady-states obtained in Section 4 in the case with or without maintenance and the results of stability obtained in Section 5 in the case without maintenance. These results permit also to elucidate the problem of the local stability of SS2 and SS3, which was left open in Section 4.5.

Operating diagram: case (a)
We endeavoured to find numerically the operating conditions under which SS3 is unstable, previously unreported by [10]. Given that we have determined analytically in Proposition 2 that when SS3 is close to SS2 it becomes unstable, we performed numerical simulations with the parameters defined in Table 1 over an operating region similar to that shown in Fig. 2 from [10] whilst also satisfying our conditions. In Fig. 7 we show the case when maintenance is excluded. When magnified, we observe more clearly that region J 5 does exist for the conditions described above, and also note that the region J 4 occurs in a small area between J 1 and J 5 , which corresponds to the results shown in Fig. 3(ii), and is in agreement with Proposition 4. In Fig. 8 we confirm that region J 5 does exist for the conditions described above, when maintenance is included, but could not be determined analytically, the curve Γ 3 is absent in Fig. 3(i). Furthermore, we demonstrate that a Hopf bifurcation occurs along the boundary of F 3 (D) for values of D < D 3 by selecting values of S ch,in (indicated by (α) − (δ) in Fig. 8) at a fixed dilution rate of 0.01 d −1 , and running dynamic simulations for 10000 d.
The three-dimensional phase plots, with the axes representing biomass concentrations, are shown in Fig. 9, and show that as S ch,in approaches J 3 from J 5 , emergent periodic orbits are shown to diminish to a stable limit cycle at the boundary (see Appendix B for proof). Subsequently, increasing S ch,in to J 3 results in the orbit reducing to a fixed point equilibrium at SS3. D S ch,in Fig. 8 Numerical analysis for the existence and stability of steady-states for case (a), with maintenance. This is a magnification for 0 < D < 0.1, showing the presence and extent of region J 5 undetectable by the analytical method. The coordinates labelled (α) − (δ) are subsequently used to simulate the system dynamics, as shown in the proceeding Fig. 9. and chlorophenol input, S ch,in (kgCOD/m 3 ) of α) 0.01 -the system converges to SS1, β) 0.097 -the system enters a periodic orbit of increasing amplitude, ultimately converging to SS1, γ) 0.10052 -the system is close to a stable limit cycle, δ) 0.16 -the system undergoes damped oscillations and converges to SS3.

Operating diagram: case (b)
Whilst the numerical parameters chosen for this work are taken from the original study [10], there somewhat arbitrary nature leaves room to explore the impact of the parameters on the existence and stability of the steady-states. Case (b), discussed in Sections 4.2 and 5.2, involves a small increase to the half-saturation constant (or inverse of substrate affinity), K S,H2,c , of the chlorophenol degrader on hydrogen. Following the same approach as with the preceding case, we confirm in Fig. 10(i) the Proposition 5 in the scenario without maintenance. Furthermore, the extension of this proposition with maintenance included, corresponding to the existence and stability of all three steady-states given in Table 9, is show in Fig 10(ii). It shows the region J 5 that cannot be obtained analytically (cf. Fig. 4(i)). In both cases, region J 2 has disappeared, as observed analytically. Additionally, the ideal condition J 3 , where all organisms are present and stable, diminishes.

Operating diagram: case (c)
Here, K S,H2,c , was further increased and confirm the Proposition 6, where the function SS3 never exist and SS2 never stable for the case without maintenance. The extension of this proposition to the case with mainte-nance, shown in Table 10, produce similar results as shown in the comparison of Figs. 11(i) and (ii).

Operating diagram: case (d)
With the final investigated scenario, where k m,H2 < k m,ch and K S,H2 < K S,H2,c , we observe once again the presence of all operating regions, J 1 − J 5 , without and with maintenance, as shown in Fig. 12. It can be seen that regions J 4 and J 5 increase at low dilution rates across a much larger range of S ch,in than in the default case (a), and the desired condition (stable SS3) is restricted to a much narrower set of D.
As with the previous cases, the numerical analysis for case (d) confirms the Proposition 7 without maintenance and its extension to the case with maintenance, indicated in Table 8.

The role of kinetic parameters
Finally, we give brief consideration to the characterisation of the four cases discussed in the preceding sections. The main difference between cases (a) or (b) and cases (c) or (d) is that, for small values of D, the coexistence steady-state SS3 can exist for cases (a) and (b), but cannot exist for cases (c) or (d). The cases (a) or (b) occur if and only if s 0 2 (0) < M 2 (0) holds or s 0 2 (0) = M 2 (0) and ds 0 2 dD (0) < dM2 dD (0) hold, that is to say Notice that it is easy to make the difference between case (c) and case (d): the first occurs when M 2 (D 1 ) < s 0 2 (D 1 ) and the second when M 2 (D 1 ) > s 0 2 (D 1 ). Since D 1 is the positive solution of the algebraic quadratic equation s 0 2 (D) = s 1 2 (D), it is possible to have an expression for D 1 with respect to the biological param-eters. However, this is a complicated expression involving many parameters and the preceding conditions M 2 (D 1 ) < s 0 2 (D 1 ) or M 2 (D 1 ) > s 0 2 (D 1 ) have no biological interpretation. We simply remark here that the function s 0 2 (D) has a vertical asymptote for D = m 0 −a 0 and the function M 2 (D) has a vertical asymptote for D = m 2 − a 2 . Therefore, if m 0 − a 0 < m 2 − a 2 then case (c) occurs, so that a necessary (but not sufficient) condition for case (d) to occur is m 0 − a 0 > m 2 − a 2 . If m 2 is sufficiently small then case (d) can occur.
The observations from the numerical analysis suggest that the role of the chlorophenol degrader as a secondary hydrogen scavenger is critical in maintaining full chlorophenol mineralisation and system stability, particularly at higher dilution rates, as shown by comparing cases (c) and (d) . More significantly, the results coupled with the parameter relationships shown in Eqs. 51-54, highlight the necessary conditions under which the ideal case (SS3 stable) is achieved and, in general, this is a coupling of the two key parameters describing the half-saturation constant and maximum specific growth rates between the two hydrogen competitors.

Conclusions
In this work we have generalised a simplified mechanistic model describing the anaerobic mineralisation of chlorophenol in a two-step food-web. We are able to show complete analytical solutions describing the existence and stability of the steady-states in the case that maintenance is excluded from the system, whilst with a decay term present, purely analytical determination of stability is not possible.
We confirm the findings of previous numerical analysis by [10] that with chlorophenol as the sole input substrate, three steady-states are possible. However, the analysis goes further and we determine that under certain operating conditions, two of these steady-states (SS2 and SS3) can become stable, whilst SS1 always exists and is always stable. Furthermore, without maintenance we can explicitly determine the stability of the system, and form analytical expressions of the boundaries between the different stability regions.
As the boundary of J 3 is not open to analytical determination in the case with maintenance, we determined numerically (substituting the general growth function with the classical Monod-type growth kinetics) the existence and stability of the system over a range of practical operating conditions (dilution rate and chlorophenol input). For comparison and confirmation, we also performed this for the case without maintenance and found the same regions in both cases, with variations only in their shape and extent. For example, whilst the boundary between J 1 and J 4 terminates at the origin without maintenance, with maintenance it is located at F 1 (0)/Y 3 Y 4 ≈ 0.0195. More interestingly, the addition of a decay term results in an extension of the SS3 unstable steady-state, reducing the potential for successful chlorophenol demineralisation at relatively low dilution rates and substrate input concentrations. Additionally, we show that at the boundary between J 3 and J 5 , a Hopf bifurcation occurs and a limit cycle in SS3 emerges.
Finally, we gave an example of how the model could be used to probe the system to answer specific questions regarding model parameterisation. Here we have indicated that a switch in dominance between two organisms competing for hydrogen results in the system becoming unstable and a loss in viability. This is perhaps intuitive to microbiologists, but here it has been proven using mathematical analysis, and could be used to determine critical limits of the theoretical parameter values in shifting between a stable and unstable system. Whilst parameters are not arbitrary in real organisms, the potential for microbial engineering or synthetic biology to manipulate the properties of organisms makes this observation all the more pertinent.

A Numerical methods
We consider sets of operating parameters (D and S ch,in ) for each of the three steady-states, and using Matlab, the complex polynomials for each steady-state can be solved by substitution of parameter values (see Table 1) into the explicit solution. By investigating the signs of the solutions and the eigenvalues, respectively, we determine which steady-states are meaningful and stable. By exploring a localised region of suitable operating parameters, we then generate a phase plot showing where each steady-state is stable, bistable or unstable.

B Proof for Hopf Bifurcation
In Section 6, we show the operating diagrams with the parameters given in Table 1, and determine numerically that as the parameter S ch,in increases at a fixed dilution rate (D = 0.01 d −1 ), the system bifurcates through several stability domains. We claim that as we cross the boundary between regions J 5 and J 3 , we observe a Hopf bifurcation, and, in J 5 , close to the boundary with J 3 , a limit cycle appears. In order to test this numerically, we checked the real parts of the six eigenvalues at each point along the transect shown in Fig 8 (10000 points in total), and plotted their values.
We use Eq. 36, Eq. 38 and Eq. 39 to define F 1 (D), F 2 (D) and F 3 (D) in the general case: The function F 1 (D) is defined for The function F 2 (D) and F 3 (D) are defined for For all for D ∈ I 2 , F 1 (D) ≤ F 2 (D). The equality F 1 (D) = F 2 (D) holds if, and only if, M 2 (D + a 2 ) = s 2 (D) that is, dψ ds2 (M2(D + a 2 )) = 0, that is if, and only if, F 3 (D) = 0. As it will be shown in Appendix D, the Lemmas 2 and 3, stated in Section 5 in the particular case of the Monod type growth functions (Eq. 23), are true in the general case of growth functions satisfying assumptions H1-H8.

D Proofs
In this Section we give the proofs of the results. In these proofs, we do not assume that the growth function are of Monod type (Eq. 23). We only assume that the growth functions satisfy H1-H8.
Let (x 0 , x 1 , x 2 , s in 0 , 0, 0) be a steady-state. The Jacobian matrix of Eqs. 67-72 has the block triangular form are evaluated at the steady-state. Since J is a block triangular matrix, its eigenvalues are −D (with multiplicity 3) together with the eigenvalues of the 3 × 3 upper-left matrix J 1 . Note that we have used the opposite sign for the partial derivative H = −∂µ 1 /∂s 2 , so that all constants involved in the computation become positive, which will simplify the analysis of the characteristic polynomial of J 1 .
Proof [Proposition 1] Evaluated at SS2, the matrix J 1 is Since J 1 is a block triangular matrix, its eigenvalues are simply µ 2 −D, together with the eigenvalues of the 2×2 upper-left matrix. Note that the trace of this 2 × 2 matrix is negative. Hence, its eigenvalues are of negative real part if, and only if, its determinant is positive, that is if, and only if, Therefore, the condition of stability, (Eq. 73), is equivalent to dψ ds2 > 0. Hence, we have proved that SS2 is stable if, and only if, µ 2 (s 2 ) < D and dψ ds2 > 0.
Proof [Proposition 3] Since we always have f 0 > 0 and f 2 > 0, from the previous proof it follows that SS3 is stable if, and only if, f 1 f 2 − f 0 > 0. Indeed, this condition implies that we have also f 1 > f 0 /f 2 > 0. Using Eq. 49 and Eq. 79, we see that Therefore, the condition f 1 f 2 − f 0 > 0 is equivalent to

D.3 Operating diagrams
Proof [Proposition 4] We know that SS1 always exist and is stable. We know that SS2 is unstable if it exists. Using Table  7 and Remark 2, we obtain the following results -J 1 is defined by D ≥ D 1 or 0 < D < D 1 and S ch,in < F 1 (D)/Y 3 Y 4 . Therefore, SS1 is the only existing steady state in this region. -J 2 if defined by D 3 < D < D 1 and F 1 (D)/Y 3 Y 4 < S ch,in < F 2 (D)/Y 3 Y 4 . Therefore, both steady state SS2 exist and SS2 is stable since F 3 (D) > 0. -J 3 if defined by 0 < D < D 2 and F 2 (D)/Y 3 Y 4 < S ch,in and F 4 (D, S ch,in /Y 3 Y 4) > 0 when 0 < D < D 3 . Therefore, SS3 exists and is stable, both steady state SS2 exist and SS2 is unstable since F 3 (D) < 0.. -J 4 if defined by 0 < D < D 3 and F 1 (D)/Y 3 Y 4 < S ch,in < F 2 (D)/Y 3 Y 4 . Therefore, both steady state SS2 exist and SS2 is unstable since F 3 (D) < 0. -J 5 if defined by 0 < D < D 3 , F 2 (D)/Y 3 Y 4 < S ch,in and F 4 (D, S ch,in /Y 3 Y 4) < 0. Therefore, SS3 exists and is unstable and both steady state SS2 exist and SS2 is unstable since F 3 (D) < 0.. Proof [Propositions 5, 6 and 7] The result follows from Table  7 and Remark 2. The details are as in the proof of Proposition 5.