Study of performance criteria of serial configuration of two chemostats

1.   Introduction

The chemostat device was invented concomitantly by Monod ^[1] and Novick & Szilard ^[2] in 1950. Widely used as a biochemical laboratory-pilot, it consists essentially in a continuously-fed bioreactor characterized by the equality of the input and the output flow rates. It is designed as a vessel in which different microorganisms grow, also called continuous culture of microorganisms. Its importance for the continuous culture of microorganisms has been reported in several books and publications, among them ^[3,4,5,6,7]. In other words, the classical model of the chemostat consists of a perfect mixed media at a constant temperature, a constant pH, a filtered feed and a unique flow rate. Although this model is used for industrial applications with continuously-fed bioreactors such as waste-water treatment, see for instance ^[8], in physical reality, industrial applications which use large bioreactors hardly satisfy the assumption of the perfect mixed media. Several mathematical representations of the spatial heterogeneity have been studied in the literature with partial differential equations, see for instance ^[9,10]. However, discrete spatial representations, such as the gradostat model ^[6,11], are also a way to represent spatial heterogeneity ^[12,13,14]. Serial configurations, as a simple gradostat, have received a great interest in the literature in view of optimizing bioprocesses. Indeed, it has been shown that having two tanks (or more) in series (each of them being assumed to be perfectly mixed) can produce the same substrate conversion than a single vessel, but with a significant lower total volume, and thus a lower residence time. Serial configurations have been also studied in view of ecological insight, see for instance ^[15,16]. In this paper, we propose to revisit the serial configuration of two chemostats in series with a constant total volume $V$, as shown in Figure 1. We focus on the analysis of the performance at steady-state for different criteria with the aim of drawing comparisons with the single chemostat. Notice that these different criteria of comparison are known in the literature, see for instance ^[3]. However, to our knowledge, a complete and deep analysis of all possible configurations for a general class of growth functions and the various criteria is missing in the literature, which is the aim of the present work.

It is well known ^[4,6] that, for the simple chemostat, the output concentration at steady steady $S^{out}$ is independent of the input concentration $S^{in}$, provided that there is no washout, see also (2.5). This property is no longer satisfied when there is a spatial structure, see for instance ^[15] and the references therein. Since $S^{out}$ measures the performances of the chemostat to convert the substrate $S$, our purpose is to distinguish which configuration guarantees the minimal output substrate concentration at steady state. Actually, reducing the output substrate concentration is one of the biological objectives in waste-water treatments and this minimizing problem is well known in the literature. The novelty of our work is that $S^{out}$ is considered as a function which depends on the three operating parameters: the input substrate concentration, the dilution rate and the volume of each chemostat. In fact, what has already been treated, see for instance ^[17,18,19], corresponds to the case where the input substrate concentration $S^{in}$ is fixed and the total volume $V$ can be chosen. Thus, we give conditions which involve the input substrate concentration $S^{in}$ and ensure the optimal way to slice the two serial reactors volume. These conditions can ensure a lower output substrate concentration.

Our study is somehow a generalization of the main results presented in ^[16]. The conditions that we found are necessary and sufficient to reduce the output substrate concentration in contrast of the result in ^[16] where the given conditions are only sufficient. In addition, the originality of this article consists in comparing both configurations according to two other performance indexes which are the productivity of the biomass and the biogas flow rate. The biogas flow rate represents the quantity of natural gas per unit of time produced by the decomposition of organic matter in absence of oxygen and the productivity of the biomass represents the amount of biomass per unit of time produced by the decomposition of organic matter. The productivity of the biomass of several configurations including the serial device of two interconnected chemostats has been graphically and numerically analyzed in ^[12,20]. However, these two criteria have not yet been deeply mathematically analyzed. The global analysis shows that the different performance criteria involve the same performance threshold. This threshold is explicitly defined by a function which depends on the dilution rate $D$. It defines the set of the values of $S^{in}$ and $D$ that allow or not a better performance of the serial configuration with two chemostats. Several numerical applications are given to illustrate all the results of the study.

This paper is organized as follows. Section 2 presents the model. Subsequently, the main part of the paper constituting section 3 is dedicated to the study of the equilibria and the performance analysis of the configuration. Indeed, the output substrate concentration, the productivity of the biomass and the biogas flow rate are respectively treated in sections 3.1, 3.2 and 3.3. Next, the operating diagram of the model is depicted in section 4. Afterwards, several numerical simulations illustrating the results of our analysis and using some specific growth functions are represented in section 5. Finally, section 6 contains a global conclusion. Most of the proofs corresponding of the theorems and propositions stated along the paper are proved in Appendixes A, B, C and D. Firstly, Appendix A contains the proof related to the existence and the stability of steady states. Secondly, Appendixes B and C contain respectively the proofs related to the output substrate concentration, the productivity of the biomass and the biogas flow rate. Finally, Appendix D contains proofs related to some of technical results of the paper.

2.   Mathematical model

If $S$ and $X$ denote respectively the substrate and the biomass concentration in a single chemostat of volume $V$, the input flow rate $Q$ and the input concentration of substrate $S^{in}$, their time evolution are modeled by the following system of ordinary differential equations:

where $Y$ is the yield conversion of substrate into biomass, $f(\cdot)$ the specific growth rate of the microorganisms that is assumed null at $S = 0$ and to be increasing for $S > 0$, and $D = Q/V$ is the dilution rate. Without loss of generality, one can assume $Y = 1$ in Eq (2.1) by using the change of variable $x = X/Y$. System (2.1) becomes

The detailed mathematical analysis of the model (2.2) may be found in ^[4,6]. Let us recall classical results about the asymptotic behavior of (2.2). We define

As $f$ is increasing then the break-even concentration is defined by

When $S^{in} > \lambda(D)$ (or, equivalently, $f(S^{in}) > D$), any solution of (2.2) with $S(0)\geq 0$ and $x(0) > 0$ converges toward the positive steady state $E_1 = \left(\lambda(D), S^{in}-\lambda(D)\right)$. On the contrary, when $D\geq m$ or $S^{in}\leq\lambda(D)$ (or, equivalently, $f(S^{in})\leq D$), any solution of (2.2) with $S(0)\geq 0$ and $x(0)\geq 0$ converges toward the wash-out steady state $E_0 = (S^{in}, 0)$. Thus, the output concentration at steady state $S^{out}(S^{in}, D)$ is given by

We consider now the serial interconnected chemostats, where the volume $V$ is divided into two volumes, $rV$ and $(1-r)V$ with $r\in (0, 1)$, as shown in Figure 1, with $Q$ the flow rate and $S^{in}$ the input substrate concentration in the first chemostat. The mathematical model is given by the following equations:

The dilution rates of the first and second chemostats are $D/r$ and $D/(1-r)$, respectively, where $D$ is the dilution rate, defined by $D: = Q/V$, of the single chemostat of volume $V$ and flow rate $Q$. For the limiting cases $r = 0$ and $r = 1$, these equations are not valid. Indeed, the limiting cases correspond to the single chemostat model defined by (2.2).

In ^[16], the mathematical analysis of (2.6) was performed for a linear growth function $f(S) = aS$ ($a > 0$) and numerical simulations were given for a Monod growth function $f(S) = mS/(K+S)$. The results of ^[16] were extended to Monod growth function and for increasing and concave growth function in ^[21].

Remark 1. The main result in ^[16,21], see also ^[15], predicts that there exists a threshold ${S}^{in}_1$ such that for ${S}^{in}\leq {S}^{in}_1$, the output $S_r^{out}(S^{in}, D)$, which is the output density of the substrate at steady state, satisfies $S_r^{out}(S^{in}, D) > \lambda(D)$, for all $r\in(0, 1)$ and, if ${S}^{in} > {S}^{in}_1$, there exists a threshold $r_1\in(0, 1)$, such that $S_r^{out}(S^{in}, D) < \lambda(D)$ if and only if $r_{1} < r < 1$.

As it was noticed in (2.5), for a single chemostat, one has $S^{out}(S^{in}, D) = \lambda(D)$. Therefore, if $S^{in}\leq{S}^{in}_1$, the serial configuration is always less efficient than the single chemostat of the same total volume $V$. In contrast, for $S^{in} > {S}^{in}_1$ and $r$ large enough (i.e. $r > r_1$), the serial configuration is more efficient than the single chemostat.

In this paper, we extend this result to general increasing growth functions where the concavity of $f$ is not required and we provide explicit formulas for the thresholds ${S}^{in}_1(D)$ and $r_{1}(S^{in}, D)$. Hence, we consider a growth function satisfying only the following qualitative property:

Assumption 1. The function $f$ is $\mathcal{C}^1$, with $f(0) = 0$ and $f'(S) > 0$ for all $S > 0$.

The following result is classical in the mathematical theory of the chemostat and is left to the reader.

Lemma 1. The solutions $(S_1(t), x_1(t), S_2(t), x_2(t))$ of (2.6) with nonnegative initial conditions, exist for all $t\geq 0$, are nonnegative, bounded and $\lim_{t\to +\infty}\left(S_i(t)+x_i(t)\right) = S^{in}$ for $i = 1, 2$.

The existence and stability of steady states of (2.6) are given by the following result. We use the abbreviation LES for locally exponentially stable and GAS for globally asymptotically stable in the positive orthant.

Theorem 1. Assume that Assumption 1 is satisfied. The steady states of (2.6) are:

● The washout steady state $E_0 = (S^{in}, 0, S^{in}, 0)$ which always exists. It is GAS if and only if

It is LES if and only if: $D > \max\{r, 1-r\}f(S^{in})$.

