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

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 two-tiered feeding chain, populated with a set of parameters that apply to propionate degraders and hydrogenotrophic methanogens (Xu et al., 2011). 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 behaviour 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 (Dolfing and Tiedje, 1986, Dolfing and Tiedje, 1991). 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 multi-species relationships. It should be noted that chlorophenol is the chosen compound for this study due to the availability of data for parameterisation of the model, but the analysis will hold, generally, for any similar food-web demonstrating equivalency.

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 (Batstone et al., 2002). On the other hand, simplified or reduced models have received more attention in process monitoring, control design (Bernard et al., 2001, Gaida et al., 2012), and optimisation (Vavilin et al., 2001). Simplified models have also been applied to determine the global (Benyahia et al., 2012, Shen et al., 2007) or local (Simeonov and Diop, 2010) stability of the system under investigation. Typically, these are two-species models, whereas here a case considering a three-species food web is presented.

Whilst a strictly analytical approach is not possible given the dimensions of the ensuing model, the approach taken here derives analytic expressions for all steady-states, 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.

Model and method

Mechanistic model of a three-tiered food web

The model developed here is based on Anaerobic Digestion Model No. 1 (ADM1) (Batstone et al., 2002). 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, . Chlorophenol, phenol and hydrogen are introduced with an input concentration , , , respectively, and a dilution rate D. The inhibition of hydrogen on the phenol degrader is defined as

The substrate and biomass concentrations then evolve according to the six-dimensional dynamical system of ODEsHere the functions f0, f1 and f2 are defined aswhere and 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 1, and their derivation is presented Appendix 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 phenol COD consumed as hydrogen by .

Table 1

Parameters used in the two-tiered food chain and three-tiered food web models, where d represents days, and and are the Chemical Oxygen Demand (COD) of the biomass and substrate, respectively. The derivation of these parameter values is given in Appendix A.

Parameters	Nominal values	Units
k_m,ch	29	kgCOD_S/kgCOD_X/d
K_S,ch	0.053	kgCOD/m^3
Y_ch	0.019	kgCOD_X/kgCOD_S
k_m,ph	26	kgCOD_S/kgCOD_X/d
K_S,ph	0.302	kgCOD/m^3
Y_ph	0.04	kgCOD_X/kgCOD_S
k_m,H_2	35	kgCOD_S/kgCOD_X/d
K_S,H_2	2.5 ×10^−5	kgCOD/m^3
K_S,H_2,c	1.0 ×10^−6	kgCOD/m^3
Y_H_2	0.06	kgCOD_X/kgCOD_S
k_dec,i	0.02	d^−1
K_I,H_2	3.5 ×10^−6	kgCOD/m^3

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. (2011), after Baltzis and Frederickson (1984), the following dimensionless terms are defined as

With these transformations, the dynamical system of ODEs given in Eqs. (2), (3), (4), (5), (6), (7) reduces to where, for simplicity, the following parameters are introduced

Steady-states

The steady-states of this system of ODEs are obtained by setting (for ) and solving simultaneously. This predicts eight possible steady-states (defined equilibria):

The trivial steady-state where all three populations are washed out .

Only the methanogen population is maintained .

The phenol degraders and methanogens are washed out .

The hydrogenotrophic methanogens are washed out while the chlorophenol and phenol degraders are maintained .

Only the phenol degraders are washed out .

All three populations are present .

Only the phenol degraders are present .

The phenol degraders and methanogens are maintained in the system .

For full chlorophenol mineralisation the only desired operating condition, SS6, must be stable. In a continuous system, full mineralisation occurs when all three organisms are present and are able to continually remove the phenolic compounds. This implies that setting in the two-tiered system and in the three-tiered system, at , would result in the phenolic substrate concentrations decreasing to zero over time.

Stability analysis

The challenge of analytically characterising dynamical systems of ODEs is well understood. Typically, for systems of high dimensions, it is impossible to find explicit and succinct solutions, so one must resort to numerical analysis. For example, the Routh–Hurwitz theorem allows for an explicit analysis of the stability of steady-states, but is intractable beyond five dimensions (May, 1973) 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 Appendix B.

In order that a steady-state be meaningful, all variable concentrations must be non-negative, while Eq. (5) also gives the condition that (or in dimensionless form). 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 (Pimm, 202). The Jacobian for the system of Eqs. (11), (12), (13), (14), (15), (16) corresponds to the 6×6 matrix

Note that for the two-tier model, for which , the Jacobian reduces to a 4×4 matrix. The Jacobian is then evaluated at a given steady-state, denoted , 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 3, 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 , reducing the system of ODEs to four equations. Subsequently, the full model for chlorophenol mineralisation is studied.

Results

Two-tier phenol model

Previous studies of two-tier ecological systems have obtained complete stability analysis, locally and globally, without a maintenance term () (El Hajji et al., 2010, Weedermann et al., 2013). Here, the case including this coefficient for the two-species model fed with phenol is considered by setting , and with . As expected from the model structure, the results are similar to those shown for a propionate-degrading bi-culture with a maintenance term () (Xu et al., 2011), where three stable steady-states emerge: SS1, SS7 and SS8. Fig. 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, . Although numerical analysis of the steady-states was performed here, it is possible to get explicit quadratic functions of the parameter pair (, D) for the steady-state partitions. Using the method described in Xu et al. (2011), the functions for the partitions are (in non-dimensionless form)where where relates to methanogen washout and to phenol degrader washout. These partitions are shown by the green lines in Fig. 1.

Fig. 1

Steady-state diagram for operational parameters D and , in the two-tier phenol model. Green lines between steady-states: explicit functions of the steady-state partitions, after Xu et al. (2011).

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.

Three-tier chlorophenol model

Numerical simulation of the model was performed with , using the parameters listed in Table 1 over a range of operating conditions (α,u). The stable steady-state regions are shown in Fig. 2. Solutions were found for the three steady-states identified (SS1, SS4 and SS6). Definingthen the Jacobians for SS1, SS4 and SS6, respectively, are

Fig. 2

Steady-state diagram for operational parameters D and in the three-tier chlorophenol model ().

A full analysis of the stability for each steady-state is provided in Appendix C. By inspecting the eigenvalues () of the Jacobian matrix for each steady-state, the following conditions for stability are found, in which the term ‘meaningful’ implies non-negativity throughout the proceeding analysis sections:

Meaningful and stable, with .

Meaningful when the roots of the cubic function (C.11) for s0 are positive and the conditions and are met. Additionally, it can be seen from Eq (C.18) that x0 can only be meaningful if . Solving this inequality results in the condition , such that x0 is meaningful provided and . Stability requires that and ().

Meaningful for the following conditions: ;and stable when .

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 observation not seen with the two-tiered models is the existence of SS1 at low dilution rates but increasing substrate input. It is likely that the condition 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. However, mathematically, the phenomenon can be seen as an expected condition when considering the three-tiered model as a perturbation of the two-tiered counterpart, such that the stability of is preserved for and (the conditions giving rise to only stable in the two-tiered model), as long as is small.

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 (). 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, in Eq. (7), and the dimensionless form, in Eq. (16), gives rise to a further three steady-states are defined with this model; SS2, SS3 and SS5. It will be shown that these are stable equilibria observed for the parameter values provided in Table 1.

Following the same approach as for the three-tier chlorophenol model, six Jacobian matrices were derived for each steady-state, and their eigenvalues λ (where ) were found. In this case, the conditions for stability for each steady-state are as follows:

Meaningful and stable providingGiven that, in this model, , then for , the second condition becomes the prevailing factor in determining stability, and can be reduced to , where .

Meaningful when and and stable when . Here it can be noted that the stability of SS2 does not depend on u, and as such will remain fixed for any hydrogen input concentration.

Meaningful given the roots of the quadratic function (C.3) for s0 are positive, and . For steady-state stability, it is necessary that , and for (). Here, the stability of SS3 is dependent on both u and u. The first condition dictates that at higher u, a greater amount of u is also required to maintain stability, whereas the second condition indicates that at a fixed u, increasing u will increase the dilution rate under which SS3 is stable up to a boundary described by a function .

Meaningful given the roots of the cubic function (C.11) for s0 are positive and the conditions , , and are met. Stability is assured when and for (). As u affects s0 through the highly non-linear cubic expression, it is difficult to draw any conclusion about its effect on stability by analytical means.

Meaningful when the conditions , , and are met. Stability of this steady-state is assured when and .

Meaningful when , , , , and . SS6 is stable when . 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 , where is the non-dimensionless form of ). 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:

Fig. 3

Steady-state diagram for operational parameters D and in the three-tier chlorophenol model with addition, (a) ; (b) ; (c) and (d) () (all ).

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 , 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 steady-state, 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 boundary 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 , with the additional hydrogen resource stabilising the methanogen population at higher dilution rates, limited now only by its maximum growth rate (). 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 region 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 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 (). However, unlike the classical case, here the exclusion principle dictates that consumer A () is never excluded and either consumer B () or producer A () is washed from the system. This can be seen in the system simulation shown in Fig. 4, which considers a single operational point (, ), with increasing .

Fig. 4

Evolution of three-tier system with H2 addition, where , and (green), (red), (blue). Initial conditions: , , , , , (all ).

