Mathematical model of anaerobic digestion in a chemostat: effects of syntrophy and inhibition.

Three of the four main stages of anaerobic digestion: acidogenesis, acetogenesis, and methanogenesis are described by a system of differential equations modelling the interaction of microbial populations in a chemostat. The microbes consume and/or produce simple substrates, alcohols and fatty acids, acetic acid, and hydrogen. Acetogenic bacteria and hydrogenotrophic methanogens interact through syntrophy. The model also includes the inhibition of acetoclastic and hydrogenotrophic methanogens due to sensitivity to varying pH-levels. To examine the effects of these interactions and inhibitions, we first study an inhibition-free model and obtain results for global stability using differential inequalities together with conservation laws. For the model with inhibition, we derive conditions for existence, local stability, and bistability of equilibria and present a global stability result. A case study illustrates the effects of inhibition on the regions of stability. Inhibition introduces regions of bistability and stabilizes some equilibria.

Introduction

Anaerobic digestion is a complex naturally occurring process during which organic matter is broken down to biogas and various byproducts in an oxygen-free environment. It is used for waste and wastewater treatment and for production of biogas, which is a mixture consisting of methane, carbon dioxide, and some trace gases.

First scientific experiments involving biogas date back to the eighteenth century and are attributed to Alessandro Volta, who discovered a direct correlation between the amount of organic material used and the amount of flammable gas produced. In the nineteenth century, Louis Pasteur, who coined the term anaerobic to describe the process of anaerobiosis, and his collaborators conducted experiments involving the production of biogas from cattle manure. The first known utilization of biogas generation was recorded in the late nineteenth century in Bombay, India [19], where a waste digester was installed in a leper colony to produce gas used for lighting and later electricity. The first anaerobic digester for wastewater treatment was developed in 1907 by Karl Imhoff (Imhoff tank). In the early twentieth century, several small-scale utilizations of anaerobic digestion for biogas productions were developed in France and England [6].

Spurts in the development of anaerobic digestion facilities for wastewater treatment occurred during the 1970s and in recent years due to rising oil prices and an increased interest in the commercial utilization of anaerobic digestion for its environmental and economic benefits. When anaerobic digestion is used in waste treatment facilities, especially for the treatment of sewage sludge, the biogas is captured before it escapes into the atmosphere. It can then be used as renewable energy either by combusting the gas to produce electrical energy or by extracting the methane and using it as a natural gas fuel.

In both laboratory and large-scale industrial installations, anaerobic digestion appears to be difficult to control. Reactors often experience break-down resulting in little or no biogas production. Since the mid-1970s, a number of mathematical models for anaerobic digestion in a continuous-flow stirred-tank reactor have been developed and used to predict biogas/methane production. These modelling efforts have been two-fold: quantitative and qualitative. A significant number of models exist that have been fit to experimental data and predict biogas production levels (see [1,15,17,23] and the references therein). The comprehensive Anaerobic Digestion Model 1 (ADM1) model developed in [1] is a structured model consisting of over 30 state variables. The complexity of this model makes a qualitative analysis very difficult. Studies concerning the qualitative analysis are less common [3,8,10,13,15,22]. This study contributes to these qualitative results. While the model studied here does not possess the complexity of the model introduced in [1], the general structure of the equations is maintained. The results presented here thus give insight into the behaviour of solutions of the high-dimensional model considered in [1].

Anaerobic digestion is a four-stage process comprised hydrolysis, acidogenesis, acetogenesis, and methanogenesis (Figure 1). During the first stage, complex organic molecules are broken down into simple sugars, fatty acids, and amino acids (monomers). During acidogenesis, acidogenic bacteria convert these monomers into alcohols and volatile fatty acids (VFAs) such as propionic acid, butyric acid, and acetic acid, carbon dioxide, and hydrogen. With the exception of acetic acid, VFAs and alcohols are utilized by acetogenic bacteria and converted to acetic acid as well as carbon dioxide and hydrogen. In the final stage, acetoclastic methanogens convert acetic acid into methane and carbon dioxide, while hydrogenotrophic bacteria convert hydrogen and carbon dioxide to methane. Anaerobic digestion is inhibited by a number of factors. In this study, two of these are incorporated: inhibition due to pH and inhibition due to high partial pressure of hydrogen.

Figure 1.

Schematic view of the four stages of anaerobic digestion.

While all microbes involved in acidogenesis, acetogenesis, and methanogenesis can only survive within a certain pH-range, methanogens are the most sensitive to pH and can only thrive within a pH-range of 6.7–7.2 [7]. If acid production outpaces consumption, the pH of the environment is lowered, inhibiting utilization of acetic acid and thus leading to a further accumulation of acetic acid. Acetoclastic methanogens utilize acetic acid for their growth but are limited by pH. Assuming that the pH-levels are directly related to acetic acid, we describe the growth of acetoclastic methanogens with a non-monotone growth function as proposed in [2,10,16]. The pH inhibition on hydrogenotrophic methanogens is incorporated in the form of inhibition due to high levels of acetic acid rather than of all VFAs. The same assumptions were made in [14,16,20,23,28].

The acetogenic conversion of monomers is only thermodynamically favourable under low partial pressures of hydrogen. The hydrogen transfer between acetogenic bacteria and hydrogenotrophic methanogens is referred to as syntrophy. The inhibitory effect of hydrogen on the growth of acetogenic bacteria is incorporated using a general function with the same properties as a modified Monod function first introduced in [21] and considered subsequently in [1,5,17,23,25].

The goal of this paper is to study the effects of inhibition on the microbes involved in the process. The study combines and extends recent results by Hess and Bernard [10] and Hajji et al. [8], who considered different subprocesses of anaerobic digestion. As in these studies, we assume that the temperature remains constant and consider continuously stirred tanks only.