● The steady state $E_1 = \left(S^{in}, 0, \overline{S}_2, S^{in}-\overline{S}_2\right)$ of washout in the first chemostat but not in the second one, where $\overline{S}_2$ is given by $\overline{S}_2 = \lambda\left({D}/({1-r})\right)$. This steady state exists if and only if $D < (1-r)f(S^{in})$. It is GAS if and only if

It is LES if and only if: $rf(S^{in}) < D < (1-r)f(S^{in})$.

● The steady state $E_{2} = (S_{1}^{*}, S^{in}-S_{1}^{*}, S_{2}^{*}, S^{in}-S_{2}^{*})$ of persistence of the species in both chemostats, where $S_{1}^{*}$ is given by $S^{*}_{1} = \lambda\left({D}/{r}\right)$ and $S^{*}_{2} = S_2^ {*}(S^{in}, D, r)$ is the unique solution of the equation

This steady state exists if and only if $D < rf(S^{in})$. It is GAS and LES whenever it exists.

Proof. The proof is given in Appendix A.

Remark 2. Transcritical bifurcations occur in the limit cases $D = rf(S^{in})$ and $D = (1-r)f(S^{in})$.

1. For $0 < r < 1/2$, we have a transcritical bifurcation of $E_0$ and $E_1$ when $D = (1-r)f(S^{in})$ and a transcritical bifurcation of $E_1$ and $E_2$ when $D = rf(S^{in})$.

2. For $1/2 < r < 1$, we have a transcritical bifurcation of $E_0$ and $E_2$ when $D = rf(S^{in})$ and a transcritical bifurcation of $E_0$ and $E_1$ when $D = (1-r)f(S^{in})$.

3. For $r = 1/2$ and $D = f(S^{in})/2$, we have transcritical bifurcations of $E_0$ and $E_1$, and $E_0$ and $E_2$, simultaneously.

Figure 2$(a)$ shows the functions $f$ and $h$, and the solution $S^{*}_{2} = S_2^ {*}(S^{in}, D, r)$ of the Eq (2.9), which is unique since $f$ is strictly increasing and the graph of $h$ is a hyperbola. If $S^{in, 1} > S^{in, 2}$, then $h_i(S_2) = \frac{D(S_1^*-S_2)}{(1-r)(S^{in, i}-S_2)}$, $i = 1, 2$, satisfies $h_2(S_2) > h_1(S_2)$, for all $S_2 \in (0, S_1^{*})$, as shown in Figure 2$(b)$. Therefore, we have the following result:

Proposition 1. Let $S^{in\mathit{\mbox{, }}1}$ and $S^{in\mathit{\mbox{, }}2}$ be two different input substrate concentrations. If $S^{in, 1} > S^{in, 2} > 0$ then for all $r \in \left(D/f(S^{in, 2}), 1\right)$ and $D > 0$, one has $S_2^ {*}(S^{in, 1}, D, r) < S_2^ {*}(S^{in, 2}, D, r)$.

Proof. The proof is given in Appendix B.1.

3.   The three performance criteria

In this section we give the expressions of the output substrate concentration at steady state, the productivity of the biomass and the biogas flow rate, for the serial configuration of two chemostats.

3.1.   Output substrate concentration

Let us consider the dependency of the output substrate concentration with respect to the dilution rate $D$ and the input concentration $S^{in}$. As stated in Theorem 2.1, for $0 < r < 1$, the output substrate concentration at steady state is given by the formulas:

Although $S_r^{out}(S^{in}, D)$ is defined by (3.1) only for $0 < r < 1$, we extend it, by continuity, for $r = 0$ and $r = 1$ by

The continuity follows from the facts that $\lim_{r\rightarrow 1}S_2^*(S^{in}, D, r) = \lambda(D)$ and the second case, where $S_r^{out}(S^{in}, D) = \lambda(D/(1-r))$, is possible only if $0\leq r\leq {1}/{2}$.

We have to compare $S_r^{out}(S^{in}, D)$, given by (3.1) and (3.2), with $S^{out}(S^{in}, D)$, given by (2.5). Let $r\in(0, 1)$ be fixed. Let $g_r:[0, rm) \to \mathbb{R}$, where $m$ is given by (2.3), be defined by

The following result asserts that the serial configuration of two chemostats of volumes $rV$ and $(1-r)V$ respectively, shown in Figure 1, is more efficient than the simple chemostat of volume $V$, if and only if $S^{in} > g_r(D)$.

Theorem 2. For any $r\in(0, 1)$, one has $S_r^{out}(S^{in}, D) < S^{out}(S^{in}, D)$ if and only if $S^{in} > g_r(D)$.

Proof. The proof is given in Appendix B.2

We need the following assumption, which is satisfied by any concave growth function, but also by Hill function, which is not concave, as it is shown in section 5.

Assumption 2. For every $D\in[0, m)$, the function $r\in(D/m, 1)\mapsto g_r(D)\in \mathbb{R}$ is strictly decreasing.

Let $g:[0, m) \mapsto \mathbb{R}$ be defined by

We have the following result:

Lemma 2. Assume that Assumptions 1 and 2 are satisfied. For all $(S^{in}, D)$ verifying the condition $S^{in} > g(D)$, there exists a unique $r_1 = r_1(S^{in}, D) \in$ $(0, 1)$ such that $S^{in} = g_{r_1}(D)$. One has $r > r_1(S^{in}, D)$ if and only if $S^{in} > g_r(D)$.

Proof. The proof is given in Appendix B.3.

We can state now our main result which compares $S^{out}_r(S^{in}, D)$ and $S^{out}(S^{in}, D)$.

Theorem 3. Assume that Assumptions 1 and 2 are satisfied.

● If $S^{in} \leq g(D)$ then for any $r\in(0, 1)$, $S^{out}_r(S^{in}, D) > S^{out}(S^{in}, D)$.

● If $S^{in} > g(D)$ then $S^{out}_r(S^{in}, D) < S^{out}(S^{in}, D)$ if and only if $r_1(S^{in}, D) < r < 1$ with $r_1(S^{in}, D)$ defined in Lemma 2.

The equality is fulfilled for $r = 0$, $r = r_1$ and $r = 1$.

Proof. The proof is given in Appendix B.4.

Lemma 2 and Theorem 3 give analytical expression for the thresholds $S^{in}_1$ and $r_1$ mentioned in Remark 1. Indeed, we have $S^{in}_1 = g(D)$ and $r_1$ depends on $D$ and $S^{in}$, and is given implicitly by equation $S^{in} = g_{r_1}(D)$. In section 5, we give explicit formulas for $r_1(S^{in}, D)$ in the cases of linear growth functions, see (5.1), or Monod growth functions, see (5.2). To have a better understanding of the role of the parameter $r$, we also analyze the function $r\mapsto S_r^{out}(S^{in}, D)$ when $S^{in}$ and $D$ are fixed. According to the conditions on $S^{in}$ and $D$, related to the global stability of the equilibria, several cases must be distinguished. The following result encompasses the whole possible cases.

Proposition 2. Let $D > 0$ and $S^{in} > 0$. We denote by $r_0$ the ratio $r_0 = D/f(S^{in})$.

1) If $S^{in}\leq\lambda(D)$ then for any $r\in[0, 1]$, one has $S^{out}_r(S^{in}, D) = S^{out}(S^{in}, D) = S^{in}$.

2) If $\lambda(D) < S^{in} < \lambda(2D)$ then one has $\frac{1}{2} < r_0 < 1$ and

3) If $\lambda(2D)\leq S^{in}$ then one has $0 < r_0\leq\frac{1}{2}$ and

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

For a deeper analysis, we consider the functions $D\mapsto S^{out}_r(S^{in}, D)$ and $D\mapsto S^{out}(S^{in}, D)$ where we fix the input substrate density $S^{in}$ and the parameter $r$. We add the following assumption, which is satisfied by concave growth functions and also by Hill functions as it is shown in section 5.

Assumption 3. For every $r\in(0, 1)$, the function $D\in(0, rm)\mapsto g_r(D)\in \mathbb{R}$ is strictly increasing.

We have the following result:

Proposition 3. Assume that Assumptions 1 and 3 are satisfied. For any $r\in(0, 1)$ and $S^{in} > 0$, there exists a critical value $D_r = D_r(S^{in})$, which is the unique solution of the implicit equation $S^{in} = g_r(D)$, such that the serial configuration of two interconnected chemostats is more efficient than a simple chemostat if and only if $0 < D < D_r(S^{in})$. That is to say, for any $0 < D < D_r(S^{in})$, one has $S^{out}_r(S^{in}, D) < S^{out}(S^{in}, D)$.

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

The result of Proposition 3 is illustrated by Figure 3. In this Figure, the critical value $D_r = D_r(S^{in})$ is depicted for various values of $r$ and $S^{in}$, illustrating then Proposition 1 which asserts that, for a fixed dilution rate $D$, the output substrate concentration decreases when $S^{in}$ increases.

The following Lemmas 3 and 4 provide sufficient conditions for Assumption 2 and 3 to be satisfied. These conditions are useful for the applications given in section 5. For this purpose, we consider the function $\gamma$ defined by

which consists simply in considering $g_r(D)$, given by (3.3), as a function of both variables $r$ and $D$. If $\frac{\partial \gamma}{\partial r}(r, D) < 0$ for all $(r, D)\in\mathrm{dom}(\gamma)$, then Assumption 2 is satisfied. Similarly, if $\frac{\partial \gamma}{\partial D}(r, D) > 0$ for all $(r, D)\in\mathrm{dom}(\gamma)$, then Assumption 3 is satisfied. The following lemmas give equivalent conditions, and also sufficient conditions, for partial derivatives of $\gamma$ to have their signs as indicated above.

Lemma 3. For $D\in(0, m)$, let $l_D$ be defined on $(D/m, 1]$ by $l_D(r) = \lambda(D/r)$. The following conditions are equivalent

1. For all $(r, D)\in\mathrm{dom}(\gamma)$, $\frac{\partial \gamma}{\partial r}(r, D) < 0$.

2. For all $D\in(0, m)$ and $r\in(D/m, 1)$, $l_D(1) > l_D(r)+(1-r)l_D'(r)$.

If $l_D$ is strictly convex on $(D/m, 1]$, then the condition 2 is satisfied. If, in addition, $f$ is twice derivable, then $l_D$ is twice derivable and the following conditions are equivalent

3. For all $D\in(0, m)$ and $r\in(D/m, 1]$, $l_D''(r) > 0$.

4. For all $S > 0$, $f(S)f''(S) < 2\left(f'(S)\right)^2$.

Proof. The proof is given in Appendix D.1.

Lemma 4. The following conditions are equivalent

1. For all $(r, D)\in\mathrm{dom}(\gamma)$, $\frac{\partial \gamma}{\partial D}(r, D) > 0$.

2. For all $(r, D)\in\mathrm{dom}(\gamma)$, $f'\left(\lambda({D}/{r})\right) < f'\left(\lambda(D)\right)/r^2$.

If $f'$ is decreasing, then the condition 2 is satisfied.

Proof. The proof is given in Appendix D.2.

Remark 3. If the growth function is twice derivable and satisfies $f''(S)\leq 0$ for all $S > 0$, then the condition 4 in Lemma 3 and the condition 2 in Lemma 4 are satisfied. Thus, Assumptions 2 and 3 are satisfied. Therefore, our results apply for concave growth functions. The previous Lemmas allow to consider a non-concave growth function such as the Hill function as shown in section 5.3.

3.2.   Biomass productivity

Let us consider the dependency of the productivity of the biomass with respect to the dilution rate $D$ and the input concentration $S^{in}$. Recall that for a simple chemostat the output biomass at steady state is given by $x^{out} = S^{in}-S^{out}$. Thus, the productivity of a single chemostat is defined by

Let $D^{opt}(S^{in})$ be the dilution rate which maximizes $P(S^{in}, D)$ i.e.

Assumption 4. The dilution rate $D^{opt}(S^{in})$ defined by (3.9) is unique.

Proposition 4. The dilution rate $D^{opt}(S^{in})$ defined by (3.9) is the solution of equation $S^{in} = g(D)$ where $g$ is defined by (3.4).

Proof. The proof is given in Appendix C.1.

The productivity of the two serial interconnected chemostats at steady-state is defined as

with $x^{out}_r = S^{in}-S^{out}_r$. Using the definitions (3.1) and (3.2) of $S^{out}_r(S^{in}, D)$, for $r\in (0, 1)$ we have

and $P_r(S^{in}, D) = P(S^{in}, D)$, when $r = 0$ and $r = 1$. As a consequence of Theorem 3 we obtain the following result.

Corollary 1. Assume that Assumptions 1 and 2 are satisfied.

● If $S^{in} \leq g(D)$ then for any $r\in(0, 1)$, $P_r(S^{in}, D) < P(S^{in}, D)$.

● If $S^{in} > g(D)$ then $P_r(S^{in}, D) > P(S^{in}, D)$ if and only if $r\in (r_1, 1)$, where $r_1 = r_1(S^{in}, D)$ is the unique solution of $S^{in} = g_r(D)$

and $P_r(S^{in}, D) = P(S^{in}, D)$ for $r = 0$, $r = r_1$ and $r = 1$.

Proof. The proof is given in Appendix C.2.

This Corollary ensures that if $S^{in} > g(D)$ and for any $r_1 < r < 1$, the productivity of the biomass of the serial configuration is larger than the one of the simple chemostat. These conditions, related to the productivity of the biomass, are the same conditions that arose in the case of the minimization of the output substrate concentration, see section 3.1. We illustrate this Corollary in section 4 in Figure 8.

3.3.   Biogas flow rate

Let us consider the dependency of the biogas flow rate with respect to the dilution rate $D$ and the input concentration $S^{in}$. Recall that for a simple chemostat the output biomass at steady state is given by $x^{out} = S^{in}-S^{out}$. Classically, the biogas flow rate at steady-state of the simple chemostat model is defined as

The biogas flow rate of the serial configuration of two chemostats is the sum with the same propositional coefficient kept equal to one, which is thus defined as

with $V_i$ the volume, $x^{out, i}$ the output steady-state biomass and $S^{out, i}$ the output steady-state substrate concentration, all corresponding to the tank $i = 1, 2$. In this respect, for $r = 0$ and $r = 1$ we have $G_r(S^{in}, D) = G(S^{in}, D)$ and when $r\in (0, 1)$ it is formulated by

Proposition 5. For any $D\in[0, m)$, $S^{in} > 0$ and $r\in (0, 1)$, one has $G_r(S^{in}, D) = P_r(S^{in}, D)$.

Proof. The proof is given in Appendix C.3.

We know that for a single chemostat, the biogas flow rate and the productivity of the biomass at steady state are identical. Proposition 5 asserts this same conclusion in the case of two serial interconnected chemostats. Thereby, we deduce that analyzing the productivity of the biomass or the biogas flow rate at the steady state of two interconnected chemostats are equivalent. In this respect, Corollary 1 and the following result are verified for both performance criteria.

Proposition 6. Let $S^{in} > 0$. Let $G^{max}(S^{in}): = \max_{D\in(0, f(S^{in}))}G(S^{in}, D)$. For any $D > 0$ and $r\in(0, 1)$, one has $G_r(S^{in}, D) < G^{max}(S^{in})$.

Proof. The proof is given in Appendix C.4.

The two functions $D\mapsto G_r(S^{in}, D)$ and $D\mapsto G(S^{in}, D)$ are depicted in Figure 4. It shows that, for fixed values $S^{in}$ and $r$, the biogas production of the serial configuration of two chemostats is more efficient than the one of the single chemostat if and only if $0 < D < D_r$ with $D_r$ solution of $S^{in} = g_r(D)$, as it was proved in Proposition 3. In addition, Proposition 6 guarantees that the biogas flow rate of the serial device will never exceed the maximal biogas flow rate of the single chemostat. In other words, the extrema of the blue curve of the serial configuration will never exceed the extremum of the black curve of the simple chemostat.

This result has been graphically shown in ^[12] and ^[20] for the productivity of the biomass in the particular case of the Monod growth function. The simulations depicted in these references predicted that spatialization as we proposed it, does not give a better productivity of the biomass than a simple chemostat. According to Proposition 5, we know that at steady-state, the biogas flow rate and the productivity of the biomass are the same, which explains why predictions of the authors of ^[12] and ^[20] correspond to Proposition 6.

4.   Operating diagram

The operating diagram is the bifurcation diagram for which the values of the biological parameters are fixed. The various regions of the operating diagram reflect qualitatively different dynamics. The operating parameters which are the input concentration $S^{in}$ and the dilution rate $D$ of the chemostat can be chosen by the practitioners and the behavior of the model is discussed with respect to them. In contrast, the biological parameters are the ones of the growth function since they depend on the organisms, the substrates and the conversion rate $Y$, and are usually estimated in the laboratory.

Let the curves $\Phi_r$ and $\Phi_{1-r}$ in the $(S^{in}, D)$ positive plane be defined by

The curves $\Phi_r$ and $\Phi_{1-r}$ split the positive plane $(S^{in}, D)$ in several regions denoted $I_0(r)$, $I_1(r)$, $I_2(r)$ and $I_3(r)$ defined by:

We fix $r$ in $(0, 1)$ and we depict in the plane $(S^{in}, D)$ the regions in which the solution of system (2.6), with positive initial condition, globally converges towards one of the steady states $E_0$, $E_1$ or $E_2$. In the case $0\leq r < \frac{1}{2}$ [res. $\frac{1}{2} < r\leq 1$], the regions $I_0(r)$, $I_1(r)$ and $I_2(r)$ [res. $I_0(r)$, $I_2(r)$ and $I_3(r)$] form a partition of the positive plane. The region $I_1(r)$ for $\frac{1}{2}\leq r\leq1$ [res. $I_3(r)$ for $0\leq r\leq\frac{1}{2}$] is empty. The behavior of the system in each region is given in Table 1. Let the curve $\Gamma_r$ in the positive plane $(S^{in}, D)$ be defined by

with $g_r$ defined by (3.3).

Lemma 5. For all $r \in (0, 1)$ the curve $\Phi_r$ defined by (4.1) is always above the curve $\Gamma_r$ defined by (4.2) in the plane $(S^{in}, D)$.

Proof. The proof is given in Appendix D.3.

In this respect, for any growth function $f$ verifying Assumption 1, the operating diagram of system (2.6) looks like Figure 5. Thus, according to Theorem 2, for the minimization of the output substrate concentration criterion, the serial configuration is more efficient than the simple chemostat if and only if $S^{in} > g_r(D)$ i.e. if and only if $(S^{in}, D)$ is strictly below the curve $\Gamma_r$, see Figure 5. Let use the operating plane to give a better understanding of the results of Proposition 2 on the behavior of the function $r\mapsto S_r^{out}\left(S^{in}, D\right)$, according to $(S^{in}, D)$. To this end, we consider the curves $\Phi_1$, $\Gamma$ and $\Phi_{1/2}$ in the operating plane defined by:

with $g$ defined by (3.4).

Curves $\Gamma$ and $\Phi_{1/2}$ lie below curve $\Phi_1$. Therefore, curves $\Gamma$, $\Phi_{1/2}$ and $\Phi_1$ split the operating plane into at most five regions defined by:

These regions are shown in Figure 6 which is given for an illustrative example but does not correspond to any particular growth function.

Regions $J_0$, $J_1$ and $J_3$ always exist and are connected. However, regions $J_2$ and $J_4$ do not necessarily exist and if they exist, in general, they are not necessarily connected, depending on the relative positions of curves $\Gamma$ and $\Phi_{1/2}$. For instance, for linear growth rates, $\Gamma = \Phi_{1/2}$ and regions $J_2$ and $J_4$ do not exist (see section 5.1); for Monod growth function, curve $\Gamma$ is above curve $\Phi_{1/2}$ and region $J_4$ does not exist (see section 5.2); for Hill growth function, regions $J_2$ and $J_4$ both exist (see section 5.3) and are connected. Notice that for plotting operating diagrams we must choose the growth function $f$ and the values of the biological parameters, see Figures 9 and 14. We can state now the main result on function $r\mapsto S_r^{out}\left(S^{in}, D\right)$, for $(S^{in}, D)\in J_i$, $i = 0, ..., 4$.

Proposition 7. Let $J_i$, $i = 0, 1, ...4$ be defined by (4.5). The behavior of function $r\mapsto S_r^{out}\left(S^{in}, D\right)$, according to $(S^{in}, D)$ is as follows:

● If $(S^{in}, D)\in J_0$ then, for all $r\in[0, 1]$, $S^{out}_r(S^{in}, D) = S^{out}(S^{in}, D) = S^{in}$.

● If $(S^{in}, D)\in J_1$ then when $\lambda(D) < S^{in} < \lambda(2D)$, $S^{out}_r(S^{in}, D)$ is given by (3.5) and when $S^{in} = \lambda(2D)$, $S^{out}_r(S^{in}, D)$ is given by (3.6). In addition, for all $r\in(0, 1)$, $S^{out}_r(S^{in}, D) > S^{out}(S^{in}, D)$. The equality is fulfilled for $r = 0$ and $r = 1$, see Figure 7(a).

● If $(S^{in}, D)\in J_2$ then $S^{out}_r(S^{in}, D)$ is given by (3.5) and $S^{out}_r(S^{in}, D) < S^{out}(S^{in}, D)$ if and only if $r\in(r_1, 1)$ where $r_1 = r_1(S^{in}, D)$ is the unique solution of $S^{in} = g_r(D)$. The equality is fulfilled for $r = 0$, $r = r_1$ and $r = 1$, see Figure 7(b).

● If $(S^{in}, D)\in J_3$ then $S^{out}_r(S^{in}, D)$ is given by (3.6) and $S^{out}_r(S^{in}, D) < S^{out}(S^{in}, D)$ if and only if $S^{in} > g(D)$ and $r\in(r_1, 1)$ where $r_1 = r_1(S^{in}, D)$ is the unique solution of $S^{in} = g_r(D)$. The equality is fulfilled for $r = 0$, $r = r_1$ and $r = 1$, see Figure 7(c).

● If $(S^{in}, D)\in J_4$ then $S^{out}_r(S^{in}, D)$ is given by (3.6) and for all $r\in(0, 1)$, $S^{out}_r(S^{in}, D) > S^{out}(S^{in}, D)$. The equality is fulfilled for $r = 0$ and $r = 1$, see Figure 7 (d).

Proof. The proof is given in Appendix B.7

According to the regions depicted in Figure 6, we obtain Figure 7 which covers the whole possible cases of the behavior of the function $r\mapsto S^{out}_r(S^{in}, D)$. Thusly, we can minimize the output substrate concentration at the steady state by using a serial configuration of two interconnected chemostats instead of one chemostat if $(S^{in}, D)$ is fixed in the regions $J_2$ or $J_3$ (i.e. $S^{in} > g(D)$) and for $r_1 < r < 1$.

We have previously shown that Corollary 1 is a consequence of Theorem 3 and one can see in Proof C.2 of Corollary 1 that comparing the two quantities $P_r(S^{in}, D)$ and $P(S^{in}, D)$ involves the comparison of the two quantities $S^{out}_r(S^{in}, D)$ and $S^{out}(S^{in}, D)$. That is why, the curves representing the productivity of the biomass depicted in Figure 8 are analogous to the curves of Figure 7. In Figure 8, we fix $r\in(0, 1)$ and we plot the functions $r\mapsto P_r(S^{in}, D)$ and $r\mapsto P(S^{in}, D)$ for $(S^{in}, D)$ fixed in the regions $J_1$, $J_2$, $J_3$ and $J_4$. As in the case of the output substrate concentration, it is shown that the productivity of the biomass or the biogas flow rate of the serial configuration is larger than the one of the simple chemostat if and only if $r\in (r_1, 1)$ and $(S^{in}, D)$ is fixed in one of the regions $J_2$ or $J_3$.

5.   Applications and numerical illustrations

In this section, we consider three different kinetics: the linear function, the Monod function and the Hill function. Table 2 gives the analytical expressions of most of the results previously presented. These expressions show that an analytical study of the different performance criteria is possible.

5.1.   Linear function

We consider $f$ as a linear function defined by $f(S) = aS$. According to Table 2, remark that $\lambda(2D) = g(D)$ then, the curves $\Phi_{1/2}$ and $\Gamma$ defined respectively by (4.3) and (3.7) merge and constitute only one curve. The behavior of the maps $r\mapsto S^{out}_r(S^{in}, D)$ and $r\mapsto S^{out}(S^{in}, D)$ or $r\mapsto P_r(S^{in}, D)$ and $r\mapsto P(S^{in}, D)$ depends on the position of $(S^{in}, D)$ in the three regions $J_i$, $i = 0, 1, 3$ represented in Figure 9$(a)$. These regions are defined by

For a fixed value of $D$, the passageway form the region $J_0$ to $J_1$ is defined by the critical value $S^{in}_0 = \lambda(D)$ and the passageway form the region $J_1$ to $J_3$ is defined by the critical value $S^{in}_1 = g(D) = \lambda(2D)$ as shown in Figures 9$(a)$ and 10$(b)$. As stated in Lemma 2, for any $S^{in} > S^{in}_1$ there exists a threshold $r_1 = r_1(S^{in}, D)$ solution of $S^{in} = g_r(D)$ which is explicitly given by

Then, according to the three performance criteria which are the minimization of the output substrate concentration, the maximization of the productivity of the biomass and the maximization of the biogas flow rate, the serial configuration is more efficient than the simple chemostat if and only if $S^{in} > S^{in}_1$ and $r\in (r_1, 1)$. This result is illustrated in Figure 10 for minimization of the output substrate concentration criterion. Figure 10$(a)$ should be compared with Figure 6 of ^[16], where the part of the curves represented in Figure 10$(a)$ corresponding to $r > r_0$, for which $S^{out}_r(S^{in}, D) = S_2^*(S^{in}, D, r)$, are depicted. Indeed, in ^[16], the authors were only interested in the case where the positive equilibrium $E_2$ is GAS. The threshold $S^{in}_1 = 2$ shown in Figure 6 of ^[16] is given by $S^{in}_1 = g(1)$ and for any $S^{in} > 2$ the threshold $r_1(S^{in}, D)$ is explicitly given by (5.1).

5.2.   Monod function

The Monod function is given by $f(S) = mS/(K+S)$, see the second line of Table 2.

Lemma 6. The curve $\Gamma$ is located strictly above the curve $\Phi_{1/2}$ in the $(S^{in}, D)$ plane.

Proof. The proof in given in the Appendix D.4

Thus, considering a Monod function induces four regions $J_i$, $i = 0, 1, 2, 3$ in the operating plane, that describe the behaviors of the maps $r\mapsto S^{out}_r(S^{in}, D)$ and $r\mapsto P_r(S^{in}, D)$, which depend on the position of $(S^{in}, D)$ in these regions, as depicted in Figure 9$(b)$. The behavior of the map $r\mapsto S^{out}_r(S^{in}, D)$ through these regions is depicted in Figure 11$(a)$. For a fixed dilution rate $D$, the limit curves $\Phi_1$, $\Gamma$ and $\Phi_{1/2}$ define critical values denoted $S^{in}_0 = \lambda(D)$, $S^{in}_1 = g(D)$ and $S^{in}_2 = \lambda(2D)$, that respectively characterize the passageways between the regions $J_i \mbox{, } i = 0, 1, 2, 3$, see Figures 9$(b)$ and 11$(b)$. As stated in Lemma 2, for any $S^{in} > S^{in}_1$ there exists a threshold $r_1 = r_1(S^{in}, D)$ solution of $S^{in} = g_r(D)$ which is explicitly given by

Then, according to the three studied performance criteria, the serial configuration is more efficient than the simple chemostat if and only if $S^{in} > S^{in}_1$ and $r_1 < r < 1$. Figure 11$(a)$ should be compared with Figure 9 of ^[16], where the part of the curves represented in Figure 11$(a)$ corresponding to $r > r_0$, for which $S^{out}_r(S^{in}, D) = S_2^*(S^{in}, D, r)$, are depicted. Indeed, in ^[16], the authors were only interested to the case where the positive equilibrium $E_2$ is GAS. If $D = 1$ as shown in Figure 9$(b)$, the threshold $S^{in}_1$ is given by $S^{in}_1 = g(1) = 2.2$ and for any $S^{in} > 2.2$ the threshold $r_1$ is explicitly given by (5.2). Notice that Figures 10$(a)$ and Figure 11$(a)$ illustrate Proposition 1. As stated in this Proposition, when $D$ is fixed, one can remark that when increasing $S^{in}$, the output substrate concentration at the steady state decreases. Thus, the minimum of the curve $r\mapsto S^{out}_r(S^{in}, D)$, representing the optimal point that gives the best possible serial configuration, decreases as $S^{in} > S_1^{in} = g(D)$ and $S^{in}$ increases.

For the purpose of comparing the productivity of the biomass of both configurations, for a fixed $r\in(0, 1)$, we characterize the operating parameters $(S^{in}, D)$ that allow the optimal biomass productivity of the serial configuration. Let $\Delta_r$ be the curve defined by

where $P_r$ is defined by (3.11). This curve is obtained numerically and depicted in the operating plane $(S^{in}, D)$, see Figure 12$(a)$, $(b)$. For the values of the parameters used in Figure 12$(a)$, corresponding to the case $0 < r < 1/2$, there exits a threshold $S^{in}\approx 0.84$ such that for $0 < S^{in} < 0.84$, the maximum of $P_r(S^{in}, D)$ is reached when $P_r(S^{in}, D) = VD(S^{in}-S_2^*(S^{in}, D, r))$, and for $S^{in} > 0.84$, it is reached when $P_r(S^{in}, D) = VD(S^{in}-\lambda(D/(1-r)))$, as shown in Figure 13. Therefore, for $0 < S^{in} < 0.84$, the maximum of $P_r(S^{in}, D)$ is reached when $E_2$ is stable, i.e. when $D < rf(S^{in})$, as illustrated for $S^{in} = 0.4$ in Figure 13$(a)$. That is why, for $0 < S^{in} < 0.84$, the curve $\Delta_r$ is strictly below the curve $\Phi_r$. In contrast, for $S^{in} > 0.84$ the maximum of $P_r(S^{in}, D)$ is reached when $E_1$ is stable, i.e., when $D \geq rf(S^{in})$, as illustrated for $S^{in} = 1.2$ in Figure 13$(c)$. That is why, for $S^{in} > 0.84$, the curve $\Delta_r$ is strictly above the curve $\Phi_r$. In the limit case $S^{in} = 0.84$, both maxima of $P_r(S^{in}, D)$ are equal, as shown in Figure 13$(b)$. This corresponds to the leap of the curve $\Delta_r$, shown in Figure 12$(a)$. On the other hand, for $1/2 < r < 1$, the equilibrium $E_1$ cannot be stable and $P_r(S^{in}, D) = VD(S^{in}-S_2^*(S^{in}, D, r))$, whenever it is positive. Therefore, its maximum is reached when the positive equilibrium $E_2$ is stable, that is why, the curve $\Delta_r$ is strictly below the curve $\Phi_r$, see Figure 12$(b)$.

According to Proposition 4, $\Gamma$ is the curve of equation $D = D^{opt}(S^{in})$, where $D^{opt}(S^{in})$ is defined in (3.9). In other words, $D^{opt}(S^{in})$ is the optimal dilution rate corresponding to the maximal productivity of the biomass, of the simple chemostat. We observe on Figure 12 that $\Delta_r$ is strictly below the curve $\Gamma$. Hence $D^{opt}(S^{in}) > D^{opt}_r(S^{in})$, as it was also depicted in Figure 4. We conjecture that this property is always verified.

5.3.   Hill function

For all $p > 1$, the non-concave Hill function is given by $f(S) = mS^{p}/(K^p+S^{p})$.

Proposition 8. The Hill function verifies Assumptions 2 and 3.

Proof. The proof is given in Appendix D.5

Proposition 8 shows that we can use effectively a non-concave growth function in our analysis. In the following, we consider the case where $p = 2$, see third line of Table 2.

Lemma 7. Let us denote $D_1: = m(3-\sqrt{5})/4$.

If $0 < D < D_1$ then the curve $\Phi_{1/2}$ defined by (4.3) is strictly above the curve $\Gamma$ defined in (4.4). In contrast, if $D_1 < D < \frac{m}{2}$ then the curve $\Phi_{1/2}$ is strictly below the curve $\Gamma$.

Proof. The proof is given in Appendix D.6

According to Lemma 7, considering an Hill function with $p = 2$ induces five regions $J_i$, $i = 0, 1, 2, 3, 4$ in the operating plane, defined in (4.5), that describe the behavior of the maps $r\mapsto S^{out}_r(S^{in}, D)$ and $r\mapsto S^{out}(S^{in}, D)$ or $r\mapsto P_r(S^{in}, D)$ and $r\mapsto P(S^{in}, D)$, which depends on the position of $(S^{in}, D)$ in these regions (see Figure 14). For a fixed dilution rate $D$, the limit curves $\Phi_1$, $\Gamma$ and $\Phi_{1/2}$ define critical values denoted $S^{in}_0 = \lambda(D)$, $S^{in}_1 = g(D)$ and $S^{in}_2 = \lambda(2D)$ that characterize the passageways between the different regions $J_i$, $i = 0, 1, 2, 3, 4$. Notice that, if $D < D_1$, as shown in Figure 14$(b)$, where $D_1$ is defined in Lemma 7 then, we have $S^{in}_1 > S^{in}_2$ and the behavior of the maps $r\rightarrow S^{out}_r(S^{in}, D)$ is as depicted in Figure 15$(b)$. Remark that, in this case, the region where $S^{out}_r(S^{in}, D) = S^{in}$ disappears before the emergence of the threshold $r_1$ solution of $S^{in} = g_r(D)$, that is, before the emergence of the region where the serial configuration is more efficient than the simple chemostat i.e. $S^{out}_r(S^{in}, D) < S^{out}(S^{in}, D)$. On the other hand, if $D > D_1$, as shown in Figure 14$(a)$ then, we have $S^{in}_1 < S^{in}_2$ and the behavior of the maps $r\rightarrow S^{out}_r(S^{in}, D)$ is as depicted in Figure 15$(a)$.

As stated in Lemma 2, for any $S^{in} > S^{in}_1$, there exists a threshold $r_1 = r_1(S^{in}, D)$ solution of $S^{in} = g_r(D)$ such that, for $r_1 < r < 1$, the performance of the serial configuration is more efficient than the one of the simple chemostat. In other words, the output substrate concentration at steady-state of the serial configuration is lesser than the one of the simple chemostat if and only if $(S^{in}, D)\in J_2\cup J_3$ and $r\in (r_1, 1)$.

6.   Conclusions

This work presents an in-depth mathematical study of a model of two serial interconnected chemostats with one species and a monotonic growth function. We analyze, at steady-state, three different performance criteria: the minimization of the output substrate concentration, the maximization of the productivity of the biomass and the maximization of the biogas flow rate. The aim is to compare with the performance of the single chemostat. A part of this paper extends some of the results published in ^[16] and presented in the thesis ^[21]. In these both references, the concavity of the function $f$ is a required assumption but this assumption is not necessary in our analysis. The thorough study of our model reveals three main results. First, we provide an explicit expression depending on the dilution rate $D$, that represents the threshold $S^{in}_1 = g(D)$ on the input concentration for the performance. We deduce that there exists a configuration of two tanks that is better than a single tank. Actually, through the optimization of the distribution of the volume $V$ and the threshold $S^{in}_1$, we distinguish which configuration is the best. Secondly, we infer that maximizing the production of the biomass is equivalent to maximize the biogas flow rate at steady-state even in the case of a serial device of two interconnected chemostats. At the end, we obtain the same conditions for the three performance criteria. Thus, reducing the output substrate concentration, maximizing the production of the biomass or maximizing the biogas flow rate at steady state involve the same conditions and the same threshold $S^{in}_1$. These conditions are necessary and sufficient to allow the best performance, and they are characterized by the input concentration $S^{in}$, the dilution rate $D$ and the parameter $r$. Finally, for deeper understanding, we depict the corresponding operating diagram of the model which describes the behavior of the steady states. This diagram presents the conditions which induce an optimal configuration with regions characterized by the parameter $r$ and the operating parameters $S^{in}$ and $D$.

To broaden and deepen the present work, a forthcoming paper will present the analysis of performance, of an extension, of the model of two serial interconnected chemostats, with death rates. This future work will also include a comparison with the simple chemostat with death rate.

Acknowledgments

The first author thanks the Algerian Government for her PhD grant. The authors thank the Euro-Mediterranean research network TREASURE (http://www.inra.fr/treasure) for support. The authors thank Jérôme Harmand for valuable and fruitful exchanges, and the anonymous reviewers for their constructive comments which have greatly improved this work.

Appendix

A.   Proof of Theorem 1

A.1.   Existence of equilibria

System (2.6) has a cascade structure. Let us consider $z_i(t) = S_i(t)+x_i(t)$ $(i = 1, 2)$ then, we have the following system:

One can easily show that $\lim_{t\rightarrow +\infty}z_i(t) = S^{in}$ $(i = 1, 2)$. Therefore $(x_1(t), x_2(t))$ satisfies an asymptotically autonomous dynamics, whose limiting system

is defined in the square $\Sigma: = [0, S^{in}]\times[0, S^{in}]$. System (A.2) has a cascade structure. It admits at most three equilibria:

with $x_2^*\in (0, S^{in})$ a solution, if it exists, of equation

The equilibria $E_0$, $E_1$ and $E_2$ of (2.6) corresponding to $e_0$, $e_1$ and $e_2$, respectively, have the same values $x_i$, $i = 1, 2$, and their corresponding $S_i$ are given by $S_i = S^{in}-x_i$, $i = 1, 2$. Note that $e_0$, $e_1$ and $e_2$ give

where $S_2^* = S^{in}-x_2^*$. This proves that one has $\overline{S}_2 = \lambda\left({D}/({1-r})\right)$ and $S^{*}_{1} = \lambda\left({D}/{r}\right)$ as stated in the Theorem. The equilibrium $e_0$, and hence the corresponding equilibrium $E_0$, always exists. The equilibrium $e_1$, exists if and only if $S^{in}-\lambda\left({D}/({1-r})\right) > 0$, that is $D < (1-r)f(S^{in})$, which is the condition of existence of $E_1$ in the Theorem. For the existence and uniqueness of $e_2$, note that $x_2^*$ is a solution of (1.3), if and only if $S_2^* = S^{in}-x_2^*$ satisfies $f(S_2^*) = h(S_2^*)$, which proves (2.9). The function $h$ is positive, strictly decreasing and $h(S_1^*) = 0$, where $S_1^* = \lambda(D/r)$, if and only if $S^{in} > \lambda\left({D}/{r}\right)$, see Figure 2. Thus, as $f$ is strictly increasing (see Assumption 1), there exists a unique solution of $h(S_2) = f(S_2)$ denoted $S_2^*$ in $[0, S_1^*)$. Therefore, the equilibrium $e_2$ exists if and only if $S^{in} > \lambda\left({D}/{r}\right)$, that is $D < rf(S^{in})$, which is the condition of existence of $E_2$ in the statement of the Theorem.

A.2.   Local stability

For the local stability, the Jacobian matrix associated to system (A.2) is defined by

The eigenvalues of this triangular matrix are its diagonal elements. For $e_0$, the eigenvalues are $-D/r+f(S^{in})$ and $-D/(1-r)+f(S^{in})$. Therefore ${e_0}$, and hence $E_0$, is LES if and only if $D > \max \{r, 1-r\}f(S^{in})$. For $e_1$, the eigenvalues are $-D/r+f(S^{in})$ and $-f'(\lambda(D/(1-r)))\left(S^{in}-\lambda(D/(1-r))\right)$. The second eigenvalue is negative if and only if $D < (1-r)f(S^{in})$, that is, $e_1$ exists. Therefore ${e_1}$, and hence $E_1$, is LES if and only if $rf(S^{in}) < D < (1-r)f(S^{in})$. For $e_2$, the eigenvalues are $-f'(\lambda(D/r))(S^{in}-\lambda(D/r))$ and $-D/(1-r)+f(S_2^*)-f'(S_2^*)x_2^*$, where $x_2^* = S^{in}-S_2^*$, and $S_2^*$ is the solution of Eq (2.9). As it is seen in Figure 2 $(a)$, $0 < S_2^* < S_1^* > S^{in}$. Hence $(S_1^*-S_2^*)/(S^{in}-S_2^*) < 1$. Thus, we have $h(S_2^*) < D/(1-r)$ and, according to (2.9), this is equivalent to $f(S_2^*) < D/(1-r)$. Consequently, the second eigenvalue is negative, which proves that $e_2$, and hence $E_2$, is LES, if and only if the first eigenvalue is negative. This condition is equivalent to $S^{in} > \lambda(D/r)$, i.e $D < rf(S^{in})$, which is the condition of existence of $E_2$.

A.3.   Global stability.

For the global asymptotic stability we use phase plane arguments, as in the proof of Proposition 7 in ^[22], or in section 2.1.2.3 of ^[4]. We give the details of the proof when $e_2$ exists. The case where $e_2$ does not exist but $e_1$ exists and the case where neither $e_2$ nor $e_1$ exist are similar. The isoclines $x_1 = S^{in}-\lambda(D/r)$ and $x_1 = \varphi(x_2)$, where $\varphi$ is defined by (A.3), separate the interior of $\Sigma$ into four regions defined by

Two cases must be distinguished, according to the existence, or not of $e_1$, see Figure 1. We consider the case where $e_1$ exists. The case where it does not exist is similar. The isocline $x_1 = \varphi(x_2)$ is as shown in Figure 16 $(b)$, that is, it is the graph of a strictly increasing function. Indeed, using the definition (A.3) of $\varphi$, we have

Note that $\varphi'(0) = 1-\frac{1-r}{D}f(S^{in})$. Therefore, $e_1$ exists if and only if $\varphi'(0) < 0$ as shown in the Figure. For $x_2\in(\overline{x}_2, S^{in})$, where $\overline{x}_2 = S^{in}-\lambda(D/(1-r))$ is the $x_2$ component of $e_1$, we have

which proves that $\varphi$ is strictly increasing. The vector field associated to (A.2) is horizontal if $x_1 = \varphi(x_2)$ and vertical if $x_1 = 0$ or $x_1 = S^{in}-\lambda(D/r)$. It is directed as shown in the Figure. Assume first that $\left(x_1(0), x_2(0)\right)\in I\cup III$. These regions are positively invariant. Since in $I$ [resp. $III$], $x_1(t)$ and $x_2(t)$ are strictly decreasing [resp. increasing], the following limits exist:

Therefore, $\left(x_{1\infty}, x_{2\infty}\right)$ is an equilibrium of (A.2), which belongs to the closure $\overline{I}$ or the closure $\overline{III}$. Since $e_0$, $e_1$ and $e_2$ (resp. $e_2$) are the only steady states in $\overline{III}$ (resp. $\overline{I}$) and, since $e_1$ attracts only solution with $x_1(0) = 0$ and $e_0$ attracts no solutions with positive initial conditions, it follows that

Assume now that $(x_1(0), x_2(0))\in IV$. If $(x_1(t), x_2(t))$ remains in $IV$ for all $t > 0$ then $x_1(\cdot)$ is strictly decreasing and $x_2(\cdot)$ is strictly increasing. Thus, the limits (1.4) exist. Hence, $\left(x_{1\infty}, x_{2\infty}\right)$ is an equilibrium of (A.2), which belongs to the closure $\overline{IV}$. Since $e_2$ is the only equilibrium in $\overline{IV}$, we conclude that (A.5) holds. If $(x_1(t), x_2(t))$ leaves the region $IV$, then it can only enter in the region $I$. Hence, as shown previously, it necessarily tends to $e_2$ and hence, (A.5) holds. The same argument shows that any solution starting with initial condition in $II$ always remains in $II$ and then converges to $e_2$ or leaves the region $II$, then enters necessarily in region $III$, and then, as shown previously, it tends to $e_2$. Therefore, $e_2$ is GAS in the interior of $\Sigma$. Using the theory of asymptotically autonomous systems (see Appendix F in ^[6]), we deduce that $E_2$ is GAS if and only if it exists.

B.   Output substrate concentration

B.1.   Proof of Proposition 1

Let $S_2^{*i} = S_2^{*}(S^{in, i}, D, r)$, $i = 1, 2$. Suppose that $S_2^{*1} \geq S_2^{*2}$. Since $f$ is increasing then, we have $f(S_2^{*1}) \geq f(S_2^{*2})$. Since $f(S_2^{*1}) = h_1(S_2^{*1})$ and $f(S_2^{*2}) = h_2(S_2^{*2})$ then, we have $h_1(S_2^{*1}) \geq h_2(S_2^{*2})$. Since $h_2 > h_1$ then, we have $h_2(S_2^{*2}) > h_1(S_2^{*2})$. Since $h_1$ is decreasing then, we have $h_1(S_2^{*2})\geq h_1(S_2^{*1})$. Therefore, we have $h_1(S_2^{*1}) > h_1(S_2^{*1})$ which is a contradiction. Hence $S_2^{*1} < S_2^{*2}$.

B.2.   Proof of Theorem 2

Recall that $S_2^*(S^{in}, D, r)$ is the unique solution of Eq (2.9). Let us first prove that

Since $f$ is strictly increasing and $h$ is strictly decreasing then, $S_2^*(S^{in}, D, r) < \lambda(D)$ is equivalent to $h(\lambda(D)) < f(\lambda(D)) = D$. Thus, using the definition of $h$, the condition $h(\lambda(D)) < D$ is written as

which is equivalent to $S^{in} > \lambda(D)+\left(\lambda(D/r)-\lambda(D)\right)/(1-r)$. Hence, according to the definition (3.3) of $g_r$, this is equivalent to $S^{in} > g_r(D)$. Notice also that the function $g_r$, defined by (3.3), satisfies

Therefore, one has $g_r(D) > \lambda(D/r)$.

Let us go now to the proof of the Theorem. Assume that $S^{in} > g_r(D)$. Then, $S^{in} > \lambda(D/r) > \lambda(D)$, so that, as shown by (2.5) and (3.1), we have

Therefore, using (B.1), we have $S^{out}_r(S^{in}, D) < S^{out}(S^{in}, D)$. Assume now that $S^{in}\leq g_r(D)$. When $r < 1/2$, three cases must be distinguished. First, if $\lambda(D) < \lambda(D/r) < S^{in}\leq g_r(D)$, then, by (2.5) and (3.1), we obtain (B.3). Hence, using (B.1), we have $S^{out}_r(S^{in}, D)\geq S^{out}(S^{in}, D)$. Secondly, if $\lambda(D) < \lambda(D/(1-r)) < S^{in}\leq \lambda(D/r)$ then, by (2.5) and (3.1), $S^{out}_r(S^{in}, D) = \lambda(D/(1-r))$ and $S^{out}(S^{in}, D) = \lambda(D)$. Therefore, we have $S^{out}_r(S^{in}, D) > S^{out}(S^{in}, D)$. Finally, if $S^{in}\leq\lambda(D)$, then $S^{out}_r(S^{in}, D) = S^{out}(S^{in}, D) = S^{in}$. When $r\geq 1/2$, the proof is similar, excepted that we must distinguish only two cases, $\lambda(D) < S^{in}\leq \lambda(D/r)$ and $S^{in}\leq\lambda(D)$.

In conclusion, for any $r\in(0, 1)$, $S^{out}_r(S^{in}, D) < S^{out}(S^{in}, D)$ if and only if $S^{in} > g_r(D)$.

B.3.   Proof of Lemma 2

Let $D < m$. From Assumption 2, the function $r\in(D/m, 1)\mapsto g_r(D)$ is strictly decreasing. From Assumption 1, we have $\lim_{r\to D/m}\lambda(D/r) = \lambda(m) = +\infty$. Thus, $\lim_{r\to D/m}g_r(D) = +\infty$. Using L'Hôspital's rule, one has $\lim_{r\to 1}g_r(D) = g(D)$. Then, using Intermediate Value Theorem, we deduce that for $S^{in} > g(D)$ there exists a unique $r_1 = r_1(S^{in}, D)$ in $(0, 1)$ such that $S^{in} = g_{r_1}(D)$. Since the function $r\mapsto g_r(D)$ is strictly decreasing then, $r > r_1(S^{in}, D)$ if and only if $S^{in} = g_{r_1}(D) > g_r(D)$ which ends the proof of the Lemma.

B.4.   Proof of Theorem 3

● From Assumptions 1 and 2, the function $r\in(D/m, 1)\mapsto g_r(D)$ is strictly decreasing. Thus, for any $r \in (0, 1)$, $g(D) < g_r(D)$. If $S^{in}\leq g(D)$ then $S^{in} < g_r(D)$ and according to Theorem 2 we deduce that $S^{out}_r(S^{in}, D) > S^{out}(S^{in}, D)$.

● If $S^{in} > g(D)$ then according to Lemma 2, there exists a unique $r_1 = r_1(S^{in}, D)$ in $(0, 1)$ such that $S^{in} = g_{r_1}(D)$, where for all $r > r_1$, we have $S^{in} > g_r(D)$. Thus, according to Theorem 2 we deduce that $S^{out}_r(S^{in}, D) < S^{out}(S^{in}, D)$.

The equality in the limiting cases $r = 0$ and $r = 1$ is already verified, see (3.2).

If $r = r_1$ then $S^{in} = g_{r_1}(D)$. According to (B.2), one has $\lambda(D/r_1) < g_{r_1}(D)$, that is, $\lambda(D/r_1) < S^{in}$. Thus, one has $S^{out}_{r_1}(S^{in}, D) = S_2^*(S^{in}, D, r_1)$ where $S_2^*(S^{in}, D, r_1)$ is the unique solution of $h(S_2)|_{r = r_1} = f(S_2)$. Consequently, one has $S_2^*(S^{in}, D, r_1) = \lambda(D)$ if and only if $h(\lambda(D))|_{r = r_1} = f(\lambda(D))$, which is equivalent to $D\left(\lambda(D/r_1)-\lambda(D)\right)/\left((S^{in}-\lambda(D))(1-r_1)\right) = D$, that is, $\lambda(D/r_1)-\lambda(D) = (1-r_1)(S^{in}-\lambda(D))$. Consequently, one obtains that $\lambda(D)+\left(\lambda(D/r_1)-\lambda(D)\right)/(1-r_1) = S^{in}$, which is equivalent to $g_{r_1}(D) = S^{in}$. This ends the proof of the Theorem.

B.5.   Proof of Proposition 2

Let us consider $r_0 = D/f(S^{in})$ i.e. $S^{in} = \lambda(D/r_0)$.

1) When $S^{in}\leq\lambda(D)$ one has, for all $r\in(0, 1)$, $\lambda(D)\leq\min\{\lambda(D/(1-r)), \lambda(D/r)\}$ i.e. $S^{in}\leq\min\{\lambda(D/(1-r)), \lambda(D/r)\}$. Then, according to (3.1) one has $S^{out}_r(S^{in}, D) = S^{in}$.

2) When $\lambda(D) < S^{in} < \lambda(2D)$, one has $r_0\in (1/2, 1)$. Firstly, if $0\leq r\leq 1-r_0$then, one has $\lambda(D/(1-r))\leq\lambda(D/r_0)\leq\lambda(D/r)$ i.e. $\lambda(D/(1-r))\leq S^{in}\leq\lambda(D/r)$. This is equivalent to $rf(S^{in})\leq D\leq (1-r)f(S^{in})$. According to (3.1), one has $S^{out}_r(S^{in}, D) = \lambda(D/(1-r))$. Secondly, if $1-r_0 \leq r \leq r_0$, one has $\lambda(D/r_0)\leq \min\{\lambda(D/(1-r)), \lambda(D/r)\}$ i.e. $S^{in}\leq \min\{\lambda(D/(1-r)), \lambda(D/r)\}$. According to (3.1), one has $S^{out}_r(S^{in}, D) = S^{in}$. Finally, if $r_0 < r\leq1$, one has $\lambda(D/r)\leq \lambda(D/r_0)$ i.e. $\lambda(D/r)\leq S^{in}$ then, according to (3.1), one has $S^{out}_r(S^{in}, D) = S_2^*(S^{in}, D, r)$. These all prove (3.5).

