The operating diagram for a two-step anaerobic digestion model

Abstract

The anaerobic digestion model ADM1 is a complex model which is widely accepted as a common platform for anaerobic process modeling and simulation. However, it has a large number of parameters and states that hinder its analytic study. Here, we consider the two-step simple model of anaerobic digestion named AM2, which is a four-dimensional system of ordinary differential equations. The AM2 model is able to adequately capture the main dynamical behavior of the full anaerobic digestion model ADM1 and has the advantage that a complete analysis for the existence and local stability of its steady states is available. We describe its operating diagram, which is the bifurcation diagram giving the behavior of the system with respect to the operating parameters, represented by the dilution rate and the input concentrations of the substrates. This diagram is very useful to understand the model from both the mathematical and biological points of view. It is shown that six types of behavior can be obtained for the long-term dynamics of the AM2 model, concerning the coexistence or extinction of one or both bacterial populations.

Appendices

Relationship to previous work

The two-step system (2) has been often considered in the literature. As it is usual in the mathematical theory of the chemostat, see for instance [22], in this type of models, we can use a change of variables that reduces the pseudo-stoichiometric coefficients $k_i$ to 1. Indeed, the linear change of variables

$$\begin{aligned} s_1=({k_2}/{k_1})S_1, \quad x_1=k_2X_1, \quad s_2=S_2, \quad x_2=k_3X_2, \end{aligned}$$

transforms (2) into

$$\begin{aligned} \begin{array}{lcl} \dot{s}_1 &{}=&{} D\left( s_{1\mathrm{in}}-s_1\right) -f_1\left( s_1\right) x_1,\\ \dot{x}_1 &{}=&{} \left( f_1\left( s_1\right) -\alpha D\right) x_1,\\ \dot{s}_2 &{}=&{} D\left( s_{2\mathrm{in}}-s_2\right) +f_1\left( s_1\right) x_1-f_2\left( s_2\right) x_2,\\ \dot{x}_2 &{}=&{} \left( f_2\left( s_2\right) -\alpha D\right) x_2, \end{array} \end{aligned}$$

(4)

where

$$\begin{aligned} s_{1\mathrm{in}}=({k_2}/{k_1})S_{1\mathrm{in}}, \quad s_{2\mathrm{in}}=S_{2\mathrm{in}}, \\ f_1(s_1)=\mu _1\left( ({k_1}/{k_2})s_1\right) , \quad f_2(s_2)=\mu _2(s_2). \end{aligned}$$

However, since the stoichiometric coefficients have their own importance for the biologist and since we are interested in giving these later a useful tool for the understanding of the role of the operating parameters, following [6], we do not make this reduction and we present the results in the original model (2). This model can have at most six steady states, labeled below as in [6]:

• $E_1^0$, where $X_1=0$ and $X_2=0$: the washout steady state where acidogenic and methanogenic bacteria are extinct.

• $E_1^i$ ($i=1,2$), where $X_1=0$ and $X_2>0$: acidogenic bacteria are washed out, while methanogenic bacteria are maintained.

• $E_2^0$, where $X_1>0$ and $X_2=0$: methanogenic bacteria are washed out, while acidogenic bacteria are maintained.

• $E_2^i$ ($i=1,2$), where $X_1>0$ and $X_2>0$: both acidogenic and methanogenic bacteria are maintained.

Table 10 Necessary and sufficient conditions of existence and local stability of the steady states of (2). obtained in [6]. $S_1^{*}(D)$, $S_2^{i*}(D)$ and ${S}^*_{2\mathrm{in}}\left( D,S_{1\mathrm{in}},S_{2\mathrm{in}}\right)$ are defined in Table 1

Table 11 The 9 cases of existence and stability of steady states of (2) obtained in [6], and the corresponding regions defined in Table 5

As shown in Proposition 1 of [6], the components of the steady states $E_1^0$, $E_1^i$ ($i=1,2$), $E_2^0$ and $E_2^i$ ($i=1,2$) are given in Table 2, where $S_1^{*}$, $S_2^{i*}$, $S_{2\mathrm{in}}^{*}$, $X_1^{*}$, $X_2^{i}$, and $X_2^{i*}$, for $i=1,2$, are defined in Table1. The necessary and sufficient conditions of existence of the steady states, given in Proposition 1 of [6], are summarized in the second column of Table 10. The necessary and sufficient conditions of local stability of these steady states, obtained in Table A.1 of [6], are summarized in the third column of Table 10.

Remark 3

In Table 10, since the function $S_1^*$ is defined on $(0,D_1)$, the condition $S_{1\mathrm{in}}>S_1^*(D)$ means $0<D<D_1$ and $S_{1\mathrm{in}}>S_1^*(D)$. Conversely, since by convention $S_1^*(D)=+\infty$ for $D\ge D_1$, the condition $S_{1\mathrm{in}}<S_1^*(D)$ means $D\ge D_1$ and $S_{1\mathrm{in}}>0$ or $0<D<D_1$ and $0<S_{1\mathrm{in}}<S_1^*(D)$. On the other hand, since the function $S_2^{i*}$ is defined on $(0,D_2)$, the condition $S_{2\mathrm{in}}>S_2^{i*}(D)$ means $0<D<D_2$ and $S_{2\mathrm{in}}>S_2^{i*}(D)$ and, conversely, since by convention $S_2^{1*}(D)=+\infty$ for $D> D_2$, the condition $S_{2\mathrm{in}}\notin \left[ S_2^{1*}(D),S_2^{2*}(D)\right]$ means $D\ge D_2$ and $S_{2\mathrm{in}}>0$ or $0<D< D_2$ and

$$\begin{aligned} S_{2\mathrm{in}}\notin \left[ S_2^{1*}(D),S_2^{2*}(D)\right] . \end{aligned}$$

Table 12 Auxiliary functions in the case given by (3)

The existence and stability conditions of the steady states of (2) given in Table 10 depend only on the relative positions of the values of $S_{1in}$ and $S_1^*(D)$ and of the values of

$$\begin{aligned} S_2^{1*}(D),~ S_2^{2*}(D),~ S_{2in} \text{ and } {S}^*_{2in}\left( D,S_{1in},S_{2in}\right) . \end{aligned}$$

Actually, as stated in Theorem 1 of [6], we can distinguish nine cases, according to the relative positions of these numbers. These cases are summarized in Table 11, together with the corresponding regions $\mathcal {I}_k$, $k=0,\ldots ,8$ of Table 5.

Remark 4

Note that Table 10 is identical to Table 3, except for the stability condition of $E_2^0$, and the existence conditions of $E_2^{i*}$, $i=1,2$, which are expressed in Table 3 using the $H_i$, $i=1,2$, functions, defined in Table 1.

Let us prove the following lemma which shows that the existence conditions of $E_2^{i*}$, $i=1,2$, given in Table 10, can be stated using the functions $H_i(D)$, defined in Table 1.

Lemma 1

The conditions ${S}^*_{2\mathrm{in}}\left( D,S_{1\mathrm{in}},S_{2\mathrm{in}}\right) =S_2^{i*}(D)$ and ${S}^*_{2\mathrm{in}}\left( D,S_{1\mathrm{in}},S_{2\mathrm{in}}\right) <S_2^{i*}(D)$, for $i=1,2$, are equivalent to the conditions $S_{2\mathrm{in}}+\frac{k_2}{k_1}S_{1\mathrm{in}}=H_{i}(D)$ and $S_{2\mathrm{in}}+\frac{k_2}{k_1}S_{1\mathrm{in}}<H_{i}(D)$, for $i=1,2$, respectively.

Proof

The result follows from the definitions of the functions ${S}^*_{2\mathrm{in}}\left( D,S_{1\mathrm{in}},S_{2\mathrm{in}}\right)$ and $H_i(D)$, given in Table 1. Indeed, the condition ${S}^*_{2\mathrm{in}}\left( D,S_{1\mathrm{in}},S_{2\mathrm{in}}\right) =S_2^{i*}(D)$ is equivalent to:

$$\begin{aligned} S_{2\mathrm{in}}+\frac{k_2}{k_1}\left( S_{1\mathrm{in}}-S_1^*(D)\right) =S_2^{i*}(D), \end{aligned}$$

which is itself equivalent to :

$$\begin{aligned} S_{2\mathrm{in}}+\frac{k_2}{k_1}S_{1\mathrm{in}}= S_2^{i*}(D)+\frac{k_2}{k_1}S_1^*(D). \end{aligned}$$

That is to say $S_{2\mathrm{in}}+\frac{k_2}{k_1}S_{1\mathrm{in}}= H_i(D)$. The proof for the inequality is the same. $\square$

The role of $H_i$-functions, in the description of the operating diagram, has already been highlighted, see Fig. 4 in [25], where cases $D_2<D_1$ and $D_1<D_2$ are distinguished.

Fig. 12

