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

The complexity of the anaerobic digestion process has motivated the development of complex models, such as the widely used Anaerobic Digestion Model No. 1. However, this complexity makes it intractable to identify the stability profile coupled to the asymptotic behaviour of existing steady-states as a function of conventional chemostat operating parameters (substrate inflow concentration and dilution rate). In a previous study this model was simplified and reduced to its very backbone to describe a three-tiered chlorophenol mineralising food-web, with its stability analysed numerically using consensus values for the various biological parameters of the Monod growth functions. Steady-states where all organisms exist were always stable and non-oscillatory. Here we investigate a generalised form of this three-tiered food-web, whose kinetics do not rely on the specific kinetics of Monod form. The results are valid for a large class of growth kinetics as long as they keep the signs of their derivatives. We examine the existence and stability of the identified steady-states and find that, without a maintenance term, the stability of the system may be characterised analytically. These findings permit a better understanding of the operating region of the bifurcation diagram where all organisms exist, and its dependence on the biological parameters of the model. For the previously studied Monod kinetics, we identify four interesting cases that show this dependence of the operating diagram with respect to the biological parameters. When maintenance is included, it is necessary to perform numerical analysis. In both cases we verify the discovery of two important phenomena; i) the washout steady-state is always stable, and ii) a switch in dominance between two organisms competing for hydrogen results in the system becoming unstable and a loss in viability. We show that our approach results in the discovery of an unstable operating region in its positive steady-state, where all three organisms exist, a fact that has not been reported in a previous numerical study. This type of analysis can be used to determine critical behaviour in microbial communities in response to changing operating conditions.

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 [23]. In engineered biological systems, mechanistic modelling reached consensus with the development of the Activated Sludge Models [9], [10] for wastewater treatment processes, followed by the Anaerobic Digestion Model No. 1 (ADM1) [12] a few years later. The development of ADM1 was enabled largely due to the possibilities for better identification and characterisation of functional microbial groups responsible for the chemical transformations within anaerobic digesters. It describes a set of stoichiometric and kinetic functions representing the standard anaerobic processes, remaining the scientific benchmark to the present day. However, there has been a growing awareness that the model should take advantage of improved empirical understanding and extension of biochemical processes included in its structure, to acquire a better trade-off between model realism and complexity [13].

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

The most common models used to describe microbial systems are two-tiered models, which take the form of a cascade of two biological reactions where one substrate is consumed by one microorganism to produce a product that serves as the main limiting substrate for a second microorganism. When the second organism has no feedback on the first organism, the system is known as commensalistic [16], [20]. The system has a cascade structure and the number of steady-states and their (mathematical) stability as a function of model inputs and parameters may be investigated [1], [2], [19]. When the growth of the first organism is affected by the substrate produced by the second organism the system is known as syntrophic. For instance, if the first organism is inhibited by high concentrations of the product, the extent to which the first substrate is degraded by the first organism depends on the efficiency of the removal of the product by the second organism. The mathematical analysis of such a model is more delicate than for commensalistic models, (see for instance early work by [4], [14], [15], [27] and the more recent papers [5], [11], [17], [18], [22]). An important and interesting extension should be mentioned here: [25], [26] analysed an 8-dimensional mathematical model, which includes syntrophy and inhibition, both mechanisms considered by [2] and [6].

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 [28]. Maintenance is defined as the energy consumed by an organism that is used for all biological processes other than growth. In [28] and here, it is analogous to a first-order decay rate constant, or biomass death term. Although the authors were not able to determine the general conditions under which this four dimensional syntrophic consortium was stable, further work has shown that a model with generality can be used to answer the question posed, determining that the two-tiered food-chain is always stable when maintenance is included [18].

More recently, the model described by [28] was extended by the addition of a third organism and substrate to create a three-tiered ‘food-web’ [24], as shown in Fig. 1. In this paper, the existence and stability of the steady-states were determined only numerically. Although the results were important in revealing emergent properties of this extended model, the motivation of this work is to give an analytical study of the model. Moreover, our analysis does not require that growth functions are of the specific form considered and are valid for a large class of growth functions. This is critical as it provides the means by which microbiologists can theoretically test the influence of the growth characteristics of organisms on the properties of the system and the interactions between multiple species.

Fig. 1

Schematic of the three-tier chlorophenol mineralising food-web indicating the flow and conversion of chemical substrates and products in the system.

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

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

The model

The model developed in [24] has six components, three substrate (chlorophenol, phenol and hydrogen) and three biomass (chlorophenol and phenol degraders, and a hydrogenotrophic methanogen) variables. The substrate and biomass concentrations evolve according to the six-dimensional dynamical of ODEs: where Sch and Xch are the chlorophenol substrate and degrader concentrations, Sph and Xph for phenol and and for hydrogen; Ych, Yph and YH2 are the yield coefficients, represents the part of chlorophenol degraded to phenol, and 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 f1 (see (Eq. 7)) as a product inhibition term. In [24], this model is given in dimensionless form that significantly reduces the number of parameters describing the dynamics. In the analysis of the generalised model (Sections 5 and 7), we do not assume that the growth functions f0, f1 and f2 have the specific analytical expression shown. We will only assume that the growth functions satisfy properties that are listed in Section 4. Therefore, we cannot benefit from the dimensionless rescaling used by [24], because it 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. For the numerical analyses given in Section 8, we return to an assumed set of Monod growth functions given by (Eq. 7), in order to directly compare with the results found by [24].

Here, apart from the four operating (or control) parameters, which are the inflowing concentrations Sch, in, Sph, in, and the dilution rate D, that can vary, all others have biological meaning and are fixed depending on the organisms and substrate considered. We use the following notations in ((1), (2), (3), (4), (5), (6)): With these notations ((1), (2), (3), (4), (5), (6)) can be written as follows: Furthermore, we restrict our analysis to the case where we only have one substrate added to the system, such that: and . As shown in [24], the general case with and present many important and interesting behaviours and deserves future work.

Model reduction

To ease the mathematical analysis, we can rescale the system ((8), (9), (10), (11), (12), (13)) using the following change of variables adapted from [18]: 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 ((14), (15), (16), (17), (18), (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 ((14), (15), (16), (17), (18), (19)) with the following Monod-type growth functions: where: For the numerical simulations we will use the nominal values in Table 1 given in [24].

Table 1

Nominal parameter values. We use units expressed in Chemical Oxygen Demand (COD) as used by [12], [24].

Parameters	Nominal values	Units
km, ch	29	kgCODS/kgCODX/d
KS, ch	0.053	kgCOD/m3
Ych	0.019	kgCODX/kgCODS
km, ph	26	kgCODS/kgCODX/d
KS, ph	0.302	kgCOD/m3
Yph	0.04	kgCODX/kgCODS
k_m,H_2	35	kgCODS/kgCODX/d
K_S,H_2	2.5×10^−5	kgCOD/m3
K_S,H_2,c	1.0×10^−6	kgCOD/m3
Y_H_2	0.06	kgCODX/kgCODS
kdec, i	0.02	d^−1
K_I,H_2	3.5×10^−6	kgCOD/m3

Assumptions on the model and notations

Our study does not require that growth functions are of Monod type (Eq. 23). Actually, the results are valid for a more general class of growth functions satisfying the following conditions, which concur with those given by (Eq. 23):

For all s0 > 0 and s2 > 0 then and .

For all s1 > 0 and s2 ≥ 0 then and .

For all s2 > 0 then and .

For all s0 > 0 and s2 > 0,

For all s1 > 0 and s2 > 0,

For all s2 > 0, .

The function is monotonically increasing and the function is monotonically decreasing.

The proof of the following result is standard and hence omitted.

For non-negative initial conditions, all solutions of the system (Eqs. 14–19) are bounded and remain non-negative for all t > 0.

In the following lemma we define the functions M0(y, s2), M1(y, s2) and M2(y), which are the inverse functions of the functions s0↦μ0(s0, s2), s1↦μ1(s1, s2) and s2↦μ2(s2), respectively.

Let s2 ≥ 0 be fixed. There exists a unique function such that for s0 ≥ 0, s2 ≥ 0 and we have: Let s2 ≥ 0 be fixed. There exists a unique function such that for s1 ≥ 0, s2 ≥ 0 and we have: There exists a unique function such that, for s2 ≥ 0 and we have:

Using the functions M0, M1 and M2 we define now the function ψ(s2, D), F1(D), F2(D) and F3(D). Let us observe first that by H7, for and there exist unique values and see Fig. 2(a): Assume that ω < 1. Let D be fixed such that and exist. We consider the function defined on by: It should be noted that ψ(s2, D) > 0 for . The function ψ together with the values all depend on D. However, to avoid cumbersome notation we will use the more precise form ψ(s2, D), and only if necessary. From (Eq. 25), (Eq. 26) and (Eq. 28) we deduce that: Therefore, we have, see Fig. 2(b): Hence, the function ψ(s2), which is positive and tends to at the extremities of the interval has a minimum value on this interval. We add the following assumption:

Fig. 2

Graphical definitions. (a): and . (b) : ψ(s2), F1(D) and F2(D). (c): and .

The function ψ has a unique minimum on the interval and is negative on and positive on respectively.

Table 2

Parameter values of the maintenance terms a, for cases (a), (b) and (c) of Fig. 3. Unspecified parameter values are as in Table 1. The table gives the values of D1, D2 and D3 where and .

	K_S,H_2,c	ai	D1	D2	D3
(a)	1.0×10^−6	0.02	0.432	0.373	0.058
		0	0.452	0.393	0.078
(b)	4.0×10^−6	0.02	0.329	0.236	I_3=I_2
		0	0.349	0.256	I_3=I_2
(c)	7.0×10^−6	0.02	0.287	I_2=∅
		0	0.303	I_2=∅

Table 3

Parameter values of the maintenance terms a, for case (d) of Fig. 3: and . Unspecified parameter values are as in Table 1. The table gives the values of D1, D2, D2 and D3 where and .

	ai	D1	D2min	D2max	D3
(d)	0.02	0.238	0.101	0.198	0.161
	0	0.258	0.121	0.218	0.181

The value depends on D. However, to avoid cumbersome notation we will use the more precise form only if necessary.

We consider the function F1(D), which is the infimum of ψ(s2, D), with respect to the variable s2, see Fig. 2(b): We also define the functions F2(D) and F3(D): The function F1(D) is defined for D ∈ I1 where: The functions F2(D) and F3(D) are defined for D ∈ I2 where: For all for D ∈ I2, F1(D) ≤ F2(D). The equality holds if, and only if, that is, that is if, and only if, . We define For the Monod-type growth functions (Eq. 23), straightforward computations show that the functions M0, M1 and M2 are given respectively by:

For For For Moreover, we have: and Hence, for all so that the function s2↦ψ(s2, D) is convex and, thus, it satisfies assumption H8, see Fig. 2(b). The minimum is a solution of an algebraic equation of degree 4 in s2. Although mathematical software, such as Maple, cannot give solutions explicitly with respect to the parameters, 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 and, hence, can be obtained numerically. Since M2 and ψ are given explicitly the functions F2(D) and F3(D) are given explicitly with respect to the biological parameters in (Eq. 23). Since is increasing and is decreasing, and assuming the domain of definition I1 of F1(D) is an interval where D1 is the solution of see Fig. 3. For the domain of definition I2 of F2(D), several cases can be distinguished. I2 is an interval where D2 is the solution of see Fig. 3(a), or the solution of equation see Fig. 3(b). I2 is empty, see Fig. 3(c). I2 is an interval where D2 and D2 are the solutions of and respectively, see Fig. 3(d).

Fig. 3

Graphs of and (in black) and (in red) and graphical depiction of where D1 is the solution of and I2. (a): where D2 is the solution of . (b) : where D2 is the solution of . (c): I2 is empty. (d) : where D2 and D2 are the solutions of and respectively. Cases (a)–(d) are obtained with the numerical parameter values listed in Tables 2 and 3. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Existence of steady-states

In this paper, as we are restricted only to local stability analysis, any reference to steady-state stability should be considered as local stability. A steady-state of ((14), (15), (16), (17), (18), (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.

The only steady-state of (Eqs. 14–19), for which or is the washout steady-state where and. This steady-state always exists and is always locally stable.

Besides the steady-state SS1, the system can have at most two other type of steady-states.

SS2: x0 > 0, x1 > 0 and where species x2 is washed out, while species x0 and and x1 exist. We show below that, generally, the system can have two steady-states SS2♭ and SS2♯.

SS3: x0 > 0, x1 > 0, and x2 > 0, where all populations are maintained. We show below that the steady-state SS3 is unique if it exists.

We can state now the necessary and sufficient conditions of existence of SS2 and SS3.

If ω ≥ 1 then SS2 does not exist. If ω < 1 then SS2 exists if, and only if,. Therefore, a necessary condition for the existence of SS2 is that D ∈ I1, where I1 is defined by (Eq. 33). If then each solution s2 of equation gives a steady-state where

If ω ≥ 1 then SS3 does not exist. If ω < 1 then SS3 exists if, and only if,. Therefore, a necessary condition of existence of SS3 is that D ∈ I2, where I2 is defined by (Eq. 34). If then the steady-state is given by called the break-even concentrations, and

If then (Eq. 41) has exactly two solutions denoted by and and such that, see Fig. 2(c), If then .

To these solutions, and correspond two steady-states of SS2, which are denoted by SS2♭ and SS2♯. These steady-states coalesce when .

Since F1(D) ≤ F2(D), the condition for the existence of the positive steady-state SS3 implies that the condition 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 then SS3 coalesces with SS2♭ if F3(D) < 0, and with SS2♯ if F3(D) > 0, respectively.

Using (Eq. 20), the conditions and of existence of the steady-state SS2 and SS3, respectively are equivalent to the conditions respectively expressed with respect to the inflowing concentration Sch, in. Moreover, for the nominal parameter values given in Table 1, we have ω ≈ 0.53. Therefore, the steady-states SS2 and SS3 exist if the conditions on the operating parameters stated in Lemmas 3 and 4, respectively are satisfied, see Table 7.

Table 7

Existence (with or without maintenance) and stability (without maintenance) of steady-states. The conditions for existence and stability are satisfied given that the inequalities are observed.

	Existence	Stability
SS1	Always exists	Always stable
SS2♭	s0in>F_1(D)	Always unstable
SS2♯	s0in>F_1(D)	F3(D) > 0 and s0in<F_2(D)
SS3	s0in>F_2(D)	F3(D) ≥ 0 or
		F3(D) < 0 and F_4(D,s0in)>0

Operating diagram

The Operating diagrams show how the system behaves when we vary the two control parameters Sch, in and D (i.e., substrate inflow rate and medium flow-rate per culture volume) in ((1), (2), (3), (4), (5), (6)). This conventional bifurcation plot is used to visualise the existence and stability of steady-states with respect to these operating parameters, as they are the parameters most easily manipulated in a chemostat by which the system under scrutiny can be examined.

According to Remark 2, 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. 4. If we want to plot the operating diagram we must fix the values of the biological parameters. In the remainder of the Section we plot the operating diagrams corresponding to cases (a)–(d) depicted in Fig. 3.

Fig. 4

The curves Γ1 (black), Γ2 (red) and Γ3 (green) for case (a). (i) : regions of steady-state existence, with maintenance. On the right, a magnification for showing the region . (ii) : regions of steady-state existence and their stability, without maintenance. On the right, a magnification for showing the regions and . (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Operating diagram: case (a)

This case corresponds to the parameter values used by [24]. We have seen in Table 2 that the curves Γ1 and Γ2 are defined for D < D1 and D < D2, respectively and that they are tangent for where and . Therefore, they separate the operating plane (Sch, in, D) into four regions, as shown in Fig. 4(i), labelled and and .

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. 4(i).

Table 4

Existence of steady-states in the regions of the operating diagrams of Fig. 4(i) and Fig. 7(i).

Region	Steady-states
J_1	SS1
J_2∪J_4	SS1, SS2♭, SS2♯
J_3	SS1, SS2♭, SS2♯, SS3

Operating diagram: case (b)

This case corresponds to the parameter values used by [24], except that is changed from to . We have seen in Table 2 that the curves Γ1 and Γ2 are defined for D < D1 and D < D2, respectively and F1(D) < F2(D) for all D < D2, where and . Therefore, they separate the operating plane (Sch, in, D) in three regions, as shown in Fig. 5(i), labelled and .

Fig. 5

The curves Γ1 (black), Γ2 (red) and Γ3 (green) for case (b). (i) : regions of steady-state existence, with maintenance. (ii) : regions of steady-state existence and their stability, without maintenance. On the right, a magnification for 0 < D < 0.1. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

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. 5(i). Note that the region has disappeared.

Table 5

Existence of steady-states in the regions of the operating diagram of Fig 5(i).

Region	Steady-states
J_1	SS1
J_4	SS1, SS2♭, SS2♯
J_3	SS1, SS2♭, SS2♯, SS3

Operating diagram: case (c)

This case corresponds to the parameter values used by [24], except that is changed from to . We have seen in Table 2 that the curve Γ1 is defined for and that I2 is empty so that SS3 does not exist. Therefore, Γ1 separates the operating plane (Sch, in, D) in two regions, as shown in Fig. 6(i), labelled and .

Fig. 6

The curve Γ1 for case (c). (i) : regions of steady-state existence, with maintenance. (ii) : regions of steady-state existence and their stability, without maintenance. On the right, a magnification for 0 < D < 0.1.

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. 6(i). Note that the region of existence of SS3 has disappeared.

Table 6

Existence of steady-states in the regions of the operating diagram of Fig 6(i).

Region	Steady-states
J_1	SS1
J_4	SS1, SS2♭, SS2♯

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 D2 such that for D < D2 the system cannot have a positive steady-state SS3. This case corresponds to the parameter values used by [24], 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 < D1 and D2 < D < D2 and that they are tangent for where and . Therefore, Γ1 and Γ2 separate the operating plane (Sch, in, D) in four regions, as shown in Fig. 7(i), labelled and . 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. 7(i).

Fig. 7

The curves Γ1 (black), Γ2 (red) and Γ3 (green) for case (d). (i) : regions of steady-state existence, with maintenance. (ii) : regions of steady-state existence and their stability, without maintenance. On the right, a magnification for 0 < D < 0.1. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Local stability without maintenance

We know that SS1 is always stable. However, the analytical study of the stability of SS2 and SS3 is very difficult because the conditions for Routh-Hurwitz in the 6-dimensional case are intractable. For this reason we will consider in this section the question of the stability only for the case without maintenance (), since the system reduces to a three-dimensional system. The general case will be considered only numerically in Section 8. When maintenance is not considered in the model, the steady-states SS1, SS2 and SS3 are given by:

where s2 a solution of equation and

where

Let be a steady-state. Then SS2 is stable if, and only if, μ2(s2) < D and.

Therefore, SS2♭ is always unstable and SS2♯ is stable if, and only if, μ2(s2) < D. This last condition is equivalent to which implies that F3(D) > 0. Hence, if SS3 exists then SS2♯ is necessarily unstable. Therefore, SS2♯ is stable if, and only if, F3(D) > 0 and SS3 does not exist.

Let be a steady-state. If F3(D) ≥ 0 then SS3 is stable as long as it exists. If F3(D) < 0 then SS3 can be unstable. The instability of SS3 occurs inparticular when s2 is sufficiently close to i.e. when SS3 is sufficiently close to SS2♭.

The condition F3(D) ≥ 0 is equivalent to namely . If viz. then SS3 can be unstable.

When D is such that F3(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 be a steady-state and evaluated at the steady-state SS3 defined by (48), i.e. for For D ∈ I3 and we define: where . Notice that to compute we must replace x0, x1, x2, s0, s1 and s2 by their values at SS3, given by (48). Hence, this function depends on the operating parameters D and . For each fixed D ∈ I3, is polynomial in of degree 3 and tends to when tends to . Therefore, it is necessarily positive for large enough . The values of the operating parameters D and for which is positive correspond to the stability of SS3 as shown in the following proposition.

Let be a steady-state. If F3(D) < 0 then SS3 is stable if, and only if,.

The results on the existence of steady-states (with or without maintenance) of Lemmas 2–4, and their stability (without maintenance) of Props 2–4, are summarised in Table 7.

This case corresponds to the parameter values used by [24] 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 < D1 and D < D2, respectively and that they are tangent for where and . Beside these curves, we plot also on the operating diagram of Fig. 4(ii), the curve Γ3 of equation According to Prop. 4, this curve is defined for and it separates the region of existence of SS3 into two subregions labelled and such that SS3 is stable in and unstable in . The other regions and are defined as in the previous section. The operating diagram is shown Fig. 4(ii). It looks very similar to Fig. 4(i), except near the origin, as it is indicated in the magnification for . From Table 7, we deduce the following result

Table 8 shows the existence and stability of the steady-states SS1, SS2 and SS3 in the regions of the operating diagram in Fig. 4(ii).

Table 8

Existence and stability of steady-states in the regions of the operating diagrams of Fig. 4(ii) and Fig. 7(ii)., where S and U indicate stable and unstable steady-states, respectively.

Region	SS1	SS2♭	SS2♯	SS3
J_1	S
J_2	S	U	S
J_3	S	U	U	S
J_4	S	U	U
J_5	S	U	U	U

We see from Table 2 that the curves Γ1 and Γ2 are defined now for and respectively and that F1(D) < F2(D) for all D. Beside these curves, we plot also on the operating diagram of Fig. 5(ii), the curve Γ3 of equation (Eq. 50) which separates the region of existence of SS3 into two subregions labelled and such that SS3 is stable in and unstable in . Therefore, the curves Γ1, Γ2 and Γ3 separate the operating plane (Sch, in, D) into four regions, as shown in Fig. 5(ii), labelled and .

Table 9 shows the existence and stability of the steady-states SS1, SS2 and SS3 in the regions of the operating diagram in Fig. 5(ii).

Table 9

Existence and stability of steady-states in the regions of the operating diagram of Fig. 5(ii).

Region	SS1	SS2♭	SS2♯	SS3
J_1	S
J_3	S	U	U	S
J_4	S	U	U
J_5	S	U	U	U

We see from Table 2 that Γ1 is defined for and that I2 is empty so that SS3 does not exist. Therefore, Γ1 separates the operating plane (Sch, in, D) into two regions, as shown in Fig. 6(ii), labelled and .

Table 10 shows the existence and stability of the steady-states SS1, SS2 and SS3 in the regions of the operating diagram in Fig. 6(ii).

Table 10

Existence and stability of steady-states in the regions of the operating diagram of Fig. 6(ii).

Region	SS1	SS2♭	SS2♯	SS3
J_1	S
J_4	S	U	U

We see in Table 3 that the curves Γ1 and Γ2 are defined for D < D1 and D2 < D < D2 and that they are tangent for where and and . Beside these curves, we plot also on the operating diagram of Fig. 7(ii), the curve Γ3 defined by (Eq. 50), which separates the region of existence of SS3 into two subregions labelled and such that SS3 is stable in and unstable in . Therefore, the curves Γ1, Γ2 and Γ3 separate the operating plane (Sch, in, D) into five regions, as shown in Fig. 7(ii), labelled and .

Table 8 shows the existence and stability of the steady-states SS1, SS2 and SS3 in the regions of the operating diagram in Fig. 7(ii).

Numerical analysis to illustrate 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 5 and 7. The results in this section confirm the results on existence of the steady-states obtained in Section 5 in the cases with or without maintenance and the results of stability obtained in Section 7 in the case without maintenance. These results permit also to extend the analytical results and elucidate the problem of the local stability of SS2 and SS3, which was left open in Section 7.

We endeavoured to find numerically the operating conditions under which SS3 is unstable, previously unreported by [24]. Given that we have determined analytically in Proposition 3 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 [24] whilst also satisfying our conditions. In Fig. 8 we show the case when maintenance is excluded. When magnified, we observe more clearly that region does exist for the conditions described above, and also note that the region occurs in a small area between and which corresponds to the results shown in Fig. 4(ii), and is in agreement with Proposition 5. In Fig. 9 we confirm that region does exist for the conditions described above, when maintenance is included, but could not be determined analytically, the curve Γ3 is absent in Fig. 4(i). Furthermore, we demonstrate that there is evidence of Hopf bifurcation, which occurs along the boundary of F3(D) for values of D < D3 by selecting values of Sch, in (indicated by in Fig. 9) at a fixed dilution rate of and running dynamic simulations for 10000 d. The three-dimensional phase plots, with the axes representing biomass concentrations, are shown in Fig. 10, and show that as Sch, in approaches from emergent periodic orbits are shown to diminish to a stable limit cycle at the boundary (see  Appendix B for details). Subsequently, increasing Sch, in to results in the orbit reducing to a fixed point equilibrium at SS3.

Fig. 8

Numerical analysis for the existence and stability of steady-states for case (a), without maintenance. On the right, a magnification for 0 < D < 0.16.

Fig. 9

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 undetectable by the analytical method. The coordinates labelled are subsequently used to simulate the system dynamics, as shown in the proceeding Fig. 10.

Fig. 10

Three-dimensional phase plane diagrams of the biomass dynamics for showing initial (green dot) and final (red dot) conditions for a dilution rate, and chlorophenol input, Sch, (kgCOD/m3), of (α): - the system converges to SS1, (β): - the system enters a periodic orbit of increasing amplitude, ultimately converging to SS1, (γ): - the system is close to a stable limit cycle, (δ): - the system undergoes damped oscillations and converges to SS3. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Whilst the numerical parameters chosen for this work are taken from the original study [24], there somewhat arbitrary nature leaves room to explore the impact of the parameters on the existence and stability of the steady-states. Case (b), discussed in Sections 6.2 and 7.2, involves a small increase to the half-saturation constant (or inverse of substrate affinity), of the chlorophenol degrader on hydrogen. Following the same approach as with the preceding case, we confirm in Fig. 11(i) the Proposition 6 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 shown in Fig 11(ii). It shows the region that cannot be obtained analytically (cf. Fig. 5(i)). In both cases, region has disappeared, as observed analytically. Additionally, the ideal region where all organisms are present and stable, diminishes.

Fig. 11

Numerical analysis for the existence and stability of steady-states for case (b). (i) : without maintenance. (ii) : with maintenance.

Here, was further increased and confirm the Proposition 7, where the function SS3 never exist and SS2 never stable for the case without maintenance. The extension of this proposition to the case with maintenance, shown in Table 10, produce similar results as shown in the comparison of Fig. 12(i) and (ii).

Fig. 12

Numerical analysis for the existence and stability of steady-states for case (c). (i) : without maintenance. (ii) : with maintenance. On the right, a magnification for 0 < D < 0.1.

With the final investigated scenario, where and we observe once again the presence of all operating regions, without and with maintenance, as shown in Fig. 13. It can be seen that regions and increase at low dilution rates across a much larger range of Sch, in than in the default case (a), and the desired condition (stable SS3) is restricted to a much narrower set of D.

Fig. 13

Numerical analysis for the existence and stability of steady-states for case (d). (i) : without maintenance. (ii) : with maintenance.

As with the previous cases, the numerical analysis for case (d) confirms the Proposition 8 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 holds or and hold, viz. The cases (c) or (d) occur if and only if holds or and hold, viz.

Notice that it is easy to make the difference between case (c) and case (d): the first occurs when and the second when . Since D1 is the positive solution of the algebraic quadratic equation it is possible to have an expression for D1 with respect to the biological parameters. However, this is a complicated expression involving many parameters and the preceding conditions or have no biological interpretation. We simply remark here that the function has a vertical asymptote for and the function M2(D) has a vertical asymptote for . Therefore, if then case (c) occurs, so that a necessary (but not sufficient) condition for case (d) to occur is . If m2 is sufficiently small then case (d) can occur.

The observations from the numerical analysis suggest that the role of the chlorophenol degrader as a secondary hydrogen scavenger is critical in maintaining full chlorophenol mineralisation and system stability, particularly at higher dilution rates, as shown by comparing cases (c) and (d) . More significantly, the results coupled with the parameter relationships shown in ((51), (52), (53), (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 give conditions for the existence and stability of the steady-states in the case that maintenance is excluded from the system. However, with a decay term present, purely analytical determination of stability was not achievable.

We confirm the findings of previous numerical analysis by [24] 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 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 and terminates at the origin without maintenance, with maintenance it is located at F1(0)/Y3Y4 ≈ 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 and there is numerical evidence of Hopf bifurcation occurring and that a limit cycle in SS3 emerges.

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