Global analysis of a simplified model of anaerobic digestion and a new result for the chemostat

TL;DR

This paper provides a comprehensive global analysis of a simplified anaerobic digestion model, proving convergence to equilibrium and exploring the impact of environmental perturbations through stochastic algorithms.

Contribution

It introduces a global stability analysis for a simplified ADM1 model and develops algorithms to simulate environmental effects on anaerobic digestion.

Findings

No periodic orbits exist even with bistability.

Solutions always converge to an equilibrium.

Environmental perturbations significantly affect transient dynamics.

Abstract

A. Bornh\"oft, R. Hanke-Rauschenbach, and K. Sundmacher, [Nonlinear Dyn., 73 (2013), pp. 535-549] introduced a qualitative simplification to the ADM1 model for anaerobic digestion. We obtain global results for this model by first analyzing the limiting system, a model of single species growth in the chemostat in which the response function is non-monotone and the species decay rate is included. Using a Lyapunov function argument and the theory of asymptotically autonomous systems, we prove that even in the parameter regime where there is bistability, no periodic orbits exist and every solution converges to one of the equilibrium points. We then describe two algorithms for stochastically perturbing the parameters of the model. Simulations done with these two algorithms are compared with simulations done using the Gillespie and tau-leaping algorithms. They illustrate the severe impact environmental factors may have on anaerobic digestion in the transient phase.

Tables

Table 1. Table 1: The parameter values used in the following bifurcation diagrams are the ones used in [ 4 ] , except m I ​ I , m I ​ I ​ I , r , subscript 𝑚 𝐼 𝐼 subscript 𝑚 𝐼 𝐼 𝐼 𝑟 m_{II},m_{III},r, and a 𝑎 a , which were chosen so that the response functions μ 2 , I ​ I subscript 𝜇 2 𝐼 𝐼 \mu_{2,II} and μ 2 , I ​ I ​ I subscript 𝜇 2 𝐼 𝐼 𝐼 \mu_{2,III} closely resemble μ 2 , I subscript 𝜇 2 𝐼 \mu_{2,I} . The parameter D 𝐷 D is the bifurcation parameter in Figs.   4(a) , 4(c) and  4(e) , and S ( 0 ) superscript 𝑆 0 S^{(0)} is the bifurcation parameter in Figs.   4(b) , 4(d) and  4(f) .

Parameter	$S^{\left(0\right)}$	$D$	$D_i,i=1,2$	$\kappa$	$K$	$k_1$	$k_2$	$r$	$r_1$
Value	50	0.15	0.16	1.2	9.28	0.05	0.5	0.1	7.1

Table 2. Table 2: Equilibria for system Eq.   1 with parameters given in Table   1 , with μ 2 ​ ( S 2 , S 3 ) = μ 2 , III ​ ( S 2 , S 3 ) subscript 𝜇 2 subscript 𝑆 2 subscript 𝑆 3 subscript 𝜇 2 III subscript 𝑆 2 subscript 𝑆 3 \mu_{2}(S_{2},S_{3})=\mu_{2,\text{III}}(S_{2},S_{3}) .

Equilibria
$E$	(50, 0, 0, 0, 0)
$E_0$	(1.092, 1.088, 135.2, 1.352, 0)
$E_1$	(1.092, 1.088, 3.304, 1.352, 0.4614)
$E_2$	(1.092, 1.088, 28.09, 1.352, 0.3747)

Equations

$$\dot{S_1}$$

$$\dot{X_1}$$

$$\dot{S_2}$$

$$\dot{S_3}$$

$$\dot{X_2}$$

$${\mu }_1({\lambda }_1)=D_1.$$

$${\lambda }_2$$

$${\lambda }_3$$

$${\mu }_2({\sigma }_i,{\lambda }_3)=D_2.$$

$$E$$

$$E_0$$

$$E_1$$

$$E_2$$

$$X_1^{\ast }=\frac{D(S^{(0)}-{\lambda }_1)}{y_1{D}_1}\text{mbox}and{X}_{2,{\sigma }_i}^{\ast }=\frac{D({\lambda }_2-{\sigma }_i)}{y_3{D}_2}.$$

$$\dot{S_2}$$

$$\dot{X_2}$$

$$X(t)=y_3{X}_2(t),S(t)=S_2(t),\mu (S(t))={\mu }_2(S_2(t),{\lambda }_3),S^0={\lambda }_2,$$

$$\dot{S}(t)$$

$$\dot{X}(t)$$

$${\mathcal{E}}_0=(S^0,0),{\mathcal{E}}_1=({\sigma }_1,X_{{\sigma }_1}^{\ast }),{\mathcal{E}}_2=({\sigma }_2,X_{{\sigma }_2}^{\ast }).$$

$$\mathcal{G}=\{(S,X):0<S<{\sigma }_2,X>0\}.$$

$$V(S,X)={\int }_{{\sigma }_1}^S\frac{(\mu (\xi )-D_2)(S^0-{\sigma }_1)}{D_2(S^0-\xi )}d\xi +\left[X-X_{{\sigma }_1}^{\ast }-X_{{\sigma }_1}^{\ast }\ln \left(\frac{X}{X_{{\sigma }_1}^{\ast }}\right)\right],$$

$$\dot{V}(S,X)$$

$$\nabla V(S,X)$$

$$V({\sigma }_2,X_{{\sigma }_1}^{\ast })={\int }_{{\sigma }_1}^{{\sigma }_2}\frac{(\mu (\xi )-D_2)(S^0-{\sigma }_1)}{D_2(S^0-\xi )}d\xi .$$

$$\mathcal{U}=\{(S,X)\in {\mathbb{R}}_{+}^2:0\le S\le {\sigma }_2,V(S,X)\le V({\sigma }_2,X_{{\sigma }_1}^{\ast })\}$$

$${\mu }_{2,\text{I}}(S_2,S_3)=\frac{m_I{S}_2{S}_3}{(K+S_2)(S_3+k_1{S}_2)(1+k_2{S}_3)},$$

$${\mu }_{2,\text{II}}(S_2,S_3)$$

$${\mu }_{2,\text{III}}(S_2,S_3)$$

$${\mu }_1(S_1)=\frac{\kappa S_1}{r_1+S_1},$$

$${\mu }_1(S_1)=\frac{\kappa S_1}{r_1+S_1}$$

$$(S_1(0),X_1(0),S_2(0),S_3(0),X_2(0))=(50,0.4,0,0,1.16),$$

$$(S_1(0),X_1(0),S_2(0),S_3(0),X_2(0))=(50,0.4,0,0,1.14),$$

$$P_{{\tau }_{i+1}}=P_0{N}_{{\tau }_i}.$$

$$P_0=[S_0,D,y_1,y_2,y_3,y_4,K,k_1,m_{II},r]$$

$$P_{{\tau }_{i+1}}=P_{{\tau }_i}{N}_{{\tau }_i}=P_0\prod _{n=1}^i{N}_{{\tau }_n},$$

Full text

{AMS}

34C60, 34D20, 34D23, 70K05, 93D30, 92D25, 93D30

1 Introduction

Anaerobic digestion is a biochemical process where microorganisms or multicellular organisms break down organic material in the absence of oxygen. Anaerobic digestion is an important part of many industrial practices, including the treatment of wastewater and the production of biogas. The role of anaerobic digestion in such applications has been an active area of recent research [1, 2, 3, 4, 10, 12, 14, 15, 19, 21]. This paper focuses on a particular model of anaerobic digestion in biogas production.

