# Effect of a new variable integration on steady states of a two-step Anaerobic Digestion Model

Boumediene Benyahia, Tewfik Sari

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

## To cite this version:

Boumediene Benyahia, Tewfik Sari. Effect of a new variable integration on steady states of a two-step Anaerobic Digestion Model. Mathematical Biosciences and Engineering, 2020, 17 (5), pp.5504-5533. 10.3934/mbe.2020296. hal-02657259v2

HAL Id: hal-02657259

https://hal.inrae.fr/hal-02657259v2

Submitted on 31 May 2022

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_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_3.png" /></figure>

Research article

# Effect of a new variable integration on steady states of a two-step Anaerobic Digestion Model

Boumediene Benyahia<fnref n="1" />,* and Tewfik Sari<fnref n="2" />

* Correspondence: Email: boumediene.benyahia@univ-tlemcen.dz; Tel: +213555267958; Fax: +21343410011.

**Abstract:** This paper deals with a mathematical analysis of two-steps model of anaerobic digestion process, including dynamics of soluble microbial products (SMP). We propose to investigate effects of the new variable SMP on qualitative properties of the process in different generic cases. Equilibria of the model are graphically established considering qualitative properties of the kinetics and, their stability are proved theoretically and/or verified by numerical simulations. It will shown that the model has a rich qualitative behavior as equilibria bifurcation and multi-stability according to the considered bifurcation parameter.

**Keywords:** anaerobic digestion; bioprocess; equilibria bifurcation; modeling; stability; steady state analysis

## 1. Introduction

Mathematical modeling of bioprocesses is a powerful tool to (i) explain observed phenomena, (ii) understand some mechanisms of the system and predict its evolution, (iii) better control the process operations and (iv) build the "roots" for dialogue and discussion with biologists. In recent years, many studies were carried out on the mathematical models analysis of biological ecosystems using chemostat. A number of mathematical modeling methods that are relevant to the field of microbial ecology and bioprocesses was presented in [1]. Di and Yang [2], evaluated how structures and parametrization of synthetic microbial communities with two or three species could affect their productivity and stability. Qualitative analysis of local and global stability of steady states of a syntrophic relationship between two consortium of bacteria in a chemostat is detailed in [3, 4].

The technology of Anaerobic Digestion is highly promising with the potential to substantially improve efficiency in wastewater treatment, digestate handling and bioenergy production. Anaerobic

<footnote n="1">¹ Laboratoire d'Automatique de Tlemcen, Université de Tlemcen, Algeria</footnote>

<footnote n="2">² ITAP, University of Montpellier, INRAE, Institut Agro, Montpellier, France</footnote>
digestion is a complex process, which is widely described by the most complete model ADM1 (Anaerobic Digestion Model n.1) [5]. Because of its high complexity and strong non-linearity, ADM1 cannot be used for analytical analysis of the steady states of the system. In the literature, a number of studies have been made on equilibria and the nature of their stability of reduced and simplified models of anaerobic digestion processes using operating diagram analysis, which allows to describe the behavior of the system with respect to the control parameters. Khedim et al. [6], investigated how operating parameters (dilution rate and substrate inflow concentration) could ensure an optimal production of biogas in a Microalgae Anaerobic Digestion process. As regards [7], authors showed that the stability of the positive equilibrium of a two-tiered microbial food-chain is not affected when maintenance is included in the model and for a large class of kinetics. A generalised form of a three-tiered microbial food-web was proposed in [8]; when maintenance is not considered in the model, it was shown that one can explicitly determine the stability of the system and, boundaries between the different stability regions are characterized by analytical expressions.

A review of mathematical modeling of anaerobic digestion with respect to the theory, applications and technologies is given in [9], where it is argued that mathematical analysis tools can be appropriately applied to reduced-order models of anaerobic digestion to investigate the qualitative behavior of the system. Even if modeling of anaerobic digestion is increasing in complexity and new challenges should be addressed [10], for a simplified modeling, the biological process may be described mainly by two-steps reactional framework as given in [11]: in the first step (acidogenesis), the acidogenic bacteria consume the organic substrate and produce Volatile Fatty Acids (VFA) and CO2, while in the second step (methanogenesis), the methanogenic population consumes VFA and produces methane and CO2. A well known model for such process is the AM2 model [11] which has four main variables (two substrates and two microbial populations). In [12], it is shown that this model of two reactions represents 97.8% of biological variability, which justifies its choice to describe the main mass transfer within the bioreactor. An extended version of the AM2 model was used in [13] to predict biogas and methane production rates. Also, AM2 was compared to the ADM1 and it was shown that a tradeoff has to be made between model complexity and tractability. The AM2 model can successfully support on-line control and supervision strategies, based on state observers and feedback control [14, 15]. These literature examples of some applications of the AM2 model, show that this simple model is able to predict the main dynamical behavior of ADM1, which would be considered as a virtual anaerobic bioreactor for simulation.

Many mathematical studies were carried out on the qualitative behavior of the AM2 model in generic cases [16, 17], or in particular cases [18–20]. It is shown in [16] that the AM2 model can have at most six equilibria and it can have a monostability or a bistability behavior, according to the functioning conditions. A comparison of performance of one-reactor vs two-reactors configurations for a two-reaction (acidogenesis and methanogenesis) anaerobic digestion model were discussed in [21]. Using the AM2 model, authors have proven that separation of the reactions in two bioreactors does not improve the stability of the process nor the soluble organic matter removal capacity. Weedermann studied the effects of an external toxin on the behavior of a two-step model of anaerobic digestion [22]. He showed under what conditions the toxin can alter the steady states of the system (wash-out of bacteria, fluctuations (limit cycles) or bistabilities).

Even if the AM2 model has proven its usefulness for the control and supervision of anaerobic digesters, it remains a very simple model, which would not be able to explain certain biological
phenomena as the dynamics of lower concentrations or small bacteria populations [23] and, which has a limited applicability as typically the case for anaerobic digestion of waste-activated sludge [24]. This is why more or less extented versions of AM2 have been proposed in the literature in order to better describe anaerobic digestion processes with the integration of new main variables, while remaining simple from a mathematical modeling point of view. For instance in [24], one proposed the AM2HN model, which is a modification of AM2 by adding one additional state variable $X_T$ (total particulate substrate), i.e., one additional differential equation in order to include the disintegration/hydrolysis step and, initial differential equations of AM2 was accordingly modified. The model proposed in [22] is exactly a perturbation of the AM2 model to study the effects of an externally introduced toxin. Author added one differential equation of the dynamic of the toxin $T$, which inhibits the growth of bacteria $X_1$ while it is broken by bacteria $X_2$.

A model for anaerobic membrane digesters has been proposed in [25] for control design purposes. This model named AM2b is based on the modification of the two step model AM2 and integrates the dynamic of a new variable (SMP: Soluble Microbial Products) in the system. Recently, the AM2b model was combined with a simple fouling models to describe both biological and membrane dynamics in an Anaerobic Membrane BioReactors (AnMBR) [26, 27] and to assess system performance and membrane fouling [28]. A state-of-the-art on coupling of membrane fouling models with biological dynamics is provided in [29]. Authors reviewed modeling and control aspects of AnMBR and, focused on existing challenges and future perspectives to improve them. Stochastic versions of the AM2b model was proposed in [30] and [31], to provide a deeper description of the process when modeling lower concentrations or small bacteria populations, which can be seen as uncertain and noisy dynamics.

It is shown in [25] that the AM2b model is highly sensitive to the maximum growth rate of acidogenic bacteria on SMP (which is considered as a bifurcation parameter). This paper is complementary to [25] and it proposes a detailed mathematical analysis of the qualitative behavior of the model AM2, especially with respect to the bifurcation parameter. It reports briefly some results which presented in [25] and, gives pertinently their mathematical backgrounds. The paper is structured as follows: first, we recall the AM2b model and we prove positivity and boundedness of its variables. Then, we characterize equilibria in some generic cases and we explain the background of their graphical determination. Finally, we investigate through numerical simulation equilibria and their stability of the system, before conclusions and perspectives are formulated.

## 2. Mathematical model

### 2.1. Mass balance equations

In Figure 1, we give a schematic representation of the anaerobic membrane bioreactor for which the model (2.5–2.9) is proposed below and, where the membrane retention of soluble and particulate components is illustrated.

We consider the anaerobic mathematical model AM2b presented in [25], where we have four reaction networks:

$$k_1 S_1 \xrightarrow{\mu_1(.)X_1} X_1 + k_2 S_2 + b_3 S + k_4 CO_2 \quad (2.1)$$

$$k_3 S_2 \xrightarrow{\mu_2(.)X_2} X_2 + b_4 S + k_5 CO_2 + k_6 CH_4 \quad (2.2)$$
$$b_1 S \xrightarrow{\mu(.)X_1} X_1 + b_2 S_2 + k_7 CO_2 \qquad (2.3)$$

$$D_0 X_1 \rightarrow D_0 S, \quad D_0 X_2 \rightarrow D_0 S \qquad (2.4)$$

<figure><img src="image_4.png" /><figcaption>**Figure 1.** Schematic representation of the compartment bioreactor-membrane.</figcaption></figure>

In the first reaction, the substrate $S_1$ (organic matter) is degraded into substrates $S_2$ (Volatile Fatty Acids) and $S$ (SMP) by acidogenic bacteria $X_1$ and then in the second reaction, $S_2$ is converted into $S$ by methanogenic bacteria $X_2$. The third reaction network consist in degrading $S$ into $S_2$ by the consortium $X_1$. A part of $S$ is produced from biomasses decay. During reactions (2.1), (2.2) and (2.3), there is a production of biogas.

Mass balance equations are given by:

$$\dot{S}_1 = D(S_{1\text{in}} - S_1) - k_1\mu_1(S_1)X_1, \qquad (2.5)$$

