Mathematical Study of a Two-Stage Anaerobic Model When the Hydrolysis Is the Limiting Step

Abstract

A two-step model of the anaerobic digestion process is mathematically and numerically studied. The focus of the paper is put on the hydrolysis and methanogenesis phases when applied to the digestion of waste with a high content of solid matter: existence and stability properties of the equilibrium points are investigated. The hydrolysis step is considered a limiting step in this process using the Contois growth function for the bacteria responsible for the first degradation step. The methanogenesis step being inhibited by the product of the first reaction (which is also the substrate for the second one), and the Haldane growth rate is used for the second reaction step. The operating diagrams with respect to the dilution rate and the input substrate concentrations are established and discussed.

1. Introduction

Two-step models are very common in environmental engineering literature to describe engineered biological systems. They are of particular interest to design feedback control laws since they are usually complex enough to capture the most important process dynamics while mathematically tractable.

The most common two-step model is used to describe the so-called ‘commensal ecological relationship’: it takes the form of a cascade of two biological reactions where one limiting substrate $S_1$ is consumed by one microorganism/ecosystem $X_1$ to produce $S_2$, which serves as the main limiting substrate for a second microorganism/ecosystem $X_2$, as schematically represented by the following reaction:

$$S_1\stackrel{{\mu }_1(.)}{⟶}{X}_1+S_2,S_2\stackrel{{\mu }_2(.)}{⟶}{X}_2+C{O}_2+C{H}_4$$

This model is given by the following general dynamical system:

$$\left\{\begin{array}{c}\hfill \begin{array}{c}\dot{S_1}=D(S_1^{in}-S_1)-{\mu }_1(.)\frac{X_1}{Y_1},\hfill \\ \dot{X_1}=[{\mu }_1(.)-\alpha D-k_1]X_1,\hfill \\ \dot{S_2}=D(S_2^{in}-S_2)+{\mu }_1(.)\frac{X_1}{Y_3}-{\mu }_2(.)\frac{X_2}{Y_2},\hfill \\ \dot{X_2}=[{\mu }_2(.)-\alpha D-k_2]X_2\hfill \end{array}\end{array}\right.$$

(1)

where D is the dilution rate, while $S_1^{in}$ and $S_2^{in}$ are the input substrate concentrations, respectively. Parameter $Y_i$ is yield the coefficient associated to the bio-reactions, $k_i$ is the mortality term, while $\alpha \in [0,1]$ is a term that allows decoupling the retention time applied to substrates (supposed to be soluble) and biomass (supposed to be particulate). The kinetics ${\mu }_1$ and ${\mu }_2$ are the growth rate functions associated to $X_1$ and $X_2$, respectively.

The different analyses of the class of models (1) available in the literature essentially differ on the growth rate functions used and whether a specific input for $S_2$ is considered or not (i.e., if there is a source term $S_2^{in}$ in the dynamic equation of $S_2$ or not). They differ also on the values (i) of the coefficient $\alpha$ (allowing to decouple the solid and liquid retention times) and (ii) of the mortality terms $k_i$. For details on the various models of this type considered in the literature, the reader can refer to Tables 2 and 3 of the review paper [1].

Among the numerous biological systems used in environmental engineering, the anaerobic systems are among those that are the most studied. Anaerobic digestion (AD) is the biological degradation of organic matter in the absence of oxygen. The final product is methane, a renewable energy source. This is why this process is used more and more for the treatment of liquid and solid wastes. Because of its relative instability due to the possible accumulation of intermediate products, notably the volatile fatty acids (VFA), the modeling of this process has been extensively studied over the last few years. Their complexity highly depends on the objectives pursued by the modeler. On the one hand, when the objective is to develop models for integrating and formalizing the available knowledge—typically to better understand bioprocesses—models are generally of high order and not tractable from a mathematical viewpoint, for instance, the ADM1 developed by the International Water Association or IWA [2]. On the other hand, as already mentioned above, when the aim of the modeling is to develop decision tools or control systems, low order models are better suited, as for instance, the well-known ‘AM2’, cf. [3]. Considering a system where the second step was limiting and subject to inhibition, the authors considered a Monod function for ${\mu }_1$ and a Haldane function for ${\mu }_2$. In [4], this model was studied for $\alpha =1$, while the most interesting case, where $0<\alpha <1$ and growth functions were characterized by qualitative properties, was studied in [5].

On the one hand, it is particularly important to underline here that such a model has been validated many times on experimental data and that it has been used in a large number of studies related to the AD of waste, notably as a basis for control design, cf. [3,6] or [7] and related references. On the other hand, using real data, several authors have shown that such low-order models are appropriate to describe the AD processes (cf. [8,9]). Following this idea, a method has even been proposed to easily, directly and systematically calibrate AM2 parameters from real data or data obtained by simulating a complex model, such as the ADM1, cf. [10].

Most of the waste treated in the anaerobic systems studied in the above-mentioned studies are liquid wastes. Depending on the nature of the wastes and, in particular, whether they are liquid or solid, the limiting step of the AD is not the same. If the treated waste is liquid, the main limiting step is usually considered to be the methanogenesis: in such a case, simple models including only acidogenesis and methanogenesis can be used as in the AM2. However, if the waste contains a high proportion of solid matter, it is the rule rather than the exception to consider that the hydrolysis is the main limiting step of the AD. In such a case, a model including only hydrolysis and methanogenesis can be used. As proposed in [11], when hydrolysis is the limiting step, rates depending only on substrate concentrations, such as the Monod function, are not the most appropriate. It is better to describe such complex phenomena by density-dependent models, such as density-dependent kinetics, a family in which the Contois models falls, cf. [12]. Using this model, the rate of the hydrolysis step is modeled as

$${\mu }_1(S_1,X_1)=\frac{m_1{S}_1}{K_1{X}_1+S_1}=\frac{m_1\frac{S_1}{X_1}}{K_1+\frac{S_1}{X_1}}$$

which exhibits the following properties specific to hydrolysis, cf. [12]:

$$\left\{\begin{array}{c}\hfill \begin{array}{c}\frac{S_1}{X_1}\gg K_1\Rightarrow {\mu }_1(S_1,X_1)X_1\approx m_1{X}_1\propto X_1,\hfill \\ \frac{S_1}{X_1}\ll K_1\Rightarrow {\mu }_1(S_1,X_1)X_1\approx \frac{m_1}{K_1}{S}_1\propto S_1\hfill \end{array}\end{array}\right.$$

(2)

While the analysis of the general model of AD initially proposed in [3] (representing acidogenesis and methanogenesis steps) has been realized in [5], to the best of authors knowledge, a two-step model where the kinetic of the first step is modeled by generic density-dependent kinetics and the second step exhibits a Haldane-type function has never been studied in the literature. It is the aim of this paper to study such a generic model. This analysis takes advantage of the fact that the system has a cascade structure: known results are then applied to study the whole fourth-order system as the coupling of two second-order chemostat models. The main contribution of the paper is the set of operating diagrams of the fourth-order system that is provided in Section 4.

The paper is organized as follows. In Section 2, the two-step model with two input substrate concentrations is presented, and the general hypotheses on the growth functions are given. In Section 3, the expressions of the steady states are given, and their stability properties are established. In Section 4, the effect of the second input substrate concentration on the steady states is illustrated in designing the operating diagrams first with respect to the first input substrate concentration and the dilution rate and second with respect to the second input substrate concentration and the dilution rate.

2. Mathematical Model

The two-step model reads:

$$\left\{\begin{array}{c}\hfill \begin{array}{c}\dot{S_1}=D(S_1^{in}-S_1)-{\mu }_1(S_1,X_1)\frac{X_1}{Y_1},\hfill \\ \dot{X_1}=[{\mu }_1(S_1,X_1)-D_1]X_1,\hfill \\ \dot{S_2}=D(S_2^{in}-S_2)+{\mu }_1(S_1,X_1)\frac{X_1}{Y_3}-{\mu }_2\left(S_2\right)\frac{X_2}{Y_2},\hfill \\ \dot{X_2}=[{\mu }_2\left(S_2\right)-D_2]X_2\hfill \end{array}\end{array}\right.$$

(3)

where $S_1$ and $S_2$ are the substrate concentrations introduced in the chemostat with input concentrations $S_1^{in}$ and $S_2^{in}$. $D_1=\alpha D+k_1$ and $D_2=\alpha D+k_2$ are the sink terms of biomass dynamics, where D is the dilution rate, $k_1$ and $k_2$ represent maintenance terms and parameter $\alpha \in [0,1]$ represents the fraction of the biomass affected by the dilution rate, while $Y_i$ is the yield coefficient. $X_1$ and $X_2$ are the hydrolytic bacteria and methanogenic bacteria concentrations, respectively. The functions ${\mu }_1:(S_1,X_1)\to {\mu }_1(S_1,X_1)$ and ${\mu }_2:\left(S_2\right)\to {\mu }_2\left(S_2\right)$ are the specific growth rates of the bacteria.

To ease the mathematical analysis of the system, it is rescaled. Notice that it is simply equivalent to changing the units of the variables:

$$s_1=S_1,x_1=\frac{1}{Y_1}{X}_1,s_2=\frac{Y_3}{Y_1}{S}_2,x_2=\frac{Y_3}{Y_1{Y}_2}{X}_2$$

The following system is obtained:

$$\left\{\begin{array}{c}\hfill \begin{array}{c}\dot{s_1}=D(s_1^{in}-s_1)-f_1(s_1,x_1)x_1,\hfill \\ \dot{x_1}=[f_1(s_1,x_1)-D_1]x_1,\hfill \\ \dot{s_2}=D(s_2^{in}-s_2)+f_1(s_1,x_1)x_1-f_2\left(s_2\right)x_2,\hfill \\ \dot{x_2}=[f_2\left(s_2\right)-D_2]x_2\hfill \end{array}\end{array}\right.$$

(4)

where $s_2^{in}=\frac{Y_3}{Y_1}{S}_2^{in}$, and $f_1$ and $f_2$ are defined by

$$f_1(s_1,x_1)={\mu }_1(s_1,Y_1{x}_1)and{f}_2\left(s_2\right)={\mu }_2\left(\frac{Y_1}{Y_3}{s}_2\right)$$

It is assumed that the functions ${\mu }_1(.,.)$ and ${\mu }_2(.)$ satisfy the following hypotheses.

Hypothesis 1

(H1). ${\mu }_1(s_1,x_1)$ is positive for $s_1>0$, $x_1>0$ and satisfies ${\mu }_1(0,x_1)=0$ and ${\mu }_1(+\infty ,x_1)=m_1\left(x_1\right)$. Moreover, ${\mu }_1(s_1,x_1)$ is strictly increasing in $s_1$ and decreasing in $x_1$, that is to say $\frac{\partial {\mu }_1}{\partial s_1}>0$ and $\frac{\partial {\mu }_1}{\partial x_1}\le 0$ for $s_1>0$, $x_1>0$.

Hypothesis 2

(H2). ${\mu }_2\left(s_2\right)$ is positive for $s_2>0$ and satisfies ${\mu }_2\left(0\right)=0$ and ${\mu }_2(+\infty )=0$. Moreover, ${\mu }_2\left(s_2\right)$ increases until a concentration $s_2^M$ and then decreases; thus, ${\mu }_2^{{}^{\prime }}\left(s_2\right)>0$ for $0\le s_2<s_2^M$, and ${\mu }_2^{{}^{\prime }}\left(s_2\right)<0$ for $s_2>s_2^M$.

As underlined in the introduction, particular kinetics models, such as the Contois function, verify H1, while the Haldane function verifies H2. Since the functions ${\mu }_1$ and ${\mu }_2$ satisfy the hypotheses H1 and H2, it follows from the above that functions $f_1$ and $f_2$ satisfy the following assumptions.

Assumption 1

(A1). $f_1(s_1,x_1)$ is positive for $s_1>0$, $x_1>0$ and satisfies $f_1(0,x_1)=0$ and $f_1(+\infty ,x_1)=m_1\left(x_1\right)$. Moreover, $\frac{\partial f_1}{\partial s_1}>0$ and $\frac{\partial f_1}{\partial x_1}\le 0$ for $s_1>0$, $x_1>0$.

Assumption 2

(A2). $f_2\left(s_2\right)$ is positive for $s_2>0$ and satisfies $f_2\left(0\right)=0$ and $f_2(+\infty )=0$. Moreover, $f_2\left(s_2\right)$ increases until a concentration of $s_2^M$ and then decreases, with $f_2^{{}^{\prime }}\left(s_2\right)>0$ for $0<s_2<s_2^M$ and $f_2^{{}^{\prime }}\left(s_2\right)<0$ for $s_2>s_2^M$.

3. Analysis of the Model

3.1. The Dynamics of $s_1$ and $x_1$

3.1.1. Study of the Steady States of System (5)

Model (4) has a cascade structure that renders its analysis easier. In other terms, $s_1$ and $x_1$ are not influenced by variables $s_2$ and $x_2$, and their dynamics are given by:

$$\left\{\begin{array}{c}\hfill \begin{array}{c}\dot{s_1}=D(s_1^{in}-s_1)-f_1(s_1,x_1)x_1,\hfill \\ \dot{x_1}=[f_1(s_1,x_1)-D_1]x_1.\hfill \end{array}\end{array}\right.$$

(5)

The behaviour of this system is well-known, cf. [13]. A steady state $(s_1^{*},x_1^{*})$ must be the solution of the system

$$\left\{\begin{array}{c}\hfill \begin{array}{c}0=D(s_1^{in}-s_1)-f_1(s_1,x_1)x_1,\hfill \\ 0=[f_1(s_1,x_1)-D_1]x_1\hfill \end{array}\end{array}\right.$$

(6)

From the second equation, it is deduced that $x_1^{*}=0$, which corresponds to the washout $E_0=(s_1^{in},0)$, or $s_1^{*}$ and $x_1^{*}$ must satisfy both equations

$$f_1(s_1^{*},x_1^{*})=D_{1}and{x}_1^{*}=\frac{D}{D_1}(s_1^{in}-s_1^{*}).$$

(7)

Let $\gamma$ a function defined by:

$$\gamma \left(s_1\right)=f_1\left(s_1,\frac{D}{D_1}(s_1^{in}-s_1)\right),$$

so $s_1^{*}$ is a solution of $\gamma \left(s_1\right)=D_1$, and it is noticed that ${\displaystyle {\gamma }^{{}^{\prime }}\left(s_1\right)=\frac{\partial f_1}{\partial s_1}-\frac{D}{D_1}\frac{\partial f_1}{\partial x_1}}$.

According to the hypothesis A1, $\gamma \left(s_1\right)$ is strictly increasing over the interval $[0,s_1^{in}]$, with $\gamma \left(0\right)=0$ and $\gamma \left(s_1^{in}\right)=f_1(s_1^{in},0)$. According to the theorem of intermediate values, the equation $\gamma \left(s_1\right)=D_1$ has a solution between 0 and $s_1^{in}$ if and only if $D_1<\gamma \left(s_1^{in}\right)$—that is, if $D_1<f_1(s_1^{in},0)$, cf. Figure 1.

Hence, for $x_1^{*}\ne 0$, the equilibrium $E_1(s_1^{*},x_1^{*})$ exists if and only if $D_1<f_1(s_1^{in},0)$.

The local stability of the steady state is given by the sign of the real part of eigenvalues of the Jacobian matrix evaluated at this steady state. In the following, the abbreviation LES for locally exponentially stable is used.

Proposition 1.

Assume that Assumptions A1 and A2 hold. Then, the local stability of steady states of System (5) is given by:

1. 

$E_0=(s_1^{in},0)$ is LES if and only if $f_1(s_1^{in},0)<D_1$ (i.e., $s_1^{in}<s_1^{*}$);

2. 

$E_1=(s_1^{*},x_1^{*})$ is LES if and only if $f_1(s_1^{in},0)>D_1$ (i.e., $s_1^{in}>s_1^{*}$), ($E_1$ is stable if it exists).

The reader may refer to [13] for the proof of this proposition. In the same book, notice that global stability results for System (5) are also provided. When $E_0$ and $E_1$ coincide, the equilibrium is attractive (the eigenvalues are equal to zero). The results of Proposition 1 are summarized in the following

Table 1. Summary of the results of Proposition 1.

Steady State	Existence Condition	Stability Condition
$E_0$	Always exists	$f_1(s_1^{in},0)<D_1$
$E_1$	$f_1(s_1^{in},0)>D_1$	Stable when it exists

.

3.1.2. Operating Diagram of the System (5)

Apart from the two operating (or control) parameters, which are the input substrate concentration $s_1^{in}$ and the dilution rate D, which can vary, all others parameters ($\alpha$, $k_1$ and the parameters of the growth function $f_1(s_1,x_1)$) have biological meaning and are fixed depending on the organisms and substrate considered. The operating diagram shows how the steady states of the system behave when the two control parameters $s_1^{in}$ and D are varied. The operating diagram for System (5) is shown in Figure 2. The condition $f_1(s_1^{in},0)>D_1$ of the existence of $E_1$ is equivalent to $D<\frac{1}{\alpha }[f_1(s_1^{in},0)-k_1]$. Therefore, the curve

$$\Gamma :\left\{(s_1^{in},D):D=\frac{1}{\alpha }[f_1(s_1^{in},0)-k_1]\right\}$$

separates the operating plan in two regions as defined in Figure 2.

The curve $\Gamma$ is the border that makes $E_0$ unstable, and at the same time, $E_1$ exists.

Table 2. Stability properties of steady states of system (5) in each region.

Region	Equ. ${\mathit{E}}_0$	Equ. ${\mathit{E}}_1$
$(s_1^{in},D)\in {\mathcal{R}}_0$	S
$(s_1^{in},D)\in {\mathcal{R}}_1$	U	S

indicates the stability properties of steady states of System (5) in each region where S and U read for LES and unstable, respectively, and no letter means that the steady state does not exist.

Except for small values of D and $s_1^{in}$, notice that the operating diagram of this first part of the two-step system under study is qualitatively similar to that one of the first part of the AM2 model, that is when a Monod-like growth function is considered, cf. [3].

3.2. The Dynamics of $s_2$ and $x_2$

3.2.1. Study of the Steady States of System (8)

Due to the cascade structure of System (4), the dynamics of the state variables $s_2$ and $x_2$ are given by

$$\left\{\begin{array}{c}\hfill \begin{array}{c}\dot{s_2}=D(F\left(t\right)-s_2)-f_2\left(s_2\right)x_2,\hfill \\ \dot{x_2}=[f_2\left(s_2\right)-D_2]x_2,\hfill \end{array}\end{array}\right.$$

(8)

with

$$F\left(t\right)=s_2^{in}+\frac{1}{D}{f}_1(s_1\left(t\right),x_1\left(t\right))x_1\left(t\right)$$

where $s_1\left(t\right),x_1\left(t\right)$ are a solution of System (5). A steady state $(s_1^{*},x_1^{*},s_2^{*},x_2^{*})$ of System (4) corresponds to a steady state $(s_2^{*},x_2^{*})$ of System (8) where either $(s_1\left(t\right),x_1\left(t\right))=E_0$ or $(s_1\left(t\right),x_1\left(t\right))=E_1$. Therefore, $(s_2^{*},x_2^{*})$ must be a steady state of the system

$$\left\{\begin{array}{c}\hfill \begin{array}{c}\dot{s_2}=D(s_2^{in*}-s_2)-f_2\left(s_2\right)x_2,\hfill \\ \dot{x_2}=[f_2\left(s_2\right)-D_2]x_2\hfill \end{array}\end{array}\right.$$

(9)

where

$$s_2^{in*}=s_2^{in}or{s}_2^{in*}=s_2^{in}+\frac{D_1}{D}{x}_1^{*}.$$

(10)

The first case corresponds to $(s_1^{*},x_1^{*})=E_0$ and the second to $(s_1^{*},x_1^{*})=E_1$.

System (9) corresponds to a classical chemostat model with Haldane-type kinetics, including a mortality term for $x_2$ and an input substrate concentration. Notice that $s_2^{in*}$, given by Equation (10), depends explicitly on the input flow rate. For a given D, the long-term behavior of such a system is well-known, cf. [13].

A steady state $(s_2^{*},x_2^{*})$ must be a solution of the system

$$\left\{\begin{array}{c}\hfill \begin{array}{c}0=D(s_2^{in*}-s_2)-f_2\left(s_2\right)x_2,\hfill \\ 0=[f_2\left(s_2\right)-D_2]x_2.\hfill \end{array}\end{array}\right.$$

(11)

From the second equation, it is deduced that $x_2^{*}=0$, which corresponds to the washout $F_0=(s_2^{in*},0)$, or $s_2^{*}$ must satisfy the equation

$$f_2\left(s_2\right)=D_2.$$

(12)

Under hypothesis A2, and if

$$D_2<f_2\left(s_2^M\right)$$

(13)

this equation has two solutions that satisfy $s_2^1<s_2^2$. Therefore, the system has two positive steady states $F_1=(s_2^1,x_2^{1*})$ and $F_2=(s_2^2,x_2^{2*})$, where

$$x_2^{i*}=\frac{D}{D_2}(s_2^{in*}-s_2^i),i=1,2.$$

(14)

For $i=1,2$, the steady states $F_i$ exist if and only if $s_2^{in*}>s_2^i$.

Proposition 2.

Assume that Assumptions A1, A2 and Condition (13) hold. Then, the local stability of steady states of System (9) is given by:

1. 

$F_0$ is LES if and only if $s_2^{in*}<s_2^1$ or $s_2^{in*}>s_2^2$;

2. 

$F_1$ is LES if and only if $s_2^{in*}>s_2^1$ (stable if it exists);

3. 

$F_2$ is unstable if it exists (unstable if $s_2^{in*}>s_2^2$).

The reader may refer to [13] for the proof of this proposition.

The results of Proposition 2 are summarized in the following

Table 3. Summary of the results of Proposition 2.

Steady-State	Existence Condition	Stability Condition
$F_0$	Always exists	$s_2^{in*}<s_2^1$ or $s_2^{in*}>s_2^2$
$F_1$	$s_2^{in*}>s_2^1$	Stable if it exists
$F_2$	$s_2^{in*}>s_2^2$	Unstable if it exists

.

3.2.2. Operating Diagram of the System (8)

Now, the operating diagram shows how the system behaves when the two control parameters $s_2^{in*}$ and D are varied. The operating diagram is shown in Figure 3. The conditions $s_2^{in*}=s_2^1$ or $s_2^{in*}=s_2^2$ are equivalent to $f_2\left(s_2^{in*}\right)=D_2$, that is to say $D=\frac{1}{\alpha }(f_2\left(s_2^{in*}\right)-k_2)$. Therefore, the horizontal line

$${\Gamma }_1:\left\{(s_2^{in*},D):D=\frac{1}{\alpha }(f_2\left(s_2^M\right)-k_2),s_2^{in*}>s_2^M\right\}$$

together with the curve

$${\Gamma }_2:\left\{(s_2^{in*},D):D=\frac{1}{\alpha }(f_2\left(s_2^{in*}\right)-k_2)\right\}$$

separate the operating diagram plane in three regions, as defined in Figure 3.

The following

Table 4. Stability properties of the steady states of System (8).

Region	Equ. ${\mathit{F}}_0$	Equ. ${\mathit{F}}_1$	Equ. ${\mathit{F}}_2$
$(s_2^{in*},D)\in {\mathcal{R}}_2$	S
$(s_2^{in*},D)\in {\mathcal{R}}_3$	U	S
$(s_2^{in*},D)\in {\mathcal{R}}_4$	S	S	U

summarizes the stability properties of steady states of System (8).

3.3. Analysis of the Whole System (4)

3.3.1. Steady States

The aim of this section is to study the dependence of the steady state of System (4) with respect to the operating parameters D, $s_1^{in}$ and $s_2^{in}$. Let $(s_1^{*},x_1^{*},s_2^{*},x_2^{*})$ be a steady state of System (4), then $(s_1^{*},x_1^{*})$ is a steady state of System (5) and $(s_2^{*},x_2^{*})$ is a steady state of System (9) where $s_2^{in*}$ is given by Equation (10).

If $(s_1^{*},x_1^{*})=E_0=(s_1^{in},0)$ then $s_2^{in*}=s_2^{in}$ and three possibilities can occur

• $(s_2^{*},x_2^{*})=(s_2^{in},0)$, and $E_1^0:=(s_1^{in},0,s_2^{in},0)$;
• $(s_2^{*},x_2^{*})=(s_2^1,x_2^1)$, and $E_1^1:=(s_1^{in},0,s_2^1,x_2^1)$;
• $(s_2^{*},x_2^{*})=(s_2^2,x_2^2)$, and $E_1^2:=(s_1^{in},0,s_2^2,x_2^2)$.

If $(s_1^{*},x_1^{*})=E_1=(s_1^{*},x_1^{*})$ then three others possibilities can occur

• $(s_2^{*},x_2^{*})=(s_2^{in*},0)$, and $E_2^0:=(s_1^{*},x_1^{*},s_2^{in*},0)$;
• $(s_2^{*},x_2^{*})=(s_2^1,x_2^{1*})$, and $E_2^1:=(s_1^{*},x_1^{*},s_2^1,x_2^{1*})$;
• $(s_2^{*},x_2^{*})=(s_2^2,x_2^{2*})$, and $E_2^2:=(s_1^{*},x_1^{*},s_2^2,x_2^{2*})$.

These results are summarized in the following proposition.

Proposition 3.

System (4) has, at most, six steady states:

• $E_1^0=(s_1^{in},0,s_2^{in},0)$ always exists;
• $E_1^1=(s_1^{in},0,s_2^1,x_2^1)$ exists if and only if $s_2^{in}>s_2^1$;
• $E_1^2=(s_1^{in},0,s_2^2,x_2^2)$ exists if and only if $s_2^{in}>s_2^2$;
• $E_2^0=(s_1^{*},x_1^{*},s_2^{in*},0)$ exists if and only if $f_1(s_1^{in},0)>D_1$;
• $E_2^1=(s_1^{*},x_1^{*},s_2^1,x_2^{1*})$ exists if and only if $f_1(s_1^{in},0)>D_1$ and $s_2^{in*}>s_2^1$;
• $E_2^2=(s_1^{*},x_1^{*},s_2^2,x_2^{2*})$ exists if and only if $f_1(s_1^{in},0)>D_1$ and $s_2^{in*}>s_2^2$.

3.3.2. Steady States Stability of the System (4)

In this section, the stability of the steady states given in Proposition 3 is studied. For this, the following Jacobian matrix is considered:

$$J=\left[\begin{array}{cc}{J}_{11}& J_{12}\\ 0& J_{22}\end{array}\right],$$

where $J_{11}$ and $J_{22}$ are defined by Equations (15) and (16), respectively, given by:

$$J=\left[\begin{array}{cc}-D-M{x}_1& -N{x}_1-f_1(s_1,x_1)\\ M{x}_1& N{x}_1+[f_1(s_1,x_1)-D_1]\end{array}\right],$$

(15)

and

$$J=\left[\begin{array}{cc}-D-f_2^{{}^{\prime }}\left(s_2\right)x_2& -f_2\left(s_2\right)\\ f_2^{{}^{\prime }}\left(s_2\right)x_2& f_2\left(s_2\right)-D_2\end{array}\right].$$

(16)

This matrix has a block-triangular structure. Hence, the eigenvalues of J are the eigenvalues of $J_{11}$ and the eigenvalues of $J_{22}$. The existence of steady states depends only on the relative positions of the two numbers $s_1^{in}$ and $s_1^{*}$ defined by Equation (7) with respect to the four numbers $s_2^1$ and $s_2^2$, defined by Equation (12) on the one hand, and $s_2^{in}$ and $s_2^{in*}$, defined by Equation (10) on the other hand. Equilibrium stability is summarized in

Table 5. The stability conditions for the steady states of System (4).

Equ.	Matrices ${\mathit{J}}_{11}$ and ${\mathit{J}}_{22}$	Conditions of Stability
$E_1^0$		$Tr\left(J_{11}\right)<0$ if $f_1(s_1^{in},0)<D_1$,
	$J_{11}=\left[\begin{array}{cc}-D& -f_1(s_1^{in},0)\\ 0& f_1(s_1^{in},0)-D_1\end{array}\right]$	$Tr\left(J_{22}\right)<0$ if $s_2^{in}<s_2^1$ or $s_2^{in}>s_2^2$,
		$det\left(J_{11}\right)>0$ and $det\left(J_{22}\right)>0$
	$J_{22}=\left[\begin{array}{cc}-D& -f_2\left(s_2^{in}\right)\\ 0& f_2\left(s_2^{in}\right)-D_2\end{array}\right]$	$\Rightarrow \left\{\begin{array}{c}{E}_1^{0}isstableif{f}_1(s_1^{in},0)\le D_{1}and\hfill \\ (s_2^{in}<s_2^{1}or{s}_2^{in}>s_2^2),\hfill & \hfill \\ E_1^{0}isunstableif{f}_1\left(s_1^{in}\right)>D_{1}or\hfill \\ s_2^1<s_2^{in}<s_2^2).\hfill & \hfill \end{array}\right.$
$E_1^i$	($i=1,2$)	$Tr\left(J_{11}\right)<0$ and $det\left(J_{11}\right)>0$ if $f_1(s_1^{in},0)<D_1$
	$J_{11}=\left[\begin{array}{cc}-D& -f_1(s_1^{in},0)\\ 0& f_1(s_1^{in},0)-D_1\end{array}\right]$	$Tr\left(J_{22}\right)<0$ and $det\left(J_{22}\right)>0$ at $E_1^1$,
		$det\left(J_{22}\right)<0$ at $E_1^2$
	$J_{22}=\left[\begin{array}{cc}-[D+f_2^{{}^{\prime }}\left(s_2^i\right)x_2^i]& -D_2\\ f_2^{{}^{\prime }}\left(s_2^i\right)x_2^i& 0\end{array}\right]$	$\Rightarrow \left\{\begin{array}{cc}{E}_1^{1}isstable,\hfill & \hfill \\ E_1^{2}isunstable\hfill & \hfill \\ E_1^{i}arebothunstableif{s}_1^{in}>s_1^{*}.\hfill & \hfill \end{array}\right.$
$E_2^0$		$Tr\left(J_{11}\right)<0$ and $det\left(J_{11}\right)>0$ by A1
	$J_{11}=\left[\begin{array}{cc}-[D+\left(\frac{\partial f_1}{\partial s_1}\right)x_1^{*}]& -[D_1+\left(\frac{\partial f_1}{\partial x_1}\right)x_1^{*}]\\ \left(\frac{\partial f_1}{\partial s_1}\right)x_1^{*}& \left(\frac{\partial f_1}{\partial x_1}\right)x_1^{*}\end{array}\right]$	$Tr\left(J_{22}\right)<0$ and $det\left(J_{22}\right)>0$
		if $s_2^{in*}<s_2^1$ or $s_2^{in*}>s_2^2$
	$J_{22}=\left[\begin{array}{cc}-D& -f_2\left(s_2^{in*}\right)\\ 0& [f_2\left(s_2^{in*}\right)-D_2]\end{array}\right]$	$\Rightarrow \left\{\begin{array}{c}{E}_2^{0}isstableif{s}_1^{in}\ge s_1^{*}and\hfill \\ (s_2^{in*}<s_2^{1}or{s}_2^{in*}>s_2^2),\hfill & \hfill \\ E_2^{0}isunstableif{s}_1^{in}>s_1^{*}and\hfill \\ s_2^1<s_2^{in*}<s_2^2.\hfill & \hfill \end{array}\right.$
$E_2^i$	($i=1,2$)	$Tr\left(J_{11}\right)<0$ and $det\left(J_{11}\right)>0$
	$J_{11}=\left[\begin{array}{cc}-[D+\left(\frac{\partial f_1}{\partial s_1}\right)x_1^{*}]& -[D_1+\left(\frac{\partial f_1}{\partial x_1}\right)x_1^{*}]\\ \left(\frac{\partial f_1}{\partial s_1}\right)x_1^{*}& \left(\frac{\partial f_1}{\partial x_1}\right)x_1^{*}\end{array}\right]$	$Tr\left(J_{22}\right)<0$ and $det\left(J_{22}\right)>0$ at $E_2^1$,
		$det\left(J_{22}\right)<0$ at $E_2^2$
	$J_{22}=\left[\begin{array}{cc}-[D+f_2^{{}^{\prime }}\left(s_2^i\right)x_2^{i*}]& -D_2\\ f_2^{{}^{\prime }}\left(s_2^i\right)x_2^{i*}& 0\end{array}\right]$	$\Rightarrow \left\{\begin{array}{cc}{E}_2^{1}isstable,\hfill & \hfill \\ E_2^{2}isunstable\hfill & \hfill \end{array}\right.$

, while the different regions of the operating diagram are synthesized in

Table 6. The three cases when $s_1^{in}<s_1^{*}$.

Case	Area	Condition	${\mathit{E}}_1^0$	${\mathit{E}}_1^1$	${\mathit{E}}_1^2$	${\mathit{E}}_2^0$	${\mathit{E}}_2^1$	${\mathit{E}}_2^2$
1.1	${\mathcal{A}}_1$	$s_2^{in}<s_2^1<s_2^2$	S
1.2	${\mathcal{A}}_2$	$s_2^1<s_2^{in}<s_2^2$	U	S
1.3	${\mathcal{A}}_3$	$s_2^1<s_2^2<s_2^{in}$	S	S	U

and

Table 7. The six cases when $s_1^{in}>s_1^{*}$.

Case	Area	Condition	${\mathit{E}}_1^0$	${\mathit{E}}_1^1$	${\mathit{E}}_1^2$	${\mathit{E}}_2^0$	${\mathit{E}}_2^1$	${\mathit{E}}_2^2$
2.1	${\mathcal{A}}_4$	$s_2^{in}<s_2^{in*}<s_2^1<s_2^2$	U			S
2.2	${\mathcal{A}}_5$	$s_2^{in}<s_2^1<s_2^{in*}<s_2^2$	U			U	S
2.3	${\mathcal{A}}_6$	$s_2^{in}<s_2^1<s_2^2<s_2^{in*}$	U			S	S	U
2.4	${\mathcal{A}}_7$	$s_2^1<s_2^{in}<s_2^{in*}<s_2^2$	U	U		U	S
2.5	$A_8$	$s_2^1<s_2^{in}<s_2^2<s_2^{in*}$	U	U		S	S	U
2.6	${\mathcal{A}}_9$	$s_2^1<s_2^2<s_2^{in}<s_2^{in*}$	U	U	U	S	S	U

.

Remark 1.

Here, the limit values in the stability condition (Ex: $s_2^{in}=s_2^1$ or $s_2^{in}=s_2^2$) are excluded. These limit cases are related to the eigenvalues of the Jacobian matrix with a real part equal to 0, in which case the corresponding states are named non-hyperbolic stationary states. Otherwise, they are named hyperbolic steady states.

The different possible cases of non-hyperbolic (NH) equilibrium are summarized in the Theorem 1.

Theorem 1.

If $s_1^{in}<s_1^{*}$, then there are three sub-cases:

Case	Condition	NH	S	U
1.4	$s_2^{in}=s_2^1<s_2^2$	$E_1^0=E_1^1$
1.5	$s_2^1<s_2^{in}=s_2^2$	$E_1^0=E_1^2$	$E_1^1$
1.6	$s_2^1=s_2^2<s_2^{in}$	$E_1^1=E_1^2$	$E_1^0$

If $s_1^{in}>s_1^{*}$, then there are nine sub-cases:

Case	Condition	NH	S	U
2.7	$s_2^{in}<s_2^{in*}=s_2^1<s_2^2$	$E_2^0=E_2^1$		$E_1^0$
2.8	$s_2^{in}<s_2^1<s_2^2=s_2^{in*}$	$E_2^0=E_2^2$	$E_2^1$	$E_1^0$
2.9	$s_2^{in}=s_2^1<s_2^{in*}<s_2^2$	$E_1^0=E_1^1$	$E_2^1$	$E_2^0$
2.10	$s_2^{in}=s_2^1<s_2^2<s_2^{in*}$	$E_1^0=E_1^1$	$E_2^1$,$E_2^0$	$E_2^2$
2.11	$s_2^{in}=s_2^1<s_2^2=s_2^{in*}$	$E_1^0=E_1^1$,$E_2^0=E_2^1$	$E_2^1$
2.12	$s_2^1<s_2^{in}=s_2^2<s_2^{in*}$	$E_2^1=E_2^2$	$E_2^0$	$E_1^0$
2.13	$s_2^1<s_2^{in}<s_2^2=s_2^{in*}$	$E_2^0=E_2^2$	$E_2^1$	$E_1^0$,$E_1^1$
2.14	$s_2^1<s_2^{in}=s_2^2<s_2^{in*}$	$E_1^0=E_1^2$	$E_2^0$,$E_2^1$	$E_1^0$,$E_2^2$
2.15	$s_2^1=s_2^2<s_2^{in}<s_2^{in*}$	$E_1^1=E_1^2$,$E_2^1=E_2^2$	$E_2^0$	$E_1^0$

Proof. 

Let us give the details of the proof for case 2.9. The other cases can be studied similarly. Assume that $s_2^{in}=s_2^1<s_2^{in*}<s_2^2$, then $x_2^1=0$, $x_1^{*}>0$ and $x_2^{1*}>0$ by Equation (14). Therefore (cf. Proposition 3), the system has three equilibria $E_1^0=E_1^1$, $E_2^0$ and $E_2^1$. Using the linearization, it is established that $E_2^0$ and $E_2^1$ are hyperbolic. □

Remark 2.

In each case, among those cited in the Table 6 and Table 7 and in the previous Theorem 1, a corresponding figure that represents the relative position of $s_2^1$, $s_2^2$, $s_2^{in}$ and $s_2^{in*}$ [5] can be associated.

4. Simulations

To illustrate these results, the operating diagrams of System (3) under Hypotheses H1 and H2 in a number of situations are plotted. Recall that the operating diagram summarizes the existence and the nature of the the steady states of a dynamical system as a function of its input variables. Here, the control inputs are D, $S_1^{in}$ and $S_2^{in}$. More particularly, either the operating diagrams in the plan $\left\{S_1^{in},D\right\}$ for a fixed value of $S_2^{in}$ or in the plan $\left\{S_2^{in},D\right\}$ for a fixed value of $S_1^{in}$ are plotted. All simulations are performed with the following growth functions:

$${\mu }_1(S_1,X_1)=\frac{m_1{S}_1}{K_1{X}_1+S_1},{\mu }_2\left(S_2\right)=\frac{m_2{S}_2}{\frac{S_2^2}{I}+S_2+K_2}.$$

Of course, these operating diagrams depend on the model parameters. The choice of their values is a difficult task. Indeed, our objective here is not to match a set of data but rather to highlight interesting qualitative properties of the model under interest. To do so, most parameters are taken from [3], while others are changed significantly as the inhibition coefficient I of the Haldane function. Indeed, as underlined in [14], the inhibition of the second reaction is not visible if the original parameters proposed in [3] are used considering reasonable ranges of variations for $S_1$ and $S_2$. With respect to this later, the Haldane parameter I was thus significantly decreased to willingly increase the inhibition effect of $S_2$ on the growth of $X_2$. Finally, the parameter values used are summarized in the following

Table 8. Parameters values for the simulations.

Parameter	Unit	Nominal Value
$m_1$	d${}^{-1}$	0.5
$K_1$	g/L	2.1
$m_2$	d${}^{-1}$	1
I	mmol/L	60
$K_2$	mmol/L	24
$k_1$	d${}^{-1}$	0.1
$k_2$	d${}^{-1}$	0.06
$\alpha$	in $[0,1)$	0.5
$Y_1$	g/g	1/25
$Y_2$	g/mmol	1/250
$Y_3$	g/mmol	1/268

:

To compute the different regions of the operating diagrams, the numerical method reported in [14] is used. The algorithm is recalled hereafter (emphcf. Algorithm 1).

4.1. Algorithm for the Determination of the Operating Diagrams

The algorithm is as follows: for each value of input variables chosen on a grid, the equilibria are computed. The eigenvalues of the Jacobian matrix are then calculated for each equilibrium. Finally, according to the conditions of existence and the sign of the real parts of the eigenvalues, a ‘flag’ is assigned to each of the 6 equilibria: ‘S’ for stable, ‘U’ for unstable or nothing if the steady state does not exist. This procedure stops when all the values of the grid {$S^{in},D$} have been scanned. As a result, a number of ‘signatures’ composed of sequences of ‘S’, ‘U’ or ‘nothing’ are obtained. These cases code for the existence and stability of the equilibria that are grouped into regions, as summarized in the tables at the end of Section 3.1 and Section 3.2, respectively. This algorithm may be formalized as follows: let $N_1,N_2$ be two integers in ${\mathbb{N}}^{*}$ and $h_1=\frac{D}{N_1}$ and $h_2=\frac{S_{in}}{N_2}$ the two iteration steps:

Algorithm 1 Operating diagram

for i varying from 1 to $N_1$ do;

for j varying from 1 to $N_2$;

determine 6 equilibria of the model $E_1...E_6$

for k varying from 1 to 6 do

calculate the Jacobian matrix at $E_k$ $\left(J_{E_k}\right)$

calculate the eigenvalues of $\left(J_{E_k}\right)$

if all the conditions of existence of $E_k$ are fulfilled and all real parts of

the eigenvalues of $\left(J_{E_k}\right)$ are non-positive then $E_k$ is stable

else if all conditions of existence of $E_k$ are fulfilled and at least one

real part of eigenvalue of $\left(J_{E_k}\right)$ is positive then $E_k$ is unstable

else $E_k$ does not exist

end if

end for (k)

end for (j)

end for (i)

4.2. Operating Diagrams

In this part, the results are illustrated by plotting the operating diagrams and are discussed.

Figure 4 represents the operating diagram of Model (1) in the plan $\left\{S_1^{in},D\right\}$ for $S_2^{in}=1.5$ mmol/L. The regions are defined as follows. ${\mathcal{A}}_1$ (in green) is the stability region of the washout $E_1^0$, ${\mathcal{A}}_5$ (in blue) is the stability region of steady-state $E_2^1$, ${\mathcal{A}}_6$ (in red) is the bi-stability region of the steady-states $E_2^0$ and $E_2^1$, and ${\mathcal{A}}_7$ (in dark blue) is the stability region of steady-state $E_2^1$, the difference between the regions ${\mathcal{A}}_5$ and ${\mathcal{A}}_7$ being that the equilibrium $E_1^1$ does not exist in the region ${\mathcal{A}}_5$ but exists and is unstable in ${\mathcal{A}}_7$, cf. Table 6 and Table 7.

Figure 5 represents the operating diagram of Model (1) in the plan $\left\{S_2^{in},D\right\}$ for $S_1^{in}=0.8$ g/L. The regions are defined as follows. ${\mathcal{A}}_2$ (in yellow) is the stability region of steady-state $E_1^1$, ${\mathcal{A}}_3$ (in orange) is the bi-stability region of the washout $E_1^0$ and the steady-state $E_1^1$, ${\mathcal{A}}_8$ and ${\mathcal{A}}_9$ (in pink) and (in dark pink), respectively, are the bi-stability regions of the steady-states $E_2^0$ and $E_2^1$, the difference between ${\mathcal{A}}_8$ and ${\mathcal{A}}_9$ being that the equilibrium $E_1^2$ does not exist in the region ${\mathcal{A}}_8$ but exists and is unstable in ${\mathcal{A}}_9$, cf. Table 6 and Table 7.

Figure 6 represents the operating diagram of Model (1) in the plan $\left\{S_2^{in},D\right\}$ for $S_1^{in}=0.03$ g/L, which is a smaller value of $S_1^{in}$ than before. The differences with the previous case are (i) the appearance of a little region ${\mathcal{A}}_4$ (in dark orange) that is the stability region of steady-state $E_2^0$ and (ii) a sharp decrease in the size of region ${\mathcal{A}}_8$ that almost disappears (it is reduced to a very narrow surface along the frontier with the region ${\mathcal{A}}_9$, as can seen in Figure 6).

The region ${\mathcal{A}}_8$ becomes very small and narrow compared to Figure 5 because it is numerically linked to the value of $S_1^{in}$. Indeed, when reducing the value of $S_1^{in}$, the region ${\mathcal{A}}_4$ appears, while the size of the region ${\mathcal{A}}_8$ is becoming smaller and smaller (compare Figure 5 and Figure 6). In other words, in decreasing $S_1^{in}$, the attraction basin of the positive stable equilibrium of ${\mathcal{A}}_7$ ($E_2^1$) increases. It is equivalent to say that, given two values of $S_1^{in}$, say $S_1^{in1}$ and $S_1^{in2}$, where $S_1^{in2}>S_1^{in1}$, a greater dilution rate is needed to destabilize the process if $S_1^{in}=S_1^{in2}$ than if $S_1^{in}=S_1^{in1}$; thus, the shard decrease in ${\mathcal{A}}_8$ observed in Figure 6.

4.3. Practical Interpretations of the Operating Diagrams

Below, it is explained how the operating diagrams may be used in practice. First, notice that the operator must avoid the process to operate in regions where either $E_1^0$ and/or $E_2^0$ would be the only stable equilibrium points. Indeed, in such regions, $x_2^{*}=0$ and the process does not produce any methane, then ${\mathcal{A}}_1$ and ${\mathcal{A}}_4$ must be avoided. In addition, a particular attention must be paid to operating conditions in which there is bi-stability with one of the stable equilibrium point is $E_1^0$ or $E_2^0$, which are ${\mathcal{A}}_3$, ${\mathcal{A}}_6$, ${\mathcal{A}}_8$ and ${\mathcal{A}}_9$. In the following, what may happen within regions ${\mathcal{A}}_6$ and ${\mathcal{A}}_7$ that are methane-producing areas is analyzed. To illustrate their practical interest, let us consider the operating diagram pictured in Figure 4 (that is for $S_2^{in}=1.5$ g/L) and let us browse it for increasing values of D at a fixed value of $S_1^{in}$.

Example 1.

Let $S_1^{in}=18g/L$.

In such a situation, the following regions are browsed in considering successive equilibrium when increasing D: ${\mathcal{A}}_7⤏{\mathcal{A}}_6⤏{\mathcal{A}}_5⤏{\mathcal{A}}_1$, cf. Figure 7. To better interpret whose steady states the system passes through, the bifurcation diagram is plotted in Figure 8.

This last diagram allows us to see the appearance/disappearance of steady states as a function of the input variable D (recall that $S_1^{in}$ and $S_2^{in}$ are fixed). As long as D is small enough (i.e., such that $\alpha D+k_2<d_1$), the quantity of substrates entering the second step of the reaction is very important: the system is in the region ${\mathcal{A}}_7$ where the positive equilibrium is the only stable equilibrium. As D increases, the size of the attraction basin of this equilibrium decreases until D reaches a critical value (corresponding to the point $P_1$ in Figure 7, which is situated on the frontier between regions ${\mathcal{A}}_7$ and ${\mathcal{A}}_6$). This critical value corresponds to the one for which the term $S_2^{in*}$ becomes exactly the largest solution of the equation ${\mu }_2\left(S_2\right)=D_2$ (equivalent to Equation (12) for System (4)): the system enters the region ${\mathcal{A}}_6$. From a biological point of view, the interpretation is as follows: as D increases, $X_1^{*}$ decreases, and thus, $S_2^{in*}$ decreases, as can be seen from Equation (10). When $D_2=d_2$, the quantity of available resources necessary for $X_2$ to grow may become limiting for some initial conditions, leading the system to enter a bi-stability zone. With the values of the parameters chosen, further increasing D leads $X_2$ to the washout: the system enters ${\mathcal{A}}_5$ when crossing the point $P_2$ of Figure 7. Finally, if D is such that $D_1=d_2$ (the critical value corresponding to the maximum growth rate of $X_1$), $X_1$ also goes extinct, and the system enters into ${\mathcal{A}}_1$.

Example 2.

Let $S_1^{in}=14$ g/L.

This case is even more interesting since, when D increases, the system goes back to ${\mathcal{A}}_5$ once before leaving it indefinitely to browse the following regions: ${\mathcal{A}}_7⤏{\mathcal{A}}_6⤏{\mathcal{A}}_7⤏{\mathcal{A}}_5⤏{\mathcal{A}}_1$. While D is small enough (i.e., such that $\alpha D+k_2<d_3$, cf. Figure 9), the reasoning remains the same as before. The only difference is that the value of D leading the system to enter into ${\mathcal{A}}_6$ through $P_3$ is a little bit higher than in the previous case ($d_3>d_1$). It is due to the fact that the second step of the reaction receives less input from the first step when compared to the case where $S_1^{in}=18g/L$, thus enlarging the attraction basin of the stable positive equilibrium. Then, when D is further increased, an interesting phenomenon may happen: the system enters back into ${\mathcal{A}}_7$ through point $P_4$ instead of entering ${\mathcal{A}}_5$ as it did in the case before. In fact, this strongly depends on model parameters and, in particular, on the relative rate at which $S_2^{in*}$ and the largest solution of the Equation (10) vary as functions of D, cf. Figure 10. In other words, it depends on how the input concentration of the second step $S_2^{in*}$, which includes the part of $S_1$ transformed into $S_2$ during the first reaction, is affected by D. On the one hand, if the system is in a ‘flat’ zone of the growth rate ${\mu }_2$ assuming concentrations at the right of the maximum of ${\mu }_2$ are considered, a small variation of D (and thus of $D_2$) will change the largest solution of Equation (10) a lot while $S_2^{in*}$ may almost remain constant. On the other hand, if D is such that $D_2$ crosses ${\mu }_2$ in a sharper zone of the Haldane function (typically around the inflexion point, still considering that $S_2$ evolves at concentrations such as the system is on the right of the maximum of ${\mu }_2$), a small change on D (and thus on $D_2$) will affect the solution of Equation (10) much more than before. In any of these situations, the relative positions of the largest solution of Equation (10) and of $S_2^{in*}$ determine whether the system will evolve in the region ${\mathcal{A}}_7$ or ${\mathcal{A}}_6$, and as D increases to a value, such as $\alpha D+k_2=d_4$, a return of the system from ${\mathcal{A}}_7$ into ${\mathcal{A}}_6$ can be observed. It goes through $P_3$ and then goes back into ${\mathcal{A}}_7$ through $P_4$. It is well illustrated in the bifurcation diagram plotted in Figure 10.

5. Conclusions

In this paper, a model of the AD has been studied. A two-step mass-balance model, corresponding to hydrolysis and methanogenesis phases, is considered. An unusual growth function for the hydrolysis step that is a generic density-dependent growth rate has been used while a Haldane-type function is considered for the methanogenesis step. To the best of the authors’ knowledge, it is the first time such a model with the association of Contois–Haldane-type growth functions in a two-step model of the AD has been considered. In this analysis, it has been shown that this model has six steady states $(E_1^0,E_1^1,E_1^2,E_2^0,E_2^1,E_2^2)$. Conditions under which they exist and are stable or unstable have been highlighted. The regions of stability of these equilibria have been established with the help of the operating diagrams and their practical use in a number of situations has been discussed. Finally, it should be underlined that the proposed change in the growth rate function of the first reaction step when considering the AD of solid waste—using a density function instead of a substrate-dependent function—changes the qualitative properties of the system when compared to systems treating liquid wastes.