The foundation of previous work on the mathematical analysis of the production of biogas is the Anaerobic Digestion Model 1 (ADM1) [1] introduced in 2002. If implemented as a system of differential equations, this model has 32 state variables, including seven different species of microorganisms. Understandably, anything other than numerical analysis has not been feasible.

In an effort to formally analyze the system, several groups [4, 10, 12, 24] have studied various subsystems of ADM1. Recently, Weedermann et al. [24, 25] combined two previous models [10, 12] to get a reasonably complete picture using only eight state variables. Due to the inclusion of two pathways to biogas production in [24] and because the model captures the ADM1’s sensitivity to the accumulation of acetic acid, [24] illustrates some of the complexity of ADM1, which must exhibit the same or an even richer dynamics than the model in [24].

Bornhöft et al. [4] introduced a model with five state variables based on their observations from a numerical steady-state analysis of the ADM1 model, and conjectured that their model undergoes the same bifurcations as the ADM1 model with the substrate inlet concentration as bifurcation parameter. The model in [4] is the first simplified model to consider the effects of ammonia. It is demonstrated that the proposed model is able to capture the same effects of ammonia on anaerobic digestion that are displayed by the ADM1 model. However, the analysis in this paper shows that the model does not possess all of the dynamics of ADM1, even if a broader class of growth functions is considered than the ones that were initially proposed. The model is missing some of the dynamics shown in [24], namely the possible bistability between two equilibria that both correspond to biogas production, a behaviour of the full ADM1 model that is also noted in [4].

The model in [4] considers two stages of anaerobic digestion, acidogenesis and methanogenesis. In the first stage, simple substrates are broken down by acidogenic microorganisms. The microorganisms use the energy from the simple substrates to grow, and produce volatile fatty acids (VFAs) and ammonia as byproducts. The VFAs and ammonia have opposing effects on the pH of the system; an increase in the concentration of VFAs will decrease the pH, while an increase in the concentration of ammonia will increase the pH. In the second stage, methanogenic microorganisms convert the VFAs to biogas. The methanogenic microorganisms are very sensitive to the environment, and can only tolerate a relatively narrow pH range. Furthermore, ammonia is toxic to the methanogenic microorganisms in large quantities and will restrict their growth. The flow chart in Fig. 1 summarizes the process.

In this paper we provide a formal mathematical analysis of the model proposed in [4], allowing a more general class of response functions. In Section 2, we describe the model and assumptions, and give properties of the solutions of the system. If the substrate input concentration is too low, the system converges to a state where no microorganisms are present. We show that if the substrate input concentration is high enough to allow the acidogenic microbial population to survive, the system reduces to a limiting system that is a two-dimensional basic model of growth in the chemostat that includes the decay rates and allows for any non-monotone response function.

In Section 3, we study the dynamics of this limiting system. We obtain a new global result in the case that the parameters allow bistability by proving that no nontrivial periodic orbits exist.

In Section 4 we use the theory of asymptotically autonomous systems and the results for the limiting system from Section 3 to provide a complete global analysis of the anaerobic digestion model in [4] for a more general class of response functions.

In Section 5, we propose two alternative prototype functions to model the growth of the methanogenic archaea and capture the inhibition by ammonia. These prototypes complement the one used in [4], which has the property that there is no growth in the absence of ammonia. The prototypes we introduce allow growth in the absence of ammonia, but are either unimodal or decreasing in ammonia. We provide bifurcation diagrams for all three prototypes, and compare how they influence the outcome.

In Section 6 we further investigate the model when the parameters are selected so that there are two stable steady states. In industrial applications of processes such as anaerobic digestion, operators must be aware of how physical and environmental processes, and changes in the biology of the species can affect the long term health of the reactor. One way to address these challenges is to include the effects of stochasticity in simulations of the model. Stochasticity can be a result of random births and deaths, or of fluctuations in the model parameters, possibly due to mutations or changes in the environment. Models of chemostats that include stochasticity have been considered in the literature (e.g., [6, 13]), but our approach differs from the ones presented in those papers. We consider stochasticity in the case where there are fluctuations in the parameters, and compare the results to two well known methods for simulating stochasticity in the case where there are random births and deaths. Our studies give new insight into why seemingly identical reactor setups can lead to very different reactor performances. In Section 7, we summarize our results and discuss their implications for biogas production using anaerobic digestion.

2 The Model

Let $X_{1}$, $X_{2}$, $S_{1}$, $S_{2}$ and $S_{3}$ denote the concentrations of the acidogenic microorganisms, methanogenic microorganisms, simple substrates, acetic acid and ammonia, respectively. The model is described by the system

[TABLE]

where $D$ is the dilution rate, $S^{(0)}$ is the input concentration of simple substrates, $D_{i}=D+k_{i}$, where $k_{i}\geq 0$ are the respective decay rates of $X_{i}$, and $y_{i}$ are yield constants.

Let $\mathbb{R}_{+}$ and $\mathbb{R}^{2}_{+}$ denote the set of non-negative real numbers and the non-negative plane, respectively. We make the following assumptions concerning $\mu_{1}$ and $\mu_{2}$:

1. (H1)

$\mu_{1}(S_{1})\in C^{1}(\mathbb{R}_{+})$, and $\mu_{1}^{\prime}(S_{1})>0$ for all $S_{1}>0$. 2. (H2)

$\mu_{1}(0)=0$, $\mu_{1}(S_{1})>0$ for all $S_{1}>0$. 3. (H3)

$\mu_{2}(S_{2},S_{3})\in C^{1}(\mathbb{R}^{2}_{+})$, and $\mu_{2}(S_{2},S_{3})>0$ if $S_{2}>0$ and $S_{3}>0$. 4. (H4)

$\lim_{S_{3}\to\infty}\mu_{2}(S_{2},S_{3})=0$ for all $S_{2}\geq 0$. 5. (H5)

$\lim_{S_{2}\to\infty}\mu_{2}(S_{2},S_{3})=0$ for all $S_{3}\geq 0$. 6. (H6)

$\mu_{2}(0,S_{3})=0$ for all $S_{3}\geq 0$ and $\mu_{2}(S_{2},0)\geq 0$ for all $S_{2}>0$ 7. (H7)

There exists $\Gamma(S_{3})\in C(\mathbb{R}_{+})$ such that for $S_{2}<\Gamma(S_{3}),~{}\partial_{S_{2}}\mu_{2}(S_{2},S_{3})>0$ and for $S_{2}>\Gamma(S_{3}),$ $\partial_{S_{2}}\mu_{2}(S_{2},S_{3})<0$.

Unlike in [4], we do not assume that both $X_{i}$ have identical decay rates. (H1) and (H2) are satisfied by any of the Holling type I, II or III growth functions, which are standard to chemostat models. (H3), (H4), and (H5) capture the inhibitory nature of $S_{2}$ and $S_{3}$, guaranteeing that large quantities of either $S_{2}$ or $S_{3}$ will be detrimental to the growth of $X_{2}$. (H6) ensures that an absence of acetic acid will result in no growth of the methanogenic microorganisms, while an absence of ammonia does not necessarily have this effect. We would like to note that the prototype describing the growth of methanogens proposed in [4] has the property that $\lim_{S_{3}\to 0}\mu_{2}(S_{2},S_{3})=0$ for all $S_{2}\geq 0$. We decided to relax this condition. (H7) intends to capture the nature of the inhibition mechanisms outlined in [4], whereby small concentrations of $S_{2}$ are limiting on the growth of $X_{2}$, while large concentrations of $S_{2}$ will increase the pH and hence be inhibitory. The curve $\Gamma(S_{3})$ indicates that for each fixed $S_{3}$ there is at most one $S_{2}$ such that $\partial_{S_{2}}\mu_{2}(S_{2},S_{3})=0$. The curve $\Gamma(S_{3})$ therefore gives the optimal concentration of $S_{2}$ for growth of $X_{2}$ as a function of $S_{3}$. We make no further assumptions about how $\mu_{2}(S_{2},S_{3})$ changes with $S_{3}$. In many cases, including ADM1, $\mu_{2}(S_{2},S_{3})$ will have a unimodal shape for fixed $S_{2}$, but we do not want to rule out the possibility of other profiles that may be useful. For examples of functions satisfying these hypotheses, refer to Section 5.