$$\dot{X}_1 = (\mu_1(S_1) + \mu(S) - D_0 - D_1)X_1, \qquad (2.6)$$

$$\dot{S}_2 = D(S_{2\text{in}} - S_2) - k_3\mu_2(S_2)X_2 + (k_2\mu_1(S_1) + b_2\mu(S))X_1, \qquad (2.7)$$

$$\dot{X}_2 = (\mu_2(S_2) - D_0 - D_1)X_2, \qquad (2.8)$$

$$\dot{S} = (b_3\mu_1(S_1) + D_0 - b_1\mu(S))X_1 + (b_4\mu_2(S_2) + D_0)X_2 - MS, \qquad (2.9)$$

where $S_{1\text{in}}$ and $S_{2\text{in}}$ are input substrate concentrations, $D$, $D_0$ and $D_1$ are the dilution rate, the decay rate of biomass and the withdraw rate respectively. $M = [\beta D + (1 - \beta)D_1]$, where $\beta \in [0, 1]$ represents the fraction of $S$ leaving the bioreactor (see [25] for more detail on the model development). Parameters $k_i$ and $b_i$ are pseudo-stoichiometric coefficients associated to the bioreactions, which represent degradation and production rates of different substrates. The identifiability and the estimation of such parameters are discussed in [11] and [12].

We make the following matter conservation principles:

* over a given period of time, the quantity of biomass (or products) produced is always smaller than the quantity of substrate degraded. Thus, from (2.1–2.3) one has:

$$k_1 \geq 1 + b_3 + k_2, \qquad (2.10)$$

$$k_3 \geq 1 + b_4, \qquad (2.11)$$

$$b_1 \geq 1 + b_2. \qquad (2.12)$$

* the quantity $S_2$ produced from $S_1$ is higher than the quantity produced from the SMP (see (2.1) and (2.3)):

$$k_2 > b_2. \qquad (2.13)$$
The kinetics $\mu_1, \mu_2$ and $\mu$ are assumed to be dependent on $S_1, S_2$ and $S$ respectively, satisfying the following hypotheses:

**Hypothesis 2.1.** $\mu_1(S_1)$ and $\mu(S)$ are of class $C^1$ and satisfy the following properties:

* $\mu_1(0) = \mu(0) = 0$,

* $\mu'_1(S_1) > 0$ and $\mu'(S) > 0$ for $S_1 > 0$ and $S > 0$ respectively,

* $\mu_1(\infty) = m_1$ and $\mu(\infty) = m$.

**Hypothesis 2.2.** $\mu_2(S_2)$ is of class $C^1$ and satisfies the following properties:

* $\mu_2(0) = \mu_2(\infty) = 0$,

* $\mu_2(S_2)$ has a maximum $\mu_2(S_2^M) > 0$ for $S_2 = S_2^M$,

* $\mu'_2(S_2) > 0$ for $0 < S_2 < S_2^M$,

* $\mu'_2(S_2) < 0$ for $S_2 > S_2^M$.

The model analysis given in this paper, is valid for all functions verifying the hypotheses (2.1) and (2.2). Examples of functions satisfying these assumptions are (see appendix 1 of [32]):

* The Monod kinetics $\mu(\xi) = m\frac{\xi}{\xi+K}$, the Tessier kinetics $\mu(\xi) = m(1 - e^{-\frac{\xi}{K}})$, the Moser or the Ming et al. kinetics $\mu(\xi) = \frac{m\xi^2}{K+\xi^2}$ (with $m$ and $K$ are constants), which all satisfy hypothesis 2.1.

* The Haldane kinetics $\mu(\xi) = m\frac{\xi}{\frac{\xi^2}{K_i}+\xi+K}$, or the function $\mu(\xi) = K(e^{-\alpha_1\xi} - e^{-\alpha_2\xi})$ (with $m, K, K_i$ and $\alpha_2 > \alpha_1$ are constants), which satisfy hypothesis 2.2.

Positivity and boundedness are very important properties for biological systems. We have to check that for zero or positive initial conditions, all variables of system (2.5–2.9) are non-negative and bounded for all time.

**Proposition 2.3.** *The variables $(S_1, X_1, S_2, X_2, S)$ of system (2.5–2.9) are positive and bounded.*

*Proof.* The proof is given in Appendix A.1. □

### 3. Equilibria of model

The equilibria of system are solutions of the following nonlinear algebraic system:

$$0 = D(S_{1\text{in}} - S_1) - k_1\mu_1(S_1)X_1 \qquad (3.1)$$

$$0 = [\mu_1(S_1) + \mu(S) - D_0 - D_1]X_1 \qquad (3.2)$$

$$0 = D(S_{2\text{in}} - S_2) - k_3\mu_2(S_2)X_2 + [k_2\mu_1(S_1) + b_2\mu(S)]X_1 \qquad (3.3)$$

$$0 = [\mu_2(S_2) - D_0 - D_1]X_2 \qquad (3.4)$$

$$0 = [b_3\mu_1(S_1) + D_0 - b_1\mu(S)]X_1 + [b_4\mu_2(S_2) + D_0]X_2 - MS \qquad (3.5)$$

We use the following notations:

$$A = \frac{b_4(D_0 + D_1) + D_0}{k_3(D_0 + D_1)}, \qquad B = \frac{M}{D} = [\beta + (1 - \beta)\frac{D_1}{D}]. \qquad (3.6)$$
If $D_0 + D_1 < \mu_2(S_2^M)$ then $S_2^{1*} < S_2^M < S_2^{2*}$ are the roots of equation $\mu_2(S_2) = D_0 + D_1$ and, we note:

$$ \alpha_i := \frac{A}{B} (S_{2\text{in}} - S_2^{i*}), \quad \beta_i = \frac{D}{k_3(D_0 + D_1)} (S_{2\text{in}} - S_2^{i*}), \quad i = 1, 2 \qquad (3.7) $$

From Eq (3.2) one deduce that $X_1 = 0$ or $\mu_1(S_1) + \mu(S) = D_0 + D_1$. The following lemma describes the equilibria points for which $X_1 = 0$, that is to say, there is a washout of $X_1$.

**Lemma 3.1.** The equilibria $(S_1^*, 0, S_2^*, X_2^*, S^*)$ of the system (2.5–2.9) are given by:

* the washout equilibrium of $X_1$ and $X_2$, $E_0^0 = (S_{1\text{in}}, 0, S_{2\text{in}}, 0, 0)$, which always exists,

* the washout equilibrium of $X_1$ but not of $X_2$,

$$ E_1^i = (S_{1\text{in}}, 0, S_2^{i*}, X_2^{i*}, S^{i*}), \quad i = 1, 2 $$

where $S_2^{i*}$ are the roots of equation $\mu_2(S_2) = D_0 + D_1$, $X_2^{i*}$ and $S^{i*}$ are given by the formulas:

$$ X_2^{i*} = \beta_i, \quad S^{i*} = \alpha_i, \quad i = 1, 2. $$

The equilibrium $E_1^i$ exists if and only if:

$$ S_{2\text{in}} > S_2^{i*}. \qquad (3.8) $$

*Proof.* The proof is given in Appendix A.2.1. □

**Remark 1.** In this paper, we present a detailed mathematical analysis of the model equilibria. From biological realism point of view, equilibria $E_1^i$ would not occur except in certain cases. Indeed, $S_2$ (VFA) available for the reaction (2.2) is produced in reactions (2.1) and (2.3) when bacteria $X_1$ degrade $S_1$ (organic substrate) and $S$ (SMP). Also, it can come from outside the bioreactor in $S_{2\text{in}}$ (see Figure 1). Often, this is not possible in the biological realism, unless we consider a third acetogenic microorganisms which produce $S_{2\text{in}}$ from external organic matter or, if we carry out a bench-scale study, by introducing $S_{2\text{in}}$ into the bioreactor.

Now, we consider equilibria for which there is no washout of $X_1$ but washout of $X_2$. We introduce the following notations:

$$ F(S) := \mu_1^{-1} (D_0 + D_1 - \mu(S)), \qquad (3.9) $$

$$ G(S_1) := (S_{1\text{in}} - S_1) \left( B_1 - \frac{B_2}{\mu_1(S_1)} \right), \qquad (3.10) $$

where:

$$ B_1 = \frac{b_1 + b_3}{k_1\beta}, \quad B_2 = \frac{b_1(D_0 + D_1) - D_0}{k_1\beta}, \qquad (3.11) $$

**Lemma 3.2.** Let $E_2^0 = (S_1^*, X_1^*, S_2^*, 0, S^*)$ an equilibrium point of the system (2.5–2.9), such that $X_1^* > 0$. Then $S_1^*$ and $S^*$ are solutions of the system of equations:

$$ \begin{cases} S_1 = F(S) \\ S = G(S_1) \end{cases} \qquad (3.12) $$
and $X_1^*$ and $S_2^*$ are given by the formulas:

$$X_1^* = \frac{D}{k_1\mu_1(S_1^*)} (S_{1\text{in}} - S_1^*), \quad S_2^* = S_{2\text{in}} + \frac{k_2\mu_1(S_1^*) + b_2\mu(S^*)}{k_1\mu_1(S_1^*)} (S_{1\text{in}} - S_1^*).$$

The equilibrium $E^*$ exists if and only if:

$$S_{1\text{in}} > S_1^*. \qquad (3.13)$$

*Proof.* The proof is given in Appendix A.2.2. □

Now, we consider equilibria for which there is no washout of $X_1$ nor $X_2$. We introduce the following notations:

$$H(S_1) := (S_{1\text{in}} - S_1) \left( C_1 - \frac{C_2}{\mu_1(S_1)} \right), \qquad (3.14)$$

$$H_i(S_1) := \alpha_i + H(S_1), \quad i = 1, 2. \qquad (3.15)$$

where:

$$C_1 = B_1 + \frac{A(k_2 - b_2)}{k_1\beta}, \quad C_2 = B_2 - \frac{Ab_2}{k_1\beta}, \qquad (3.16)$$

**Lemma 3.3.** Let $E_2^i = (S_1^*, X_1^*, S_2^*, X_2^*, S^*)$, $i = 1, 2$ an equilibrium point of the system (2.5–2.9) such that $X_1^* > 0$ and $X_2^* > 0$. Then one has $S_2^i$, $i = 1, 2$ are the roots of equation $\mu_2(S_2) = D_0 + D_1$, and $S_1^*$ and $S^*$ are solutions of the system of equations:

$$\begin{cases} S_1 = F(S), \\ S = H_i(S_1), \quad i = 1, 2. \end{cases} \qquad (3.17)$$

and $X_1^*$ and $X_2^*$ are given by the formulas:

$$X_1^* = \frac{D}{k_1\mu_1(S_1^*)} (S_{1\text{in}} - S_1^*), \quad X_2^* = \beta_i + \frac{D}{k_3(D_0 + D_1)} \frac{k_2\mu_1(S_1^*) + b_2\mu(S^*)}{k_1\mu_1(S_1^*)} (S_{1\text{in}} - S_1^*)$$

The equilibrium $E^*$ exists if and only if the following conditions hold:

$$S_{1\text{in}} > S_1^* \text{ and } H_i(S_1^*) > G(S_1^*), \quad i = 1, 2. \qquad (3.18)$$

*Proof.* The proof is given in Appendix A.2.3. □

*Remark 2.* When the system (3.12) or the system (3.17) has several solutions $(S_{1j}^*, S_j^*)$, one notes $E_{2j}^0$ (respectively $E_{2j}^1$ and $E_{2j}^2$), $j = 1, 2$ the corresponding equilibria (see section 6.3).

# 4. Graphical determination of equilibria

Equilibria of system (2.5–2.9) are determined, by finding graphically solutions of system (3.12) and (3.17). Values of $S_1^*$ and $S^*$ should be positive and satisfy conditions (3.13) and (3.18). Thus, we should study sign of functions $G(S_1)$ and $H_i(S_1)$ and, specify the domain where they are positive.

First, let us give the following lemma:
**Lemma 4.1.** We have $\lambda_H < \lambda_G < \lambda_1$ where $\lambda_H$, $\lambda_G$ and $\lambda_1$ are defined by:

$$ \lambda_1 = \mu_1^{-1}(D_0 + D_1), \quad \lambda_G = \mu_1^{-1}(D_G), \quad \lambda_H = \mu_1^{-1}(D_H), $$

with $D_G = B_2/B_1$ and $D_H = C_2/C_1$.

*Proof.* The proof is given in Appendix A.3. □

In Figure 2 on the left, we illustrate positions of $\lambda_H$, $\lambda_G$ and $\lambda_1$. On the right, we show solutions $S_2^{i*}$, $i = 1, 2$ of equation $\mu_2(S_2) = D_0 + D_1$.

<figure><img src="image_5.png" /><figcaption>**Figure 2.** Positions of values $\lambda_1$, $\lambda_G$ and $\lambda_H$ (left). Solutions $S_2^{i*}$, $i = 1, 2$ of $\mu_2(S_2) = D_0 + D_1$ (right).</figcaption></figure>

The function $G(S_1)$ defined by (3.10) is positive for $S_1$ between $S_{1\text{in}}$ and $\lambda_G$, the root of:

$$ g(S_1) = B_1 - \frac{B_2}{\mu_1(S_1)} $$

We have two cases (see Figure 3):

* $D_G > m_1$, where $g(S_1)$ is always negative for $S_1 < S_{1\text{in}}$ (Figure 3, left) and values of $S_1 > S_{1\text{in}}$ do not satisfy the condition (3.13), or

* $D_G < m_1$, where $g(S_1) > 0$ if and only if $S_1 > \lambda_G$ (Figure 3, right).

The case corresponding to Figure 3, center, is not considered since it does not satisfy the condition (3.13).

<figure><img src="image_6.png" /><figcaption>**Figure 3.** Graphical representation of $G(S_1)$, left: if $D_G > m_1$, center: if $D_G < m_1$ and $S_{1\text{in}} < \lambda_G$, right: if $D_G < m_1$ and $S_{1\text{in}} > \lambda_G$. Gray area represents zone where (3.13) is not satisfied.</figcaption></figure>
**Proposition 4.2.** A necessary condition for (3.12) to have positive solutions is $\lambda_G < S_{1\text{in}}$ that is to say $\mu_1(S_{1\text{in}}) > D_G$.

*Proof.* The proof is given in Appendix A.4. □

The function $H(S_1)$ defined by (3.14) is positive for $S_1$ between $S_{1\text{in}}$ and $\lambda_H$, the root of:

$$H(S_1) = C_1 - \frac{C_2}{\mu_1(S_1)}$$

Two cases can be distinguished:

* $D_H > m_1$, where $H(S_1)$ is always negative for $S_1 < S_{1\text{in}}$, or

* $D_H < m_1$, where $H(S_1) > 0$ if and only if $S_{1\text{in}} > S_1 > \lambda_H$ (see Figure 4).

In the following, we assume that $\mu_1(S_{1\text{in}}) > D_G$. Let us notice that:

* $H(S_1)$ is positive if and only if $\lambda_H < S_1 < S_{1\text{in}}$,

* $G(S_1)$ is positive if and only if $\lambda_G < S_1 < S_{1\text{in}}$,

* $H(S_1) > G(S_1)$ for all $\lambda_H < S_1 < S_{1\text{in}}$.

**Proposition 4.3.** The equilibrium $E_2^i$, $i = 1,2$ exists if and only if the graph of $H_i(S_1)$ intersects the axis of $S_1$ on the right of $S_{1\text{in}}$.

*Proof.* The proof is given in Appendix A.5. □

<figure><img src="image_7.png" /><figcaption>**Figure 4.** Graphical representation of $H(S_1)$, $H_i(S_1)$, $i = 1,2$ and $G(S_1)$, left: $H(S_1)$ and $G(S_1)$, center: $\alpha_i > 0$ thus $H_i(S_1) > H(S_1)$, right: $\alpha_i < 0$ thus $H_i(S_1) < H(S_1)$. Gray area represents zone where (3.13) is not satisfied. Character $T$ is to say *Trivial Equilibria* given by Lemma 3.1.</figcaption></figure>

## 5. Necessary conditions of existence of equilibria

The existence of equilibria depend on the relative positions of the value of $S_{1\text{in}}$ and the values of $\lambda_H$, $\lambda_G$ and $\lambda_1$ (see Figure 2). We have four cases:

* $S_{1\text{in}} < \lambda_H < \lambda_G < \lambda_1$

* $\lambda_H < S_{1\text{in}} < \lambda_G < \lambda_1$
* $\lambda_H < \lambda_G < S_{\text{lin}} < \lambda_1$

* $\lambda_H < \lambda_G < \lambda_1 < S_{\text{lin}}$

Recall that an equilibrium exists if and only if the conditions (3.18) are satisfied. We list in the Table 1 the possible existence of equilibria in the four above cases. In all the following figures, the dotted vertical line represents $\lambda_1 = F(0)$, the blue graph represents $G(S_1)$, the green one represents $H(S_1)$ and those in red represent $H_1(S_1)$ (top red graph) and $H_2(S_1)$ (bottom red graph).

*Remark 3.* The function $F(S)$ depends on $\mu(S)$, but functions $G(S_1)$ and $H_i(S_1)$, $i = 1, 2$ do not depend on it. For $\mu(0) = 0$, intersections of $F(0) = \lambda_1$ with $G(S)$ and $H_i(S)$ correspond to cases of [16] and [25] as mentioned in the last column of Table 1 and seen on Figure 8.

<table><caption>**Table 1.** Existence of equilibria in the four cases. The symbol 'X' indicates that the equilibrium can exist. If there is no symbol, that indicates that equilibrium does not exist.</caption><thead><tr><th>Case</th><th>Figure</th><th>F ∩ G</th><th>F ∩ H<sub>1</sub></th><th>F ∩ H<sub>2</sub></th><th>Corresponding cases in [16] and/or in [25]</th></tr></thead><tbody><tr><td rowspan="3">S<sub>lin</sub> &lt; λ<sub>H</sub> &lt; λ<sub>G</sub> &lt; λ<sub>1</sub></td><td>Figure 5, left</td><td></td><td></td><td></td><td>1.1 of [16]</td></tr><tr><td>Figure 5, center</td><td></td><td>X</td><td></td><td>1.2 of [16]</td></tr><tr><td>Figure 5, right</td><td></td><td>X</td><td>X</td><td>1.3 of [16]</td></tr><tr><td rowspan="3">λ<sub>H</sub> &lt; S<sub>lin</sub> &lt; λ<sub>G</sub> &lt; λ<sub>1</sub></td><td>Figure 6, left</td><td></td><td>X</td><td>X</td><td>1.1 of [16]</td></tr><tr><td>Figure 6, center</td><td></td><td>X</td><td>X</td><td>2.1 of [16]</td></tr><tr><td>Figure 6, right</td><td></td><td>X</td><td>X</td><td>1.3 of [16]</td></tr><tr><td rowspan="3">λ<sub>H</sub> &lt; λ<sub>G</sub> &lt; S<sub>lin</sub> &lt; λ<sub>1</sub></td><td>Figure 7, left</td><td>X</td><td>X</td><td>X</td><td>1.1 of [16]</td></tr><tr><td>Figure 7, center</td><td>X</td><td>X</td><td>X</td><td>1.2 of [16]<br>C of [25]</td></tr><tr><td>Figure 7, right</td><td>X</td><td>X</td><td>X</td><td>1.3 of [16]<br>B of [25]</td></tr><tr><td rowspan="6">λ<sub>H</sub> &lt; λ<sub>G</sub> &lt; λ<sub>1</sub> &lt; S<sub>lin</sub></td><td>Figure 8, top left</td><td>X</td><td>X</td><td>X</td><td>2.1 of [16]</td></tr><tr><td>Figure 8, top center</td><td>X</td><td>X</td><td>X</td><td>2.2 of [16]</td></tr><tr><td>Figure 8, top right</td><td>X</td><td>X</td><td>X</td><td>2.3 of [16]</td></tr><tr><td>Figure 8, bottom left</td><td>X</td><td>X</td><td>X</td><td>2.4 of [16]</td></tr><tr><td>Figure 8, bottom center</td><td>X</td><td>X</td><td>X</td><td>2.5 of [16]</td></tr><tr><td>Figure 8, bottom right</td><td>X</td><td>X</td><td>X</td><td>2.6 of [16]<br>A of [25]</td></tr></tbody></table>

Case: $S_{\text{lin}} < \lambda_H < \lambda_G < \lambda_1$

In the case 1, illustrated by Figure 5 on the left, we have $H > H_1 > H_2$ and, intersections $F \cap H_i$ and $F \cap G$ do not give any positive equilibria, because it does not satisfy the condition (3.13). In the case 2, illustrated by Figure 5 on the center, we have $H_1 > H > H_2$. The equilibrium of $F \cap H_1$ can exist, but there are no equilibria of $F \cap H_2$ and $F \cap G$. The last case 3, represented by Figure 5 on the right, we have $H_1 > H_2 > H$. There is a possibility of existence of equilibria $F \cap H_1$ and $F \cap H_2$ for values of $m$ enough high, but equilibria of $F \cap G$ do not exist.
<figure><img src="image_8.png" /><figcaption>**Figure 5.** Different graphical representations corresponding to the case $S_{\text{lin}} < \lambda_H$, blue: $G$, top red: $H_1$, bottom red: $H_2$, green: $H$.</figcaption></figure>

*Remark 4.* When we have intersection of the function $F(S)$ with functions $H_1(S_1)$ and $H_2(S_1)$ at $S_1 = S_{\text{lin}}$, then we obtain equilibria $E_1^i = (S_{\text{lin}}, 0, S_2^{i*}, X_2^{i*}, S^{i*})$ where $S^{i*} = \alpha_i = H_i(S_{\text{lin}})$, $i = 1, 2$ as it can be seen on Figure 5, center, for $E_1^i$ and on Figure 5, right, for $E_1^i$, $i = 1, 2$.

*Case :* $\lambda_H < S_{\text{lin}} < \lambda_G < \lambda_1$

This case is illustrated by Figure 6, left, center and right for $H > H_1 > H_2$, $H_1 > H > H_2$ and $H_1 > H_2 > H$ respectively. Equilibria of $F \cap H_i$, $i = 1, 2$ can exist for higher values of $m$, but not those of $F \cap G$ for all values of $m$.

<figure><img src="image_9.png" /><figcaption>**Figure 6.** Different graphical representations corresponding to the case $\lambda_H < S_{\text{lin}} < \lambda_G$, blue: $G$, top red: $H_1$, bottom red: $H_2$, green: $H$.</figcaption></figure>

*Case :* $\lambda_H < \lambda_G < S_{\text{lin}} < \lambda_1$

We represent this case by Figure 7, where all equilibria of $F \cap H_i$, $i = 1, 2$ and $F \cap G$ can exist since condition (3.13) is satisfied. Also, some equilibria bifurcations can occur for higher values of $m$.

*Case :* $\lambda_H < \lambda_G < \lambda_1 < S_{\text{lin}}$

Here we have rich situations, equilibria for $F \cap G$ exist always, while $F \cap H_1$ and $F \cap H_2$ may give both equilibria for all $m$ (see Figure 8, top-right, bottom-center and bottom-right), only $F \cap H_1$ gives always equilibria (see Figure 8, top-center and bottom-left) or there is equilibria bifurcations for large values of $m$ for $F \cap H_1$ and/or $F \cap H_2$ (see Figure 8, top-left for $F \cap H_1$, top-center and bottom-left for $F \cap H_2$).
<figure><img src="image_10.png" /><figcaption>**Figure 7.** Different graphical representations corresponding to the case $\lambda_G < S_{1in} < \lambda_1$, blue: G, top red: $H_1$, bottom red: $H_2$, green: $H$.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Figure 8.** Different graphical representation corresponding to the case $\lambda_1 < S_{1in}$, blue: G, top red: $H_1$, bottom red: $H_2$, green: $H$.</figcaption></figure>

# 6. Existence and stability of equilibria

## 6.1. Stability of equilibria with washout of $X_1$

For trivial equilibria given by lemma 3.1, the results on their stability are summarized in Theorem 6.1.

**Theorem 6.1.** Existence and stability of washout equilibria of $X_1$ are as follows:

1. The equilibrium $E_0^0$ exists always and it is stable if and only if:

$$ \mu_1(S_{1in}) < D_0 + D_1 \text{ and, } \mu_2(S_{2in}) < D_0 + D_1 \quad (6.1) $$

2. The equilibrium $E_1^2$ exist if and only if $S_{2in} > S_{2in}^{2*}$ and it is always unstable.
3. The equilibrium $E_1^1$ exist if and only if $S_{2\text{in}} > S_2^{1*}$ and it is stable if and only if:

$$ \mu_1(S_{1\text{in}}) + \mu(S^{1*}) < D_0 + D_1 \quad (6.2) $$

*Proof.* The proof is given in Appendix A.6. □

The condition (6.2) may be graphically explained on the Figure 9. The graph $f(S_1, S) = \mu_1(S_1) + \mu(S) - D_0 - D_1 = 0$ separates the plane $(S_1, S)$ into two zones:

* Zone $Z_0$: where $f(S_1, S) > 0$,

* Zone $Z_1$: where $f(S_1, S) < 0$.

According to the condition (6.2), the equilibrium $E_1^1$ is stable if and only if: $E_1^1 \in Z_1$ (the case represented on left in Figure 9).

<figure><img src="image_12.png" /><figcaption>**Figure 9.** Condition of stability of the equilibrium $E_1$ (Left: $E_1^1 \in Z_1$, stable. Right: $E_1^1 \in Z_0$, unstable).</figcaption></figure>

## 6.2. Stability of equilibria with washout of $X_2$ and stability of equilibria with $X_1 > 0$ and $X_2 > 0$

Here, we improve numerical simulations to check the system stability. Values of model parameters are chosen as in Tables 3 and 4 except the parameter $m$. According to the considered generic case 1, 2 or 3 represented by Figure 10, 11 and 12 respectively, the value of the parameter $m$ is varying in specific intervals for which we could have all possible equilibria bifurcations. Stability nature does not depend on values of $m$ in those intervals (see the column *Condition* in Table 5 for values of $m$). Then, we proceed as follows:

* Develop the Jacobian matrix $J$ of system (2.5–2.9) as given by (6.3),

* Evaluate this matrix for each equilibrium characterized by lemma 3.2 or 3.3,

* Develop the characteristic equation of the evaluated matrix,

* Use Routh-Herwitz criterion to analyze the system stability by numerical simulations (plot the coefficients of the first column of the Routh Table 2 according to $m$).
The jacobian matrix of system (2.5–2.9) evaluated at equilibria is given by (6.3).

$$
J = \begin{bmatrix}
-D - k_1 \mu_1'(S_1^*) X_1^* & -k_1 \mu_1(S_1^*) & 0 & 0 & 0 \\
\mu_1'(S_1^*) X_1^* & \mu_1(S_1^*) + \mu(S^*) - D_0 - D_1 & 0 & 0 & 0 \\
k_2 \mu_1'(S_1^*) X_1^* & k_2 \mu_1(S_1^*) + b_2 \mu(S^*) & -D - k_3 \mu_2'(S_2^*) X_2^* & -k_3 \mu_2(S_2^*) & b_2 \mu'(S^*) X_1^* \\
0 & 0 & \mu_2'(S_2^*) X_2^* & \mu_2(S_2^*) - D_0 - D_1 & 0 \\
b_3 \mu_1'(S_1^*) X_1^* & b_3 \mu_1(S_1^*) + D_0 - b_1 \mu(S^*) & b_4 \mu_2'(S_2^*) X_2^* & b_4 \mu_2(S_2^*) + D_0 & -M - b_1 \mu'(S^*) X_1^*
\end{bmatrix} \quad (6.3)
$$

Which can be symbolized as follows:

$$
J = \begin{bmatrix}
j_{11} & j_{12} & 0 & 0 & 0 \\
j_{21} & j_{22} & 0 & 0 & j_{25} \\
j_{31} & j_{32} & j_{33} & j_{34} & j_{35} \\
0 & 0 & j_{43} & j_{44} & 0 \\
j_{51} & j_{52} & j_{53} & j_{54} & j_{55}
\end{bmatrix} \quad (6.4)
$$

We can distinguish two cases according to lemma 3.2 where $X_1 > 0$ and $X_2 = 0$ or, lemma 3.3 where $X_1 > 0$ and $X_2 > 0$.

• In the case $X_2 = 0$, one has: $j_{33} = -D$ and $j_{43} = j_{53} = 0$.

• In the case $X_1 > 0$ and $X_2 > 0$, one has: $j_{22} = j_{44} = 0$ (from (3.2) and (3.4)).

The characteristic equation of the linearized system of (2.5–2.9) is:

$$
|\lambda I - J| = 0 \Leftrightarrow \lambda^5 + a_1 \lambda^4 + a_2 \lambda^3 + a_3 \lambda^2 + a_4 \lambda + a_5 = 0. \quad (6.5)
$$

where $a_i$ are coefficients depending on $j_{ik}$, $(i, k = 1..5)$ given by (6.4). Now, one establishes the
following Routh table:

<caption>**Table 2.** Table of Routh for the linearized system of (2.5–2.9).</caption>

$$
\begin{vmatrix}
\lambda^5 & 1 & a_2 & a_4 & 0 \\
\lambda^4 & a_1 & a_3 & a_5 & 0 \\
\lambda^3 & n_1 & n_2 & 0 & 0 \\
\lambda^2 & l_1 & a_5 & 0 & 0 \\
\lambda^1 & r_1 & 0 & 0 & 0 \\
\lambda^0 & a_5 & 0 & 0 & 0 \\
0 & 0 & 0 & 0 & 0
\end{vmatrix}
$$

with:

$$
n_1 = \frac{a_1 a_2 - a_3}{a_1}, \quad n_2 = \frac{a_1 a_4 - a_5}{a_1}, \quad l_1 = \frac{n_1 a_3 - n_2 a_1}{n_1}, \quad r_1 = \frac{l_1 n_2 - a_5 n_1}{l_1}
$$

The Routh-Herwitz criterion imposes that all coefficients of the first column of the Table 2 must
have the same sign, i.e., they must be positive (because the first element of the column is positive).

$$
a_1 > 0, \quad n_1 > 0, \quad l_1 > 0, \quad r_1 > 0, \quad a_5 > 0 \qquad (6.6)
$$
## 6.3. Numerical simulations

To illustrate our approach, we improve numerical simulations. We present three generic cases illustrated by Figures 10, 11 and 12, which are obtained for the biological parameters values given in Tables 3 and 4 and, kinetics functions (6.7), satisfying hypotheses 2.1 and 2.2.

$$ \mu_1(S_1) = m_1 \frac{S_1}{S_1 + K_1}, \quad \mu(S) = m \frac{S}{S + K}, \quad \mu_2(S_2) = m_2 \frac{S_2}{\frac{S_2^2}{K_i} + S_2 + K_2}. \qquad (6.7) $$

If they exist, equilibria are noted on figures by:

* $E_{2j}^1$: equilibria given by the intersection of $F(S)$ with $H_1(S_1)$, $j = 1, 2$,

* $E_{2j}^2$: equilibria given by the intersection of $F(S)$ with $H_2(S_1)$, $j = 1, 2$,

* $E_{2j}^0$: equilibria given by the intersection of $F(S)$ with $G(S_1)$, $j = 1, 2$.

<table><caption>**Table 3.** Nominal values for the parameters of the AM2b model [25].</caption><thead><tr><th>Parameter</th><th>Value</th><th>Parameter</th><th>Value</th><th>Parameter</th><th>Value</th><th>Parameter</th><th>Value</th></tr></thead><tbody><tr><td>$m_1$</td><td>1.2</td><td>$\beta$</td><td>0.6</td><td>$b_1$</td><td>5</td><td>$m$</td><td>varying</td></tr><tr><td>$m_2$</td><td>1.5</td><td>$k_1$</td><td>25</td><td>$b_2$</td><td>0.6</td><td>$K$</td><td>3</td></tr><tr><td>$K_2$</td><td>0.3</td><td>$k_2$</td><td>15</td><td>$b_3$</td><td>7</td><td>$D$</td><td>1</td></tr><tr><td>$K_i$</td><td>0.9</td><td>$k_3$</td><td>16.08</td><td>$b_4$</td><td>5</td><td>$D_0$</td><td>0.25</td></tr></tbody></table>

<table><caption>**Table 4.** Values for adjustable parameters $K_1$, $S_{1in}$, $S_{2in}$ and $D_1$ for each figure.</caption><thead><tr><th>Parameter</th><th>Generic case 1 (Figure 10)</th><th>Generic case 2 (Figure 11)</th><th>Generic case 3 (Figure 12)</th></tr></thead><tbody><tr><td>$K_1$</td><td>10</td><td>16</td><td>18</td></tr><tr><td>$S_{1in}$</td><td>15</td><td>15</td><td>10</td></tr><tr><td>$S_{2in}$</td><td>1</td><td>1</td><td>0.6</td></tr><tr><td>$D_1$</td><td>0.4</td><td>0.4</td><td>0.25</td></tr></tbody></table>

The form of $F(S)$ changes according to the value of the parameter $m$, the maximum growth rate of $\mu(S)$ (see (3.9)). Consequently, $F(S)$ can have one or two intersections with each one of functions $H_1(S_1)$, $H_2(S_1)$ or $G(S_1)$ as illustrated in Figures 10, 11 and 12.

In Figure 10 corresponding to the generic case 1, we have only one equilibrium noted $E_1^1$, $E_1^2$ and $E_1^0$ for each intersection of $F$ with $H_1$, $H_2$ and $G$ respectively.
<figure><img src="image_13.png" /><figcaption>**Figure 10.** Intersection of the graph of $F(S)$ with the graphs $G(S_1)$ and $H_i(S_1)$ in the Case 1.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>**Figure 11.** Intersection of the graph of $F(S)$ with the graphs $G(S_1)$ and $H_i(S_1)$ in the Case 2.</figcaption></figure>
<figure><img src="image_15.png" /><figcaption>**Figure 12.** Intersection of the graph of $F(S)$ with the graphs $G(S_1)$ and $H_i(S_1)$ in the Case 3.</figcaption></figure>

In generic cases 2 and 3, we have equilibria bifurcation when $m$ varies. For some values $c_i$, $i = 1, \dots, 5$ of $m$, (of course, they are different between cases 2 and 3), the graph of $F(S)$ intersects graphs of $H_i(S_1)$ and $G(S_1)$ (see Figures 13 and 14), leading to the apparition of new equilibria. The reader can refer to [25, 33] for more details.

<figure><img src="image_16.png" /><figcaption>**Figure 13.** Values $c_i$, $i = 1, \dots, 5$ of $m$ giving equilibria bifurcation in the generic case 2.</figcaption></figure>
<figure><img src="image_17.png" /><figcaption>**Figure 14.** Values $c_i$, $i = 1..5$ of $m$ giving equilibria bifurcation in the generic case 3.</figcaption></figure>

<table><caption>**Table 5.** Equilibria and their nature in generic cases represented in Figures 10, 11 and 12. Values of $c_i$, $i = 1..5$ of the case 2 are different from those of the case 3.</caption><thead><tr><th rowspan="3">Cases</th><th rowspan="3">Condition</th><th colspan="8">Equilibria and nature</th></tr><tr><th colspan="3">F ∩ T</th><th colspan="2">F ∩ G</th><th colspan="2">F ∩ H<sub>1</sub></th><th colspan="2">F ∩ H<sub>2</sub></th></tr><tr><th>E<sub>0</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>21</sub><sup>0</sup></th><th>E<sub>22</sub><sup>0</sup></th><th>E<sub>21</sub><sup>1</sup></th><th>E<sub>22</sub><sup>1</sup></th><th>E<sub>21</sub><sup>2</sup></th><th>E<sub>22</sub><sup>2</sup></th></tr></thead><tbody><tr><td colspan="2"><strong>Case 1 (Figure 10)</strong></td><td>m ≥ 0</td><td>U</td><td>U</td><td>U</td><td>S</td><td></td><td>U</td><td></td><td>S</td></tr><tr><td colspan="2"><strong>Case 2 (Figure 11)</strong></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>2.1</td><td>0 ≤ m &lt; c<sub>1</sub></td><td>S</td><td>S</td><td>U</td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>2.2</td><td>c<sub>1</sub> &lt; m &lt; c<sub>2</sub></td><td>S</td><td>S</td><td>U</td><td></td><td></td><td>S</td><td>U</td><td></td><td></td></tr><tr><td>2.3</td><td>c<sub>2</sub> &lt; m &lt; c<sub>3</sub></td><td>S</td><td>S</td><td>U</td><td></td><td></td><td>S</td><td>U</td><td>U</td><td>U</td></tr><tr><td>2.4</td><td>c<sub>3</sub> &lt; m &lt; c<sub>4</sub></td><td>S</td><td>S</td><td>U</td><td>S</td><td>U</td><td>S</td><td>U</td><td>U</td><td>U</td></tr><tr><td>2.5</td><td>c<sub>4</sub> &lt; m &lt; c<sub>5</sub></td><td>S</td><td>U</td><td>U</td><td>S</td><td>U</td><td>S</td><td></td><td>U</td><td>U</td></tr><tr><td>2.6</td><td>c<sub>5</sub> &lt; m</td><td>S</td><td>U</td><td>U</td><td>S</td><td>U</td><td>S</td><td></td><td>U</td><td>U</td></tr><tr><td colspan="2"><strong>Case 3 (Figure 12)</strong></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>3.1</td><td>0 ≤ m &lt; c<sub>1</sub></td><td>U</td><td>S</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>3.2</td><td>c<sub>1</sub> &lt; m &lt; c<sub>2</sub></td><td>U</td><td>S</td><td></td><td></td><td></td><td>S</td><td>U</td><td></td><td></td></tr><tr><td>3.3</td><td>c<sub>2</sub> &lt; m &lt; c<sub>3</sub></td><td>U</td><td>S</td><td></td><td></td><td></td><td>S</td><td>U</td><td>U</td><td>U</td></tr><tr><td>3.4</td><td>c<sub>3</sub> &lt; m &lt; c<sub>4</sub></td><td>U</td><td>S</td><td></td><td>S</td><td>U</td><td>S</td><td>U</td><td>U</td><td>U</td></tr><tr><td>3.5</td><td>c<sub>4</sub> &lt; m &lt; c<sub>5</sub></td><td>U</td><td>S</td><td></td><td>S</td><td>U</td><td>S</td><td>U</td><td>U</td><td>U</td></tr><tr><td>3.6</td><td>c<sub>5</sub> &lt; m</td><td>U</td><td>U</td><td></td><td>S</td><td>U</td><td>S</td><td></td><td>U</td><td>U</td></tr></tbody></table>

