# Mathematical study of a two-stage anaerobic model when the hydrolysis is the limiting step

Mohammed Hanaki, Jérôme Harmand, Zoubida Mghazli, Alain Rapaport,
Tewfik Sari, Pablo Ugalde

► To cite this version:

Mohammed Hanaki, Jérôme Harmand, Zoubida Mghazli, Alain Rapaport, Tewfik Sari, et al.. Mathematical study of a two-stage anaerobic model when the hydrolysis is the limiting step. Processes, 2021, 9, 10.3390/pr9112050. [https://doi.org/10.34254/9112050](https://doi.org/10.34254/9112050)

HAL Id: hal-03425496

https://hal.inrae.fr/hal-03425496

Submitted on 10 Nov 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Mathematical study of a two-stage anaerobic model when the hydrolysis is the limiting step

Mohammed Hanaki<sup>a</sup>, Jérôme Harmand<sup>b</sup>, Zoubida Mghazli<sup>a</sup>, Alain Rapaport<sup>c</sup>, Tewfik Sari<sup>d</sup>, Pablo Ugalde<sup>b</sup>

<sup>a</sup>LIRNE, Ibn Tofail University, Kenitra, Morocco;

<sup>b</sup>LBE, University of Montpellier, INRAE, Narbonne, France;

<sup>c</sup>Mistea, University of Montpellier, INRAE, Institut Agro, Montpellier, France;

<sup>d</sup>ITAP, University of Montpellier, INRAE, Institut Agro, Montpellier, France

## ARTICLE HISTORY

Compiled October 28, 2021

## 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 in solid matter: existence and stability properties of the equilibrium points are investigated. The hydrolysis step is considered as 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), 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.

## KEYWORDS

Commensal system; Anaerobic digestion; Steady state; Mortality; Stability; Operating diagrams.

## 1. Introduction

Two-step models are very common in the environmental engineering literature to describe engineered biological systems. They are of particular interest to design feedback control laws since they are usually enough complex 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 \xrightarrow{\mu_1(.)} X_1 + S_2, \quad S_2 \xrightarrow{\mu_2(.)} X_2 + CO_2 + CH_4$$

CONTACT Jérôme Harmand. Email: jerome.harmand@inrae.fr
This model is given by the following general dynamical system:

$$
\left\{
\begin{aligned}
\dot{S}_1 &= D(S_1^{in} - S_1) - \mu_1(\cdot) \frac{X_1}{Y_1}, \\
\dot{X}_1 &= [\mu_1(\cdot) - \alpha D - k_1] X_1, \\
\dot{S}_2 &= D(S_2^{in} - S_2) + \mu_1(\cdot) \frac{X_1}{Y_3} - \mu_2(\cdot) \frac{X_2}{Y_2}, \\
\dot{X}_2 &= [\mu_2(\cdot) - \alpha D - k_2] X_2
\end{aligned}
\right.
\qquad (1)
$$

where D is the dilution rate, while $S_1^{in}$ and $S_2^{in}$ are the input substrate concentrations, respectively. Parameters $Y_i$ are yield coefficients associated to the bio-reactions, $k_i$ are the mortality terms while $\alpha \in [0, 1]$ is a term allowing to decouple 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 Wade et al. (2016).

Among the numerous biological systems used in environmental engineering, the anaerobic systems are among those which 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 more and more used for the treatment of liquid and solid wastes. Because of its relative instability due to the possible accumulation of intermediate products however, notably the volatile fatty acids (VFA), the modeling of this process has been extensively studied over these last 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 formalize the available knowledge - typically to better understand bioprocesses - models are generally of high order and not tractable from a mathematical viewpoint, cf. for instance the ADM1 developed by the International Water Association or IWA Batstone et al. (2002). On the other hand, as already mentioned hereabove, 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. Bernard, Hadj-Sadock, Dochain, Genovesi, and Steyer (2001). 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 Sbarciog, Loccufier, and Nolus (2010), this model was studied for $\alpha = 1$ while the most interesting case where $0 < \alpha < 1$ and where growth functions were characterized by qualitative properties was studied in Benyahia, Sari, Cherki, and Harmand (2012).

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 wastes notably as a basis for control design, cf. Bernard et al. (2001), Lopez and Borzacconi (2009) or still Arzate et al. (2017) 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. Garcia-Dieguez, Bernard, and Roca (2013) or Weedermann, Seo, and Wolkowicz (2013)). 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 like the ADM1, cf. Hassam, Ficara, Leva, and Harmand (2015).

Most of the wastes 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 methano-genesis: in such a case, simple models including only acidogenesis and methanogenesis can be used as in the AM2. If the waste contains a high proportion of solid matter however, 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 Vavilin, Fernandez, Palatsi, and Flotats (2008), when hydrolysis is the limiting step, rates depending only on substrate con-centrations, 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 Contois models falls, cf. Ramirez et al. (2009). 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. Ramirez et al. (2009):

$$ \left\{ \begin{array}{l} \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, \\ \\ \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 \end{array} \right. \qquad (2) $$

While the analysis of the general model of AD initially purposed in Bernard et al. (2001) (representing acidogenesis and methanogenesis steps) has been realized in Benyahia et al. (2012), from the best of authors knowledge, a two-step model where the kinetic of the first step is modeled by a generic density-dependent kinetics while the second step exhibits a Haldane-type function has never been studied in the literature. It is the aim of the actual paper to study such a generic model. This analysis in realized in taking advantage of the fact that the system has a cascade structure: known results are then applied to study the whole 4th order system as the coupling of two 2nd order chemostat models. The main contribution of the paper is the set of operating diagrams of the 4th order system that are provided in section 4.

The paper is organized as follows. In section 2 the two-step model with two in-put 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, firstly, with respect to the first input substrate concentration and the dilu-tion rate and secondly, with respect to the second input substrate concentration and the dilution rate.
## 2. Mathematical model

The two-step model reads:

$$ 
\left\{
\begin{aligned}
\dot{S}_1 &= D(S_1^{in} - S_1) - \mu_1(S_1, X_1) \frac{X_1}{Y_1}, \\
\dot{X}_1 &= [\mu_1(S_1, X_1) - D_1]X_1, \\
\dot{S}_2 &= D(S_2^{in} - S_2) + \mu_1(S_1, X_1) \frac{X_1}{Y_3} - \mu_2(S_2) \frac{X_2}{Y_2}, \\
\dot{X}_2 &= [\mu_2(S_2) - D_2]X_2
\end{aligned}
\right.
\qquad (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$ are the yield coefficients. $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 : (S_2) \to \mu_2(S_2)$ 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 units of variables:

$$ s_1 = S_1, \quad x_1 = \frac{1}{Y_1}X_1, \quad s_2 = \frac{Y_3}{Y_1}S_2, \quad x_2 = \frac{Y_3}{Y_1Y_2}X_2 $$

The following system is obtained:

$$ 
\left\{
\begin{aligned}
\dot{s}_1 &= D(s_1^{in} - s_1) - f_1(s_1, x_1)x_1, \\
\dot{x}_1 &= [f_1(s_1, x_1) - D_1]x_1, \\
\dot{s}_2 &= D(s_2^{in} - s_2) + f_1(s_1, x_1)x_1 - f_2(s_2)x_2, \\
\dot{x}_2 &= [f_2(s_2) - D_2]x_2
\end{aligned}
\right.
\qquad (4) $$

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

$$ f_1(s_1, x_1) = \mu_1(s_1, Y_1x_1) \quad \text{and} \quad f_2(s_2) = \mu_2\left(\frac{Y_1}{Y_3}s_2\right) $$

It is assumed that the functions $\mu_1(\cdot, \cdot)$ and $\mu_2(\cdot)$ satisfy the following hypotheses.

**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(x_1)$. 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$.

**H2.** $\mu_2(s_2)$ is positive for $s_2 > 0$, and satisfies $\mu_2(0) = 0$ and $\mu_2(+\infty) = 0$. Moreover $\mu_2(s_2)$ increases until a concentration $s_2^M$ and then decreases; thus $\mu_2'(s_2) > 0$ for $0 \le s_2 < s_2^M$, and $\mu_2'(s_2) < 0$ for $s_2 > s_2^M$.

As underlined in the introduction, particular kinetics models as 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.

**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(x_1)$. 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$.

**A2.** $f_2(s_2)$ is positive for $s_2 > 0$, and satisfies $f_2(0) = 0$ and $f_2(+\infty) = 0$. Moreover $f_2(s_2)$ increases until a concentration $s_2^M$ and then decreases, with $f_2'(s_2) > 0$ for $0 < s_2 < s_2^M$, and $f_2'(s_2) < 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 which 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:

$$ \begin{cases} \dot{s}_1 = D(s_1^{in} - s_1) - f_1(s_1, x_1)x_1, \\ \dot{x}_1 = [f_1(s_1, x_1) - D_1]x_1. \end{cases} \tag{5} $$

The behaviour of this system is well-known, cf. Harmand, J. and C. Lobry and A. Rapaport and T. Sari (2016). A steady state $(s_1^*, x_1^*)$ must be solution of the system

$$ \begin{cases} 0 = D(s_1^{in} - s_1) - f_1(s_1, x_1)x_1, \\ 0 = [f_1(s_1, x_1) - D_1]x_1 \end{cases} \tag{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 \quad \text{and} \quad x_1^* = \frac{D}{D_1}(s_1^{in} - s_1^*). \tag{7} $$

Let $\gamma$ a function defined by :

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

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

According to the hypothesis **A1**, $\gamma(s_1)$ is strictly increasing over the interval $]0, s_1^{in}[$, with $\gamma(0) = 0$ and $\gamma(s_1^{in}) = f_1(s_1^{in}, 0)$. According to the theorem of intermediate values, the equation $\gamma(s_1) = D_1$ has a solution between 0 and $s_1^{in}$ if and only if $D_1 < \gamma(s_1^{in})$, that is to say 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
<figure><img src="image_2.png" /><figcaption>Figure 1. The existence of the solution of $\gamma(s_1) = D_1$.</figcaption></figure>

abbreviations LES for locally exponentially stable is used.

**Proposition 3.1.** Assume that Assumptions A1 and A2 hold. Then, the local stability of steady states of (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 Harmand, J. and C. Lobry and A. Rapaport and T. Sari (2016) 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 3.1 are summarized in the following table.

<table><caption>**Table 1.** Summary of the results of Proposition 3.1</caption><thead><tr><th>Steady state</th><th>Existence condition</th><th>Stability condition</th></tr></thead><tbody><tr><td>$E_0$</td><td>Always exists</td><td>$f_1(s_1^{in}, 0) < D_1$</td></tr><tr><td>$E_1$</td><td>$f_1(s_1^{in}, 0) > D_1$</td><td>Stable when it exists</td></tr></tbody></table>

### 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$ - that 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 the system (5) is shown in Figure 2. The condition $f_1(s_1^{in}, 0) > D_1$ of existence of $E_1$ is equivalent to $D < \frac{1}{\alpha}[f_1(s_1^{in}, 0) - k_1]$. Therefore, the curve

$$ \Gamma : \{ (s_1^{in}, D) : D = \frac{1}{\alpha}[f_1(s_1^{in}, 0) - k_1] \} $$
separates the operating plan in two regions as defined in Figure 2.

<figure><img src="image_3.png" /><figcaption>Figure 2. Operating diagram of (5)</figcaption></figure>

The curve $\Gamma$ is the border which makes $E_0$ unstable and at the same time $E_1$ exists. Table 2 indicates the stability properties of steady states of the 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.

<table><caption>**Table 2.** Stability properties of steady states of system (5) in each region</caption><thead><tr><th>Region</th><th>Equ. E<sub>0</sub></th><th>Equ. E<sub>1</sub></th></tr></thead><tbody><tr><td>(s<sub>1</sub><sup>in</sup>, D) ∈ R<sub>0</sub></td><td>S</td><td></td></tr><tr><td>(s<sub>1</sub><sup>in</sup>, D) ∈ R<sub>1</sub></td><td>U</td><td>S</td></tr></tbody></table>

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. Bernard et al. (2001).

## 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 (4), the dynamics of the state variables $s_2$ and $x_2$ are given by

$$ \begin{cases} \dot{s}_2 = D(F(t) - s_2) - f_2(s_2)x_2, \\ \dot{x}_2 = [f_2(s_2) - D_2]x_2, \end{cases} \qquad (8) $$

with

$$ F(t) = s_2^{in} + \frac{1}{D} f_1(s_1(t), x_1(t)) x_1(t) $$

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

$$ \left\{ \begin{array}{l} \dot{s}_2 = D(s_2^{in*} - s_2) - f_2(s_2)x_2, \\ \dot{x}_2 = [f_2(s_2) - D_2]x_2 \end{array} \right. \qquad (9) $$

where

$$ s_2^{in*} = s_2^{in} \quad \text{or} \quad s_2^{in*} = s_2^{in} + \frac{D_1}{D} x_1^*. \qquad (10) $$

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

The 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 (10), depends explicitly on the input flow rate. For a given $D$, the long term behaviour of such a system is well-known, cf. Harmand, J. and C. Lobry and A. Rapaport and T. Sari (2016).

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

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

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

$$ f_2(s_2) = D_2. \qquad (12) $$

Under hypothesis **A2**, and if

$$ D_2 < f_2(s_2^M) \qquad (13) $$

this equation has two solutions which 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), \quad i = 1, 2. \qquad (14) $$

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

**Proposition 3.2.** Assume that Assumptions **A1**, **A2** and condition (13) hold. Then, the local stability of steady states of (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 Harmand, J. and C. Lobry and A. Rapaport and T. Sari (2016) for the proof of this proposition.

The results of Proposition 3.2 are summarized in the following table.
<table><caption>**Table 3.** Summary of the results of Proposition 3.2</caption><thead><tr><th>Steady-state</th><th>Existence condition</th><th>Stability condition</th></tr></thead><tbody><tr><td>F<sub>0</sub></td><td>Always exists</td><td>s<sub>2</sub><sup>in*</sup> &lt; s<sub>2</sub><sup>1</sup> or s<sub>2</sub><sup>in*</sup> &gt; s<sub>2</sub><sup>2</sup></td></tr><tr><td>F<sub>1</sub></td><td>s<sub>2</sub><sup>in*</sup> &gt; s<sub>2</sub><sup>1</sup></td><td>Stable if it exists</td></tr><tr><td>F<sub>2</sub></td><td>s<sub>2</sub><sup>in*</sup> &gt; s<sub>2</sub><sup>2</sup></td><td>Unstable if it exists</td></tr></tbody></table>

### 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(s_2^{in*}) = D_2$, that is to say $D = \frac{1}{\alpha}(f_2(s_2^{in*}) - k_2)$. Therefore, the horizontal line

$$ \Gamma_1 : \left\{ (s_2^{in*}, D) : D = \frac{1}{\alpha}(f_2(s_2^M) - 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(s_2^{in*}) - k_2) \right\} $$

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

<figure><img src="image_4.png" /><figcaption>**Figure 3.** Operating diagram of (9)</figcaption></figure>

The following table indicates the stability properties of steady states of system (8).

<table><caption>**Table 4.** Stability properties of the steady states of system (8)</caption><thead><tr><th>Region</th><th>Equ. F<sub>0</sub></th><th>Equ. F<sub>1</sub></th><th>Equ. F<sub>2</sub></th></tr></thead><tbody><tr><td>(s<sub>2</sub><sup>in*</sup>, D) ∈ R<sub>2</sub></td><td>S</td><td></td><td></td></tr><tr><td>(s<sub>2</sub><sup>in*</sup>, D) ∈ R<sub>3</sub></td><td>U</td><td>S</td><td></td></tr><tr><td>(s<sub>2</sub><sup>in*</sup>, D) ∈ R<sub>4</sub></td><td>S</td><td>S</td><td>U</td></tr></tbody></table>

### 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 (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 (4), then $(s_1^*, x_1^*)$ is a steady state of (5) and $(s_2^*, x_2^*)$ is a steady state of (9) where $s_2^{in*}$ is given by (10).

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

$$ (1) \quad (s_2^*, x_2^*) = (s_2^{in}, 0), \text{ and } E_1^0 := (s_1^{in}, 0, s_2^{in}, 0). $$

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

$$ (3) \quad (s_2^*, x_2^*) = (s_2^2, x_2^2), \text{ 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

$$ (1) \quad (s_2^*, x_2^*) = (s_2^{in*}, 0), \text{ and } E_2^0 := (s_1^*, x_1^*, s_2^{in*}, 0). $$

$$ (2) \quad (s_2^*, x_2^*) = (s_1^1, x_2^{1*}), \text{ and } E_2^1 := (s_1^*, x_1^*, s_2^1, x_2^{1*}). $$

$$ (3) \quad (s_2^*, x_2^*) = (s_2^2, x_2^{2*}), \text{ and } E_2^2 := (s_1^*, x_1^*, s_2^2, x_2^{2*}). $$

These results are summarized in the following proposition.

**Proposition 3.3.** The 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_1^1$.

* $E_1^2 = (s_1^{in}, 0, s_2^2, x_2^2)$, exists if and only if $s_2^{in} > s_1^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_1^2$.

* $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_1^2$.

### 3.3.2. Steady states stability of the system (4)

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

$$ J = \begin{bmatrix} J_{11} & J_{12} \\ 0 & J_{22} \end{bmatrix}, $$

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

$$ J = \begin{bmatrix} -D - Mx_1 & -Nx_1 - f_1(s_1, x_1) \\ Mx_1 & Nx_1 + [f_1(s_1, x_1) - D_1] \end{bmatrix}, \quad (15) $$

and

$$ J = \begin{bmatrix} -D - f_2'(s_2)x_2 & -f_2(s_2) \\ f_2'(s_2)x_2 & f_2(s_2) - D_2 \end{bmatrix}. \quad (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 (7) with respect to the four numbers $s_2^1$ and $s_2^2$, defined by (12) on the one hand, and $s_2^{in}$ and $s_2^{in*}$, solutions of (10) on the other hand. Equilibrium stability is summarized in Table 5 while the different regions of the operating diagram are synthesized in Tables 6 and 7.
<table><caption>**Table 5.** The stability conditions for the steady states of the system (4)</caption><thead><tr><th>Equ.</th><th>Matrices J<sub>11</sub> and J<sub>22</sub></th><th>Conditions of stability</th></tr></thead><tbody><tr><td>E<sub>1</sub><sup>0</sup></td><td>J<sub>11</sub> = [-D -f<sub>1</sub>(s<sub>1</sub><sup>in</sup>, 0) 0 f<sub>1</sub>(s<sub>1</sub><sup>in</sup>, 0) -D<sub>1</sub>]</td><td>Tr(J<sub>11</sub>) &lt; 0 if f<sub>1</sub>(s<sub>1</sub><sup>in</sup>, 0) &lt; D<sub>1</sub>,<br>Tr(J<sub>22</sub>) &lt; 0 if s<sub>2</sub><sup>in</sup> &lt; s<sub>2</sub><sup>1</sup> or s<sub>2</sub><sup>in</sup> &gt; s<sub>2</sub><sup>2</sup>,<br>det(J<sub>11</sub>) &gt; 0 and det(J<sub>22</sub>) &gt; 0</td></tr><tr><td>J<sub>22</sub></td><td>J<sub>22</sub> = [-D 0 -f<sub>2</sub>(s<sub>2</sub><sup>in</sup>) f<sub>2</sub>(s<sub>2</sub><sup>in</sup>) -D<sub>2</sub>]</td><td>⇒ {<br>E<sub>1</sub><sup>0</sup> is stable if f<sub>1</sub>(s<sub>1</sub><sup>in</sup>, 0) ≤ D<sub>1</sub> and<br>(s<sub>2</sub><sup>in</sup> &lt; s<sub>2</sub><sup>1</sup> or s<sub>2</sub><sup>in</sup> &gt; s<sub>2</sub><sup>2</sup>),<br>E<sub>1</sub><sup>0</sup> is unstable if f<sub>1</sub>(s<sub>1</sub><sup>in</sup>) &gt; D<sub>1</sub> or<br>s<sub>2</sub><sup>1</sup> &lt; s<sub>2</sub><sup>in</sup> &lt; s<sub>2</sub><sup>2</sup>.)</td></tr><tr><td>E<sub>1</sub><sup>i</sup></td><td>(i = 1, 2)<br>J<sub>11</sub> = [-D 0 -f<sub>1</sub>(s<sub>1</sub><sup>in</sup>, 0) f<sub>1</sub>(s<sub>1</sub><sup>in</sup>, 0) -D<sub>1</sub>]</td><td>Tr(J<sub>11</sub>) &lt; 0 and det(J<sub>11</sub>) &gt; 0 if f<sub>1</sub>(s<sub>1</sub><sup>in</sup>, 0) &lt; D<sub>1</sub><br>Tr(J<sub>22</sub>) &lt; 0 and det(J<sub>22</sub>) &gt; 0 at E<sub>1</sub><sup>1</sup>,<br>det(J<sub>22</sub>) &lt; 0 at E<sub>1</sub><sup>2</sup></td></tr><tr><td>J<sub>22</sub></td><td>J<sub>22</sub> = [-[D + f'<sub>2</sub>(s<sub>2</sub><sup>i</sup>)x<sub>2</sub><sup>i</sup>] -D<sub>2</sub> f'<sub>2</sub>(s<sub>2</sub><sup>i</sup>)x<sub>2</sub><sup>i</sup> 0]</td><td>⇒ {<br>E<sub>1</sub><sup>1</sup> is stable,<br>E<sub>1</sub><sup>2</sup> is unstable<br>E<sub>1</sub><sup>i</sup> are both unstable if s<sub>1</sub><sup>in</sup> &gt; s<sub>1</sub><sup>*</sup>.</td></tr><tr><td>E<sub>2</sub><sup>0</sup></td><td>J<sub>11</sub> = [-[D + (∂f<sub>1</sub>/∂s<sub>1</sub>)]x<sub>1</sub><sup>*</sup> -[D<sub>1</sub> + (∂f<sub>1</sub>/∂x<sub>1</sub>)]x<sub>1</sub><sup>*</sup> (∂f<sub>1</sub>/∂s<sub>1</sub>)x<sub>1</sub><sup>*</sup> (∂f<sub>1</sub>/∂x<sub>1</sub>)x<sub>1</sub><sup>*</sup></td><td>Tr(J<sub>11</sub>) &lt; 0 and det(J<sub>11</sub>) &gt; 0 by A1<br>Tr(J<sub>22</sub>) &lt; 0 and det(J<sub>22</sub>) &gt; 0<br>if s<sub>2</sub><sup>in*</sup> &lt; s<sub>2</sub><sup>1</sup> or s<sub>2</sub><sup>in*</sup> &gt; s<sub>2</sub><sup>2</sup></td></tr><tr><td>J<sub>22</sub></td><td>J<sub>22</sub> = [-D 0 -f<sub>2</sub>(s<sub>2</sub><sup>in*</sup>) [f<sub>2</sub>(s<sub>2</sub><sup>in*</sup>) - D<sub>2</sub>]</td><td>⇒ {<br>E<sub>2</sub><sup>0</sup> is stable if s<sub>1</sub><sup>in</sup> ≥ s<sub>1</sub><sup>*</sup> and<br>(s<sub>2</sub><sup>in*</sup> &lt; s<sub>2</sub><sup>1</sup> or s<sub>2</sub><sup>in*</sup> &gt; s<sub>2</sub><sup>2</sup>),<br>E<sub>2</sub><sup>0</sup> is unstable if s<sub>1</sub><sup>in</sup> &gt; s<sub>1</sub><sup>*</sup> and<br>s<sub>2</sub><sup>1</sup> &lt; s<sub>2</sub><sup>in*</sup> &lt; s<sub>2</sub><sup>2</sup>.</td></tr><tr><td>E<sub>2</sub><sup>i</sup></td><td>(i = 1, 2)<br>J<sub>11</sub> = [-[D + (∂f<sub>1</sub>/∂s<sub>1</sub>)]x<sub>1</sub><sup>*</sup> -[D<sub>1</sub> + (∂f<sub>1</sub>/∂x<sub>1</sub>)]x<sub>1</sub><sup>*</sup> (∂f<sub>1</sub>/∂s<sub>1</sub>)x<sub>1</sub><sup>*</sup> (∂f<sub>1</sub>/∂x<sub>1</sub>)x<sub>1</sub><sup>*</sup></td><td>Tr(J<sub>11</sub>) &lt; 0 and det(J<sub>11</sub>) &gt; 0<br>Tr(J<sub>22</sub>) &lt; 0 and det(J<sub>22</sub>) &gt; 0 at E<sub>2</sub><sup>1</sup>,<br>det(J<sub>22</sub>) &lt; 0 at E<sub>2</sub><sup>2</sup></td></tr><tr><td>J<sub>22</sub></td><td>J<sub>22</sub> = [-[D + f'<sub>2</sub>(s<sub>2</sub><sup>i</sup>)x<sub>2</sub><sup>i*</sup>] -D<sub>2</sub> f'<sub>2</sub>(s<sub>2</sub><sup>i</sup>)x<sub>2</sub><sup>i*</sup> 0]</td><td>⇒ {<br>E<sub>2</sub><sup>1</sup> is stable,<br>E<sub>2</sub><sup>2</sup> is unstable</td></tr></tbody></table>

<table><caption>**Table 6.** The three cases when s<sub>1</sub><sup>in</sup> &lt; s<sub>1</sub><sup>*</sup></caption><thead><tr><th>Case</th><th>Area</th><th>Condition</th><th>E<sub>1</sub><sup>0</sup></th><th>E<sub>1</sub><sup>1</sup></th><th>E<sub>1</sub><sup>2</sup></th><th>E<sub>2</sub><sup>0</sup></th><th>E<sub>2</sub><sup>1</sup></th><th>E<sub>2</sub><sup>2</sup></th></tr></thead><tbody><tr><td>1.1</td><td>A<sub>1</sub></td><td>s<sub>2</sub><sup>in</sup> &lt; s<sub>2</sub><sup>1</sup> &lt; s<sub>2</sub><sup>2</sup></td><td>S</td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>1.2</td><td>A<sub>2</sub></td><td>s<sub>2</sub><sup>1</sup> &lt; s<sub>2</sub><sup>in</sup> &lt; s<sub>2</sub><sup>2</sup></td><td>U</td><td>S</td><td></td><td></td><td></td><td></td></tr><tr><td>1.3</td><td>A<sub>3</sub></td><td>s<sub>2</sub><sup>1</sup> &lt; s<sub>2</sub><sup>2</sup> &lt; s<sub>2</sub><sup>in</sup></td><td>S</td><td>S</td><td>U</td><td></td><td></td><td></td></tr></tbody></table>

<table><caption>**Table 7.** The six cases when s<sub>1</sub><sup>in</sup> &gt; s<sub>1</sub><sup>*</sup></caption><thead><tr><th>Case</th><th>Area</th><th>Condition</th><th>E<sub>1</sub><sup>0</sup></th><th>E<sub>1</sub><sup>1</sup></th><th>E<sub>1</sub><sup>2</sup></th><th>E<sub>2</sub><sup>0</sup></th><th>E<sub>2</sub><sup>1</sup></th><th>E<sub>2</sub><sup>2</sup></th></tr></thead><tbody><tr><td>2.1</td><td>A<sub>4</sub></td><td>s<sub>2</sub><sup>in</sup> &lt; s<sub>2</sub><sup>in*</sup> &lt; s<sub>2</sub><sup>1</sup> &lt; s<sub>2</sub><sup>2</sup></td><td>U</td><td></td><td></td><td>S</td><td></td><td></td></tr><tr><td>2.2</td><td>A<sub>5</sub></td><td>s<sub>2</sub><sup>in</sup> &lt; s<sub>2</sub><sup>1</sup> &lt; s<sub>2</sub><sup>in*</sup> &lt; s<sub>2</sub><sup>2</sup></td><td>U</td><td></td><td></td><td>U</td><td>S</td><td></td></tr><tr><td>2.3</td><td>A<sub>6</sub></td><td>s<sub>2</sub><sup>in</sup> &lt; s<sub>2</sub><sup>1</sup> &lt; s<sub>2</sub><sup>in*</sup> &lt; s<sub>2</sub><sup>2</sup></td><td>U</td><td></td><td></td><td>S</td><td>S</td><td>U</td></tr><tr><td>2.4</td><td>A<sub>7</sub></td><td>s<sub>2</sub><sup>1</sup> &lt; s<sub>2</sub><sup>in</sup> &lt; s<sub>2</sub><sup>in*</sup> &lt; s<sub>2</sub><sup>2</sup></td><td>U</td><td>U</td><td></td><td>U</td><td>S</td><td></td></tr><tr><td>2.5</td><td>A<sub>8</sub></td><td>s<sub>2</sub><sup>1</sup> &lt; s<sub>2</sub><sup>in</sup> &lt; s<sub>2</sub><sup>2</sup> &lt; s<sub>2</sub><sup>in*</sup></td><td>U</td><td>U</td><td></td><td>S</td><td>S</td><td>U</td></tr><tr><td>2.6</td><td>A<sub>9</sub></td><td>s<sub>2</sub><sup>1</sup> &lt; s<sub>2</sub><sup>2</sup> &lt; s<sub>2</sub><sup>in</sup> &lt; s<sub>2</sub><sup>in*</sup></td><td>U</td><td>U</td><td>U</td><td>S</td><td>S</td><td>U</td></tr></tbody></table>

**Remark 1.** Here the limit values in the stability condition (Ex: s<sub>2</sub><sup>in</sup> = s<sub>2</sub><sup>1</sup> or s<sub>2</sub><sup>in</sup> = s<sub>2</sub><sup>2</sup>)
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 3.4.

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

<table><thead><tr><th>Case</th><th>Condition</th><th>NH</th><th>S</th><th>U</th></tr></thead><tbody><tr><td>1.4</td><td>$s_2^{in} = s_2^1 < s_2^2$</td><td>$E_1^0 = E_1^1$</td><td></td><td></td></tr><tr><td>1.5</td><td>$s_2^1 < s_2^{in} = s_2^2$</td><td>$E_1^0 = E_1^2$</td><td>$E_1^1$</td><td></td></tr><tr><td>1.6</td><td>$s_2^1 = s_2^2 < s_2^{in}$</td><td>$E_1^1 = E_1^2$</td><td>$E_1^0$</td><td></td></tr></tbody></table>

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

<table><thead><tr><th>Case</th><th>Condition</th><th>NH</th><th>S</th><th>U</th></tr></thead><tbody><tr><td>2.7</td><td>$s_2^{in} < s_2^{in*} = s_2^1 < s_2^2$</td><td>$E_2^0 = E_2^1$</td><td></td><td>$E_1^0$</td></tr><tr><td>2.8</td><td>$s_2^{in} < s_2^1 < s_2^2 = s_2^{in*}$</td><td>$E_2^0 = E_2^2$</td><td>$E_2^1$</td><td>$E_1^0$</td></tr><tr><td>2.9</td><td>$s_2^{in} = s_2^1 < s_2^{in*} < s_2^2$</td><td>$E_1^0 = E_1^1$</td><td>$E_2^1$</td><td>$E_2^0$</td></tr><tr><td>2.10</td><td>$s_2^{in} = s_2^1 < s_2^2 < s_2^{in*}$</td><td>$E_1^0 = E_1^1$</td><td>$E_2^1, E_2^0$</td><td>$E_2^2$</td></tr><tr><td>2.11</td><td>$s_2^{in} = s_2^1 < s_2^2 = s_2^{in*}$</td><td>$E_1^0 = E_1^1, E_2^0 = E_2^1$</td><td>$E_2^1$</td><td></td></tr><tr><td>2.12</td><td>$s_2^1 < s_2^{in} = s_2^2 < s_2^{in*}$</td><td>$E_2^1 = E_2^2$</td><td>$E_2^0$</td><td>$E_1^0$</td></tr><tr><td>2.13</td><td>$s_2^1 < s_2^{in} < s_2^2 = s_2^{in*}$</td><td>$E_2^0 = E_2^2$</td><td>$E_2^1$</td><td>$E_1^0, E_1^1$</td></tr><tr><td>2.14</td><td>$s_2^1 < s_2^{in} = s_2^2 < s_2^{in*}$</td><td>$E_1^0 = E_1^2$</td><td>$E_2^0, E_2^1$</td><td>$E_1^0, E_2^2$</td></tr><tr><td>2.15</td><td>$s_2^1 = s_2^2 < s_2^{in} < s_2^{in*}$</td><td>$E_1^1 = E_1^2, E_2^1 = E_2^2$</td><td>$E_2^0$</td><td>$E_1^0$</td></tr></tbody></table>

**Proof.** Let us give the details of the proof in the 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 (14). Therefore (cf. Proposition 3.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 Tables 6, 7 and in the previous Theorem 3.4, a corresponding figure which represents the relative position of $s_2^1$, $s_2^2$, $s_2^{in}$ and $s_2^{in*}$ Benyahia et al. (2012) can be associated.

## 4. Simulations

To illustrate these results, the operating diagrams of the system (3) under hypothesis **H1** and **H2** in a number of situations are plotted. Recall that operating diagram summarizes the existence and the nature of 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 $\{S_1^{in}, D\}$ for a fixed value of $S_2^{in}$ or in the plan $\{S_2^{in}, D\}$ 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}, \qquad \mu_2(S_2) = \frac{m_2 S_2}{\frac{S_2^2}{T} + 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 Bernard et al. (2001), while others are changed significantly as the inhibition coefficient *I* of the Haldane function. Indeed, as underlined in Khedim, Benyahia, Cherki, Sari, and Harmand (2018), inhibition of the second reaction is not visible if original parameters proposed in Bernard et al. (2001) are used considering reasonable ranges of variations for *S*₁ and *S*₂. With respect to this later, the Haldane parameter *I* was thus significantly decreased to willingly increase the inhibition effect of *S*₂ on the growth of *X*₂. Finally, parameter values used are summarized in the following table To compute the different regions

<table><caption>**Table 8.** Parameters values for the simulations</caption><thead><tr><th>Parameter</th><th>Unit</th><th>Nominal value</th></tr></thead><tbody><tr><td>m<sub>1</sub></td><td>d<sup>-1</sup></td><td>0.5</td></tr><tr><td>K<sub>1</sub></td><td>g/L</td><td>2.1</td></tr><tr><td>m<sub>2</sub></td><td>d<sup>-1</sup></td><td>1</td></tr><tr><td>I</td><td>mmol/L</td><td>60</td></tr><tr><td>K<sub>2</sub></td><td>mmol/L</td><td>24</td></tr><tr><td>k<sub>1</sub></td><td>d<sup>-1</sup></td><td>0.1</td></tr><tr><td>k<sub>2</sub></td><td>d<sup>-1</sup></td><td>0.06</td></tr><tr><td>α</td><td>in [0, 1)</td><td>0.5</td></tr><tr><td>Y<sub>1</sub></td><td>g/g</td><td>1/25</td></tr><tr><td>Y<sub>2</sub></td><td>g/mmol</td><td>1/250</td></tr><tr><td>Y<sub>3</sub></td><td>g/mmol</td><td>1/268</td></tr></tbody></table>

of the operating diagrams, the numerical method reported in Khedim et al. (2018) is used. The algorithm is recalled hereafter.

### 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 sections 3.1 and 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:

### 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 {$S_1^{in}, D$} for $S_2^{in} = 1.5$ mmol/L. The regions are defined as follows. $A_1$ (in green) is the stability region of the washout $E_1^0$, $A_5$ (in blue) is the stability region of steady-state $E_2^1$, $A_6$ (in red) is the bi-stability region of the steady-states $E_2^0$ and $E_2^1$, and $A_7$ (in dark blue) is the stability region of steady-state $E_2^1$, the difference between the regions $A_5$ and $A_7$ being that the equilibrium $E_1^1$ does not exist in the region $A_5$ but exists and is unstable in $A_7$, cf. Tables 6 and 7.
<table><thead><tr><th colspan="2">Algorithm 1 Operating diagram</th></tr></thead><tbody><tr><td>for</td><td>i varying from 1 to N<sub>1</sub> do;</td></tr><tr><td>for</td><td>j varying from 1 to N<sub>2</sub>;</td></tr><tr><td></td><td>determine 6 equilibria of the model E<sub>1</sub>...E<sub>6</sub></td></tr><tr><td>for</td><td>k varying from 1 to 6 do</td></tr><tr><td></td><td>calculate the Jacobian matrix at E<sub>k</sub> (J<sub>E<sub>k</sub></sub>)</td></tr><tr><td></td><td>calculate the eigenvalues of (J<sub>E<sub>k</sub></sub>)</td></tr><tr><td></td><td>if all the conditions of existence of E<sub>k</sub> are fulfilled and all real parts of the eigenvalues of (J<sub>E<sub>k</sub></sub>) are non-positive then E<sub>k</sub> is stable</td></tr><tr><td></td><td>else if all conditions of existence of E<sub>k</sub> are fulfilled and at least one real part of eigenvalue of (J<sub>E<sub>k</sub></sub>) is positive then E<sub>k</sub> is unstable</td></tr><tr><td></td><td>else E<sub>k</sub> does not exist</td></tr><tr><td></td><td>end if</td></tr><tr><td></td><td>end for (k)</td></tr><tr><td></td><td>end for (j)</td></tr><tr><td></td><td>end for (i)</td></tr></tbody></table>

<figure><img src="image_5.png" /><figcaption>**Figure 4.** The operating diagram for $S_2^{in} = 1.5 \text{ mmol/L}$.</figcaption></figure>

Figure 5 represents the operating diagram of model (1) in the plan {$S_2^{in}, D$} for $S_1^{in} = 0.8 \text{ 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. Tables 6 and 7.

Figure 6 represents the operating diagram of model (1) in the plan {$S_2^{in}, D$} for $S_1^{in} = 0.03 \text{ g/L}$, that 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) which is the stability region of steady-state $E_2^0$ and ii) a sharp decrease of the size of region $\mathcal{A}_8$ which 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
<figure><img src="image_6.png" /><figcaption>Figure 5. The operating diagram for $S_1^{\text{in}} = 0.8 \text{ g/L}$.</figcaption></figure>

is numerically linked to the value of $S_1^{\text{in}}$. Indeed, when reducing the value of $S_1^{\text{in}}$, the region $\mathcal{A}_4$ appears while the size of the region $\mathcal{A}_8$ is getting smaller and smaller (compare Figures 5 and 6). In other words, in decreasing $S_1^{\text{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^{\text{in}}$, say $S_1^{\text{in}1}$ and $S_1^{\text{in}2}$ where $S_1^{\text{in}2} > S_1^{\text{in}1}$, a greater dilution rate is needed to destabilize the process if $S_1^{\text{in}} = S_1^{\text{in}2}$ than if $S_1^{\text{in}} = S_1^{\text{in}1}$, thus the shard decrease of $\mathcal{A}_8$ observed in Figure 6.

<figure><img src="image_7.png" /><figcaption>Figure 6. The operating diagram for $S_1^{\text{in}} = 0.03 \text{ g/L}$.</figcaption></figure>

### 4.3. *Practical interpretations of the operating diagrams*

Here 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 points is $E_1^0$ or $E_2^0$, that 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 \text{ 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} = 18 \text{ g/L}$. In such a situation, the following regions are browsed

<figure><img src="image_8.png" /><figcaption>**Figure 7.** Biological interpretation for $S_1^{in} = 18 \text{ g/L}$</figcaption></figure>

in considering successive equilibrium when increasing $D$: $A_7 \rightarrow A_6 \rightarrow A_5 \rightarrow A_1$, cf. Figure 7. To better interpret whose 'steady states the system passes through', the bifurcation diagram is plotted in Figure 8.

<figure><img src="image_9.png" /><figcaption>**Figure 8.** The bifurcation diagram for the input control $D$ for $S_1^{in} = 18 \text{ g/L}$.</figcaption></figure>

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 substrate 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 that one for which the term $S_2^{\text{in}*}$ becomes exactly the largest solution of the equation $\mu_2(S_2) = D_2$ (equivalent to equation (12) for the system (4)): the system enters then in 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^{\text{in}*}$ decreases as can be seen from equation (10). When $D_2 = d_2$, the quantity of available resource necessary to $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 definitely $X_2$ to the washout: the system enters into $\mathcal{A}_5$ in 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$ goes also to extinction and the system enters into $\mathcal{A}_1$.

**Example 2:** Let $S_1^{\text{in}} = 14 \text{ g/L}$. This case is even more interesting since, when $D$

<figure><img src="image_10.png" /><figcaption>**Figure 9.** Biological interpretation for $S_1^{\text{in}} = 14 \text{ g/L}$</figcaption></figure>

increases, the system goes back to $\mathcal{A}_5$ once before leaving it definitely in browsing the following regions: $\mathcal{A}_7 \dashrightarrow \mathcal{A}_6 \dashrightarrow \mathcal{A}_7 \dashrightarrow \mathcal{A}_5 \dashrightarrow \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 than 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^{\text{in}} = 18 \text{ g/L}$, thus enlarging the attraction basin of the stable positive equilibrium. Then, when $D$ is further increased, it may happen an interesting phenomenon: the system enters back into $\mathcal{A}_7$ through point $P_4$ instead of entering $\mathcal{A}_5$ as it was the case before. In fact, this strongly depends on model parameters and in particular on the relative rate at which $S_2^{\text{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^{\text{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 very much the largest solution of the equation (10) while $S_2^{\text{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 $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 much more the solution of the equation (10) than before. In any of these situations, the relative positions of the largest solution of the equation (10) and of $S_2^{in*}$ determine whether the system will evolve in the region $A_7$ or $A_6$ and, as $D$ increases to a value such $\alpha D + k_2 = d_4$ a return of the system from $A_7$ into $A_6$ can be observed. It goes through $P_3$ and then goes back into $A_7$ through $P_4$. It is well illustrated in the bifurcation diagram plotted in Figure 10.

<figure><img src="image_11.png" /><figcaption>Figure 10. The bifurcation diagram of the operating diagram shown in Figure 4 for $S_1^{in} = 14g/L$.</figcaption></figure>

## 5. Conclusions

In this paper, a model of the AD has been studied. A two stages mass-balance model, corresponding to hydrolysis and methanogenesis phases, is considered. A non usual 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. From the best of authors knowledge, it is the first time such a model including the association of Contois-Haldane-type growth functions in a two-step model of the AD is 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 wastes - using a density function instead of a substrate-dependent function - changes the qualitative properties of the system when compared to systems treating liquid wastes.
Acknowledgements

The authors thank the projects PHC TOUBKAL TBK 17/47-36773WE and the TREASURE Euro-Mediterranean research network (www.inra.fr/treasure) for their financial supports.

References

Arzate, J., Kirstein, M., Ertem, F., Kielhorn, E., Malule, H., Neubauer, P., ... Junne, S. (2017). Anaerobic digestion model (am2) for the description of biogas processes at dynamic feedstock loading rates. *Chem. Ing. Tech.*, **89**(5), 686–695.

Batstone, D., Keller, J., Angelidaki, I., Kalyzhnyi, S., Pavlostathis, S., Rozzi, A., ... Vavilin, V. (2002). The iwa anaerobic digestion model no 1 (adm1). *Water Sci. Technol.*, **45**, 65–73.

Benyahia, B., Sari, T., Cherki, B., & Harmand, J. (2012). Bifurcation and stability analysis of a two step model for monitoring anaerobic digestion prosesses. *J. Process Control*, **22**, 1008–1019.

Bernard, O., Hadj-Sadock, Z., Dochain, D., Genovesi, A., & Steyer, J.-P. (2001). Dynamical model development and parameter identification for an anaerobic wastewater treatment process. *J. Biotechnol. Bioeng.*, **75**, 424–438.

Garcia-Dieguez, C., Bernard, O., & Roca, E. (2013). Reducing the anaerobic digestion model no. 1 for its application to an industrial wastewater treatment plant treatingwinery effluent wastewater. *Bioresour. Technol.*, **132**, 244–253.

Harmand, J. and C. Lobry and A. Rapaport and T. Sari. (2016). *Chemostat: Theory mathematics of continuous culture of microorganisms*. ISTE-Wiley.

Hassam, S., Ficara, E., Leva, A., & Harmand, J. (2015). A generic and systematic procedure to derive a simplified model from the anaerobic digestion model no. 1 (adm1). *Biochemical Engineering Journal*, **99**, 193–203.

Khedim, Z., Benyahia, B., Cherki, B., Sari, T., & Harmand, J. (2018). Effect of control parameters on biogas production during the anaerobic digestion of protein-rich substrates. *Applied Mathematical Modelling*, **61**, 351–376.

Lopez, I., & Borzacconi, L. (2009). Modelling a full scale uasb reactor using a cod global balance approach and state observers. *Chemical Engineering Journal*, **146**, 1–5.

Ramirez, I., Mottet, A., Carrère, H., Déléris, S., Vedrenne, F., & Steyer, J.-P. (2009). Modified adm1 disintegration/hydrolysis structures for modeling batch thermophilic anaerobic digestion of thermally pretreated waste activated sludge. *Water Research*, **43**, 3479–3492.

Sbarciog, M., Loccfufer, M., & Nolus, E. (2010). Determination of appropriate operating strategies for anaerobic digestion systems. *Biochemical engineering journal*, **51**(3), 180–188.

Vavilin, V., Fernandez, B., Palatsi, J., & Flotats, X. (2008). Hydrolysis kinetics in anaerobic degradation of particulate organic material: An overview. *Waste Management*, **28**, 939–951.

Wade, M., Harmand, J., Benyahia, B., Bouchez, T., Chailou, S., Cloez, B., ... Lobry, C. (2016). Perspectives in mathematical modelling for microbial ecology. *Ecological Modelling*, **321**, 64–74.

Weedermann, M., Seo, G., & Wolkowicz, G. (2013). Mathematical model of anaerobic digestion in a chemostat: Effects of syntrophy and inhibition. *J. Biol. Dyn*, **7**, 59–85.