Here, we introduce some notation. The break-even concentration of $S_{1}$, $\lambda_{1}$, is the unique positive extended real number that satisfies

[TABLE]

If no such number exists, we take $\lambda_{1}=+\infty$. When $S_{1}=\lambda_{1}$ and $X_{2}=0$, the equilibrium concentrations of $S_{2}$ and $S_{3}$, $\lambda_{2}$ and $\lambda_{3}$, respectively, are given by

[TABLE]

The break-even concentrations of $S_{2}$, $\sigma_{1}$ and $\sigma_{2}$, are the extended real numbers $\sigma_{2}\geq\sigma_{1}$ that solve

[TABLE]

If no such numbers exist, which is the case when $\mu_{2}(\Gamma(S_{3}),S_{3})<D_{2}$, then we write $\sigma_{1}=\sigma_{2}=+\infty$. Note that by (H7), there is at most one turning point of $\mu_{2}(S_{2},\lambda_{3})$, and therefore at most two solutions to $\mu_{2}(S_{2},\lambda_{3})=D_{2}$.

System Eq. 1 has a total of four possible equilibria,

[TABLE]

where

[TABLE]

These equilibria are only biologically meaningful if each of the components is non-negative. $E_{1}$ and $E_{2}$, when they exist, are called interior equilibria, since they lie in the interior of the positive cone $\mathbb{R}_{+}^{5}$. $E$ and $E_{0}$ are called boundary equilibria, since they lie on the boundary of the positive cone $\mathbb{R}_{+}^{5}$.

The following propositions give well-posedness results for system Eq. 1, provide conditions for the washout of the microorganisms in the reactor when the substrate input concentration is too low, and introduce the limiting system when the input concentration is high enough so that the acidogens survive. The proofs are given in Appendix A.

Proposition 1**.**

Assume that $X_{i}(0)\geq 0$ and $S_{i}(0)\geq 0$.

1. i)

If $X_{1}(0)=0$, then solutions converge to $E$ as $t\to\infty$. 2. ii)

If $X_{1}(0)>0$ and $X_{2}(0)=0$, then $X_{1}(t)>0$ and $S_{i}(t)>0$ for all $t>0$ while $X_{2}(t)=0$ for all $t\geq 0$. 3. iii)

If $X_{1}(0)>0$ and $X_{2}(0)>0$, then $X_{i}(t)>0$ and $S_{i}(t)>0$ are positive for all $t>0$. 4. iv)

All solutions are bounded for $t\geq 0$.

Proposition 2**.**

*If $\lambda_{1}\geq S^{(0)}$, then $E$ is a globally asymptotically stable equilibrium of Eq. 1. *

Proposition 3**.**

If $\lambda_{1}<S^{(0)}$, then Eq. 1 is a quasi-autonomous system with limiting system

[TABLE]

By Theorem 1.4 in [22], it will be enough to study the dynamics of this limiting system.

3 Global Analysis of Growth in the Chemostat

After the change of variables

[TABLE]

system Eq. 7 becomes a model of the chemostat:

[TABLE]

Recall that $\mu_{2}(S(t),\lambda_{3})$ satisfies (H3) and (H4), and hence $\mu(S(t))$ is a non-monotone response function with break-even concentrations $0<\sigma_{1}<\sigma_{2}$, the extended real numbers that solve $\mu(\sigma_{i})=D_{2}$.

We define the equilibria of Eq. 8 that correspond to $E_{0},E_{1},$ and $E_{2}$, respectively, for system Eq. 1 defined in Eq. 6b-Eq. 6d:

[TABLE]

where $X^{*}_{\sigma_{i}}=\displaystyle\frac{D(S^{0}-\sigma_{i})}{D_{2}},\ i=1,2$.

Models of the chemostat have been well studied (e.g., see [20, 11] and the references therein). Model Eq. 8 is a model of growth of a single species in the chemostat with non-monotone response function that includes the species decay rate, i.e. $D_{2}=D+\epsilon$ where $\epsilon>0$ is the species decay rate.

In Wolkowicz and Lu [26], model Eq. 8 extended to the $n$ species case was analyzed. The results of that paper, if applied to the single species growth model, completely determine the dynamics of Eq. 8 when $\mu(S(t))$ is any monotone increasing function or it is non-monotone and $\sigma_{1}<S^{0}\leq\sigma_{2}$. However, the case that $\mu(S(t))$ is a non-monotone response function and $\sigma_{1}<\sigma_{2}<S^{0}$, remained open. Here, we provide a proof in this case and thus complete the global analysis of system Eq. 8. In particular, we prove that there are no periodic orbits, and hence although the outcome is initial condition dependent, either the species dies out or it approaches an equilibrium.

In the following theorem, we summarize what is known for the dynamics of Eq. 8, and provide a proof in the case that had remained open.

Theorem 4**.**

Consider model Eq. 8. Assume $\mu(S)$ is continuously differentiable, $\mu(0)=0$, $\mu(S)\geq 0$ for all $S>0$, and there exist positive numbers $\sigma_{1}\leq\sigma_{2}$ (possibly infinite) such that $\mu(S)<D_{2}$ if $0<S<\sigma_{1}$, $\mu(S)>D_{2}$ if $\sigma_{1}<S<\sigma_{2}$, and $\mu(S)<D_{2}$ if $S>\sigma_{2}$. Let $S(0)\geq 0$ and $X(0)>0$.

1. i)

If $S^{0}\leq\sigma_{1}\leq\sigma_{2}$, then $\mathcal{E}_{0}$ is globally asymptotically stable. 2. ii)

If $\sigma_{1}<S^{0}\leq\sigma_{2}$, then $\mathcal{E}_{1}$ is globally asymptotically stable. 3. iii)

If $\sigma_{1}=\sigma_{2}<S^{0}$, then $\mathcal{E}_{0}$ is locally asymptotically stable and attracts all solutions except the solutions in the stable manifold of $\mathcal{E}_{1}=\mathcal{E}_{2}$. 4. iv)

If $\sigma_{1}<\sigma_{2}<S^{0}$, then $\mathcal{E}_{1}$ and $\mathcal{E}_{0}$ are locally asymptotically stable and $\mathcal{E}_{2}$ is a saddle. Furthermore, any orbit that is not in the stable manifold of $\mathcal{E}_{2}$ converges to either $\mathcal{E}_{1}$ or $\mathcal{E}_{0}$.

Proof 3.1**.**

i)-ii) See [26].

iii) This result follows from standard phase plane analysis. When $\sigma_{1}=\sigma_{2}$, $\mathcal{E}_{1}$ and $\mathcal{E}_{2}$ coalesce and are unstable. All orbits converge to $\mathcal{E}_{0}$ except those in the stable manifold of the degenerate saddle $\mathcal{E}_{1}=\mathcal{E}_{2}$.