Equilibria of system and their nature according to $m$ in the three generic cases are summarized in Table 5, where $T$ stands for *Trivial Equilibria* $E_0^0$, $E_1^1$ and $E_1^2$, $F ∩ H_1$, $F ∩ H_2$ and $F ∩ G$ stand for *Equilibria obtained by the intersections of the graph $F$ with graphs $H_1$, $H_2$ and $G$, respectively*, $S$ and
U stand for *Stable* and *Unstable* Equilibrium, respectively. If there is no symbol, then it means that the equilibrium does not exist.

Stability nature of equilibria corresponding to the washout of $X_2$ ($F \cap G$) and the existence of both $X_1$ and $X_2$ ($F \cap H_1$ and $F \cap H_2$) is checked by using the Routh-Herwitz criterion as detailed in the section 6.2.

On Figures 15, 16, and 17, we represent the coefficients of the first column of Table 2 with different colors: $a_1$ in black, $n_1$ in blue, $l_1$ in red, $r_1$ in magenta and $a_5$ in green. On Figures 16, and 17, vertical lines represent bifurcation values $c_i$, $i = 1..5$ of the parameter $m$ (they are different between the two figures). According to the considered case, coefficients are represented only for values of $m$, for which equilibria may exist. For instance, Routh coefficients for the equilibrium $E_{22}^1$ are represented on Figure 16, bottom-left, only for $c_1 \le m \le c_4$. If equilibrium is stable, then all coefficients must be positive in the corresponding interval of $m$.

**Remark 5.** Stability nature of equilibria $E_0^0$, $E_1^1$ and $E_1^2$ of the case 1 in Table 5, can be seen on Figure 10 as follows:

• $\lambda_1 < S_{\text{lin}} \Rightarrow \mu_1(\lambda_1) < \mu_1(S_{\text{lin}})$, that is to say $D_0 + D_1 < \mu_1(S_{\text{lin}}) \Rightarrow E_0^0$ is unstable according to theorem 6.1.1

• $(S_{\text{lin}}, S^{1*}) \in Z_0 \Rightarrow E_1^1$ is unstable according to condition (6.2) of theorem 6.1.3 and, Figure 9.

• $E_1^2$ does exist according to proposition 4.3 and, is unstable thanks to theorem 6.1.2.

**Remark 6.** Stability nature of the equilibrium $E_0^0$ of the cases 2 and 3 in Table 5, can be analyzed as follows:

• $\mu_1(S_{\text{lin}}) < D_0 + D_1$ as it is seen on Figure 11 and 12 for both cases.

• From parameters values in Tables 3 and 4 and, according to condition (6.1) of theorem 6.1 we have:

$\mu_2(S_{2\text{in}}) < D_0 + D_1$ for the case 2, thus the equilibrium $E_0^0$ is Stable.

$\mu_2(S_{2\text{in}}) > D_0 + D_1$ for the case 3, thus the equilibrium $E_0^0$ is Unstable.

As can be seen on figures, if they exist:

• The first equilibrium $E_{21}^1$ of $F \cap H_1$ and the first equilibrium $E_{21}^0$ of $F \cap G$ are always stable. Coefficients given by (6.6) of the first column of Routh table are always positive.

• The second equilibrium $E_{22}^1$ of $F \cap H_1$ and the second equilibrium $E_{22}^0$ of $F \cap G$ are always unstable. Some coefficients given by (6.6) are (or become) negative (for instance $a_5$ on Figure 16, second sub-Fig from top, is always negative, or $r_1$ in magenta on Figure 17, last sub-Fig, becomes negative).

• both equilibria $E_{21}^2$ and $E_{22}^2$ of $F \cap H_2$ are unstable. Some coefficients of (6.6) are always negative (for instance $r_1$ in magenta and $a_5$ in green on Figure 15).

At this stage of discussion about stability nature of possible equilibria, we give a conjecture on positive ones which are obtained for $F \cap H_i$, $i = 1, 2$.

**Conjecture 6.2.**

• Equilibria $E_{2j}^2$, $j = 1, 2$, resulting from $F \cap H_2$ are unstable if they exist.

• The only stable equilibrium $E_{21}^1$ resulting from $F \cap H_1$, is the one which corresponds to the smallest value of $S_1^*$.
<figure><img src="image_18.png" /><figcaption>**Figure 15.** The coefficients of the first column of the Routh Table 2 in the generic case 1: $a_1$ (—), $n_1$ (—), $l_1$ (—), $r_1$ (—) and $a_5$ (—).</figcaption></figure>
<figure><img src="image_19.png" /><figcaption>**Figure 16.** The coefficients of the first column of the Routh Table 2 in the generic case 2, (coordinates $m_{ci}$, $i = 1..5$ on the X-axis are the same for all the sub-figures): $a_1$ (—), $n_1$ (—), $l_1$ (—), $r_1$ (—) and $a_5$ (—)</figcaption></figure>
<figure><img src="image_20.png" /><figcaption>**Figure 17.** The coefficients of the first column of the Routh Table 2 in the generic case 3, (coordinates $m_{ci}$, $i = 1..5$ on the X-axis are the same for all the sub-figures): $a_1$ (—), $n_1$ (—), $l_1$ (—), $r_1$ (—) and $a_5$ (—).</figcaption></figure>
# 7. Conclusion

The simple model AM2 is widely used in the literature to describe anaerobic digestion in two-step of biological reactions. This model proved its ability to adequately predict dynamics of the main variables of the anaerobic digestion and, it was used with efficiency for control and supervision purposes. It has been shown in [16] that the AM2 model can have at most six steady states depending on its operating parameters. Nevertheless, the AM2 model was not able to simulate some phenomena in many practical biological experiments. There is why extended versions of the AM2 model were proposed in the literature, by integrating some few new variables.

In this paper we investigated the effect of a new variable $S$ (SMP: Soluble Microbial Product) integration on steady states of a two-step anaerobic digestion model. Indeed, this model initially proposed in [25] for control purposes, is an extension of the AM2 model for anaerobic membrane bioreactors. We consider the dynamics of five variables: two bacteria populations ($X_1, X_2$) and three substrates including the new variable ($S_1, S_2, S$), where one microorganism $X_1$ can growth on one substrate $S_1$ to produce both the second and the third substrates $S_2$ and $S$ and also on $S$ to produce $S_2$, while the second microorganism $X_2$ can only growth on $S_2$ and, could be inhibited by an excess quantity of this substrate. $S$ is produced from degradation of $S_1$ and $S_2$ and, death of bacteria $X_1$ and $X_2$. One important parameter which could considerably alter the system behavior is the maximum growth rate $m$ of the first bacteria population $X_1$ on the new substrate $S$. Indeed, this biological parameter is considered as a bifurcation parameter in addition to the conventional operating parameters which are the dilution rate $D$ and the inlet substrate concentrations $S_{1in}$ and $S_{2in}$.

In this paper, the model equilibria and their stability were analyzed analytically and using numerical simulations according to this bifurcation parameter $m$. We distinguished three generic cases accordingly to the system parameters values (Tables 3 and 4), where the system can exhibit rich qualitative behavior in terms of equilibria bifurcation and multistability. We have highlighted that in the first generic case (Figure 10), the behavior of the extended model is exactly similar to the AM2 one (i.e., six equilibria with bistability). While in the second generic case (Figure 11), for a set of parameters values, especially the maximum growth rate $m$ of $X_1$ on $S$ (the new variable), the system can have until nine equilibria where four of them are stable (multistability). In the third generic case (Figure 12), we can have eight equilibria with trystability.

Our study shows how the behavior of a two-step model of anaerobic digestion can be altered by the integration of the new variable in some generic cases and, how the model equilibria and their stability would be sensible to the bifurcation parameter $m$. Our results would be useful for both mathematicians and biologists communities and, could build the roots for dialogue between them as noted in the introduction. If a mathematical model as the one used in this paper is fitted accurately with experimental data of microorganisms and substrates, then biologists can use trustfully this model to predict future main behaviors of their anaerobic digesters. Also, they can explain and understand some observed phenomena using results of our analysis. For instance, depending on the value of the parameter $m$ of the growth kinetics of $X_1$ on $S$, biologists can interpret why they find different concentrations for bacteria and substrates at steady state, when doing the experiment starting from different initial concentrations. This is exactly the multistability which is predicted by the model when the value of the bifurcation parameter $m$ varies. On the other hand, experimenters can act on the operating parameters $S_{1in}$ and $S_{2in}$ and the bifurcation parameter $m$ in order to force the behavior of
the biological process towards a desired steady state.

In light of these results, our main perspective consists of establishing of a complete operating diagram of the considered model with respect of operating parameters which are $D, S_{1in}, S_{2in}$ and, especially the maximum growth rate $m$ of the first bacteria $X_1$ on the new variable $S$. In other terms, we wish to explore the different asymptotic behaviors of the system in 2 dimensional planes where one of the plane coordinates is the maximum growth rate $m$. Such operating diagrams if well established and discussed, can be really useful to interpret experimental results and, to help biologists to best choose values of operating parameters for controlling their experiments.

## Acknowledgments