3) When $\lambda(2D)\leq S^{in}$ one has $r_0\in (0, 1/2]$. If $0\leq r\leq r_0$ then $\lambda(D/(1-r))\leq\lambda(D/r_0)\leq\lambda(D/r)$ i.e. $\lambda(D/(1-r))\leq S^{in}\leq\lambda(D/r)$. According to (3.1), one has $S^{out}_r(S^{in}, D) = \lambda(D/(1-r))$. If $r_0\leq r\leq1$ then $\lambda(D/r)\leq\lambda(D/r_0)$ i.e. $\lambda(D/r)\leq S^{in}$. According to (3.1), one has $S^{out}_r(S^{in}, D) = S_2^*(S^{in}, D, r)$. These all prove (3.6).

B.6.   Proof of Proposition 3

Let $r\in(0, 1)$. Form Assumption 3, the function $D\in[0, rm)\mapsto g_r(D)$ is strictly increasing. From Assumption 1, we have $\lim_{D\rightarrow rm}\lambda(D/r) = \lambda(m) = +\infty$. Thus, $\lim_{D\rightarrow rm}g_r(D) = +\infty$ and $g_r(0) = 0$. Then, using Intermediate Value Theorem, we deduce that for $S^{in} > 0$ there exists a unique $D_r = D_r(S^{in})$ in $[0, rm)$ such that $S^{in} = g_r(D_r)$. Since the function $D\mapsto g_r(D)$ is strictly increasing then, $0 < D < D_r(S^{in})$ if and only if $0 < g_r(D) < g_r(D_r) = S^{in}$. Consequently, according to Theorem 2, $g_r(D) < S^{in}$ if and only if $S^{out}_r(S^{in}, D) < S^{out}(S^{in}, D)$, which ends the proof of the Proposition.