iv) If $\sigma_{1}<\sigma_{2}<S^{0}$, then $\mathcal{E}_{0}$, $\mathcal{E}_{1}$, and $\mathcal{E}_{2}$ all lie in $\mathbb{R}_{+}^{2}$. From standard local stability analysis, it follows that $\mathcal{E}_{0}$ and $\mathcal{E}_{1}$ are both locally asymptotically stable and $\mathcal{E}_{2}$ is a saddle.

Next we show that no nontrivial periodic solutions are possible. We proceed using proof by contradiction. Suppose that there exists a nontrivial periodic solution, $\Phi$. By 1 all solutions are bounded in forward time and the first quadrant is invariant. Orbits above the $\dot{S}=0$ nullcline (shown in Fig. 2 by the dashed curve) move from right to left, and so if they cross the $\dot{X}=0$ nullclines (Shown in Fig. 2 by dashed vertical lines) cross them from right to left. Orbits below the $\dot{S}=0$ nullcline move from left to right, and so if they cross the $\dot{X}=0$ nullclines, cross them from left to right. Orbits thatmare to the right of the line $S=\sigma_{2}$ that meet the $S^{\prime}=0$ nullcline cross it downward. Therefore, if a periodic orbit exists, it must lie entirely to the left of $S=\sigma_{2}$, since if an orbit enters the region below the $S^{\prime}=0$ nullcline and to the right of $S=\sigma_{2}$, it is trapped in that region and must converge to $\mathcal{E}_{0}$ by the Poincaré-Bendixson Theorem. If it is to the left of $S=\sigma_{2}$ and above $\dot{S}=0$ it moves to the left. Therefore, by the Poincaré-Bendixson Theorem and standard phase-plane analysis, $\Phi$ must surround $\mathcal{E}_{1}$, and must lie entirely in the set

[TABLE]

Define the Lyapunov function,

[TABLE]

as in [26]. See Fig. 2 for phase portraits of system Eq. 8 with typical level sets of the Lyapunov function. Note that Eq. 9 is a valid Lyapunov function for $\mathcal{E}_{1}$ in $\mathcal{G}$, and

[TABLE]

is non-positive, for all $S\in[0,\sigma_{2}]$ and $X\geq 0$, i.e., for all $S$ in the closure of $\mathcal{G}$. By examining

[TABLE]

we see that $\mathcal{E}_{1}$ and $(\sigma_{2},X_{\sigma_{1}}^{*})$ are the only critical points of $V(S,X)$ for which $S\leq\sigma_{2}$. The point $(\sigma_{2},X_{\sigma_{1}}^{*})$ is directly above $\mathcal{E}_{2}$ in phase space, since by definition $X_{\sigma_{1}}^{*}>X_{\sigma_{2}}^{*}$. Notice that $\partial^{2}V(S,X)/\partial X^{2}=X_{\sigma_{1}}^{*}/X^{2}>0$ for all $X>0$, $\partial V(S,X)/\partial S<0$ for $0<S<\sigma_{1}$, $\partial V(S,X)/\partial S>0$ for $\sigma_{1}<S<\sigma_{2}$, and $\partial V(S,X)/\partial S<0$ for $\sigma_{2}<S<S^{0}$. It follows that $\mathcal{E}_{1}$ is a local minimum of $V(S,X)$, and $(\sigma_{2},X_{\sigma_{1}}^{*})$ is a saddle point of $V(S,X)$. The level set $V(S,X)=V(\sigma_{2},X_{\sigma_{1}}^{*})$ is given by

[TABLE]

For $S\leq\sigma_{2}$, this level set is a closed curve surrounding $\mathcal{E}_{1}$ that passes through the point $(\sigma_{2},X_{\sigma_{1}}^{*})$. (See the bold level set in Fig. 2 where two possible configurations are shown.) Since $\dot{V}(S(t),X(t))\leq 0$, the set

[TABLE]

is a positively invariant set. Since any periodic orbit must lie entirely in $\mathcal{G}$ and must surround an equilibrium point, it must enter $\mathcal{U}$. Since $\mathcal{U}$ positively invariant, it follows that $\Phi$ is contained entirely in $\mathcal{U}$. By the minor variation of LaSalle’s invariance principle [26], any trajectory in $\mathcal{U}$ converges to the largest invariant set in $\mathcal{U}\cap\{(S,X):\dot{V}(S,X)=0\}$. The only such invariant set is $\mathcal{E}_{1}$, and hence no nontrivial periodic orbit exists, a contradiction.

*Now, from standard phase plane analysis and the Poincaré-Bendixson Theorem, it follows that all orbits converge to one of the three equilibria as $t$ tends to infinity. The one-dimensional stable manifold of $\mathcal{E}_{2}$ acts as a separatrix, defining the basins of attraction for $\mathcal{E}_{1}$ and $\mathcal{E}_{0}$. *

4 Global Analysis of the Full System

Eq. 1

Theorem 5**.**

Consider model Eq. 1.

1. i)

If $\lambda_{1}\geq S_{1}^{(0)}$, then $E$ is globally asymptotically stable. 2. ii)

If $\lambda_{1}<S_{1}^{(0)}$ and $\lambda_{2}\leq\sigma_{1}\leq\sigma_{2}$, then $E_{0}$ is a globally asymptotically stable equilibrium. 3. iii)

If $\lambda_{1}<S_{1}^{(0)}$ and $\sigma_{1}<\lambda_{2}\leq\sigma_{2}$, then $E_{1}$ is a globally asymptotically stable equilibrium. 4. iv)

If $\lambda_{1}<S_{1}^{(0)}$ and $\sigma_{1}<\sigma_{2}<\lambda_{2}$, then $E_{0}$ and $E_{1}$ are locally asymptotically stable, and $E_{2}$ is a saddle. Furthermore any orbit that does not lie on the stable manifold of $E_{2}$ converges to one of $E_{0}$ or $E_{1}$.

Proof 4.1**.**

i) was proved in 2.

*ii) - iv) Since each of the $E_{i},~{}i=0,1,2$ for model Eq. 1 corresponds to $\mathcal{E}_{i}$ for system Eq. 8, the results follow from the results for the limiting system given in 4, followed by an application of the theory for asymptotically autonomous systems, either by using Theorem 1.4 in [22], or by a direct proof using the Butler-McGehee Lemma (as stated in Lemma 5.2 in [5] and applied there). *

5 Bifurcation Analysis of Full System Eq. 1

As a result of the analysis of the previous two sections, the only possible bifurcations that can occur in Eq. 1 are transcritical bifurcations and saddle-node bifurcations.

In [4], a prototype growth function was introduced to capture the inhibition caused by ammonia. This prototype

[TABLE]

has the property that when there is no ammonia, which is toxic to the methanogenic microorganisms, the methanogenic microorganisms are unable to grow. We introduce two additional prototype functions

[TABLE]