The paper is organized as follows. In Section 2, we derive a model with nutrient dependencies and inhibitory effects as described above. The resulting system of differential equations describes the acidogenesis, acetogenesis, and methanogenesis stages of anaerobic digestion. We chose not to include hydrolysis as it can be considered as a preliminary phase of the process. An alternative interpretation of the model is that hydrolysis and acidogenesis are considered as a single stage. We first give some elementary results and consider the connection to two submodels which have been considered in the literature [8,10]. In Section 3, we provide a bifurcation result for the preservation of global stability of equilibria. This result is applied to an inhibition-free version of our model, which serves as a baseline for studying the effects of inhibition. We then return to the model introduced in Section 2. After establishing conditions for the existence and local asymptotic stability of each equilibrium, we exclude the possibility of the existence of a periodic orbit due to a Hopf bifurcation. We also show that the system can have several coexisting locally stable equilibria (bistability). Under these conditions global stability of an equilibrium is not possible. We then consider the problem of global stability. Using differential equations and ideas similar to the inhibition-free case, we are able to give a result for the global stability of a positive equilibrium, that is, an equilibrium at which all microorganisms under consideration coexist. Similar ideas can be used to develop global stability results for several of the other equilibria. We conclude with a case study of the impact of inhibition on the system and a comparison of the predictions of our model with the predictions made in [8,10]. Most parameters for this case study are based on [28]. We illustrate our findings with several diagrams corresponding to various inhibition levels. These diagrams show how regions of stability and bistability of equilibria vary with inhibition and demonstrate that inhibition has a significant impact on the long-term survival of the microorganisms and thus biogas production.

Mathematical model

We consider the acidogenesis, acetogenesis, and methanogenesis stages of anaerobic digestion with focus on the concentrations of these nutrient groups: simple substrates S (monomers such as simple sugars, long chain fatty acids, and amino acids), alcohols, and VFAs such as propionic and butyric acids V, acetic acid A, and hydrogen H. Each stage is facilitated by a distinct group of microorganisms as shown in Figure 1: acidogenic bacteria X, acetogenic bacteria X, and two groups of methanogenic archaea: acetoclastic methanogens X and hydrogenotrophic methanogens X, where X is the concentration of biomass in the reactor.

In general, the stoichiometric equation of a biochemical reaction where a nutrient N is transformed into several compounds P and biomass X is of the form

where v, and v are the respective stoichiometric coefficients. These coefficients can be normalized using carbon oxygen demand (COD) units instead, which leads to

where Y are the respective yield coefficients, YMW. Here, MW denotes the molecular weight and σ is the COD equivalent of N or P or biomass. Furthermore, Σ + Y = 1.

Specifically, acidogenesis is described by the equation −S + Y + Y + Y + Y = 0, acetogenesis by −V + Y + Y + Y = 0, and acetoclastic and hydrogenotrophic methanogenesis by −A + YBGBG + Y = 0 and −H + YBGBG + Y = 0, where BG denotes biogas.

The reaction rates are dependent on the concentration of the nutrient that is being decomposed as well as other factors that may inhibit the reaction. The following inhibitions, which were described in the previous section, will be included in the model:

Acetogenesis is inhibited by high partial pressure or abundance of hydrogen.

Acetoclastic methanogenesis can only take place within a narrow pH-range or within a certain range of concentration of acetic acid A.

Hydrogenotrophic methanogenesis is inhibited by high pH-levels or high concentrations of acetic acid A.

Assuming a continuous-flow constant-volume reactor and uniformly mixed substrate, and ignoring biomass decay rates, we obtain the following system of differential equations:

Letting g = Y for α = S, V, A, H, system (1) is equivalent to

where γ = Y/Y and c = Y for α, β represent S, V, A, H, or s, v, a, h, respectively. The meaning of the constants and functions included in Equation (2) is summarized in Table 1.

Table 1.

Interpretation of constants and functions used in Equation (2).

Name	Meaning
D	Dilution rate
S(o)	Concentration of monomers in inflow
cα	Yield coefficients for biomass
gα(·)	Bacterial growth rate
γαβ	Ratio of product yield relative to biomass yield in conversion of α to β due to Xα

Assumptions and basic results

Throughout this paper, we denote the space of vectors x = (x) with x ≥ 0 by ℝo. We write ℝ+ whenever x > 0 for all 1 ≤ i ≤ n. For most of the paper with the exception of Section 3, the bacterial growth rates satisfy the following conditions:

Condition (HS) is satisfied by linear as well as Holling type II functions (Monod growth) [11,12], which will be used for numerical simulations. Condition (HA) is satisfied by Holling type IV functions (Haldane growth). Typical representations of g and g are modified Monod functions of the form g(α, β) = g(α)I(β), where g is a Monod function and I(β) is a positive, decreasing function, or g(α, β) = m/(k + α + μ).