B.7.   Proof of Proposition 7

The result is a direct consequence of Proposition 2 and Theorem 3. We give the details for regions $J_1$ and $J_2$. The proof for other regions is similar.

If $(S^{in}, D)\in J_1$ then, according to (4.5), $\lambda(D) < S^{in}\leq\min\left(g(D), \lambda(2D)\right)$. Therefore, $\lambda(D) < S^{in}\leq\lambda(2D)$. When $\lambda(D) < S^{in} < \lambda(2D)$, from Proposition 2, $S^{out}_r(S^{in}, D)$ is given by (3.5) and if $S^{in} = \lambda(2D)$ then $S^{out}_r(S^{in}, D)$ is given by (3.6). Now, using $S^{in}\leq g(D)$, from Theorem 3, we have for all $r\in(0, 1)$, $S^{out}_r(S^{in}, D) > S^{out}(S^{in}, D)$.

If $(S^{in}, D)\in J_2$ then, according to (4.5), $g(D) < S^{in} < \lambda(2D)$. Therefore, $\lambda(D) < S^{in} < \lambda(2D)$ and from Proposition 2, $S^{out}_r(S^{in}, D)$ is given by (3.5). Now, using $g(D) < S^{in}$, from Lemma 2, there exists a threshold $r_1$, defined as the unique solution of $S^{in} = g_r(D)$. Therefore, from Theorem 3, $S^{out}_r(S^{in}, D) < S^{out}(S^{in}, D)$ if and only if $r\in(r_1, 1)$ and equality holds if and only if $r = 0$, $r = r_1$ or $r = 1$.

C.   Productivity and biogas production

C.1.   Proof of Proposition 4

The equation $P(S^{in}, D) = 0$ admits the two roots $D = 0$ and $D = f(S^{in})$. For all $D > 0$, $P(S^{in}, D)$ is positive if and only if $S^{in} > \lambda(D)$. In addition, for all $D < f(S^{in})$ we have

with $g$ defined by (3.4). Thus, $\frac{\partial P}{\partial D}(S^{in}, D) = 0$ is verified if and only if $S^{in} = g(D)$. Consequently, using Assumption 4, $D^{opt}$ defined in (3.9) is the unique solution of $S^{in} = g(D)$.

C.2.   Proof of Corollary 1

One knows that $x^{out}_r(S^{in}, D) = S^{in}-S^{out}_r(S^{in}, D)$ and $x^{out}(S^{in}, D) = S^{in}-S^{out}(S^{in}, D)$. Firstly, if $S^{in}\leq g(D)$ then according to Theorem 3, for any $r\in(0, 1)$, one has $S^{out}_r(S^{in}, D) > S^{out}(S^{in}, D)$. Thus, for any $r\in(0, 1)$, one has $x^{out}_r(S^{in}, D) < x^{out}(S^{in}, D)$. Consequently, for any $r\in(0, 1)$, one has $P_r(S^{in}, D) < P(S^{in}, D)$. Secondly, if $S^{in} > g(D)$ then according to Theorem 3, one has $S^{out}_r(S^{in}, D) < S^{out}(S^{in}, D)$ if and only if $r_1 < r < 1$ with $r_1$ defined in Lemma 2. Then, one has $x^{out}_r(S^{in}, D) > x^{out}(S^{in}, D)$ if and only if $r_1 < r < 1$. Consequently, one has $P_r(S^{in}, D) > P(S^{in}, D)$ if and only if $r_1 < r < 1$. Finally, if $r = 0$, $r = r_1$ or $r = 1$, then one has $S^{out}_r(S^{in}, D) = S^{out}(S^{in}, D)$. Thus, for $r = 0, r_1, 1$ one has $x^{out}_r(S^{in}, D) = x^{out}(S^{in}, D)$. Consequently, if $r = 0$, $r = r_1$ or $r = 1$ then, $P_r(S^{in}, D) = P(S^{in}, D)$, which ends the proof of the Corollary.

C.3.   Proof of Proposition 5

Let $V$ be a fixed volume. In the following, we use the respective definitions (3.11) and (3.14) of $P_r$ and $G_r$. In both cases: $\max\{r, 1-r\}f(S^{in})\leq D$ and $rf(S^{in})\leq D \leq(1-r)f(S^{in})$, it is clear that $G_r(S^{in}, D) = P_r(S^{in}, D)$. In addition, if $D < rf(S^{in})$ then

with $S_2^*$ the unique solution of (2.9). According to this equation, $G_r$ can be written as

Thus, we deduce that $G_r(S^{in}, D) = P_r(S^{in}, D) = VD(S^{in}-S_2^*)$ and consequently, for any $r\in(0, 1)$, we have $G_r(S^{in}, D) = P_r(S^{in}, D)$.

C.4.   Proof of Proposition 6

Let $V$ be a fixed volume and $S^{in} > 0$. Let us consider the function $\varphi(S) = f(S)(S^{in}-S)$. Considering the change of variable $S = \lambda(D)$, one can easily verify that $\varphi'(S) = 0$ is equivalent to $S^{in}-g(D) = 0$. According to Assumption 4, $\varphi$ admits a unique maximum. We maximize the biogas flow rate at steady-state with respect to $D$. On the one hand, the biogas flow rate of the simple chemostat is defined by $G(S^{in}, D) = V\varphi(S^{out}(D))$ with $S^{out}$ defined by (2.5). Then, the maximal biogas flow rate of the simple chemostat is $G^{max}(S^{in}) = V\max_{D\in(0, f(S^{in}))} \varphi(S^{out}(D))$. Since the map $\lambda$ defines a homeomorphism from $[0, f(S^{in})]$ to $[0, S^{in}]$ then $\max_{D\in(0, f(S^{in}))}\varphi(S^{out}(D)) = \max_{S\in(0, S^{in})}\varphi(S)$. On the other hand, as $S^*_1 > S^*_2$ and using the definition (3.13) of $G_r$, the biogas flow rate of the two serial interconnected chemostats at steady-state is defined by $G_r(S^{in}, D) = rV\varphi(S_1^*)+(1-r)V\varphi(S_2^*)$ with $S^*_1 = \lambda\left({D}/{r}\right)$ and $S_2^*$ the unique solution of (2.9). In addition, as for all $D < f(S^{in})$ we have $\varphi(S^*_i(D)) < \max_{S\in(0, S^{in})}\varphi(S)$, $i = 1, 2$ then, we have