that satisfy (H3)-(H7). Both $\mu_{2,\text{II}}(S_{2},S_{3})$ and $\mu_{2,\text{III}}(S_{2},S_{3})$ satisfy the additional property, $\mu_{2,\text{}}(S_{2},0)\geq 0$ with equality only when $S_{2}=0$ or in the limit as $S_{2}\to\infty$. For the parameters given in Table 1, $\mu_{2,\text{II}}(S_{2},S_{3})$ is strictly decreasing in $S_{3}$ and can be thought of as the opposite extreme of $\mu_{2,\text{I}}(S_{2},S_{3})$. It describes the scenario where ammonia is strictly inhibitory and the methanogenic microorganisms do best without any ammonia present. With a different set of parameters this response function can be unimodal in $S_{3}$. The difference term in the denominator of $\mu_{2,\text{II}}(S_{2},S_{3})$ acts as a proxy for the influence of pH in the system. Since $S_{2}$ is acidic, and $S_{3}$ is basic, we assume that a large difference between the two concentrations would cause the pH to be outside of the acceptable range for the growth of $X_{2}$. The third prototype, $\mu_{2,\text{III}}(S_{2},S_{3})$ covers the middle ground between $\mu_{2,\text{I}}(S_{2},S_{3})$ and $\mu_{2,\text{II}}(S_{2},S_{3})$; it is unimodal in $S_{3}$, like $\mu_{2,\text{I}}(S_{2},S_{3})$, but is non-zero when $S_{2}>0$ and $S_{3}=0$, like $\mu_{2,\text{II}}(S_{2},S_{3})$.

The substrate input concentration, $S^{(0)}$, and dilution rate, $D$, are the two parameters that the operator of a reactor has the ability to control. In our bifurcation analysis, we focus on how the dynamics of the full system Eq. 1 change when these parameters vary. We note that $\lambda_{2}$ and $\lambda_{3}$ depend on $S^{(0)}$ (see Eqs. 3 and 4), and hence, $\max_{S_{2}>0}\mu_{2}(S_{2},\lambda_{3})$ changes when $S^{(0)}$ changes. From the stability analysis in Section 3, two scenarios are possible. In the first scenario (see Fig. 4), there is a transcritical bifurcation when $\lambda_{2}=\sigma_{1}$, a transcritical bifurcation when $\lambda_{2}=\sigma_{2}$, and a saddle-node bifurcation when $\max_{S_{2}>0}\mu_{2}(S_{2},\lambda_{3})=D_{2}$. In the second scenario (see Fig. 5) there are two saddle-node bifurcations as $\lambda_{3}$ increases. This sequence of bifurcations occurs because $\mu_{2,\text{I}}(S_{2},S_{3})$ and $\mu_{2,\text{III}}(S_{2},S_{3})$ are unimodal in $S_{3}$. With the parameters listed in Table 1, the second prototype, $\mu_{2,\text{II}}(S_{2},S_{3})$, is strictly decreasing in $S_{3}$, and so only the first scenario is possible. The other two prototypes, $\mu_{2,\text{I}}(S_{2},S_{3})$ and $\mu_{2,\text{III}}(S_{2},S_{3})$, are unimodal in $S_{3}$, and either scenario is possible.

In the bifurcation diagrams shown in Figs. 4 and 5,

[TABLE]

and the parameters are the ones used in [4]. Any parameters not given in [4] (e.g., $m_{II}$, $m_{III}$, $r$, and $a$), were chosen so that the functions, $\mu_{2,\text{II}}$ and $\mu_{2,\text{III}}$, closely resemble the function $\mu_{2,I}$ given in [4]. See Table 1 for the parameter values used. A plot of each function is shown in Fig. 3. The bifurcation diagrams in Fig. 4 are qualitatively similar for each uptake function. The bifurcation diagrams corresponding to $\mu_{2,\text{II}}(S_{2},S_{3})$ and $\mu_{2,\text{III}}(S_{2},S_{3})$ resemble the diagram for ADM1 in [4] more closely than the diagram for $\mu_{2,\text{I}}(S_{2},S_{3})$.

In the diagrams where $D$ was used as the bifurcation parameter (Figs. 4(a), 4(c) and 4(e)), there are three clear regions. In the first region when $0<D<D^{*}_{1}$, only the equilibria $\mathcal{E}_{1}$ and $\mathcal{E}_{0}$ lie in the positive cone, $\mathcal{E}_{1}$ is globally asymptotically stable and therefore all non-stationary solutions converge to $\mathcal{E}_{1}$. When $D=D^{*}_{1}$ the washout equilibrium $\mathcal{E}_{0}$ undergoes a transcritical bifurcation. In the second region, where $D^{*}_{1}<D<D_{2}^{*}$ all three equilibria lie in the positive cone. $\mathcal{E}_{1}$ and $\mathcal{E}_{0}$ are locally asymptotically stable and $\mathcal{E}_{2}$ is a saddle. All solutions (except the stable manifold of $\mathcal{E}_{2}$) converge to one of $\mathcal{E}_{1}$ or $\mathcal{E}_{0}$, depending on initial conditions. When $D=D_{2}^{*}$, the two interior equilibria $\mathcal{E}_{1}$ and $\mathcal{E}_{2}$ undergo a saddle-node bifurcation. In the third region, where $D_{2}^{*}<D$ only $\mathcal{E}_{0}$ exists, and it is globally asymptotically stable. Therefore all solutions tend to $\mathcal{E}_{0}$.

6 Stochastic Simulations of the Full

System Eq. 1

We describe two stochastic algorithms to capture stochasticity in the parameters. For comparison we also include simulations done with Gillespie’s stochastic simulation algorithm [8] and the adaptive tau-leaping algorithm [9].

The simulations in this section are all done for the full system Eq. 1 with

[TABLE]

and $\mu_{2}(S_{2},S_{3})=\mu_{2,\text{III}}(S_{2},S_{3})$. The parameters are listed in Table 1. With these parameters, the deterministic system has two stable equilibria, $E_{0}$ and $E_{1}$, and so the long-term behaviour of the solutions is initial condition dependent. If

[TABLE]

the solution of the deterministic system converges to $E_{1}$ (see Table 2), and if

[TABLE]

the solution of the deterministic system converges to $E_{0}$ (see Table 2). Thus, for one set of initial conditions, the deterministic system Eq. 1 predicts that the methanogens survive and produce biogas, and for the other it predicts that they do not. These initial conditions simulate the start up and inoculation of the reactor. The only difference between the initial conditions in Eqs. 13 and 14 is the value of $X_{2}(0)$. Both initial conditions are close to the separatrix. We only include figures that show the population of methanogens, $X_{2}(t)$, to compare the effect of stochasticity on biogas production, which only occurs if $X_{2}$ is positive. In simulations (not shown) with initial conditions farther from the separatrix, solutions converged to the same equilibrium predicted by the deterministic model every time. The figures were produced using Matlab [17].

We use two different approaches to study the behavior of Eq. 1 under stochastic perturbations. The first method is meant to model fluctuations in the parameters due, for example, to fluctuations in the environment. The second method captures the effect of potential mutations in members of the populations. In both schemes, multiple parameters are perturbed at randomly chosen times. Because we are varying many parameters, some of which appear in the non-linearities of the system, we are unable to write the resulting stochastic equations as a linear stochastic perturbation of the original system as was done in [23, 27] for chemostat models. In [23], the dilution rate and in [27], the dilution rate and the decay rates are assumed to vary stochastically. In one algorithm the perturbations are from the mean and in the other the perturbations are accumulative. Between perturbations the system is treated as a deterministic system that is solved numerically.

Let $\tau_{0}=0$ and $\tau_{i+1}=\tau_{i}-\ln(T_{i})$, where $T_{i}\in(0,1)$ is a uniformly distributed random variable. Therefore, $\{\tau_{i}\}$ describes a monotone increasing sequence of times. By applying the inverse sampling transform, we see that the difference $\tau_{i+1}-\tau_{i}$ is exponentially distributed with unit mean and variance. Let $P_{0}$ be a row vector containing the parameter values present in the deterministic system that are affected by stochasticity. At each randomly chosen time $\tau_{i}$, these parameters values are updated to obtain a sequence of vectors $\{P_{\tau_{i}}\}_{i=1}^{\infty},$ and we set the parameters equal to $P_{t}=P_{\tau_{i}}$, for $t\in[\tau_{i},\tau_{i+1})$.