The authors thank the Euro-Mediterranean research network Treasure (https://www6.inrae.fr/treasure) who partially financed this research. Part of the work was completed during the mission of the first author in Narbonne. This mission was publicly funded through ANR (the French National Research Agency) under the "Investissements d'avenir" program with the reference ANR-16-IDEX-0006. The first author would like to thank: Direction Générale de la Recherche Scientifique et du Développement Technologique (DG RSDT), Algeria, for support. We thank the reviewers for their valuable comments, which helped to improve the quality of the paper.

## Conflict of interest

The authors declare that they have no competing interests.

## References

1. M. J. Wade, J. Harmand, B. Benyahia, T. Bouchez, S. Chaillou, B. Cloez, et al., Perspectives in mathematical modelling for microbial ecology, *Ecol. Model.*, **321** (2016), 64–74.

2. S. Di, A. Yang, Analysis of productivity and stability of synthetic microbial communities, *J. R. Soc. Interface*, **16** (2019), 20180859.

3. T. Sari, M. El Hajji, J. Harmand, The mathematical analysis of a syntrophic relationship between two microbial species in a chemostat, *Math. Biosci. Eng.*, **9** (2012), 627–645.

4. M. El Hajji, F. Mazenc, J. Harmand, A mathematical study of a syntrophic relationship of a model of anaerobic digestion process, *Math. Biosci. Eng.*, **7** (2010), 641–656.

5. D. J. Batstone, J. Keller, I. Angelidaki, S. V. Kalyuzhnyi, S. G. Pavlostathis, A. Rozzi, et al., The IWA Anaerobic Digestion Model No 1 (ADM1), *Water Sci. Technol.*, **45** (2002), 65–73.

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

7. T. Sari, J. Harmand, A model of a syntrophic relationship between two microbial species in a chemostat including maintenance, *Math. Biosci.*, **275** (2016), 1–9.

8. T. Sari, M. J. Wade, Generalised approach to modelling a three-tiered microbial food-web, *Math. Biosci.*, **291** (2017), 21–37.
9. M. J. Wade, Not Just Numbers: Mathematical Modelling and Its Contribution to Anaerobic Digestion Processes, *Processes*, **8** (2020), 888.

10. D. J. Batstone, D. Puyol, X. Flores-Alsina, R. Jorge, Mathematical modelling of anaerobic digestion processes: applications and future needs, *Rev. Environ. Sci. Biotechnol.*, **14** (2015), 595–613.

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

12. O. Bernard, G. Bastin, On the estimation of the pseudo-stoichiometric matrix for macroscopic mass balance modelling of biotechnological processes, *Math. Biosci.*, **193** (2005), 51–77.

13. J. A. Arzate, M. Kirstein, F. C. Ertem, E. Kielhorn, H. R. Malule, P. Neubauer, et al., Anaerobic Digestion Model (AM2) for the Description of Biogas Processes at Dynamic Feedstock Loading Rates, *Chem. Ing. Tech.*, **89** (2017), 686–695.

14. V. Alcaraz-González, J. Harmand, A. Rapaport, J. P. Steyer, V. González-Alvarez, C. Pelayo-Ortiz, Software sensors for highly uncertain WWTPs: a new approach based on interval observers, *Water Res.*, **36** (2002), 2515–2524.

15. V. Alcaraz-González, J. P. Steyer, J. Harmand, A. Rapaport, V. González-Alvarez, C. Pelayo-Ortiz, Application of a Robust Interval Observer to an Anaerobic Digestion Process, *Dev. Chem. Eng. Miner. Process.*, **13** (2008), 267–278.

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

17. T. Sari, B. Benyahia, The operating diagram for a two-step anaerobic digestion model, preprint, https://hal.inrae.fr/hal-02557464, (2020).

18. M. Sbarciog, M. Loccfufer, E. Noldus, Determination of appropriate operating strategies for anaerobic digestion systems, *Biochem. Eng. J.*, **51** (2010), 180–188.

19. J. Hess, O. Bernard, Design and study of a risk management criterion for an unstable anaerobic wastewater treatment process, *J. Process Control*, **18** (2008), 71–79.

20. N. Dimitrova, M. Krastanov, Nonlinear stabilizing control of an uncertain bioprocess model, *Int. J. Appl. Math. Comput. Sci.*, **19** (2009), 441–454.

21. A. Donoso-Bravo, P. Gajardo, M. Sebbah, D. Vicencio, Comparison of performance in an anaerobic digestion process: one-reactor vs two-reactor configurations, *Math. Biosci. Eng.*, **16** (2019), 2447–2465.

22. M. Weedermann, Analysis of a model for the effects of an external toxin on anaerobic digestion, *Math. Biosci. Eng.*, **9** (2012), 445–459.

23. O. H. Abdelkader, A. H. Abdelkader, Modeling Anaerobic Digestion Using Stochastic Approaches, *Trends in Biomathematics: Mathematical Modeling for Health, Harvesting, and Population Dynamics*. Springer International Publishing, **99** (2019), 373–396.

24. S. Hassam, E. Ficara, A. Leva, J. Harmand, A generic and systematic procedure to derive a simplified model from the Anaerobic Digestion Model No. 1 (ADM1), *Biochem. Eng. J.*, **99** (2015), 193–203.
25. B. Benyahia, T. Sari, B. Cherki, J. Harmand, Anaerobic membrane bioreactor modeling in the presence of Soluble Microbial Products (SMP) - the Anaerobic Model AM2b, *Chem. Eng. J.*, **228** (2013), 1011–1022.

26. A. Charfi, N. Thongmak, B. Benyahia, M. Aslam, J. Harmand, N. Ben Amar, et al., A modelling approach to study the fouling of an anaerobic membrane bioreactor for industrial wastewater treatment, *Bioresour. Technol.*, **245** (2017), 207–215.

27. B. Benyahia, A. Charfi, M. Heran, B. Cherki, N. Kalboussi, J. Harmand, Coupling a simple and generic membrane fouling model with biological dynamics: application to the modeling of an Anaerobic Membrane BioReactor (AnMBR), preprint, https://hal.inrae.fr/hal-02558408, (2020).

28. A. Charfi, E. Parka, M. Aslam, J. Kim, Particle-sparged anaerobic membrane bioreactor with fluidized polyethylene terephthalate beads for domestic wastewater treatment: Modelling approach and fouling control, *Bioresour. Technol.*, **258** (2018), 263–269.

29. Á. Robles, M. V. Ruano, A. Charfi, G. Lesage, M. Heran, J. Harmand, et al., A review on anaerobic membrane bioreactors (AnMBRs) focused on modelling and control aspects, *Bioresour. Technol.*, **270** (2018), 612–626.

30. F. Campillo, M. Chebbi, S. Toumi, Stochastic modeling for biotechnologies: Anaerobic model AM2b, *Math. Biol. Environ.*, **28** (2019), 13–23.

31. O. Hadj Abdelkader, *Méthodes stochastiques pour la modélisation et l'identification des bioprocédés de digestion anaérobie*, Ph.D Thesis, Tlemcen University, 2019.

32. G. Bastin, D. Dochain, *On-line Estimation and Adaptive Control of Bioreactors*, Elsevier, 1990.

33. B. Benyahia, *Modélisation et Observation des Bioprocédés à Membranes: Application à la digestion anaérobie*, Ph.D Thesis, Tlemcen University, 2012.

## A. Proofs

### A.1. *Proof of propositions 2.3*

#### A.1.1. Positivity

We have the following solutions for equations (2.6) and (2.8):

$$X_1(t) = X_1(0)e^{\int_0^t [\mu_1(S_1(\tau)) + \mu(S(\tau)) - D_0 - D_1]d\tau}, \\ X_2(t) = X_2(0)e^{\int_0^t [\mu_2(S_2(\tau)) - D_0 - D_1]d\tau}.$$

Thus, we have:

* $X_1(0) = 0 \Rightarrow X_1(t) = 0$ and $X_2(0) = 0 \Rightarrow X_2(t) = 0$,

* $X_1(0) > 0 \Rightarrow X_1(t) > 0$ and $X_2(0) > 0 \Rightarrow X_2(t) > 0$.

To prove the positivity of $S_1$, $S_2$ and $S$, we set these variables equal to zero in (2.5), (2.7) and (2.9) respectively and, we verify if their derivatives are positives:

* $\dot{S}_1 = DS_{1in} > 0$,
• $\dot{S}_2 = DS_{2\text{in}} > 0$,

• $\dot{S} = D_0X_1 + D_0X_2 > 0$, if $X_1 > 0$ and $X_2 > 0$.

Notice that $\dot{S}_1$, $\dot{S}_2$ and $\dot{S}$ are positives. All vector fields at bounds are inside directed. Consequently, the variables $S_1$, $S_2$ and $S$ remain positives for positive initial conditions.

### A.1.2. Boundedness

Let us define the quantity:

$$ \Sigma = S_1 + S_2 + X_1 + X_2 + S. $$

The dynamic of $\Sigma$ is written as follows:

$$ \dot{\Sigma} = D(S_{1\text{in}} + S_{2\text{in}}) - D(S_1 + S_2) - D_1(X_1 + X_2) - MS \\ -(k_1 - 1 - b_3 - k_2)\mu_1(S_1)X_1 - (k_3 - 1 - b_4)\mu_2(S_2)X_2 - (b_1 - 1 - b_2)\mu(S)X_1. $$

We have three dilution rates: $D$, $D_1$ and $M$ which is a combination of $D$ and $D_1$. Let us set $D_{\min} = \min(D, D_1)$, which allows to write:

$$ \dot{\Sigma} \le D(S_{1\text{in}} + S_{2\text{in}}) - D_{\min}\Sigma - (k_1 - 1 - b_3 - k_2)\mu_1 X_1 - (k_3 - 1 - b_4)\mu_2 X_2 - (b_1 - 1 - b_2)\mu X_1. \quad (A.1) $$

By using inequalities (2.10), (2.11) and (2.12), we can write:

$$ \dot{\Sigma} \le D(S_{1\text{in}} + S_{2\text{in}}) - D_{\min}\Sigma. $$

Since the solution of the equation $\dot{\Sigma}_0 = D(S_{1\text{in}} + S_{2\text{in}}) - D_{\min}\Sigma_0$ is:

$$ \Sigma_0(t) = \frac{D(S_{1\text{in}} + S_{2\text{in}})}{D_{\min}} + C e^{-D_{\min}t}, \quad \text{with } C \text{ is constant,} $$

then, we have $\Sigma(t) \le \Sigma_0(t)$, i.e.:

$$ \Sigma(t) \le \frac{D(S_{1\text{in}} + S_{2\text{in}})}{D_{\min}} + C e^{-D_{\min}t} \implies \lim_{t \to +\infty} \Sigma(t) \le \frac{D(S_{1\text{in}} + S_{2\text{in}})}{D_{\min}}. $$

Consequently, the variables of system (2.5–2.9) remain bounded.

### A.2. Proof of Lemma 3.1, 3.2 and 3.3

The equilibrium points are solutions of the nonlinear algebraic system obtained from (2.5–2.9) by setting the right-hand sides equal to zero.

#### A.2.1. Proof of Lemma 3.1 ($X_1^* = 0$)

From (3.2), we can have a trivial solution $X_1^* = 0$, which if replaced in (3.1), then we obtain $S_1^* = S_{1\text{in}}$. From Eq (3.4), we can have two cases:

• A trivial solution $X_2^* = 0$: which if replaced in (3.3) and (3.5), then we have $S_2^* = S_{2\text{in}}$ and $S^* = 0$ respectively. This is the equilibrium $E_0^0$.

• A nontrivial solution $S_2^* = \mu_2^{-1}(D_0 + D_1) = S_2^{i*}$, $i = 1, 2$: which if replaced in (3.3) and (3.5), then we deduce corresponding values of $X_2^{i*}$ and $S^{i*}$ respectively. These are equilibria $E_1^i$, $i = 1, 2$.
### A.2.2. Proof of Lemma 3.2 ($X_1^* > 0$ and $X_2^* = 0$)

Let $(S_1^*, X_1^*, S_2^*, X_2^*, S^*)$ a solution of system (3.1)–(3.5).

Since $X_1^* > 0$, from (3.2) we have $\mu_1(S_1^*) + \mu(S^*) = D_0 + D_1$, i.e.:

$$S_1^* = \mu_1^{-1}(D_0 + D_1 - \mu(S^*)) = F(S^*).$$

From (3.1), we deduce:

$$X_1^* = D \frac{S_{1\text{in}} - S_1^*}{k_1 \mu_1(S_1^*)},$$

which is positive and bounded if $S_1^* < S_{1\text{in}}$. By replacing $X_2^* = 0$ and $X_1^*$ in (3.3) we obtain:

$$S_2^* = S_{2\text{in}} + [k_2\mu_1(S_1^*) + b_2\mu(S^*)]\frac{S_{1\text{in}} - S_1^*}{k_1\mu_1(S_1^*)}.$$

Finally, if we replace $X_2^*, X_1^*$ and $\mu(S^*) = D_0 + D_1 - \mu_1(S_1^*)$ in (3.5), then we have after simplification:

$$S^* = (S_{1\text{in}} - S_1^*) \left( B_1 + \frac{B_2}{\mu_1(S_1^*)} \right) = G(S_1^*),$$

with:

$$B_1 = \frac{b_3 + b_1}{k_1 B}, \quad B_2 = \frac{D_0 - b_1(D_0 + D_1)}{k_1 B}, \quad B = \beta + (1 - \beta) \frac{D_1}{D}.$$

Then $S_1^*$ and $S^*$ are solutions of the system of Eqs (3.12).

### A.2.3. Proof of Lemma 3.3 ($X_1^* > 0$ and $X_2^* > 0$)

Let $(S_1^*, X_1^*, S_2^*, X_2^*, S^*)$ a solution of system (3.1)–(3.5).

Since $X_2^* > 0$, from (3.4) we have the nontrivial solution:

$$S_2^* = \mu_2^{-1}(D_0 + D_1) = S_2^{i*}, \quad i = 1, 2.$$

Since $X_1^* > 0$, from (3.2) we have:

$$S_1^* = \mu_1^{-1}(D_0 + D_1 - \mu(S^*)) = F(S^*),$$

and thus, if $0 < S_1^* < S_{1\text{in}}$ then, we deduce from (3.1):

$$X_1^* = \frac{D[S_{1\text{in}} - S_1^*]}{k_1 \mu_1(S_1^*)}.$$

By replacing $X_1^*$ in (3.3), we obtain:

$$X_2^{i*} = \beta_i + \frac{D}{k_3(D_0 + D_1)} \frac{k_2\mu_1(S_1^*) + b_2\mu(S^*)}{k_1\mu_1(S_1^*)} (S_{1\text{in}} - S_1^*),$$

with: $\beta_i = \frac{D}{k_3(D_0+D_1)} (S_{2\text{in}} - S_2^{i*})$, $i = 1, 2$

Finally, from (3.5) we have after simplification:

$$S^* = \alpha_i + (S_{1\text{in}} - S_1^*) \left( C_1 - \frac{C_2}{\mu_1(S_1^*)} \right) = H_i(S_1^*), \quad i = 1, 2,$$
with:

$$ \alpha_i = \frac{A}{B}(S_{2\text{in}} - S_2^{i*}), \quad C_1 = B_1 + \frac{A(k_2 - b_2)}{k_1\beta}, \quad C_2 = B_2 - \frac{Ab_2}{k_1\beta}. $$

$$ A = \frac{b_4(D_0 + D_1) + D_0}{k_3(D_0 + D_1)}, \quad B = \beta + (1 - \beta)\frac{D_1}{D}. $$

Then $S_1^*$ and $S^*$ are solutions of the system of Eqs (3.17).

The function $H_i(S_1^*)$ can be written as:

$$ H_i(S_1) = G(S_1) + \frac{A}{B}[S_{2\text{in}} - S_2^{i*} + (k_2\mu_1(S_1) + b_2\mu(S))\frac{S_{1\text{in}} - S_1}{k_1\mu_1(S_1)}]. $$

The condition for which $X_2^{i*} > 0$ is:

$$ S_{2\text{in}} - S_2^{i*} + (k_2\mu_1(S_1^*) + b_2\mu(S^*))\frac{S_{1\text{in}} - S_1^*}{k_1\mu_1(S_1^*)} > 0, $$

which is equivalent to: $H_i(S_1^*) > G(S_1^*)$, $i = 1, 2$. (condition of lemma 3.3).

### A.3. Proof of Lemma 4.1

Using (2.13), we have $C_1 > B_1$ and $C_2 > B_2 > 0$, consequently: $D_H < D_G < D_0 + D_1$.

Using the fact that $\mu_1$ is increasing, we deduce that:

$$ \lambda_H < \lambda_G < \lambda_1 \quad (A.2) $$

with: $\lambda_1 = \mu_1^{-1}(D_0 + D_1)$, $\lambda_G = \mu_1^{-1}(D_G)$ and $\lambda_H = \mu_1^{-1}(D_H)$, where $D_G = B_2/B_1$ and $D_H = C_2/C_1$.

### A.4. Proof of the proposition 4.2

The function $G(S_1)$ given by (3.10) is positive between $S_{1\text{in}} > S_1 > \lambda_G$ and, solutions of the system (3.12) must satisfy $S_{1\text{in}} > S_1^*$, then $S_{1\text{in}} > \lambda_G$ that is to say $\mu_1(S_{1\text{in}}) > D_G$.

### A.5. Proof of the proposition 4.3

Functions $H_i(S_1)$, $i = 1, 2$ given by (3.15) are translations of the function $H(S_1)$ with quantities $\alpha_i$ given by (3.7). The sign of this later indicates if the equilibrium $E_2^i$, $i = 1, 2$ does exist or not (see Figure 4).

### A.6. Proof of Theorem 6.1

The study of the local stability of trivial equilibria follows easily from the study of the Jacobian matrix of system (2.5–2.9), which has a block-diagonal structure:

$$ J = \begin{bmatrix} A & 0 & 0 \\ C & B & 0 \\ M_1 & M_2 & -M \end{bmatrix}, $$
Hence, the eigenvalues of J are the eigenvalues of A, the eigenvalues of B and $-M$ (which is always negative).

For the $X_1$ and $X_2$ washout equilibrium $E_0^0 = [S_{1\text{in}}, 0, S_{2\text{in}}, 0, 0]$:
one has:

$$A = \begin{bmatrix} -D & -k_1\mu_1(S_{1\text{in}}) \\ 0 & \mu_1(S_{1\text{in}}) - D_0 - D_1 \end{bmatrix},$$

$$B = \begin{bmatrix} -D & -k_3\mu_2(S_{2\text{in}}) \\ 0 & \mu_2(S_{2\text{in}}) - D_0 - D_1 \end{bmatrix}.$$

Conditions of stability are: $tr(A) < 0$, $tr(B) < 0$, $det(A) > 0$ and $det(B) > 0$. Thus, $E_0^0$ is stable if and only if:

$$\mu_1(S_{1\text{in}}) < D_0 + D_1 \text{ and, } \mu_2(S_{2\text{in}}) < D_0 + D_1$$

For the $X_1$ washout equilibria $E_1^i = [S_{1\text{in}}, 0, S_2^{i*}, X_2^{i*}, S^{i*}]$, $i = 1, 2$:
one has:

$$A = \begin{bmatrix} -D & -k_1\mu_1(S_{1\text{in}}) \\ 0 & \mu_1(S_{1\text{in}}) + \mu(S^{i*}) - D_0 - D_1 \end{bmatrix},$$

$$B = \begin{bmatrix} -D - k_3\mu_2(S_2^{i*})X_2^{i*} & -k_3\mu_2(S_2^{i*}) \\ \mu_2(S_2^{i*})X_2^{i*} & 0 \end{bmatrix}.$$

One can easily deduce that if $E_1^2$ exists, then it is unstable because $det(B) < 0$, since $\mu_2'(S_2^{2*}) < 0$.

On the other hand, stability of $E_1^1$ depends on $\mu(S^{1*})$. Indeed, $E_1^1$ is stable if and only if:

$$\mu_1(S_{1\text{in}}) + \mu(S^{1*}) < D_0 + D_1.$$

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

AIMS Press

©2020 the Author(s), licensee AIMS Press. This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0)