At low hydrogen addition, the concentration is not enough to sustain the methanogen population whilst a small population of dechlorinators are present (). 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 (steady-state 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.

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

When phenol is added to the system the phenol inflow terms, from Eq. (6) and u from the dimensionless equation (15), become positive. An additional two steady-states 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:

Stable for . For high concentrations of phenol, this is equivalent to , where . However, at lower phenol concentrations, the effect of phenol lowers the minimum dilution rate under which SS1 is stable.

Never stable as a contradiction occurs in the terms describing s2 and x2. For these variables to be meaningful, the conditions and must be satisfied, which is invalid (see Eqs. (C.1), (C.2)).

Never stable as a contradiction occurs in the terms describing s2 and x0. With for s2 to be positive. This condition, however, results in x0 always being negative and meaningless (see Eqs. (C.4), (C.5)).

Meaningful given the roots of the cubic function (C.11) for s0 are positive and the conditions , , and are met. Stability is assured when and the roots of the characteristic polynomial (Eq. (C.31)) have negative real parts. In the special case that hydrogen is also added to the system, it should be noted that if , then the condition must be strictly observed. Similarly, if , then u must be greater than s1.

Never stable as without u then x2 can never be positive (see Eq. (C.29)).

Meaningful when , , , , and . With these conditions met, stability of SS6 is guaranteed when the roots of the characteristic polynomial (C.42) have negative real parts. In this case, it can be shown that SS6 is influenced by u.

Meaningful given the roots of the quadratic function for s1 (Eq. (C.50)) are positive and, subsequently, , and . Therefore, stability is assured when , and the roots of the characteristic polynomial (Eq. (C.53)) have negative real parts. It can be seen that the stability of SS7 relies on u, but is also influenced by hydrogen if added to the system.

Stable when the following conditions are satisfied. Firstly, s2 is meaningful given . Stability is assured when s2 is meaningful, , and the roots of the characteristic polynomial (Eq. (C.62)) have negative real parts. In general, for SS8 to be stable, α must be smaller than and, therefore, u should be small.

Again, the operational parameter plots are used to visualise the outputs from the stability analysis, as shown in Fig. 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 three-tiered 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 ) and SS4 (high ), and the boundary of SS1. Of greater interest is the observation that the extent of the desired condition, SS6, is greater for between approximately 0 and , replacing both SS4 and part of SS7.

Fig. 5

Steady-state diagram for operational parameters D and in the three-tier chlorophenol model with phenol addition, (a) and ; (b) and ; (c) and ; (d) and .

Three-tier chlorophenol model with bi-substrate addition

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

Fig. 6

Steady-state diagram for operational parameters and in the three-tier chlorophenol model with (a) and ; (b) and .

An alternative situation under low chlorophenol addition is shown in Fig. 6(b). Here, it is clear that phenol addition at concentrations above provides the necessary conditions for full chlorophenol mineralisation, regardless 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 , 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.

Bistability analysis

Throughout Section 3, numerous regions of bistability were observed between SS3 and SS4 (Fig. 3(a)–(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 understanding. Through analysis of these results, the steady-state 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. Fig. 7 shows the final stable steady-states (SS3 or SS4) achieved for a range of initial conditions. Note that for SS3 or SS4, and does not deviate from this value throughout, which simplifies the analysis required here. When is small, the system is likely to head towards SS3 where washout of occurs. Otherwise, the system achieves SS4. Very similar results were obtained for biologically meaningful substrate initial conditions between and for , and for and and for (). These results suggest that as long as sufficient dechlorinating biomass is present in the system, the most likely outcome is for both and 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, , or within the bistable region once again leads to similar results.

Fig. 7

Phase plot showing final steady-state with initial conditions , , and and parameters and (all .

A similar analysis of the bistable regions for SS3–SS6 and SS1–SS6 give a more dramatic outcome. If any of the three organisms, , or , are initially zero and invariant, SS6 cannot be obtained. If all organisms are simultaneously 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. The analysis highlights not only how the main parameters of the system can impact the final stable state, but also the nature and importance of initial conditions.

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 (Xu et al., 2011, Dubey and Hussain, 2000, Wu et al., 2004). More recently, attempts at analysing more complete systems, such as the anaerobic digestion process, have been performed with some degree of rigour (Bornhöft et al., 2013, Weedermann et al., 2013), 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 (Xu et al., 2011); 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 (), 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. This case confirms that in specific operating regions hypothesis (i), which states that the dechlorinator can (partially) replace the methanogen as the syntroph, may be achieved. Although full chlorophenol mineralisation cannot be obtained in this state (some phenolic compounds remain unconsumed), the results compare well with experimental observations for an equivalent chlorobenzoate degrading system (Dolfing and Tiedje, 1986, Dolfing, 1990). With an abundance of additional hydrogen, however, the methanogen may 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. Nevertheless, the addition of phenol under specific operational conditions can result in the emergence of SS6 at low chlorophenol input conditions. Here, a specific form of competition between the dechlorinator and methanogen occurs (as with SS6 in the standard model), in 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, a further confirmation of hypothesis (i). 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.

Fig. 8

Schematic of observed microbial interactions in the three chlorophenol degradation models at low chlorophenol input concentrations. Top: Three-tier chlorophenol model with hydrogen inhibition of the phenol degrader. Middle: 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.

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.