In the first stochastic algorithm, which we call the environmental based fluctuation algorithm, we assume that the parameter values are influenced by the environment. As such, they cannot be perfectly controlled and so at random intervals of time they undergo small random changes. However, the parameters remain near their mean values given in the row vector $P_{0}$. Following this interpretation, we let $N_{t}$ be a diagonal matrix with entries given by Gaussian random variables with mean $\mu=1$ and standard deviation, $\sigma$. We assume that $N_{t}=N_{\tau_{i}}$ for $t\in[\tau_{i},\tau_{i+1})$. Then

[TABLE]

Figs. 6 and 6 show five simulations using the environmental based algorithm with $\sigma=\frac{1}{10}$ and

[TABLE]

In Fig. 6 the initial conditions are given by Eq. 14 and the solution to the deterministic system converges to $E_{0}$. In Fig. 6, the initial conditions are given by Eq. 13 and solutions converge to $E_{1}$. The solutions for the deterministic system are shown in bold for comparison.

In the second stochastic algorithm, which we call the mutation based algorithm, we assume that the parameters are dependent on properties of the microorganisms that can mutate, and therefore are subject to changes at random times that accumulate. In this case, many of the parameters are beyond the control of the operator, however we assume that the operator has complete control of the dilution rate $D$ and the input concentration $S_{0}$. Following this interpretation, we update the parameters at random times to obtain,

[TABLE]

where again $\sigma=\frac{1}{10}$,

[TABLE]

and $N_{\tau_{i}}$ are as before. Using this algorithm, $\{P_{\tau_{i}}\}_{i=1}^{\infty}$ is a random walk with mean $P_{0}$, and the mutations accumulate. Random walks have the property that $\sigma^{2}\to\infty$ as $t\to\infty$, and therefore the system is subject to wild fluctuations as time increases. Care must be taken so that the parameters, which have interpretations as positive quantities only, do not become negative. We ensure non-negativity by taking $P_{\tau_{i+1}}=\max\{0,P_{\tau_{i}}N_{\tau_{i}}\}$, and control the wild fluctuations by limiting the difference between current parameter values $P_{\tau_{i}}$ and the initial parameter values $P_{0}$ to be less than four standard deviations. Figs. 6 and 6 shows five simulations using the mutation based algorithm.

We also include simulations using Gillespie’s stochastic simulation algorithm (SSA) [8], in Figs. 7 and 7. The SSA is an essentially exact description for systems with a finite number of interacting particles. The SSA is based on the principle of mass action, and as such the deterministic system must be converted to an equivalent system that is of the form

[TABLE]

To do so, we rescale the time variable by $dt=(r_{1}+S_{1})(K+k_{1}S_{2}+rS_{2}^{2})(a+S_{3}^{2})d\hat{t}$. The resulting system has 104 different reaction terms that must be accounted for. As such, reporting the system here would be impractical. Although we have rescaled the time variable, the dynamics of system Eq. 17 are identical to those of Eq. 1. The SSA assumes that each reaction occurs independent of the others, and occurs with rates given by the coefficients of the differential equations. The SSA determines a time until each reaction takes place using the rate coefficients and the population of individuals relevant to that reaction, and increases or decreases the population(s) of the fastest reaction by a set step size. Once we have realized the simulation, we scale time back to the original time variable before plotting in order to compare with the other stochastic algorithms. Five simulations with a step size of $\tfrac{1}{100}$ are shown in Figs. 7 and 7. In reality, the step size is meant to represent a single individual in the population, but since SSA is notoriously slow, modelling a population of trillions of microorganisms and on the order of $10^{23}$ molecules in this way is computationally impossible. It is also well known that as you decrease the step size, the SSA will approach the deterministic solution [8].

Finally, we include simulations using the adaptive tau-leaping algorithm in Figs. 7 and 7. The tau-leaping algorithm is an improvement on the SSA in terms of speed, and is generally easier to implement, although it is less accurate. One interpretation of the tau-leaping algorithm is that it is analogous to Euler’s method, but instead of the derivative, a Poisson random variable with mean proportional to the derivative is used. Here, Eq. 1 takes the form

[TABLE]

where $\delta$ is the step size (typically interpreted to be an individual particle, as with the SSA). There has been much discussion on how to choose $\tau$ appropriately [7, 9]. We chose

[TABLE]

so that the fastest reaction determines $\tau$.

The stochasticity as simulated in the environmental based fluctuation algorithm and the mutation based algorithm stems from uncertainty in the system parameters, whether due to environmental noise or from mutations. The stochasticity of the SSA and tau-leaping algorithm is derived from the fact that the populations are treated as discrete quantities. Since the populations are very large in practice, it may be more realistic to implement stochasticity using continuous hybrid algorithms that reflect the uncertainty in the parameters.

In the simulations using all four algorithms, if the stochasticity caused the system to predict a different outcome than the deterministic system, it usually happened while the system was transient. Once the system neared an equilibrium, the behaviour was usually quite stable. In rare instances, noise caused the system to destabilize after nearing an equilibrium, but this seemed only to occur for the mutation based method when the noise was quite large.

7 Conclusion

We analyze the system introduced by Bornhöft et al. [4], which was proposed as a qualitative reduction of the ADM1 model, and claimed to capture the most relevant qualitative features of the ADM1 model. We give a complete global analysis of the dynamics of the model. If the concentration of the simple substrates is too low, both the acidogenic and methanogenic populations of microogranisms are eliminated from the reactor and no biogas is produced. Even if the input concentration of simple substrates is high enough, if the equilibrium concentration of VFAs produced by the acidogenic microorganisms is too low, then the methanogenic microorganisms will be eliminated from the reactor, and the system will converge to an equilibrium where no biogas is produced. If the VFA concentration is in a proper range, the system has a single globally stable interior equilibrium. Finally, if the equilibrium concentration of VFAs is very high, then the system possesses two stable equilibria and one unstable equilibrium, and no sustained oscillatory behaviour is possible. In this case the long-term behavior is initial condition dependent. Only one of the two stable equilibria corresponds to the production of biogas meaning that it depends on the initial conditions whether the reactor will produce biogas in the long-term. The system does not allow bistability involving two or more biogas producing equilibria, previously shown to be possible for the ADM1 model [1] and for the models studied in [24, 25]

The dynamics predicted by a bifurcation analysis of the model is qualitatively similar for all three prototype functions. Ammonia inhibition is included in the ADM1 model, however, in ADM1 ammonia is not included as a dynamic variable. Ammonia concentration in ADM1 is computed as the difference of the concentration of inorganic nitrogen and $NH_{4}^{+}$. In the present model, ammonia is included as a dynamic variable and it is important to determine how to best model the effect of ammonia on the growth of the methanogens to capture the behaviour of ADM1. For all three prototype functions, inhibition of the growth of acetoclastic methanogens due to ammonia is unimodal with respect to the ammonia concentration. However, for $\mu_{2,\text{I}}(S_{2},S_{3})$, acetoclastic methanogens will not grow in the absence of ammonia, while for $\mu_{2,\text{II}}(S_{2},S_{3})$ and $\mu_{2,\text{III}}(S_{2},S_{3})$ the organisms grow even if the ammonia concentration is zero. Based on a comparison with Fig. 10 in [4], using $\mu_{2,II}$ or $\mu_{2,III}$ in model Eq. 1, the behavior resembles the behavior or the ADM1 model shown in [4] more closely than using $\mu_{2,\text{I}}(S_{2},S_{3})$. This indicates that these two functions are better suited to model the dependence of acetoclastic methanogens on ammonia.