The curves $\varGamma _1$ (in blue), $\varGamma _2$ and $\varGamma _3$ (in green), $\varGamma _4$ and $\varGamma _5$ (in red) and $\varGamma _6$ (in black), corresponding to the operating diagrams in Fig. 5. a $S_{2\mathrm{in}}=0$, b $S_{2\mathrm{in}}=15$, c $S_{2\mathrm{in}}=S_2^M\simeq 36.332$ and d $S_{2\mathrm{in}}=100$. (Color figure online)

Proofs

1.1 Proof of Proposition 1

It is seen from Proposition 1 of [6] that the steady states are given by Table 2, where $S_1^{*}$, $S_2^{i*}$, $S_{2\mathrm{in}}^{*}$, $X_1^{*}$, $X_2^{i}$ and $X_2^{i*}$ are defined in Table 1. Their conditions of existence and stability are given in Table 10. Using Lemma 1 it is seen that the results in Table 10 are equivalent to those in Table 3 which completes the proof of Proposition 1.

1.2 Proof of Proposition 2

The cases 1.1, 1.2 and 1.3 correspond to the regions $\mathcal {I}_0$, $\mathcal {I}_1$ and $\mathcal {I}_2$, respectively, defined in Table 5. Now we use Lemma 1 to show that the remaining six cases 2.1 to 2.6 correspond to the six regions $\mathcal {I}_3$ to $\mathcal {I}_8$ defined in Table 5.

Since ${S}_{2\mathrm{in}}<{S}^*_{2\mathrm{in}}$ the case 2.1 corresponds to the condition ${S}^*_{2\mathrm{in}}<S_2^{1*}$ which is equivalent, using Lemma 1, to

$$\begin{aligned} S_{2\mathrm{in}}+\frac{k_2}{k_1}S_{1\mathrm{in}}<H_1(D). \end{aligned}$$

Therefore, the case 2.1 corresponds to the region $\mathcal {I}_3$ defined in Table 5. Using again Lemma 1, the condition $S_2^{1*}<{S}^*_{2\mathrm{in}}<S_2^{2*}$ in the case 2.2 is equivalent to

$$\begin{aligned} H_1(D)<S_{2\mathrm{in}}+\frac{k_2}{k_1}S_{1\mathrm{in}}<H_2(D) \end{aligned}$$

and the condition ${S}^*_{2\mathrm{in}}>S_2^{2*}$ in the case 2.3 is equivalent to

$$\begin{aligned} S_{2\mathrm{in}}+\frac{k_2}{k_1}S_{1\mathrm{in}}>H_2(D). \end{aligned}$$

Therefore, the cases 2.2 and 2.3 correspond to the regions $\mathcal {I}_4$ and $\mathcal {I}_5$, respectively, defined in Table 5. Using similar arguments, we show that the cases 2.4, 2.5 and 2.6 correspond to the regions $\mathcal {I}_6$, $\mathcal {I}_7$ and $\mathcal {I}_8$, respectively, defined in Table 5.

Table 13 Intersections of the $\varGamma _k$ surfaces, $k=0,\ldots ,8$ with a $\left( S_{1\mathrm{in}},S_{2\mathrm{in}}\right)$ plane, where D is kept constant

Table 14 The intersections of the $\varGamma _k$ surfaces, $k=0,\ldots ,8$ with a $\left( D,S_{1\mathrm{in}}\right)$ plane, where $S_{2\mathrm{in}}$ is kept constant

Excepted for cases 1.3, 2.3, 2.5 and 2.6 of bistability, the system (2) has a unique globally asymptotically stable (GAS) steady state. Therefore, in the case 1.1, $E_1^0$ is GAS; in the case 1.2, $E_1^1$ is GAS, in the case 2.1, $E_2^0$ is GAS, and in the cases 2.2 and 2.4, $E_2^1$ is GAS. In the case 1.3, $E_1^2$ is a saddle point whose attractive manifold is a three-dimensional hyper-surface surface which separates the phase space of (2) into the basins of attractions of the stable steady states $E_1^0$ and $E_1^1$. In the cases 2.3, 2.5 and 2.6, $E_2^2$ is a saddle point whose stable manifold is a three-dimensional hyper-surface which separates the phase space of (2) into the basins of attractions of the stable steady states $E_2^0$ and $E_2^1$. For details and complements on the global behavior, see section 2.4 of [6]. This completes the proof of Proposition 2.

1.3 Proof of Proposition 3

