SIR model-based verification of effect of COVID-19 Contact-Confirming Application (COCOA) on reducing infectors in Japan

1.   Introduction

Mathematical modeling of bioprocesses is a powerful tool to (ⅰ) explain observed phenomena, (ⅱ) understand some mechanisms of the system and predict its evolution, (ⅲ) better control the process operations and (ⅳ) 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 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 $CO_2$, while in the second step (methanogenesis), the methanogenic population consumes VFA and produces methane and $CO_2$. 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,19,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 bahavior 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 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:

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:

where $S_{1{\rm in}}$ and $S_{2{\rm 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:

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

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 $\mathcal{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 $\mathcal{C}^1$ and satisfies the following properties:

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

● $\mu_2(S_2)$ has a maximum $\mu_2\left(S_2^M\right) > 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 hypothesises (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\Bigl(1 - e^{\frac{-\xi}{K}}\Bigr)$, 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\Bigl(e^{-\alpha_1\xi} - e^{-\alpha_2\xi}\Bigr)$ (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:

We use the following notations:

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:

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{\rm in}}, 0, S_{2{\rm in}}, 0, 0)$, which always exists,

● the washout equilibrium of $X_1$ but not of $X_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:

The equilibrium $E_{1}^{i}$ exists if and only if:

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{\rm in}}$ (see Figure 1). Often, this is not possible in the biological realism, unless we consider a third acetogenic microorganisms which produce $S_{2{\rm in}}$ from external organic matter or, if we carry out a bench-scale study, by introducing $S_{2{\rm 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:

where:

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:

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

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

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:

where:

Lemma 3.3. Let $E_2^i = (S_1^*, X_1^*, S_2^{i*}, X_2^{i*}, 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:

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

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

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:

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$.

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

We have two cases (see Figure 3):

● $D_G > m_1$, where $g(S_1)$ is always negative for $S_1 < S_{1{\rm in}}$ (Figure 3, left) and values of $S_1 > S_{1{\rm 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).

Proposition 4.2. A necessary condition for (3.12) to have positive solutions is $\lambda_G < S_{1{\rm in}}$ that is to say $\mu_1(S_{1{\rm 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{\rm in}}$ and $\lambda_H$, the root of:

Two cases can be distinguished:

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

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

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

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

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

● $H(S_1) > G(S_1)$ for all $\lambda_H < S_1 < S_{1{\rm 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{\rm in}}$.

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

5.   Necessary conditions of existence of equilibria

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

● $S_{1{\rm in}} < \lambda_H < \lambda_G < \lambda_1$

● $\lambda_H < S_{1{\rm in}} < \lambda_G < \lambda_1$

● $\lambda_H < \lambda_G < S_{1{\rm in}} < \lambda_1$

● $\lambda_H < \lambda_G < \lambda_1 < S_{1{\rm in}}$

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 doted 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.

Case: $S_{1{\rm in}} < \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.

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_{1{\rm in}}$, then we obtain equilibria $E_1^i = (S_{1{\rm in}}, 0, S_{2}^{i*}, X_{2}^{i*}, S^{i*})$ where $S^{i*} = \alpha_i = H_i(S_{1{\rm in}})$, $i = 1, 2$ as it can be seen on Figure 5, center, for $E_1^1$ and on Figure 5, right, for $E_1^i$, $i = 1, 2$.

Case : $\lambda_H < S_{1{\rm in}} < \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$.

Case: $\lambda_H < \lambda_G < S_{1{\rm in}} < \lambda_1$

We represents 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_{1{\rm in}}$

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$).

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 equilirium $E_0^0$ exists always and it is stable if and only if:

2. The equilibrium $E_1^2$ exist if and only if $S_{2{\rm in}} > S_{2}^{2*}$ and it is always unstable.

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

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

The condition (6.2) my be graphically explained on the Figure 9. The graph $f(S_1, S) = \mu_1\left(S_1\right) + \mu\left(S\right) - 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).

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).

\normalsize Which can be symbolized as follows:

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:

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

with:

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).

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.

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$.

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.

In generic cases 2 and 3, we have equilibria bifurcation when $m$ varies. For some values $c_{i}$, $i = 1, .., 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.

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 \cap H_1$, $F \cap H_2$ and $F \cap 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} \leq m \leq 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_{1{\rm in}} \Rightarrow \mu_1(\lambda_1) < \mu_1(S_{1{\rm in}})$, that is to say $D_0+D_1 < \mu_1(S_{1{\rm in}}) \Rightarrow E_0^0$ is unstable according to theorem 6.1.1

● $(S_{1{\rm in}}, 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_{1{\rm in}}) < 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{\rm in}}) < D_0+D_1$ for the case 2, thus the equilibrium $E_0^0$ is Stable.

$\mu_2(S_{2{\rm 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 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^*$.

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_{1{\rm in}}$ and $S_{2{\rm in}}$.

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_{1{\rm in}}$ and $S_{2{\rm in}}$ 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_{1{\rm in}}$, $S_{2{\rm in}}$ 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.

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):

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_{1{\rm in}} > 0,$

● $\dot{S}_2 = DS_{2{\rm 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:

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

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:

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

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

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

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{\rm 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{\rm 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.:

From (3.1), we deduce:

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

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:

with:

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:

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

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

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

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

Finally, from (3.5) we have after simplification:

with:

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

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

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

which is equivalent to: $H_i(S_1^*) > G(S_1^*), \quad 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:

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{\rm in}} > S_1 > \lambda_G$ and, solutions of the system (3.12) must satisfy $S_{1{\rm in}} > S_1^*$, then $S_{1{\rm in}} > \lambda_G$ that is to say $\mu_1(S_{1{\rm 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:

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{\rm in}}, 0, S_{2{\rm in}}, 0, 0]$:

one has:

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:

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

one has:

One can easily deduce that if $E_1^2$ exists, then it is unstable because $det(B) < 0$, since $\mu_2'\left(S_2^{2*}\right) < 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: