Emergent behaviour in a chlorophenol-mineralising three-tiered microbial `food web'

Anaerobic digestion enables the water industry to treat wastewater as a resource for generating energy and recovering valuable by-products. The complexity of the anaerobic digestion process has motivated the development of complex models. However, this complexity makes it intractable to pin-point stability and emergent behaviour. Here, the widely used Anaerobic Digestion Model No. 1 (ADM1) has been reduced to its very backbone, a syntrophic two-tiered microbial food chain and a slightly more complex three-tiered microbial food web, with their stability analysed as function of the inflowing substrate concentration and dilution rate. Parameterised for phenol and chlorophenol degradation, steady-states were always stable and non-oscillatory. Low input concentrations of chlorophenol were sufficient to maintain chlorophenol- and phenol-degrading populations but resulted in poor conversion and a hydrogen flux that was too low to sustain hydrogenotrophic methanogens. The addition of hydrogen and phenol boosted the populations of all three organisms, resulting in the counterintuitive phenomena that (i) the phenol degraders were stimulated by adding hydrogen, even though hydrogen inhibits phenol degradation, and (ii) the dechlorinators indirectly benefitted from measures that stimulated their hydrogenotrophic competitors; both phenomena hint at emergent behaviour.


I. INTRODUCTION
Microbial degradation of organic compounds in methanogenic environments is a sequential process catalysed by a series of different micro-organisms. Syntrophy plays a pivotal role in these feeding webs: degradation of compounds like propionate and phenol is only sustainable if their degradation products, hydrogen and acetate, are removed by methanogens. The thermodynamic rational behind syntrophy is well understood, but its kinetic framework is less established. This raises questions about the stability of these feeding chains and the factors that govern them. As a first step towards answering these questions a simple mathematical model was previously developed describing the interactions in a twotiered feeding chain, populated with a set of parameters that apply to propionate degraders and hydrogenotrophic methanogens [1]. Mathematical analysis of the model indicated that the system was always stable: there were no conditions where the populations of the two organisms oscillate or show other forms of emergent behaviour.
The objective of the present paper is to introduce an additional organism into a similar feeding chain and evaluate its effect on stability and potential emergent be- * matthew.wade@ncl.ac.uk haviour of the resulting 'food web'. The organism of choice is a chlorophenol-dechlorinating bacterium. The other two organisms are a phenol degrader and a hydrogenotrophic methanogen. The complete removal of phenolic compounds from the system is hereby referred to as chlorophenol mineralisation. The salient feature of the chlorophenol degrader here is that production of phenol is coupled to consumption of hydrogen by hydrogen cycling. Thus, as a hydrogen consumer, the dechlorinator competes with the methanogen for hydrogen [2,3]. The working hypotheses are (i) that the dechlorinator can (partially) replace the methanogen as the syntroph in a phenol-degrading consortium, and (ii) that introduction of this organism can potentially lead to unexpected emergent behaviour related to the intricacies of the multispecies relationships.
It has been shown that deterministic modelling of biological systems, typically through a system of coupled ordinary differential equations (ODEs), provides important understanding of these often complex processes, specifically in determining changes to the system behaviour given perturbations in the inputs. For anaerobic digestion, higher dimensional models are useful for capturing the phenomenological behaviour of the multi-step processes and are often the de-facto method for understanding plant operation [4]. On the other hand, simplified or reduced models have received more attention in process monitoring, control design [5,6], and optimisation [7]. Simplified models have also been applied to determine the global [8,9] or local [10] stability of the system under investigation. Typically, these are two-species models, whereas here we present a case considering a threespecies food web.
Whilst a strictly analytical approach is not possible given the dimensions of the ensuing model, the approach taken here derives analytic expressions for all steadystates, supported by numerical simulations to determine the regions of local stability within sensible operating conditions. As such, it is possible to gain a formal understanding of the emergent properties of all states.

A. Mechanistic model of a three-tiered food web
The model developed here is based on Anaerobic Digestion Model No. 1 (ADM1) [4]. The general model here has six components, three substrate and three biomass variables, from which a sub-model describing phenol degradation and the extension of the full model to include addition of extraneous substrates, are formed and described in the relevant sections. The chlorophenol degrader utilises both chlorophenol and hydrogen for growth, producing phenol as a product. Phenol is consumed by the phenol degrader forming hydrogen, which also is inhibitory to its growth. The methanogen scavenges this hydrogen and acts as the primary syntroph.
The time-dependent substrate and biomass concentrations are denoted S(t) and X(t), respectively, with subscripts ch , ph and H2 referring to chlorophenol, phenol and hydrogen components, respectively. The growth functions are of Monod form with the inclusion of a product inhibition term, K I,H2 . Chlorophenol, phenol and hydrogen are introduced with an input concentration S ch,in , S ph,in , S H2,in , respectively, and a dilution rate D. The inhibition of hydrogen on the phenol degrader is defined as . (1) The substrate and biomass concentrations then evolve according to the six-dimensional dynamical system of ODEs dX ch dt =−DX ch +Y ch f 0 X ch −k dec,ch X ch (2) dX ph dt =−DX ph +Y ph f 1 X ph −k dec,ph X ph (3) dX H2 dt =−DX H2 +Y H2 f 2 X H2 −k dec,H2 X H2 (4) dS ch dt =D(S ch,in −S ch )−f 0 X ch (5) dS ph dt =D(S ph,in −S ph )+ 224 208 Here the functions f 0 , f 1 and f 2 are defined as f 0 (S ch , S H2 ) = k m,ch S H2 K S,H2,c +S H2 S ch K S,ch +S ch (8) f 1 (S ph , S H2 ) = k m,ph S ph K S,ph +S ph I 2 (9) where K S,H2,c and K S,H2 are the half-saturation constants for hydrogen in the chlorophenol degrader and hydrogenotrophic methanogen, respectively. The values used for the various parameters are listed in Table I, and their derivation is presented A. The value 224/208 represents the fraction of chlorophenol chemical oxygen demand (COD) converted to phenol, 32/224 is the fraction of phenol converted to hydrogen, and 16/208 is the fraction of hydrogen COD consumed by X ch .

B. Dimensionless form
It is beneficial to scale unit-dependent equations to a dimensionless form; this significantly reduces the number of parameters describing the dynamics, thereby simplifying the subsequent analyses. Using the notation of Xu et al. [1], after Baltzis and Fredrickson [11], the following dimensionless terms are defined With these transformations, the dynamical system of ODEs given in Equations (2-7) reduces to where, for simplicity, the following parameters are introduced

C. Steady-states
The steady-states of this system of ODEs are obtained by setting g i = 0 (for i = 1, . . . , 6) and solving simultaneously. This predicts eight possible steady-states: SS8:: The phenol degraders and methanogens are maintained in the system (x 0 = 0, x 1 = 0, x 2 = 0).
For full chlorophenol mineralisation and, as such, the only desired operating condition, SS6 must be stable.

D. Stability analysis
The challenge of analytically characterising dynamical systems of ODEs is well understood. Typically, for systems of high dimensions, it is intractable to finding explicit solutions and one must resort to obtaining solutions numerically. For example, the Routh-Hurwitz theorem allows for an explicit analysis of the stability of steady-states, but is intractable beyond five dimensions [12] and is excluded here. Relying on numerical solutions, however, is problematic as the extent and resolution of the results are limited by the choice of model parameters and the computational resources. Nevertheless, here analytical expressions may be found for some of the steady-states that can inform general rules about their viability. To present the results of the stability analysis coherently, all steady-states are solved numerically for a realistic range of operational and kinetic parameters with a suitable resolution of 5000×5000 solution points. A brief discussion of the numerical methods used in this analysis is provided in B.
In order that a steady-state be meaningful, all variable concentrations must be non-negative, while Eq. (5) also gives the condition that S ch < S ch,in , (or s 0 < u f in dimensionless form). For model extensions with phenol and hydrogen addition, the conditions that S ph < S ph,in (s 1 < u g ) and S H2 < S H2,in (s 2 < u h ) are also necessary. It is well established that the stability of a system of autonomous ordinary differential equations (ODEs) can be determined by investigating the eigenvalues of the corresponding Jacobian matrix [13]. The Jacobian for the system of Equations (11)(12)(13)(14)(15)(16) corresponds to the 6×6 matrix Note that for the two-tier model, for which X ch = S ch = 0, the Jacobian reduces to a 4×4 matrix. The Jacobian is then evaluated at a given steady-state, denoted J SSi (i=1,..., 8), and its eigenvalues calculated. If the real parts of all the eigenvalues obtained from the Jacobian are negative, then the state is stable. If one or more of the eigenvalues have a positive real part then the steady-state is unstable.
In Section III, a thorough investigation is carried out in order to determine when each steady-state is viable. Firstly, a two-tier model describing only phenol degradation is considered, where x 0 = s 0 = 0, reducing the system of ODEs to four equations. Subsequently, the full model for chlorophenol mineralisation is studied.

III. RESULTS
A. Two-tier phenol model An initial analysis was performed on the two-species model feeding on phenol. As expected from the model structure, the results are similar to those shown for a propionate-degrading bi-culture with a maintenance term (k dec ) [1], where three stable steady-states emerge: SS1, SS7 and SS8. Figure 1 shows the steady-state diagram for the phenol model, demonstrating that the system has three mutually exclusive stable steady-states under the range of parameters and operating conditions chosen. Numerical simulations under conditions producing SS7 (not shown) indicated that the concentration of phenol degraders are comparatively low (compared to SS8 populations), resulting in low hydrogen production. This in turn results in the washout of the methanogen, whilst the phenol degraders, under significant phenol concentrations, can be maintained up to dilution rates equivalent to their theoretical maximum growth rate minus the decay constant, µ max,ph −k dec,ph = 1.02d −1 . Although numerical analysis of the steady-states was performed here, it is possible to get explicit quadratic functions of the parameter pair (S ph,in , D) for the steady-state partitions. Using the method described in Xu et al. [1], the functions for the partitions are (in non-dimensionless form) where where F 1 = 0 relates to methanogen washout and F 2 = 0 to phenol degrader washout. These partitions are shown by the green lines in Fig. 1.
Having observed that the phenol model is always stable and the desired operational state (SS8) is bounded asymptotically at higher phenol input concentrations, the next sections explore the properties of the model when extended to a three species dechlorinating tri-culture with hydrogen cycling.
Numerical simulation of the model was performed using the parameters listed in Table I, over a range of operating conditions (α,u f ). The stable steady-state regions are shown in Fig. 2. As determined from the stability analysis, SS1 is meaningful under all operating conditions and bistable steady-states with SS4 and SS6 are found in specific regions within the numerical limits of the simulation, depending on the initial conditions. An interesting phenomenon not observed in previous two-tiered models is the existence of SS1 at low dilution rates but increasing chlorophenol input. It is likely that the condition u f > s 0 is not satisfied in this region and suggests that the system does not produce enough metabolites to sustain the syntrophic populations at such low substrate concentrations and flow rates.

C. Three-tier chlorophenol model with hydrogen addition
Given that the model under analysis includes the possibility of adding extraneous substrates (phenol and hydrogen) to the system, supplementing the flux produced by the biomass, it is fairly straightforward to extend its mathematical analysis with these terms included. It is hypothesised that the addition of hydrogen to the system will result in the retention of the methanogen at dilution rates up to its maximum growth rate (µ max ). It is anticipated that this leads to the extension of full chlorophenol mineralisation by allowing the establishment of a methanogenic population that can be maintained given the availability of hydrogen in high enough concentrations.
When hydrogen addition is included, the hydrogen input term, S H2,in > 0 in Eq. (7), and the dimensionless form, u h > 0 in Eq. (16). A further three steady-states are defined with this model; SS2, SS3 and SS5.
Following the same approach as for the three-tier chlorophenol model, six Jacobian matrices were derived for each steady-state, and their eigenvalues λ i (where i = 1, . . . , 6) found. In this case, the conditions for stability for each steady-state are as follows SS1:: Meaningful and stable providing Given that, in this model, k A = k C , then for φ 2 > 1, the second condition becomes the prevailing factor in determining stability, and can be reduced to µ max, SS2:: Meaningful when α < φ 2 −k C and u h > s 2 and stable when Here it can be noted that the stability of SS2 does not depend on u h , and as such will remain fixed for any hydrogen input concentration.
SS3:: Meaningful given the roots of the quadratic function (Eq. (C3)) for s 0 are positive, u f > s 0 and u h > ω 2 (u f −s 0 ). For steady-state stability, it is necessary that α > µ 1 −k B , α > µ 2 −k C and for R(λ i ) < 0 (i = 4, 5, 6). Here, the stability of SS3 is dependent on both u f and u h . The first condition dictates that at higher u f , a greater amount of u h is also required to maintain stability, whereas the second condition indicates that at a fixed u f , increasing u h will increase the dilution rate under which SS3 is stable up to a boundary described by a function f (u f , µ max,H2 −k C ).
SS4:: Meaningful given the roots of the cubic function (Eq. (C11)) for s 0 are positive and the conditions u h > ω 1 s 1 +s 2 and x 1 > 0 are met. Stability is assured when α > µ 2 −k C and for R(λ i ) < 0 (i = 2, . . . , 6). As u h affects s 0 through the highly non-linear cubic expression, it is difficult to draw any conclusion about its effect on stability by analytical means.
SS5:: Meaningful when the conditions φ 2 > α+k C , . The analysis shows that SS6 stability is affected by hydrogen addition.
The results from the stability analysis are visualised in Fig. 3 for different concentrations of hydrogen addition. As can be seen, four steady-states (SS3-6) are affected by increasing concentrations of hydrogen. At very low hydrogen concentration addition ( Fig. 3 (a)), SS3 is stable in the region between (0,0) and SS2 (SS1 is not shown as it remains fixed at D = M max,H2 −k dec = 2.08d −1 , where M max is the non-dimensionless form of µ max ). Moreover, at lower dilution rates, the system has two observable bistable conditions: 1) between SS3 and SS6 at the lowest dilution rates; and 2) between SS3 and SS4, before SS3 becomes the single stable steady-state at higher dilution rates. This can be rationalised as follows: • The chlorophenol degrader can outcompete the methanogen for the additional hydrogen in the system. Depending on the initial biomass concentrations, this can result in complete washout of the phenol degrader and methanogen at low dilution rates. However, at higher chlorophenol concentrations, the system can be in stable equilibrium with all organisms present; • As the dilution rate increases, the methanogen is washed from the system as the chlorophenol degrader outcompetes it to become the sole hydrogen scavenger, as seen in Fig. 2; • Eventually at higher dilution rates, the phenol degrader is washed out and full chlorophenol mineralisation is no longer possible. Nevertheless, beyond the theoretical M max,ch , the methanogen is able to grow again (SS2), as there is a constant supply of substrate without an active competitor.
With increasing hydrogen addition interesting phenomena are observed ( Fig. 3 (b)-(d)). Firstly, the bistability between SS3 and SS6 is replaced by a single steadystate, in which all organisms are present. Similarly, the SS3-SS4 bistability is replaced by only SS4, although this requires greater hydrogen concentrations at higher chlorophenol input. A new steady-state also begins to emerge (SS5), starting at low chlorophenol input, which eventually subsumes the other two varying steady-states (SS3 and SS4) at high hydrogen concentrations. The understanding from this is that at dilution rates higher than the upper boundary of SS6 and below the lower bound-ary of SS2, the phenol degrader is naturally washed out, but the presence of increasing amounts of hydrogen allows the methanogen to exist whilst in competition with the chlorophenol degrader.
Most importantly, the extent of SS6 is also observed to increase below a specific S ch,in , with the additional hydrogen resource stabilising the methanogen population at higher dilution rates, limited now only by its maximum growth rate (M max,H2 ). Above this dilution rate, as with the three-tiered model, the production of phenol is limited by the stoichiometry of chlorophenol degradation. This restriction results in a limited supply of phenol under these operating conditions, which cannot produce enough growth in the phenol degrader to sustain its population. However, the increasing abundance of additional hydrogen leads to the establishment of a two-species trophic level with no intermediary organism (SS5). Nevertheless, the ability to extend the stable re-gion of SS6 with addition of hydrogen is of note, particularly for practitioners.
It has been shown that the inclusion of hydrogen addition in the three-tiered system can lead to an increase in the stable region for full chlorophenol mineralisation with regard to standard operating parameters. However, this is limited to a defined operating range. At lower S ch,in there is washout of phenol degraders and establishment of methanogens at increasing hydrogen concentrations. It appears that a form of competitive exclusion occurs at certain (D, S ch,in ). However, unlike the classical case, here the exclusion principle dictates that consumer A (X ch ) is never excluded and either consumer B (X H2 ) or producer A (X ph ) are washed from the system. This can be seen in the system simulation shown in Fig. 4, which considers a single operational point ( At low hydrogen addition, the concentration is not enough to sustain the methanogen population whilst a small population of dechlorinators are present (6×10 −5 kgCOD/m 3 ). The hydrogen not utilised by the chlorophenol degrader reaches an equilibrium that inhibits the phenol degraders and they are washed out of the system and, thus, SS3 is the dominant steady-state under these conditions for the initial conditions tested. As hydrogen addition is increased, an interesting property of the system emerges. Again, the chlorophenol degrader is able to utilise the extra hydrogen to produce more biomass and more phenol. However, the concentrations of additional hydrogen are not enough to sustain the methanogen population and the equilibrium concentration returns to a value approximate to the previous case. Whilst this exerts an inhibitory effect on the phenol degraders, the additional phenol substrate availability is enough to maintain a population in the system (steadystate not shown in Fig. 4), which switches the system to SS4. In the final case, hydrogen addition is further increased, allowing chlorophenol degraders to utilise most of the available chlorophenol and reach close to its maximum growth. As this growth is asymptotic, and because of the surplus hydrogen available, the methanogens are now able to utilise the remaining concentrations to maintain their population in the system. This results in an increased hydrogen equilibrium that exerts a greater inhibition on the phenol degraders, to such an extent that the additional phenol availability is not enough to avoid washout and the system moves to SS5.
It can be seen that this situation is fairly complex and is dependent on a number of factors. Principally, the equilibrium concentration of hydrogen is determined by the dilution rate, which in turn limits the growth of either producer A or consumer B.

D. Three-tier chlorophenol model with phenol addition
For such a non-linear system the addition of hydrogen could lead to deleterious and potentially unexpected behaviour, especially given its inhibitory effect on the phenol degraders. It is therefore useful to examine the effect of adding phenol to the system, both as a sole exogenous substrate, but also in combination with hydrogen addition. In this way, it may be determined under what conditions, if any, full chlorophenol mineralisation can be extended in relation to the two operational parameters (D and S ch,in ).
When phenol is added to the system the phenol inflow term, S ph,in from Eq. (6) and u g from the dimensionless Eq. (15), become positive. An additional two steadystates are defined with this model, SS7 and SS8.
The conditions for stability in each steady-state, determined in the same manner as previously described, are as follows: However, at lower phenol concentrations, the effect of phenol lowers the minimum dilution rate under which SS1 is stable.
SS2:: Never stable as a contradiction occurs in the terms describing s 2 and x 2 . For these variables to be meaningful, the conditions α−k C < φ 2 and α−k C > φ 2 must be satisfied, which is invalid (see Eqs. (C1) and (C2)).
SS3:: Never stable as a contradiction occurs in the terms describing s 2 and x 0 . With u h = 0, s 0 > u f for s 2 to be positive. This condition, however, results in x 0 always being negative and meaningless (see Eqs. (C4) and (C5)).
In the special case that hydrogen is also added to the system, it should be noted that if u g < s 1 , then the condition u h > s 2 must be strictly observed. Similarly, if u h < s 2 , then u g must be greater than s 1 .
SS5:: Never stable as without u h then x 2 can never be positive (see Eq. (C29)).
SS7:: Meaningful given the roots of the quadratic function for s 1 (Eq. (C50)) are positive and, subsequently, φ 1 s 1 /(1+s 1 )K I (α+k B ) > 1/K I , and u g > s 1 . Therefore, stability is assured when α > µ 0 −k A , α > µ 2 −k C and the roots of the characteristic polynomial (Eq. (C53)) have negative real parts. It can be seen that the stability of SS7 relies on u g , but is also influenced by hydrogen if added to the system.
SS8:: Stable when the following conditions are satisfied. Firstly, s 2 is meaningful given φ 2 > α−k c . Stability is assured when s 2 is meaningful, α > µ 0 −k A , and the roots of the characteristic polynomial (Eq. (C62)) have negative real parts. In general, for SS8 to be stable, α must be smaller than µ max,H2 and, therefore, u f should be small.
Again, the operational parameter plots are used to visualise the outputs from the stability analysis, as shown in Figs 5 (a)-(c). At low phenol addition, emergent steady-states (SS6-SS8) are observed at low dilution rates, with the bistability present in the standard threetiered model comprising the rest of the system below the boundary for SS1. As phenol is increased, the two bistable regions disappear to be replaced by a single stable steady-state (SS4 and SS6). Both SS7 and SS8 also appear across a greater range of dilutions, although SS8 is confined to low chlorophenol addition, whilst SS8 eventually covers a region contained by the upper bounds of SS7 (low S ch,in ) and SS4 (high S ch,in ), and the boundary of SS1. Of greater interest is the observation that the extent of the desired condition, SS6, is greater for S ch,in between approximately 0 and 1kgCOD/m 3 , replacing both SS4 and part of SS7.
E. Three-tier chlorophenol model with bi-substrate addition Figure 5 (d) shows an example when both hydrogen and phenol are added to the system. Here it can be seen that the two additional substrates each contribute to distinct emergent properties. The hydrogen addition allows for the maintenance of a methanogen population at higher dilution rates, whilst the phenol addition maintains the phenol degraders at lower dilution rates. Both additional substrates provide the conditions necessary for full chlorophenol mineralisation within a specific operational parameter space previously associated with methanogen washout.
An extension of this analysis was undertaken by examining the effect of the two supplementary substrate additions on the maintenance of SS6 at higher dilution rates than observed in the standard three-tiered model. Figure 6 shows the steady-state plot for varying concentrations of hydrogen and phenol addition for two distinct operating conditions observed previously to give rise to SS6 with addition of an extraneous substrate. Figure 6 (a) demonstrates that both hydrogen and phenol addition can transform the system from undesirable bistability to a stable process with all three organisms present. Indeed, it is shown here that addition of either substrate can lead to SS6 under relatively high chlorophenol input conditions.
An alternative situation under low chlorophenol addition is shown in Figure 6 (b). Here, it is clear that phenol addition at concentrations above 0.65kgCOD/m 3 provides the necessary conditions for full chlorophenol mineralisation, irregardless of hydrogen addition. Although it can also be seen that hydrogen addition does have some effect within a very small range of phenol concentration between about 0.63 and 0.65kgCOD/m 3 , allowing for stability of the methanogenic population and moving the steady-state from SS4 to SS6, it is relatively negligible given the influence of phenol with these operating parameters.

IV. BISTABILITY ANALYSIS
Throughout Section III, numerous regions of bistability were observed between SS3-SS4 ( Fig. 3 (a), (b) and (c)), SS3-SS6 (Fig. 3 (a) and (b)) and SS1-SS6 (Fig. 1). In this Section, these bistable regions are investigated in order to provide a more in-depth description and under- standing. Through analysis of these results, the steadystate having the larger basin of attraction is determined, i.e. the most likely outcome for a range of initial conditions. Note that the other four steady-states, even when meaningless or unstable, are still present mathematically and can impact on the overall dynamics of the system.
Firstly, the region of bistability between SS3 and SS4 was considered. Figure 7 shows the final stable steadystates (SS3 or SS4) achieved for a range of initial conditions. Note that for SS3 or SS4, X H2 = 0 and does not deviate from this value throughout, which simplifies the analysis required here. When X ch is small, the system is likely to head towards SS3 where washout of X ph occurs. Otherwise, the system achieves SS4. Very similar results were obtained for biologically meaningful substrate initial conditions between 1.3×10 −2 and 1.5×10 −1 for S ch , 9.0×10 −2 and 4.9×10 −1 for S ph and 2.2×10 −12 and 4.7×10 −12 for S H2 (kgCOD/m 3 ). These results suggest that as long as sufficient dechlorinating biomass is present in the system, the most likely outcome is for both X ch and X ph to be present over time. Thus, the basin of attraction for SS4 is larger than that for SS3. In addition, changing the dilution rate D, S ch,in , or S H2,in within the bistable region once again leads to similar results.
A similar analysis of the bistable regions for SS3-SS6 and SS1-SS6 give a more dramatic outcome. If any of the three organisms, X ph , X ch or X H2 , are initially zero, SS6 is not obtained. If all of these are present, even in tiny quantities, SS6 will be reached. It is known that SS1 is always meaningful, thus a steady-state always occurs at the origin, be it stable or unstable. In the case of bistability, SS3-SS4 or SS3-SS6, the unstable steady-state is closer to SS3 rather than SS4 or SS6, respectively. This instability thus forces the dynamics away from the total washout and towards the state with the least washout. This analysis highlights how not only the main parameters of the system can impact the final stable state, but also the importance of initial conditions.

V. DISCUSSION
In this work, a mathematical analysis of a three-tiered 'food web' comprising three organisms with hydrogen addition and inhibition has been presented. Although mechanistic models of microbial interactions are somewhat ubiquitous, analytical approaches have been limited to simple two-species systems [1,14,15]. More recently, attempts at analysing more complete systems, such as the anaerobic digestion process, have been performed with some degree of rigour [16,17], however the practicality of these approaches is limited by the need for generalisation and increasing numbers of assumptions to allow for mathematical tractability.
Although the three-tiered system is considered only a sub-process of anaerobic digestion and the analysis limited to a specific compound, the approach is general enough to provide information about the characteristics of such a process and allows for the possibility of extending the work to other systems and other compounds. Indeed the results shown here point to the fact that biological knowledge can inform mathematical approaches, whilst mathematical models can indicate or confirm biological properties of a system, succinctly and rationally.
Here, four models were analysed by a combination of analytical and numerical techniques to obtain information regarding the extent and characteristics of system stability in a three-organism process anaerobically degrading chlorophenol. In this case, only the hydrogen part of the pathway was considered, with acetate being excluded from the analysis. The simple two-tiered food chain considering only phenol degradation with two species was shown to have the same characteristics observed as a previously studied system investigating propionate [1]; three steady-states and always stable. A second model introduced chlorophenol as the primary substrate and a dechlorinating organism forming a tri-culture that resulted in the emergence of bistable conditions between the complete washout steady-state and the two other viable states. Analysis of these bistabilities has shown that the desired operating condition (SS6) has the strongest basin of attraction such that complete washout can only occur when the initial condition for one of the biomass concentrations is zero.
Although it is interesting to see the extent under which full chlorophenol mineralisation can occur using this standard model, of greater interest from an engineering perspective is the possibility of driving the system towards its limits of operational viability without compromising its function. With this in mind, the inclusion of additional input terms for hydrogen and phenol was undertaken, and corresponding stability analysis performed. In the case of hydrogen addition, it was expected that this would lead to a wider region of stability for the methanogen and, thus, extend the operational domain for SS6. Indeed, this was the case under relatively high concentrations of chlorophenol input, under which the methanogen population was maintained up to a theoretical maximum, dictated by its maximum growth rate. However at low chlorophenol input (< 0.5kgCOD/m 3 ), the nature of the resource competition between the dechlorinator and methanogen is such that a switching of behaviour occurs in localities not conducive to stability of all three organisms.
Under the standard three-tiered chlorophenol model, resource competition dictates that the dechlorinator is able to utilise the available hydrogen for growth more readily than the methanogen and under SS4 a competitive exclusion principle occurs and the latter is washed out. With hydrogen addition, however, an abundance of this resource allows the methanogen to access enough hydrogen to maintain a population at steady-state. However, below a threshold concentration and under specific dilution rates, an excess of hydrogen leads to inhibition of the phenol degrader and washout. In effect, the addition of hydrogen at this threshold results in a switching between SS4 and SS5, rather than a stabilisation at SS6. The addition of phenol, however, under specific operational conditions can result in the emergence of SS6 under low chlorophenol input conditions. Here, a specific form of competition between the dechlorinator and methanogen occurs (as with SS6 in the standard model), by which both organisms benefit from the production of hydrogen by the intermediate phenol degrader, whilst the phenol degrader is stabilised by the syntrophic hydrogen removal. In other words, the presence of the chlorophenol degrader allows for the production of phenol, which benefits both the phenol degrader and, indirectly, the methanogen. The presence of the methanogen reduces the inhibition on the phenol degrader, thus allowing the hydrogen resource to be maintained for both hydrogenotrophic populations. This can be seen as a form of mutualism rather than competition. The phenol addition plays a significant role under conditions where the chlorophenol degrader cannot produce enough phenol to maintain the phenol degraders, and under such conditions it is possible to achieve full chlorophenol mineralisation beyond the standard model system. A schematic of these models and their interactions for the low chlorophenol input condition is shown in Fig. 8.
Although higher dimensional systems are analytically restricted and lack generality, the work reported here for a three-tiered 'food web' underlines the potential for applying localised stability analysis within meaningful operating and parameter ranges, to identify properties of the system that both increase fundamental understanding of the behaviour of such microbial systems, but also guide thought on ways to manipulate or control them for potential process improvement. Further work will include a thermodynamic rather than a kinetic inhibition term for the effect of hydrogen on the syntroph, and the extension of the model to polychlorinated phenols.

ACKNOWLEDGMENTS
This work was funded by the Biotechnology and Biological Sciences Research Council UK (BB/K003240/2 Engineering synthetic microbial communities for biomethane production) and by the Institute of Sustainability, Newcastle University.

Appendix A: Parameter estimation
The parameters of the model are based on the Chemical Oxygen Demand (COD) to be in consensus with the IWA Anaerobic Digestion Model No. 1 [4]. However, only those related to hydrogen can be found in Hydrogen addition model with phenol degrader washout. Bottom: Phenol addition model with extended stability. Organisms that are always present are shown as a green circle, a red circle is for when they are conditionally present, and a black circle represents biomass washout the ADM1 (for high-rate mesophilic processes). Subsequently, the parameters for chlorophenol and phenol were chosen based on a combination of literature data and deduction based on their chemical oxidation reaction.
a. Parameters for chlorophenol Adrian et al. used the Dehalococcoides mccartyi strain CBDB1 to completely convert 2,3-dichlorophenol, all six trichlorophenols, all three tetrachlorophenols, and pentachlorophenol to lower chlorinated phenols [18], and their results are used partially here, with some deductive reasoning. The observed dechlorination rates in batch cultures with cell numbers of 10 7 mL −1 was found to be 35 µMd −1 . The COD of the CBDB1 strain was first estimated using the empirical formula for bacterial cell C 5 H 7 O 2 N, with cell numbers of 10 7 mL −1 .
The oxidation reaction for C 5 H 7 O 2 N is From this equation, it is easy to show that 1 mole of C 5 H 7 O 2 N requires 5 mole of O 2 , or 160gCOD. The molar mass of C 5 H 7 O 2 N is 113g/mol. Therefore, the proportion of C 5 H 7 O 2 N is 160/113 = 1.42g. By applying this to general organic matter, a factor of 2 is used. This means that each gram of biomass (e.g. strain CBDB1 here) requires 2g of oxygen for growth. Subsequently, the cell number should be converted to grams (dry weight) using a conversion factor of 1.2×10 −14 g/cell [19], or 1.2× 10 −4 g/L.
The COD of the substrate per COD of biomass (k m,ch ) was then calculated using released chlor ions up to 35µM/Ld, which corresponds to COD consumed by biomass of It should be noted that this result does not include monochlorophenol, which have to be deduced. For monochlorophenol, the COD must be converted from the COD of Cl − . By considering the stoichiometry of its oxygenation, the COD of monochlorophenol can be calculated as 6.5× 32 = 208gCOD/mol. From this, the maximum specific uptake rate for monochlorophenol degraders can be estimated as k m,ch = 0.14×208 = 29.12gCOD S /gCOD X d, where COD S is the substrate COD.
For the half-saturation constant, K S,ch , no data was found detailing empirical studies for monochlorophenol degradation. However, from a recent study on the anaerobic microbial dechlorination of pentachlorophenol, some values were cited for dechlorinating bacteria using cosubstrates of the 5-CP with other carbon sources [20]. A K S value of 0.25µM, which converts to 0.053kgCOD/m 3 following the oxidation of monochlorophenol, was selected as a rough estimate for the chlorophenol degrader in this model.
The yield constant, Y ch was calculated using information about the growth yield of the Dehalococcoides mccartyi strain CBDB1 on 2,3-dichlorophenol (DCP), in which 7.6×10 13 cells per mole of Cl − were produced [18]. Since yield is the ratio of COD X to COD S , it can be calculated by determining these components from their stoichiometry; : For 7.6×10 13 cells, the dry weight is calculated as 7.6×10 13 ×1.2×10 −14 g/cell = 0.912g.
Knowing that 2gCOD/gbiomass, stated earlier, then the COD X = 1.824g. COD S : One mole of Cl − released should be converted to COD that DCP consumed. It is known that: DCP→ 2Cl − . That means one mole of Cl − comes from 0.5 mole of DCP. The oxidation equation of DCP is given by Therefore, 1 mole of DCP consumes 6 mole of oxygen, giving 6×32 = 192gCOD, which corresponds to 0.5 mole of DCP using 96gCOD.
Subsequently, the yield of the chlorophenol degrader can now be calculated as The assumed half-saturation constant describing the affinity of the dechlorinator for hydrogen, K S,H2,c was arbitrarily set to 1×10 −6 in the absence of supportable empirical data. The value was based on the known constant for the methanogen, and chosen to be an order of magnitude lower, accordingly.
b. Parameters for phenol There is very limited data available on the kinetics of phenol degradation in anaerobic environments. The parameters selected for the phenol degrader here are approximate values based on available literature and deduction from analogous substrate reactions.
For example, the reaction of phenol is similar to that for propionate, but consumes twice as much oxygen. The oxidation equations for the two chemicals are written as Hence, the activity of phenol is assumed to be double the activity of propionate. From the ADM1 [4], the maximum specific activity for propionate (k m,pro ) is 13gCOD S /gCOD X , so k m,ph is assumed to be 26gCOD S /gCOD X .
The half-saturation constant for phenol was taken from a study by Eismann et al. [21] using a modified Haldane kinetic model from Edwards [22], which considers the inhibitory effect of high substrate loading, for data fitting. At 35 • C, the K S,ph determined by the authors was 127mgPhenol/l, which converts to a value of 0.302kgCOD/m 3 used in this model.
In the absence of compelling literature data, the yield parameter for the phenol degraders, Y ph , was taken to be equivalent to that for another acetogenic bacteria, the propionate degrader described in previous work [1] and having a value of 0.04gCOD S /gCOD X .

Appendix B: Numerical Methods
Although it is possible to simplify the system and obtain analytical expressions for characteristic polynomials for all of the steady-states solved in C, it is impossible to gather the necessary information from these to determine the regions where each steady-state is meaningful and/or stable. Therefore, by considering sets of parameters (D, S ch,in , S ph,in and S H2,in ) in turn on each steady-state, and through the use of Matlab, the complex polynomials for each steady-state in turn can be solved to determine which are meaningful and stable. Repeating this method with numerous sets of parameters leads to results showing where each steady-state is stable.
In addition, some cases (SS3 or SS4) do not allow straightforward expressions for the steady-state (quadratic and cubic equations for s 0 are obtained, respectively). In these situations, the above procedure is followed for all possible values of s 0 at these steadystates.
Appendix C: Stability analysis for the three-tier chlorophenol model with substrate addition Here, the expressions used for each steady-states in turn when hydrogen and phenol addition are included in the model are given. To obtain similar expressions when hydrogen and/or phenol addition is not present in the model, simply set u h = 0 and/or u p = 0, respectively. For the two-tier phenol model, consider s 0 = x 0 = 0 and solve accordingly.
Equations (14), (15) and (16) result in s 0 = u f , s 1 = u g and s 2 = u h , such that the steady-state is always meaningful.
The Jacobian matrix, J SS1 , can be written as Its six eigenvalues are The condition for stability is λ 1,2,3 < 0. By substituting s 0 , s 1 and s 2 into each of µ 0 , µ 1 and µ 2 , stability is assured when Equations (14) and (15) give s 0 = u f and s 1 = u g . Solving the remaining dimensionless equations gives Here, the Jacobian matrix, J SS2 , then becomes Its eigenvalues are . For stability, it is simply required that λ 1,2 < 0, as for all meaningful values of α, I, x 2 and k c , there is always λ 3,4,5,6 < 0.
Solving Eqs (14), (11), (15) and (16) gives As mentioned previously in B, when numerically checking whether this steady-state is meaningful, each possible value for s 0 is considered in turn. The Jacobian matrix, J SS3 , then becomes Its known eigenvalues are λ 1 = −α+µ 1 −k B , λ 2 = −α+µ 2 −k C , λ 3 = −α. The other three eigenvalues are given by the characteristic polynomial For stability it is required that λ 1,2 < 0 and all three roots from Eq. (C7) have negative real parts. These are checked following the procedure described in B.
By solving Eqs. (11) to (16) for x 0 , x 1 , s 1 and s 2 , and substituting these into equation (14), the following cubic expression for s 0 is obtained Once again, when numerically checking whether this steady-state is meaningful, each possible value for s 0 in turn is checked. The remaining expressions for this steady-state are given by The Jacobian matrix, J SS4 , becomes Its known eigenvalue is λ 1 = −α+µ 2 −k C and the characteristic polynomial for the remaining five eigenvalues is given by For stability it is required that λ 1 < 0 and all five roots from equation (C20) have negative real parts.