We consider two algorithms that simulate stochastic effects in system Eq. 1. The aim of these two algorithms is to model the uncertainty and variation in environmental and biological parameters that are hard to control with numerical algorithms that are easy to implement and run relatively quickly. We compare the resulting graphs with the graphs produced using the well-known Gillespie algorithm and the the tau-leaping algorithm.The stochastic simulations from all four algorithms seem to indicate that a failure of the reactor is most likely to occur early in the reactors operating cycle, and that once the reactor has reached a steady state, it is quite resilient and less affected by minor perturbations due to mutations or small fluctuations in the environment. The one possible exception is in our mutation based stochastic algorithm that is intended to simulate the accumulation of mutations within the microbial population. Therefore, it appears to be most important to control the environment of the reactor during start up, and then to carefully monitor the characteristics of the microorganisms within the reactor after start up.

The analysis of the model of anaerobic digestion proposed by Bornhöft et al. [4] involved studying the limiting system Eq. 8, a model of growth in the chemostat in the case of a non-monotone response function with species decay rate added to the dilution rate. Armstrong and McGehee [18] considered model Eq. 8 extended to $n$ species competition in the case of monotone response functions. By ignoring the species decay rate, they were able to apply a conservation law to obtain a limiting system. They then studied the resulting limiting system, but did not apply the theory of asymptotically autonomous systems to obtain results for the full system. Butler and Wolkowicz [5] used a different method, provided a complete global analysis of this $n$ species model for both arbitrary monotone and non-monotone response functions, and applied results for asymptotically autonomous systems so that their results applied to the full system, not just the limiting system. They proved that competitive exclusion holds, i.e., all solutions approach an equilibrium that can be initial condition dependent in the non-montone case. In Wolkowicz and Lu [26], the decay rates were no longer ignored. There it was proved that for a large class of monotone and non-monotone response functions, again competitive exclusion holds and all populations approach equilibrium. However, in the case of non-monotone response they only considered the case when the species with the lowest break-even concentration also has its larger break-even concentration larger than the substrate input concentration. In the case of only one species, their method works for all monotone response functions, but for non-monotone response functions still requires the assumption that the larger break-even concentration is larger than the input concentration. In this paper we were able to eliminate this assumption, and hence complete the analysis for the model of growth in the basic chemostat.

Appendix A Proofs

Proof A.1**.**

of 1

i) Assume first that $X_{1}(0)=0$ and all other initial conditions are non-negative. It follows that $X_{1}(t)=0$, for all $t\geq 0$. Hence, Eq. 1 reduces to the system of first order differential equations

[TABLE]

Equations Eq. 20a and Eq. 20c imply that $S_{1}$ and $S_{3}$ converge exponentially to $S^{(0)}$ and [math], respectively. The hyperplane given by $S_{2}=0$ is invariant under Eq. 20b by (H6), and the hyperplane given by $X_{2}=0$ is invariant under Eq. 20d. By uniqueness of solutions to initial value problems, if $S_{2}(0)\geq 0$ and $X_{2}(0)\geq 0$, then $S_{2}(t)\geq 0$ and $X_{2}(t)\geq 0$ for all $t\geq 0$. Consider $\Sigma=S_{2}+y_{3}X_{2}$. Then $\dot{\Sigma}=-DS_{2}-y_{3}D_{2}X_{2}\leq-D\Sigma$ and thus $\Sigma(t)\to 0$ as $t\to 0$, implying $X_{2}(t)$ and $S_{2}(t)$ must each converge to 0 as $t\to\infty$.

ii) and iii) Assume that $X_{1}(0)>0$ and $X_{2}(0)\geq 0$, with all other initial conditions non-negative. Notice first that Eq. 1a and Eq. 1b decouple from the system. They describe a simple chemostat, for which it is known that if $X_{1}(0)>0$ and $S_{1}(0)\geq 0$, then $S_{1}(t)>0$ and $X_{1}(t)>0$ for all $t>0$ (e.g., see [20, 26, 11]). Note that the hyperplane $X_{2}=0$ is invariant under Eq. 20d, and so if $X_{2}(0)=0$, $X_{2}(t)=0$ for all $t\geq 0$, and if $X_{2}(0)>0$, then $X_{2}(t)>0$ for all $t\geq 0$. If $S_{3}(0)=0,$ then by Eq. 1d, $\dot{S}_{3}(0)>0$, and so there exists $\epsilon>0$ such that $S(t)>0$ for all $t\in(0,\epsilon)$. Let $S_{3}(0)\geq 0$. Suppose that there exists $\hat{t}>0$ such that $S_{3}(t)>0$ for all $t\in(0,\hat{t})$ and $S_{3}(\hat{t})=0$. Then, $\dot{S}_{3}(\hat{t})\leq 0$, but again from Eq. 1d, $\dot{S_{3}}(\hat{t})=y_{4}\mu_{1}(S_{1}(\hat{t}))X_{1}(\hat{t})>0$, a contradiction. Hence, $S_{3}(t)>0$ for all $t>0$. Using Eq. 1c, a similar argument applies to $S_{2}$.

iv) It is known (e.g., see [20, 26, 11]) that solutions to the simple chemostat (Eqs. 1a and 1b) are bounded. Hence, there exists $0<M<\infty$ such that $S_{1}(t)<M$ and $X_{1}(t)<M$ for all $t\geq 0$ . Thus, $S_{3}$ satisfies

[TABLE]

where $\tilde{M}=\mu_{1}(M)M$. This differential inequality implies that $S_{3}(t)\leq\frac{y_{4}\tilde{M}}{D}+S_{3}(0)e^{-Dt}$ for all $t\geq 0$, and thus $S_{3}(t)$ is bounded for $t>0$. Since $D_{i}\geq D$ the following differential inequality holds

[TABLE]

*Let $\Sigma=y_{3}X_{2}+S_{2}-\frac{y_{2}}{y_{4}}S_{3}$. Using Eq. 22, we see that $\dot{\Sigma}\leq-D\Sigma$, which implies that $\Sigma(t)\leq\Sigma(0)e^{-Dt}$. Since $S_{3}(t)$ is bounded above and we know that $S_{2}(t),X_{2}(t)\geq 0$ for all $t\geq 0$, they too must be bounded above. *

Proof A.2**.**

of 2

Since Eqs. 1a and 1b depend only on $S_{1}(t)$ and $X_{1}(t)$, these equations decouple from the full system Eq. 1, and it follows from known results on the basic model of the chemostat (e.g., see [20, 11] ) that if $\lambda_{1}\geq S^{(0)}$, then $(S_{1}(t),X_{1}(t))\to(S^{(0)},0)$ as $t\to\infty$.

*Therefore, for any $\epsilon>0$, there is a $T>0$ such that for $t>T$, $S_{1}(t)<S^{(0)}+\epsilon$ and $X_{1}(t)<\epsilon$. Then, for $t>T$, $\dot{S_{3}}(t)\leq-DS_{3}(t)+y_{4}\mu_{1}(S^{(0)}+\epsilon)\epsilon$, which gives $S_{3}(t)\leq S_{3}(T)e^{-Dt}+\frac{y_{4}}{D}\mu_{1}(S^{(0)}+\epsilon)\epsilon\left(1-e^{-D(t-T)}\right)$. Then, $\lim_{t\to\infty}S_{3}(t)=\frac{y_{4}}{D}\mu_{1}(S^{(0)}+\epsilon)\epsilon$. Since this holds for all $\epsilon>0$, letting $\epsilon\to 0$, gives $\lim_{t\to\infty}S_{3}(t)=0$. Next, let $\Sigma_{2}(t)=S_{2}(t)+\frac{1}{y_{3}}X_{2}(t)$. Since $D\leq D_{2}$, $\dot{\Sigma_{2}}(t)\leq-D\Sigma_{2}(t)+y_{2}\mu_{1}(S_{1}(t))X_{1}(t)$. The same argument as before proves that $\lim_{t\to\infty}\Sigma_{2}(t)=0$. Since for all $t$, $S_{2}(t)\geq 0$ and $X_{2}(t)\geq 0$, it follows that $\lim_{t\to\infty}S_{2}(t)=\lim_{t\to\infty}X_{2}(t)=0$. *

In the proof of 5 we rely on the fact that Eq. 1 is quasi-autonomous with the limiting system Eq. 8. For completeness, we include both the definition of quasi-autonomous, and Theorem 1.4 of [22].

Definition A.3**.**

Let $X$ be an open subset of $\mathbb{R}^{k}$, $k\geq 1$. A system

[TABLE]

with $x(t)\in X$ is called a quasi-autonomous system with limiting system

[TABLE]

if for any compact set $K\subset X$

[TABLE]

Theorem 6** (H. Thieme).**

Let

[TABLE]

be quasi-autonomous with limit system

[TABLE]

Let $x$ be a forward bounded solution of Eq. 26 that is defined for all forward times such that the closure of its forward orbit is contained in the open set $X\subset\mathbb{R}^{2}$. The following alternative holds for the $\omega$-limit set of $x$, $\omega$:

$\omega$* contains a periodic orbit of Eq. 27.* 2. 2.

$\omega$* contains at least one equilibrium of Eq. 27 and no periodic orbits of Eq. 27.*

To show that system Eq. 1 is a quasi-autonomous system with limiting system Eq. 7, we first prove a lemma.

Lemma 7**.**

Let $\dot{x}(t){=}f(t,x(t))$ be quasi-autonomous with limiting system $\dot{y}(t){=}$ $g(y(t))$ and assume that there exists $h(x(t))$ such that for all $K\subset X$ compact

[TABLE]

*Then $\dot{x}(t)=f(t,x(t))$ is quasi-autonomous with limiting system $\dot{y}(t)=h(y(t))$. *

Proof A.4**.**

By the triangle inequality,

[TABLE]

Proof A.5**.**

of 3

First we show that Eq. 1 is quasi-autonomous with limiting system:

[TABLE]

Since we are assuming that $\mu_{1}(S_{1})$ is a monotone response function, the results in [26] can be applied to the first two equations in Eq. 1 to prove that $(S_{1}(t),X_{1}(t))$ converge exponentially to $(\lambda_{1},X_{1}^{*})$ as $t\to\infty$. (The restriction that the results in [26] only apply to a general class of monotone response functions rather than any monotone response function does not apply to the single species growth model.)

Let $x(t)=(S_{1}(t),X_{1}(t),S_{2}(t),S_{3}(t),X_{2}(t))$ be any solution of Eq. 1, $K\subset\mathbb{R}^{5}_{+}$ be a compact set, and let $||\cdot||$ denote the Euclidean norm. For $t_{0}\geq 0$, consider

[TABLE]

where

[TABLE]

If $t_{0}=0$, then for any $0<t_{1}<\infty$, by continuity of the norm,

[TABLE]

Thus, we need only consider the case $t_{0}>0$. By the Cauchy-Schwartz inequality,

[TABLE]

The first integral, $\int_{t_{0}}^{\infty}\frac{1}{t^{2}}\,dt$, is finite. Since all of the terms of

[TABLE]

are positive, we can consider them individually. We begin with the second term,

[TABLE]

Since $\mu_{1}(S_{1})\in C^{1}$, by the Mean Value Theorem, for every $t>0$, there exists $\theta(t)$, such that $S_{1}(\theta(t))$ lies between $S_{1}(t)$ and $\lambda_{1}$. Let $M_{0}=\sup_{t\in[0,\infty)}|\mu_{1}^{\prime}(S_{1}(\theta(t))|>0$. Since $S_{1}(t)\to\lambda$ as $t\to\infty$, $\mu_{1}^{\prime}(S_{1}(\theta(t))$ remains bounded, $M_{0}$ is finite and $|-D_{1}+\mu_{1}(S_{1}(t))|=|-\mu_{1}(\lambda_{1})+\mu_{1}(S_{1}(t))|=|\mu_{1}^{\prime}(S_{1}(\theta(t)))||-\lambda_{1}+S_{1}(t)|\leq M_{0}|-\lambda_{1}+S_{1}(t)|\to 0,$ exponentially as $t\to 0$. Thus, there is a $k>0$, such that

[TABLE]

where $\overline{X_{1}}$ is the maximum value of $X_{1}(t)\in K$, and $\widetilde{M_{0}}=M_{0}|S_{1}(0)-\lambda_{1}|$.

We now consider the first term,

[TABLE]

By Young’s inequality ( i.e. that for any two real numbers, $a$ and $b$, $(a+b)^{2}\leq 2a^{2}+2b^{2}$) and using $D(S^{(0)}-\lambda_{1})=D_{1}y_{1}X_{1}^{*}$,

[TABLE]

Since this integral is a sum of positive terms we may consider each term individually. The first term is bounded above by the integral of a decaying exponential, and so is finite. We use Young’s inequality to bound the second term,

[TABLE]

where both of the terms in Eq. 33 are bounded above by a decaying exponential and so this integral is finite. For the third term in Eq. 32, write

[TABLE]

Noting that $y_{2}D_{1}X_{1}^{*}=D\lambda_{2}$, the exponential decay of $(X_{1}(t)-X_{1}^{*})^{2}$, and the same decay arguments as with the first term in Eq. 32. The finiteness of the fourth term in Eq. 32 follows from a similar idea, noting that $y_{4}D_{1}X_{1}^{*}=D\lambda_{3}$. Thus, Eq. 1 is quasi-autonomous with limiting system Eq. 29.

Now we finally show that Eq. 1 has limiting system Eq. 7. From Eq. 29b, if follows that

[TABLE]

We use this to argue that

[TABLE]

where, $Y_{5}(t)=\mu_{2}(S_{2}(t),\lambda_{3})-\mu_{2}(S_{2}(t),S_{3}(t))X_{2}(t)$, and $Y_{6}(t)=S_{3}(t)-\lambda_{3}$. The Cauchy-Schwartz inequality allows us to split the integral into more manageable pieces,

[TABLE]

By Eq. 34, the term containing $Y_{6}(t)$ is bounded above. In order to show the integral containing $Y_{5}(t)$ is bounded above we use the fact that $\mu_{2}(S_{2},S_{3})\in C^{1}$ and Eq. 34 to argue that there exists $M_{1}\geq 0$ such that

[TABLE]

Since $X_{2}(t)$ is bounded we have

[TABLE]

*Where $\overline{X_{2}}$ is the maximum value of $X(t)$ in $K$. The integral on the right is finite and therefore, by 7, Eq. 1 is quasi-autonomous with limiting system Eq. 7. *