Hence, we deduce that $G_r(S^{in}, D) < V\max_{S\in(0, S^{in})}\varphi(S)$ which is equivalent to $G_r(S^{in}, D) < G^{max}(S^{in})$. This completes the proof of the Proposition.

D.   Technical results

D.1.   Proof of Lemma 3

Using $l_D(r) = \lambda(D/r)$, $\gamma(r, D) = g_r(D)$, defined by (3.3), is given by

The partial derivative, with respect to $r$ of $\gamma$ is given then by

Therefore, $\frac{\partial\gamma}{\partial r} (r, D) < 0$ if and only if $l_D(1) > l_D(r)+(1-r)l_D'(r)$, which proves the equivalence of conditions 1 and 2 of the Lemma.

Moreover, if $l_D$ is strictly convex on $(D/m, 1]$ then for all $s$ and $r$ in $(D/m, 1]$, if $s\neq r$, then

Taking $s = 1$ and $r\in(D/m, 1)$ one obtains the condition 2.

Assume now that $f$, and hence $l_D$, are twice derivable. Using $\lambda'(D) = 1/f'(\lambda(D))$ and $\lambda''(D) = -f''(\lambda(D))/\left(f'(\lambda(D))\right)^3$, we can write

Therefore, the condition 3 is equivalent to the following condition:

Using the notation $S = \lambda(D/r)$, which is the same as $f(S) = D/r$, the condition () is equivalent to : for all $S > 0$, $f(S)f''(S) < 2\left(f'(S)\right)^2$, which is the condition 4 in the Lemma.

D.2.   Proof of Lemma 4

Using $\lambda'(D) = 1/f'(\lambda(D))$, the partial derivative, with respect to $D$ of $\gamma(r, D) = g_r(D)$, defined by (3.3), is given by

Therefore, $\frac{\partial\gamma}{\partial D} (r, D) > 0$ if and only if $f'\left(\lambda(D/r)\right) < f'\left(\lambda(D)\right)/r^2$, which proves the equivalence of conditions 1 and 2 of the Lemma.

Moreover, since $1/r > 1$ and $\lambda$ is strictly increasing, then $\lambda(D/r) > \lambda(D)$. Thus, if $f'$ is decreasing, we have $f'\left(\lambda(D/r)\right)\leq f'\left(\lambda(D)\right) < f'\left(\lambda(D)\right)/r^2$, which proves condition 2 of the Lemma.

D.3.   Proof of Lemma 5

As $0 < r < 1$ and $\lambda$ is a strictly increasing function then, we have $D/r > D$ and $\lambda(D/r) > \lambda(D)$. Consequently, using the definition (B.2) of $g_r$, we have $g_r(D) > \lambda(D/r)$. According to the respective definitions (4.1) and (4.2) of the curves $\Phi_r$ and $\Gamma_r$, we deduce that the curve $\Phi_r$ is always above the curve $\Gamma_r$.

D.4.   Proof of Lemma 6

The curves $\Phi_{1/2}$ and $\Gamma$ are respectively defined by (4.3) and (4.4). Let us define the function $H:\left[0, \frac{m}{2}\right)\mapsto \mathbb{R}$ such that $H(D) = \lambda(2D)-g(D)$. According to Table $\ref{TableApplicationsLinearAndMonod}$, for a Monod function, the function $H$ is defined explicitly by $H(D) = \frac{KmD^2}{(m-D)^2(m-2D)}$. Then, for any $D\in \left[0, \frac{m}{2}\right)$ one has $H(D) > 0$. Thus, for any $D\in \left[0, \frac{m}{2}\right)$, one has $\lambda(2D) > g(D)$ which means that $\Gamma$ is always located strictly above $\Phi_{1/2}$.

D.5.   Proof of Proposition 8

Let us prove that the Hill function satisfies Assumption 2. Straightforward computations show that

Hence, for every $p\geq 1$, $F'(S) = -\frac{p+1}{K^p}S^{p-1} < 0$ and $F(0) = \frac{p-1}{p} < 1$, which proves that $F(S) < 1$ for all $S > 0$. Therefore, condition 4 of Lemma 3 is satisfied, which is a sufficient condition for Assumption 2 to hold.

Let us prove now that the Hill function satisfies Assumption 3. It is equivalent to prove that it satisfies the condition 2 of Lemma 4. Straightforward computations show that

We have $0 < r < 1$ and $D < rm$ then, obviously, we have $0 < m-D/r < m-D$ and $0 < rm-D < m-D$. Thus, we obtain the following inequality

Straightforward calculations give

Consequently, we have

which is equivalent to $f'(\lambda(D/r)) < f'(\lambda(D))/r$ and induces $f'(\lambda(D/r)) < f'(\lambda(D))/r^2$. This completes the proof of the Proposition.

D.6.   Proof of Proposition 7

Let the function $H: \left[0, \frac{m}{2}\right) \mapsto \mathbb{R}$ be defined by $H(D): = \lambda(2D)-g(D)$. According to the analytical expressions of Table 2, we have

Thus, $H(D) > 0$ gives $4mD^2-6m^2D+m^3 < 0$. The equation $4mD^2-6m^2D+m^3 = 0$ admits the two roots $D_1 = \frac{3-\sqrt{5}}{4}m$ and $D_2 = \frac{3+\sqrt{5}}{4}m$ such that $0 < D_1 < \frac{m}{2}$ and $\frac{m}{2} < D_2$. Therefore, for any $D\in \left(D_1, \frac{m}{2}\right)$ we have, $H(D) > 0$ which means that for any $D\in \left(D_1, \frac{m}{2}\right)$, the curve $\Phi_{1/2}$ is strictly below the curve $\Gamma$.