Part of the proof follows from [6]. It is seen from Theorem 1 of [6] that non-hyperbolic steady states, that correspond to coalescence of some of the steady state, occur when two (or more) of the values of $S_2^{1*}(D)$, $S_2^{2*}(D)$, $S_{2in}$, and ${S}^*_{2in}\left( D,S_{1in},S_{2in}\right)$ are equal. Notice that the condition

$$\begin{aligned} S_2^{1*}(D)=S_2^{2*}(D), \end{aligned}$$

arising in cases 1.6, 2.11 and 2.14 of Theorem 1 of [6], corresponds of the saddle node bifurcations of $E_1^1=E_1^2$ or $E_2^1=E_2^2$. This condition holds on $\varGamma _6$,

Table 15 Nominal parameters values used in [6] and corresponding to the figures

Fig. 13

The surfaces $\varGamma _1$ (in blue), $\varGamma _2$ and $\varGamma _3$ (in green), $\varGamma _4$ and $\varGamma _5$ (in red) and $\varGamma _6$ (in Yellow), corresponding to Fig. 1a. The surfaces separate the three-dimensional operating space $\left( D,S_{1\mathrm{in}},S_{2\mathrm{in}}\right)$ in 9 regions $\mathcal {I}_k$, $k=0,\ldots ,8$. Front (a), rear (b), left (c) and right (d) view of the surfaces $\varGamma _i$. Compare with Fig. 6 of [25]. (Color figure online)

Notice the condition $S_{2\mathrm{in}}=S_2^{1*}(D)$, arising in cases 1.4, 2.8 and 2.9 of Theorem 1 of [6], corresponds of the transcritical bifurcation $E_1^0=E_1^1$. This condition holds on $\varGamma _2$. Similarly, the condition $S_{2\mathrm{in}}=S_2^{2*}(D)$, arising in cases 1.5 and 2.13 of Theorem 1 of [6], corresponds of the transcritical bifurcation $E_1^0=E_1^2$. This condition holds on $\varGamma _3$.

On the other hand the condition $S^*_{2\mathrm{in}}=S_2^{1*}(D)$, arising in cases 2.7 of Theorem 1 of [6], corresponds of the transcritical bifurcation $E_2^0=E_2^1$. Using Lemma 1, this condition holds on $\varGamma _4$. Similarly, the condition $S^*_{2\mathrm{in}}=S_2^{2*}(D)$, arising in cases 2.12 and 2.15 of Theorem 1 of [6], corresponds of the transcritical bifurcation $E_2^0=E_2^2$. Using Lemma 1, this condition holds on $\varGamma _5$.

Finally we consider the bifurcations occurring when $S_{1\mathrm{in}}=S_1^{*}(D)$. These bifurcations were not considered in Theorem 1 of [6]. The condition $S_{1\mathrm{in}}=S_1^{*}(D)$ holds on $\varGamma _1$ and corresponds to the transcritical bifurcations $E_1^0=E_2^0$, $E_1^1=E_2^1$ and $E_1^2=E_2^2$. This completes the proof of Proposition 3.

Maple code

All plots in this paper were performed with Maple. For the convenience of the reader, we give hereafter the Maple instructions to plot Figs. 5, 6 and 7. Table 12 gives explicitly the functions used in the definitions of the $\varGamma _i$ curves in 14. The plots of these curves are obtained as follows.

The $\varGamma _i$ curves and the $\mathcal {I}_k$ regions they delimit are shown in Fig. 12. This figure is identical to Fig. 5, except that in Fig. 5; the regions $\mathcal {I}_k$ have been colored using the colors of Table 6.

Tables and three-dimensional operating diagram

In this section, we give several tables that are used in the paper. In Table 12, we present the functions defined in Table 1, in the particular case of the growth functions of Monod and Haldane (3). Tables 13 and 14 give the description of the intersections with a two-dimensional operating plane, where D or $S_{2in}$ is kept constant, respectively, of the $\varGamma _i$ surfaces that separate the operating parameter space in several regions, which are defined in Table 4. In Table 15, we provide the biological parameter values of the Monod and Haldane growth functions (3) used in the figures.

For the biological parameter values given in Table 15, and $m_1=0.6$, we give in Fig. 13 front, rear, left and right views of the surfaces $\varGamma _i$, $i=1,\ldots ,6$, in the three-dimensional operating space, showing the various regions of the three-dimensional operating diagram. In this three-dimensional view, the surfaces $\varGamma _i$ are colored as in Fig. 12, except that, for clarity, $\varGamma _6$ is colored yellow, rather than black.