It is straightforward to show that any solution of Equation (2) with non-negative initial conditions has all components non-negative for all t ≥ 0. Let Σ = S + (1/c)X. Then, = D(S( − Σ) and thus Σ(t) converges to S( as t → ∞. Therefore, for all non-negative initial conditions, since both X(t) ≥ 0 and S(t) ≥ 0 for all t ≥ 0, it follows that both X(t) and S(t) are also bounded above. Next, define Σ − (1/c)X Then, = −DΣ. Therefore, Σ → 0 as t → ∞. Since V(t) = Σ(t) + γ(t) − (1/c)X(t) ≥ 0 for all t ≥ 0 and X(t) ≥ 0 is bounded above, it follows that X(t) is also bounded above. Similarly, if we define Σ − (1/c)X and Σ − (1/c)X, then = −DΣ and = −DΣ, and so Σ → 0 and Σ → 0 as t → ∞. Also, since A = −Σ − (1/c)X ≥ 0 and H = −Σ (1/c)X ≥ 0, it follows first that X is bounded above and then that X is bounded above. Since V, A, and H are defined in terms of the sum and differences of bounded quantities, they must also be bounded. This is summarized in the following proposition.

All solutions of Equation (2) with non-negative initial conditions remain non-negative and are bounded and in ℝ8o. In particular, given any ∊ ≥ 0, Δ∊ = {(S, X) ∊ ℝ8o | S + (1/c)X( + ∊, V + (1/c)X + (1/c)X + (1/c)X + ∊} ⊂ ℝ8o is abounded, globally attracting, positively invariant set of Equation (2).

Throughout this paper, conditions for the existence and local and global asymptotic stability of equilibria will be described in terms of break-even concentrations, which are defined as follows:

The assumptions (HVa), (HVb), (HHa), and (HHb) imply that λ(·) and λ(·) are increasing functions.

Relationship to previous work

Bernard et al. [2] considered anaerobic digestion as a two-stage process consisting of (i) hydrolysis and acidogenesis and (ii) methanogenesis and studied the effects of the pH restrictions on acetoclastic methanogens. While the intermediate product is labelled as VFA, during the development of the model all VFA were assumed to behave like pure acetate. Hess and Bernard [10] investigated qualitative properties of this model, which is an invariant subsystem of Equation (2) for X = X = 0:

where g(S) satisfies (HS) and g(A) satisfies (HA). The system studied in [2, 10] considered a modified outflow rate of biomass in the form of −αD with α ∊ [0, 1]. The case α = 1 considered here corresponds to an ideal scenario in which all biomass is part of the liquid flow. In addition, Hess and Bernard [10] included a constant inflow term in the equation for the A. In [2], concentrations of this inflow where significantly lower than those of S( and are ignored here. Since the equations for S and X decouple from the rest of the system, the dynamics of Equation (3) is equivalent to that of the model considered in [10] for α = 1.

System (3) may have up to four equilibria of the form (S, X),

where X = c(S( – λ) and A( = γ(S( – λ). Under the assumption that the break-even concentrations λ and λ exist, the local and global stability properties of all equilibria of Equation (3) are summarized in Table 2 and all bifurcations are shown in Figure 2 (see [10,27] for more details). For a fixed inflow rate D, we find that there is an optimal range of S( in which the unique interior equilibrium HBE1 is globally asymptotically stable, thus guaranteeing coexistence of the microbes and steady biogas generation. Should S( be too large, S( > λ2/(γ), Equation (3) exhibits the bistability of two equilibria, HBE and HBE1. HBE corresponds to acid accumulation in the system and process failure, and HBE1 corresponds to stable biogas production.

Table 2.

Conditions for existence and for local or global stability of equilibria in system (3).

Equilibria	Existence	Global stability
HBEw	Always	S(o) < λS
HBEo	S(o) > λS
HBEA1
HBEA2		Unstable

Note: When S( > λ2, equilibria HBE and HBE1 are locally asymptotically stable.

Figure 2.

Schematic view of the bifurcation diagram of A versus S( in Equation (3). The A-coordinate of stable (unstable) equilibria is represented as solid (dashed) lines. A transcritical bifurcation appears when two equilibria exchange their stability.

Haiji et al. [8] considered the acetogenic and hydrogenotrophic methanogenic phases of anaerobic digestion. In this study, parts of the overall process have been isolated in a first effort to study the effect of hydrogen inhibition on acetogenic methanogenesis. The model is somewhat problematic as acetogenesis cannot take place without prior acidogenesis, which produces the VFAs that are the growth-limiting nutrient for acetogenic bacteria. If acidogenesis takes place at a constant rate leading to a constant inflow of VFA, then hydrogen is produced along with VFA, meaning that an inflow of hydrogen should be assumed along with the inflow of VFA. The system in [8] is equivalent to a subsystem of Equation (2) with the settings described below.

In system (2), we set X = 0 and γ = 0 and assume that (S, X) approaches an interior equilibrium (λ). We also suppose that A approaches equilibrium A( = γ(S( – λ), which is attained by A in the reduced equation (λ)X. Then, Equation (2) reduces to

where V( = γ(S( − λ).

Based on the definitions of break-even concentrations given in (BV) and (BH), let λ(0), λ(A(), and λ∗(H∗), where H∗ = γ(V( – λ(H∗)). Equation (4) has three equilibria of the form (V, X):

In HMHE(V( – λ(λ)) and X[γ(V( – λ(λ)) – λ(]. The conditions for existence and stability are listed in Table 3. The existence conditions follow immediately from requiring that the equilibria lie in ℝo4. Local stability conditions can be derived from studying the eigenvalues of the equilibria. The proof that local stability implies global stability is given in Appendix 1.

Table 3.

Conditions for existence and for global stability of equilibria in system (4).

Equilibria	Existence	Global stability
HMHEo	V(o) > 0	V(o) < λV
HMHEV	V(o) > λV
HMHEVH		Whenever it exists

System (4) does not exhibit any inhibition-induced bistabilities or limit cycles. When V( > λ(λ) + λ(/γ or equivalently, S( > λ + (1/γ)(λ(λ) + λ(/γ), a single interior equilibrium HMHE is globally asymptotically stable provided that the inflow concentration is sufficiently large and the inhibition of acetogenic methanogenesis due to hydrogen is not too strong. In Figure 3(a), we see that for strong inhibition, an interior equilibrium does not exist.

Figure 3.

Bifurcation diagrams of V versus S( in Equation (4) for D = 0.2 using the growth functions g(V, H) = m(k + V + μ)−-1 and g(H, A) = m (k)−-1 with μ = 1, and (a) μ = 1 (λ(λ)) and (b) μ = 0.3 (λ(λ)); other used parameter values are provided in Appendix 2. The V-coordinate of stable (unstable) equilibria is represented as solid (dashed) lines. A transcritical bifurcation occurs when two equilibria exchange their stability.

Analysis of Equation (2) without inhibition

Without the inclusion of inhibitory effects, that is, assuming that g is monotone increasing with g(0) = 0, g(V) corresponds to g(V, 0), and g(H) to g(H, 0), Equation (2) becomes

For ℓ = S, V, A, and H, the functions g(ℓ) are such that g(0) = 0 and g′(ℓ) > 0. The inhibition-free break-even concentrations, denoted by λ, are the quantities that satisfy the equations

which correspond to the previously defined break-even concentrations defined in (BS), (BA) with λ1, (BV) with λ(0), and (BH) with λ(0). Note that here λ2 = ∞.

Remark 1 The equations for S and X evolve independently of the other variables and represent a basic chemostat. It is well known [4,24] that S( < λ implies that (S(, 0) is the only equilibrium on the boundary of ℝo2 and is globally asymptotically stable. There is no equilibrium in ℝ+2. If S( > λ, equilibrium (λ), where X(S( – λ), exists in ℝ+2. This positive equilibrium is globally asymptotically stable, while the equilibrium on the boundary of ℝo2 is unstable.

If (S, X) attains a positive equilibrium, the term g(S)X is near constant for large values of t thus leading to a near constant nutrient inflow in the equations for V, A, and H. Similarly, if (V, X) is close to a positive (stable) equilibrium, the term g(V)X in the equations for A and H can be interpreted as a near constant nutrient inflow. Motivated by this observation, we introduce the following equilibrium nutrient inflow rates:

System (5) possesses at most nine equilibria, which are listed in Table 4. Among the equilibria there is only one positive equilibrium in ℝ+8. If this interior equilibrium exists, then all other equilibria in ℝo8 are unstable and thus one might suspect that the interior equilibrium is globally stable. This is rather intuitive from the structure of Equation (5). The following lemma describes the dynamics of systems with the type of structure we encounter here.

Table 4.

Equilibria of system (5).

Equilibrium	Coordinates (S, XS, V, XV, A, XA, H, XH)
Ew	(S(o), 0, 0, 0, 0, 0, 0, 0)
Eo	(λS, XS, V(o), 0, A(o), 0, H(o), 0)
EH	(λS, XS∗, V(o), 0, A(o), 0, λH, ch(H(o) − λH))
EA	(λS, XS∗, V(o), 0, λA, ca(A(o) − λA), H(o), 0)
EAH	(λS, XS∗, V(o), 0, λA, ca(A(o) − λA), λH, ch(H(o) − λH))
EV	(λS, XS∗, λV, cv(V(o) − λV), A, 0, H, 0)
EVH	(λS, XS∗, λV, cv(V(o) − λV), A, 0, λH, ch(H − λH))
EVA	(λS, XS∗, λV, cv(V(o) − λV), λA, ca(A − λA), H, 0)
E∗	(λS, XS∗, λV, cv(V(o) − λV), λA, ca(A − λA), λH, ch(H − λH))

Notes: X(S( − λ). The break-even concentrations (λ, where ℓ = S, V, A, or H) and V(, A(, (, and H are defined in Equations (6) and (7).

Suppose F: ℝ → ℝ1 and = F(Y) has a globally asymptotically stable equilibrium Y∗. Assume g: [0, ∞) → [0, ∞) is C1, increasing, and g(0) = 0. For fixed values of D > 0 and c > 0, let λ be the solution of the equation g(s) = D. Furthermore, let f: ℝ → [0, ∞) be C1. Consider the system

with s(0) > 0 and x(0) > 0. The following statements hold:

If λ > (1/D)f (Y∗), then all solutions approach (Y∗, (1/D)f(Y∗), 0).

If 0 < λ < (1/D)f (Y∗), all solutions approach (Y∗, λ, x∗), where x∗ = c(f(Y∗)/D − λ).

Let Σ(t) = s(t) + (1/c)x(t). Then, = −DΣ + f(Y) and Σ(t) = Σ(0) e− + e− ∫ e(Y(s)) ds. Using integration by parts,

To see that , note that necessarily F(Y∗) = 0, which gives that for any ∊ > 0 there is a T > 0 so that for all t > T, |F(Y(t))| < ∊ and |f′(Y(t))| < L for some L > 0. So,

is such that . Thus,

Letting ∊ → 0 gives that .

Since , it suffices to show that under condition (1), lim x(t) = 0, while under condition (2), lim x(t) = x∗. For any ∊ > 0, there is a T∗ > 0 so that for all t > T∗,

and so, since g is increasing in s, for all t > T∗

Consider with u∊ (0) ≥ 0. This equation has two equilibria, u = 0 and u∊∗ = c(f(Y∗)/D − λ + ∊). The local stability of each equilibrium is determined by

If condition (1) holds, then for ∊ sufficiently small, λ > f(Y∗)/D + ∊, which implies u∗∊ < 0 and h(u) < 0. That is, for all solutions u∊(t), lim u∊(t) = u = 0. If condition (2) holds, then h(u) > 0, u∗∊ > 0, and h(u∗∊) < 0. Thus, lim u(t) = u∗.

Similarly, let ∊ = −v∊ (D − g(f (Y∗)/D with v∊(0) ≥ 0. The two equilibria are v = 0 and v∊∗ = c(f (Y∗)/D − λ − ∊). The local stability of each equilibrium is determined by

If condition (1) holds, then v∊∗ < 0 and h(v) < 0. That is, for all solutions v∊(t), lim u∊(t) = v0 = 0. If condition (2) holds, then for ∊ sufficiently small, λ < f(Y∗)/D − ∊ and so, h(v) > 0, v∊∗ > 0, and h(v∊∗) < 0. Thus, lim v∊(t) = v∊∗.

From Equation (8) we can now conclude that if condition (1) holds, then limx(t) = 0, while condition (2) implies that v∊∗ ≤ lim inf x(t) ≤ lim sup x(t) ≤ u∊∗. Letting ∊ → 0 gives that lim x(t) = x∗.

From Lemma 3.1, one can prove results for the global stability of all equilibria of Equation (5). Here, we outline the proof for the global stability of E∗. The conditions for global stability of other equilibria of Equation (5) can be derived similarly and are summarized in Table 5.

Table 5.

Conditions for existence and global stability of all equilibria of Equation (5) in ℝo8 (also see Table 4).

Equilibria	Conditions for existence	Conditions for global stability
Ew	Always	λS > S(o)
Eo	λS < S(o)	λV > V(o), λA > A(o), and λH > H(o)
EH	λH < H(o)	λV > V(o) and λA > A(o)
EA	λA < A(o)	λV > V(o) and λH > H(o)
EAH	λA < A(o) and λH < H(o)	λV > V(o)
EV	λV < V(o)	λA > A and λH > H
EVH	λV < V(o) and λH < H	λA > A
EVA	λV < V(o) and λA < A	λH > H
E∗	λV < V(o), λA < A, and λH < H	Whenever it exists

Suppose λ(, λ (5) has a unique interior equilibrium E+8.

If λ(, then S( > λ and thus Y = (λS, X∗) is a globally stable equilibrium of the system corresponding to Y = (S, X) (see Remark 1). Applying Lemma 3.1 with g(V) = g(V), c = c, and f(S, X) = γ(S)X implies that the solutions corresponding to Y = (S, X) approach Y∗ = (λ(V( – λ)). One can now continue to proceed in this fashion to establish the global stability of E

From a practical point of view, the main control parameters are D and S(. For parameters as in [1,28], the regions of global stability of Equation (5) are shown in Figure 4. Exact parameter values are listed in Appendix 2.

Figure 4.

Regions of global stability of equilibria E, and E∗ in system (5) varying dilution rate D and the concentration of inflowing monomers S(. All parameter values used are summarized in Table A1. For this parameter set regions exist in which the other equilibria E, and E lie in ℝo8 but are unstable.

Table A1.

Parameter values based on [28] used for all numerical calculations in this study.

Parameter (unit)	Wett et al. [28]
ms (day−1)	3
ks (g COD/L)	0.5
mv (day−1)	0.5
kv (g COD/L)	0.1
m[a] (day−1)	0.75
k[a] (g COD/L)	0.15
mh (day−1)	1.998
kh (g COD/L)	10−3a
cs[b]	0.1
cv[b]	0.05[c]
c[a][b]	0.025
ch[b]	0.06
γsv	2.0646[d]
γsa	5.2254
γsh	1.71
γva	13.3[e]
γvh	5.7[e]

This was reported as 10−6, changed to 10−3 to ensure numeric accuracy.

Yield coefficients are dimensionless due to the use of COD equivalent units. In the literature, they are usually reported as (g COD BM/g COD) or are interpreted as such.

Computed as the average of c given for propionate (0.04) and butyrate (0.06).

In [28], ppro, + f = 0.2294 was used to ensure Σ = 1. Then, γ = (1 − c)p.

p = 0.7 approximated and averaged based on f = 0.57 ≈ 0.6 and f = 0.8, and then Σp = 1 gives p = 0.3. Then, γ = (1 − c)p and γ = (1 − c)p.

Effects of inhibition

In this section, we study the effect of various inhibitions on acidogenesis, acetogenesis, and methanogenesis. We assume that the conditions (HS, HA, HV, HH) hold.

As before, (S, X) can be decoupled from the remainder of the system and forms a basic single species chemostat. Since X = 0 implies that V, A, S approach 0, and thus X → 0, we will assume that λ( In this case, limS(t) = λ and lim X(t) = c(S( – λ). For simplicity, we will assume S = λ and X(S( – λ). With this assumption, system (2) reduces to

where V(, A(, and H( were defined in Equation (7). From Proposition 2.1, we can immediately deduce that all solutions of Equation (9) corresponding to positive initial conditions are bounded.

Existence of equilibria and local stability

Based on the definition of inhibition-dependent break-even concentrations in (BA), (BV), and (BH), we introduce the following:

Since 0 < λ1 < λ2, it follows that 0 < λ1 < λ2 and λ1 < λ2.

From the equations for X, and X, we see that Equation (9) has equilibria with positive values of X, or X if and only if the equations g(V, H) = D, g(A) = D, or g(H, A) = D, respectively, have positive solutions that satisfy certain conditions. It is straightforward to establish that system (9) has 12 potential equilibria as listed in Table 6. Conditions for the existence of most of these equilibria in ℝo6 can also be derived directly (Table 7). For instance, ε lies in ℝ6o with X > 0 if λ( or, equivalently, if λ + λ/γ < S(. For two of the equilibria, ε and ε, conditions of existence in ℝo6 will be derived separately. ε is such that X = X = 0. From here, it can be seen that V = in ε must satisfy the equation

Table 6.

Possible equilibria of Equation (9) in ℝo6.

Equilibria	Coordinates (V, XV, A, XA, H, XH)
εo	(V(o), 0, A(o), 0, H(o), 0)
εH	(V(o), 0, A(o), 0, λoH, ch (H(o) − λHo)
εAi	(V(o), 0, λAi, ca (A(o) − λiA), H(o), 0)
εAHi	(V(o), 0, λiA, ca (A(o) − λiA), λiH, ch (H(o) − λHi)
εV	(, cv (V(o) − ), AV(o) − γvacv, 0, HV(o) − γvhcv, 0)
εVH	(, cv (V(o) − ), AV(o) − γvacv, 0, , ch (HV(o) − − γvhcv)
εiVA	(, cv (V(o) − ), λiA, ca (AV(o) − λAi − γvacv), HV(o) − γvhcv, 0)
εi∗	(λiV, cv (V(o) − λVi), λAi, ca (Ai − λiA), λiH, ch (Hi − λHi))

Notes: A( = A( + γ( and H( + γ(, i = 1, 2. is the solution of Equation (10) and () is the solution of Equation (11), A( − γ, and H( −γ.

Table 7.

Conditions for existence and local stability of equilibria of Equation (9) in ℝo6 (also see Table 6).

Equilibria	Conditions for existence	Conditions for local stability
εo		λVo > V(o), A(o) < λA1 or λ2A < A(o), and λHo > H(o)
εH	λHo < H(o)	λV(λHo) > V(o), A(o) < λA1, or λA2 < A(o)
εA1	λA1 < A(o)	λoV > V(o) and λ1H > H(o)
εA2	λA2 < A(o)	Always unstable
εAH1	λA1 < A(o) and λH1 < H(o)	λV1 > V(o)
εAH2	λA2 < A(o) and λH2 < H(o)	Always unstable
εV	λVo < V(o)	Â < λA1 or λA2 < Â, and Ĥ < λH(Â)
εVH	λV(λHo) < V(o), < H(o)V − γvhcv	Ă < λA1 or λ2A < Ă
εVA1	λVo < V(o) and λA1 < Â	λ1H > Ĥ
εVA2	λVo < V(o) and λA2 < Â	Always unstable
ε∗1	λV1 < V(o), λ1A A1, and λ1H < H1	Whenever it exists
ε∗2	λV < V(o), λ 2A < A2, and λH < H2	Always unstable

Notes: for i = 1, 2 with A2 < A1 and H2 < H1.

that is, is the solution to the implicit equation (H( + γ (V( − )). For X to be positive, it is necessary that < V(, which implies that the H-component in E is greater than H(. This in return implies that > λ.

Equation (10) has a unique solution ∊ (λ() if and only if λ(

If λ < V(, then g(V(, H() > D, and (H2) implies that there exists ξ > H( such that g(V(, ξ) = D. Thus, for every H in [H(, ξ], λ(H) ∊ [λ(] with λ(ξ) = V( Let . Note that (d/dH) λ (H) > 0 and thus f′(H) > 0. From f(H() < D and f(ξ) > D, we can deduce that there exists a unique Ĥ ∊ (H(, ξ) with f(Ĥ) = D. (Ĥ) is the desired solution.

On the other hand, implies λ(

In the equilibrium ε, only X = 0, while X > 0 and X > 0. Consequently, and in ε must be such that

or

Assume λ (λ), and λ(A( + γ() exist. Then, Equation (11) has a unique solution () with λ(λ) < ( and λ(A( + γ(V( – λ(λ))) if and only if λ(λ) < V(

If () is a solution of Equation (11) with (, then λ This implies that λ(λ) < thus proving λ(λ) < V( λ(λ) < also implies (A( + γ(V( – λ(λ))).

Assume λ(λ) < V( and let . Then, f is defined and differentiable for V ∊ (λ(λ), V() with f′(V) > 0. Furthermore, f(V() = g(V(, λ) > D and f(λ(λ)) < g(λ(λ), λ) = D. This implies that here is a unique value ∊ (λ(λ), V() such that is the desired value of H.

The local stability of the various equilibria is determined by the eigenvalues of the Jacobian matrix of Equation (9) evaluated at each equilibrium. The Jacobian matrix is of the form

where J11 = −D − (1/c)∂12 = −(1/c)g15 = −(1/c)∂21 = ∂22 = −D + g25 = ∂31 = γ32 = γ33 = −D − (1/c)g′34 = −(1/c) g35 = γ43 = g′44 = −D + g51 = γ52 = γ53 = −(1/c)∂55 = −D + γ − (1/c)∂56 = −(1/c)g65 = ∂, and J66 = −D + g.

The eigenvalues of most equilibria with the exception of ε and ε can be computed directly to derive conditions for local asymptotic stability as summarized in Table 7. The local stability of ε and ε will de determined using the Routh–Hurwitz criterion. We first consider ε and denote its A coordinate by Ă = AV( − γ. The characteristic polynomial corresponding to ε is (z + D)2(z – g(A) + D)p(z), where p(z) = z3 + a1z2 + a2z + a3 with coefficients

It is not hard to see that a1a2 − a3 > 0. Thus, by the Routh–Hurwitz criterion, all zeros of p(z) have negative real part. This implies that this equilibrium is locally asymptotically stable whenever g(Ă) < D, or equivalently, if Ă < λ1 or Ă > λ2.

Suppose ε (i = 1, 2) exists and lies in ℝ6∗1 is locally asymptotically stable, while ε2 is unstable.

The eigenvalues corresponding to ε are −D (double) and the zeros of the polynomial q(z) = z4 + α1z3 + α1z3 + α2z2 + α3z + α4, where

We will apply the Routh–Hurwitz criterion to q(z), which states that all solutions of q(z) = 0 have negative real part if and only if (1) α1 > 0, (2) α1α2 > α3, (3) α1α2α3 > α12α4 + α32, and (4) α4 > 0. If α3 > 0, then (3) implies (2). If α3 < 0, then necessarily α2 < 0. In this case, the Descartes change of sign rule guarantees the existence of two positive roots. Hence, we must require (5) α3 > 0 and (3) can only hold if also α2 > 0.

For i = 1, ∂(λ1) > 0 and thus α > 0 for k = 1, 2, 3, 4. To verify condition (3), let

Then, the coefficients of q(z) can be written as α1 = D + A1, α2 = DA1 + A2 + A3, α3 = DA3 + A4, and α4 = DA4. Condition (3) is equivalent to α1α2α3 − α12 α4 − α32 > 0. Direct calculations result in α1α2α3 − α12α4 − α32 = δ3D3 + δ2D2 + δ1D + δ0 = r(D), where

Thus, r(D) has no positive root and r(0) > 0. Consequently, condition (3) holds for all D ≥ 0. For i = 2, ∂(λ2) < 0 implies α4 < 0, which violates condition (4).

Under certain conditions, locally asymptotically stable equilibria of Equation (9) can coexist. The different bistabilities are summarized in the following corollary. The first statement is equivalent to the bistability found for system (3).

If λ(, λ1 < λ2 < A( and H( < λ1 ℝo6, and ε1 are locally asymptotically stable.

If λ1 > V(, λ1 < λ2 < A( and λ1 < H( < λ1 lie in ℝ6o and are locally asymptotically stable.

If λV1 < V(λV2 < A( and λ1 < H( < λ1 lie in ℝ6o and are locally asymptotically stable.

If λ1 > V(, λ1 < λ( and λ(, then εo6, and ε1 are locally asymptotically stable.

If λ1 < V( < λ(λ), λ2 < A( and λ(, then ε∗1 lie in ℝo6 and are locally asymptotically stable.

If λ(, λ1 < λ2 A and λ ℝ6o, and ε

If λ(, λ1 < λ2 < Â, and λ1 < Ĥ < λ(Â), then ε1 lie in ℝo6 and are locally asymptotically stable.

If λ(λ) < V(, λ( – γ1, and ε2 all exist in ℝ6o, and ε∗1 are locally asymptotically stable.

The functions λ(·) and λ(·) are increasing, and thus λ12( implies λ1 < λo, which then implies λ1 < λ2 < λ(λ). From here we see that (i) and (iv) follows directly from the conditions listed in Table 7. In (ii), λ1 < H( gives λ1 < λ and thus all conditions for existence and stability of ε and ε1 hold. In (iii) and (v), the conditions ensure the existence and stability of ε and ε, respectively. λ1 < V( guarantees that A( < A1 and H( < H1, which is sufficient to verify all conditions for the existence of ε1.

Similarly, λ12 implies that λ1 < λ2 < λ(Â), which together with the conditions listed in (vi) results in the bistability of ε and ε1. To see (vii), note that λ1 < Ĥ gives λ1 < λ(Ĥ) = . Therefore, Ĥ < H1 and  < A1. This means that λ1 < H1 and λ < A1 must hold. Since ∊ (λ(), we also see that λ1( and thus ε∗1 lies in ℝ+6.

Finally, to show (viii), first assume λ2 This implies λ1 < λ2 < , which gives λ1 < λ( From the definition of A and A, we now see that Ă < A2 < A1, and similarly, H( − γ2 < H1. Thus, all conditions for the existence and local stability of ε and ε∗1 hold.

No other bistabilities are possible. From Table 7, we see directly that bistabilities of ε and ε or ε or ε1 are not possible. Similarly, ε1 and ε or ε1 or ε1 cannot coexist and both be stable. The same is true for ε1 and ε1. Most other bistabilities can be excluded based on the monotonicity of the break-even concentration functions.

It is somewhat less straightforward to exclude the bistability of (1) ε1 and ε∗1, (2) ε and ε, and (3) ε and ε. To exclude these bistabilities, first recall that = λ(Ĥ) ∊ (λ() whenever λ( The stability of ε1 requires Ĥ < λ1 and thus 1. This means that Ĥ > H1 and so the condition λ1 < H1 in the existence of ε∗1 cannot hold. Therefore, (1) is not possible.

The equilibrium ε lies in ℝ6o provided that λ(λ) < ( and ( – γ So,

Both ε and ε1 only exist if λ < < V( Our assumptions (HVa) and (HVb) imply that the function f(V) = g(V, H( − γ) is increasing. Since satisfies Equation (10), we have that f() = D. Therefore, f() > f() and so . The definition of Ĥ and also gives that g() = D = g(). Thus, > implies g() < g(), or . From the definitions of  and Ă, we immediately find that > implies  < Ă. The stability of ε requires Ĥ < λ(Â). Since λ(Â) = H and λ(·) is increasing, this contradicts that < Ĥ. This excludes any bistability involving ε and ε.

Finally, to exclude (3), note that if ε1 exists, then  < Ă implies that Ă < λ1 cannot hold. Thus, if ε is stable in ℝ6o, then λ1 < λ2 < Ă. This gives λ1 < λ2 < λ(Ă). Stability of ε1 holds if Ĥ < λ1. We see that mutual stability of ε and ε1 again contradicts < Ĥ.

There are no parameter values for which equilibria of Equation (9) possess purely imaginary eigenvalues, thus eliminating the possibility of limit cycles due to a Hopf bifurcation. No periodic solutions have been observed in numerical simulations.

A global stability result

We will use differential inequalities to derive a result for the global stability of ε∗1 in system (9). For basic chemostats with a non-monotone growth function such as g, the following global behaviour is known [4].

Assume that c > 0, d > 0, and f: ℝo → ℝo is continuously differentiable. If the equation f(s) = d has two solutions λ = 1, 2, such that 0 < λ1 < λ2 and λ1 < s2, then all solutions of the system ⋅ = d(s) – (1/c)f(s)x, (s)x with positive initial conditions approach (λ1, c(s1)).

Our method of deriving global stability of ε∗1 exploits the unique structure of Equation (9), which consists of several coupled simple chemostats. We will first derive inequalities for the solutions of these simple chemostats and then use differential inequalities to obtain a global stability result.

Let g be a function satisfying (HVa) and (HVb) and assume that for 0 ≤ h1 < h2, λ(h1) < λ(h2) < s(h) is such that g(λ(h), h) = d. Assume that s(t) ≥ 0, x(t) ≥ 0, lim(s(t) + x(t)/c) = s(t) and x(t) satisfy the following differential inequalities:

Then, for any ∊ > 0, there exist functions μ(∊) and v(∊) such that lim∊→0 μ(∊) = lim∊→0 v(∊) = 0 and for t sufficiently large, λ(h1) − v(∊) ≤ s(t) ≤ λ(h2) + v(∊) and c(s(h2)) − μ(∊) ≤ x(t) ≤ c(s(h1)) + μ(∊).

Since lim (s(t) + x(t)/c) = s, for all ∊ > 0 and t sufficiently large, s(t)/c − ∊ < s(t) < s(t)/c + ∊. From (HVa), it follows that g(s(t)/c − ∊, h2) < g(s, h2) and g(s, h1) < g(s(t)/c + ∊, h1). Thus, for sufficiently large t,

If lim inf x(t) < c(s(h2) − ∊), and hence s(t)/c − ∊ > λ(h2) for all sufficiently large t, then from the left-hand inequality, x(t) → ∞, as t → ∞, a contradiction. If lim sup x(t) > c(s(h1) + ∊) > 0, and hence s(t)/c − ∊ < λ(h1) for all sufficiently large t, then from the right-hand inequality x(t) → 0, as t → ∞, a contradiction. Therefore,

for μ(∊) = c∊. The inequality for s(t) follows directly from lim (s(t) + x(t)/c) = s.

From the first two equations of Equation (9), , and hence solutions of Equation (9) satisfy lim(V(t) + X(t/c) = V(. Since we are assuming that in system (9), S = λ( and X(S( − λ), from Proposition 2.1 it follows that given any ∊ > 0, for all sufficiently large t, 0 ≤ H(t) + X(t)/c( + γ(t) + ∊. Therefore, if we define , and let Hmax = H0max, then for all sufficiently large t, 0 ≤ H(t) ≤ Hmax and Hmax → Hmax as ∊ → 0.

Define the parameters

Suppose λ(Hmax) < V(, λ1 < A2, and λ1 < H = 1, 2. Then, all solutions of Equation (9) with positive initial conditions approach ε∗1.

Since λ(Hmax) < V(, it follows that for all ∊ > 0 sufficiently small, λ(Hmax) < V( Since also for all t sufficiently large, 0 ≤ H(t) ≤ Hmax, it follows that

Since λ(0) < λ(H∊max) < V(, by Lemma 4.6,

For sufficiently small ∊ > 0, define

Then, α ≤ α2.

Note that A/D = A, and so for sufficiently small ∊, λ < A( + γ/D < λ2 holds for i = 1, 2. Thus, for t sufficiently large, we obtain the differential inequality

By Lemma 4.5 and the theory of differential inequalities (e.g. see [24], Theorem B.1), lim A(t) = λ1 and c(A + γvaα∊1/D − λ1) ≤ lim inf X(t) ≤ lim sup X(t) ≤ c(A + γ/D − λ1).

As a result, we can write inequalities for H and X. Again, for sufficiently small ∊ > 0, define H = H( + γ/D, where i = 1, 2. Then, since by hypothesis λ1 < H = 1, 2, it follows that λ1 < H. For 0 < η < λ1, recalling that α1 ≤ gv(V, H)Xv ≤ α∊1, we obtain

Specifically, for η = 0,

This implies that lim H(t) = λ1 and c(H1 − λ1) ≤ lim inf X(t) ≤ lim sup X(t) ≤ c(H2 – λ1). Thus, for η sufficiently small, there exists a function η(η) such that |H(t) – λ1| < η(η), where lim η(η) = 0. We can now write the inequalities

For η sufficiently small, λ1 + η(η) < Hmax implying that λ(λ1 ≤ V( ≤ λ(λ1 + η(ηa)) + V(ε) andc(V( − λ(λ + η(η))) − μ(∊). Letting η and ∊ approach 0 now gives lim V(t) = λ1 and lim V(t) = λ1 and lim X(t) = c(V − λ1). The convergence of X and X can be deduced in a similar fashion in conjunction with Lemmas 3.1 and 4.5.

A case study

The purpose of this section is to illustrate how inhibition affects the regions of stability and generates regions of bistability. All computations were completed assuming Monod growth for acidogenic bacteria and Haldane-type growth for acetoclastic methanogens. The growth of acetogenic bacteria and hydrogenotrophic methanogens is described by the modified Monod function specified below:

In Equation (12), m denotes the maximum growth rate and k the half-saturation constant for acidogenic bacterial growth. In Equation (13), m denotes the maximum growth rate, k the half-saturation constant, and k a coefficient describing the inhibition due to A on the growth of acetogenic bacteria X. In the growth functions g and g given in Equation (14), m and m represent the maximum growth rates, k and k the half-saturation constants, and μ and μ the inhibition factors for the growth of acetogenic bacteria and hydrogenotrophic methanogens due to hydrogen and acid, respectively. Parameter values other than k are based on [28] and are listed in Appendix 2. The paper by Wett et al. [28] is a study based on the ADM1 model developed in [1]. Most of the parameter values correspond to standard values for the conversion of manure to biogas listed in [1].

For Equations (12)–(14), the conditions for existence and local stability of each equilibrium given in Table 7 can be expressed in terms of D and S( and correspond to various regions in the (D, S()-plane as summarized in Appendix 3. The break-even concentrations can be computed explicitly and are given in Table 8.

Table 8.

Key break-even concentrations for Equation (2).

Break-even concentrations	Domain for D

Note: Note that for k1 → k(m) and λ2 → ∞.

From Equation (13) we see that g approaches a monotone function as k → ∞. If μ = 0, then for large values of k the equations in system (2) approach the equations for the inhibition-free system (5). However, even though when k is relatively large, the regions of stability for small values of S( are almost identical to the inhibition-free graph shown in Figure 4, the stability regions differ significantly from those in the inhibition-free case for relatively large values of S( (see Figure 5 where k = 100 with (a) 0 < S(0) ≤ 3 and (b) 0 < S( ≤ 150). In particular, when k > 0, it is possible to observe a region of stability of ε and regions of bistability of ε and ε∗1 and of ε and ε1 (Figures 5(b) and 6(a)). The regions of bistability of ε and ε∗1 and of ε and ε1 both persist for very large values of k, but the region of stability of ε seems to disappear. From the analytic curves for the boundaries of the regions of local stability given in Appendix 3, we can immediately derive that if μ = 0, the bistability of ε and ε1 occurs when S( > max {f91, f92}. When k → ∞, f92 → ∞ thus requiring very large values of S( to observe this bistability as illustrated in Figure 5. On the other hand, smaller values of k alter the size and prevalence of stability regions as shown, for example, in Figure 6.

Figure 5.

Regions of local stability for k = 100 and μ = 0. For small values of S(, the regions of stability are almost identical to the ones in Figure 4. The effect of the inhibition can only be seen for large values of S( where we see bistabilities of ε and ε∗1 and of ε1 and ε.

Figure 6.

Regions of local stability for (a) k = 1 and (b) k = 0.1 with μ = 0. Larger values of D give additional regions similar to those in Figure 5. Strong inhibition, that is, k = 0.1, has a significant effect on the size of the region of stability of ε1 and other equilibria.

While increasing μ but leaving μ = 0 does not introduce new regions of local stability and bistability, increasing μ has a significant effect as shown in Figure 7. For μ = 0, increasing μ decreases the region of stability of ε and leads to the introduction of regions of local stability of ε and ε1 and bistability of ε and ε∗1, ε1 and ε, and ε1 and ε as shown in Figure 7(a). Increasing μ further eliminates the bistability of ε1 and ε but introduces regions of local stability of ε1 and bistability of ε1 and ε (Figure 7(b)).

Figure 7.

Regions of local stability for k = 1 and μ = 0 with (a) μ = 2.5 and (b) μ = 8. Increasing μ introduces regions of stability of ε1 and bistability of ε and ε1. In (b), additional small regions of stability of ε1 and ε exist for 0 < S( < 0.5, and a region of bistability of ε and ε∗1 for S( > 12.

Finally, if k > 0, and μ > 0, two more bistabilities corresponding to (iii) and (v) of Corollary 4.4 can be observed. When k = 1, μ = 2.5, and μ = 0.1, we see that ε and ε∗1 can be bistable (Figure 8(a)). This means that depending on the initial conditions, one either sees the coexistence of all microorganisms or the survival of only X, the acidogenic bacteria, and no gas production will occur. Figure 8(b) illustrates the possibility bistability of ε and ε1 for k = 1, μ = 1, and μ = 1.

Figure 8.

Regions of local stability for (a) k = 1, μ = 2.5, and μ = 0.1 and (b) k = 1, μ = 1, and μ = 1. A region of stability of ε exists for larger values of S( for both sets of parameter values.

Discussion

In this model, the generation of biogas is ensured by the survival of methanogens X and X. Thus, regions that include the stability of ε or ε correspond to a complete system failure because no biogas is produced. In [7], it is pointed out that in a ‘healthy’ reactor, about 70% of the methane is produced through acetoclastic and about 30% through hydrogenotrophic methanogenesis. In our model, the equilibria ε and ε correspond to system failure. If a reactor is in a state corresponding to ε, it is often referred to as failure due to acid accumulation. All other equilibria correspond to some degree of methane production. To maintain a healthy reactor, it would appear that ideally parameter values of D and S( should be chosen in a region of stability of ε1 but not of bistability. Our case study shows that due to the effect of inhibition, reactor operators should be careful to increase S(. This could move a healthy reactor from a region of stability of ε1 into a region of bistability of ε1 and ε or even worse, bistability of ε1 and ε or ε1 and ε. Also, although when inhibition is ignored, the model predicts that a small change in D could only move the reactor from a region of stability of ε1 to one of stability of ε1. When inhibition is not ignored, besides the above scenario, increasing D could also move the reactor to a region of bistability of ε1 and ε or ε1 and ε or of stability of ε. The many possible bistabilities imply that it is very difficult to identify optimal ranges for the operating parameters and to presume stability of ε1 in experimental studies for parameter estimation.

We have shown that inhibition has a significant impact on the stability of equilibria. Even weak inhibitions can lead to dramatically different outcomes and cannot be ignored. In particular, the interaction of various inhibitions leads to changes in the size and presence of regions of stability and bistability. Some of these results were already noted in [10], which considered acidogenesis and acetoclastic methanogenesis only. In contrast, in [8], no bistabilities were observed. In [8], the focus was on studying the hydrogen transfer only without taking into account the impact of acid. This corresponds to setting μ = 0 while μ > 0, which leads to no new regions of stability in comparison to μ = 0. Therefore, the dynamics of acidogenic bacteria, acetogenic bacteria, and methanogens is driven by the interaction of all nutrient dependencies and inhibitions.

Table A2.

Inequalities corresponding to conditions of existence and local stability of the equilibria of Equation (9).

Equilibria	εo	εH	εA1	εAH1	εV	εVH	εVA1	ε∗1
f1	S(o) < f1		S(o) < f1		S(o) < f1		S(o) < f1
fi2	S(o) < f21	S(o) < f21	S(o) < f21	S(o) < f21
	or	or
	S(o) < f22	S(o) < f22
f3	S(o) < f3	S(o) < f3
f4		S(o) < f4				S(o) < f4
fi5			S(o) < f51	S(o) < f51
fi6				S(o) < f61				S(o) < f61
fi7					S(o) < f71orS(o) < f72		S(o) < f71
f8					S(o) < f8	S(o) < f8
fi9						S(o) < f91orS(o) < f92		S(o) < f91
fi10							S(o) < f101
fi11								S(o) < f111