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

HAL
open science

# Performance Study of Two Serial Interconnected Chemostats with Mortality

Manel Dali-Youcef, Alain Rapaport, Tewfik Sari

► To cite this version:

Manel Dali-Youcef, Alain Rapaport, Tewfik Sari. Performance Study of Two Serial Interconnected Chemostats with Mortality. Bulletin of Mathematical Biology, 2022, 84 (10), 10.1007/s11538-022-01068-6. http://hal-03762535

HAL Id: hal-03762535

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

Submitted on 28 Aug 2022

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
Performance study of two serial
interconnected chemostats with mortality

Manel Dali-Youcef¹,³, Alain Rapaport¹ and Tewfik Sari²*

¹MISTEA, Univ Montpellier, INRAE, Institut Agro, Montpellier, France.

²*ITAP, Univ Montpellier, INRAE, Institut Agro, Montpellier, France.

³Laboratoire de Mathématiques, Université d'Avignon, France.

*Corresponding author(s). E-mail(s): tewfik.sari@inrae.fr;

Contributing authors: daliyoucef.manel@gmail.com;

alain.rapaport@inrae.fr;

Abstract

The present work considers the model of two chemostats in series when a biomass mortality is considered in each vessel. We study the performance of the serial configuration for two different criteria which are the output substrate concentration and the biogas flow rate production, at steady state. A comparison is made with a single chemostat with the same total volume. Our techniques apply for a large class of growth functions and allow us to retrieve known results obtained when the mortality is not included in the model and the results obtained for specific growth functions in both the mathematical literature and the biological literature. In particular, we provide a complete characterization of operating conditions under which the serial configuration is more efficient than the single chemostat, i.e. the output substrate concentration of the serial configuration is smaller than that of the single chemostat or, equivalently, the biogas flow rate of the serial configuration is larger than that of the single chemostat. The study shows that the maximum biogas flow rate, relative to the dilution rate, of the series device is higher than that of the single chemostat provided that the volume of the first tank is large enough. This non-intuitive property does not occur for the model without mortality.

**Keywords:** chemostat, gradostat, mortality, bifurcations, global stability, operating diagram, biogas production
2 Performance study of two serial interconnected chemostats

33

MSC Classification: 34D20, 34H20, 65K10, 92C75

34

# 1 Introduction

35

36

37

38

39

40

The mathematical model of the chemostat has received a great attention in the literature for many years (see for instance [16] and literature cited inside). This is probably due to its relative simplicity that can explain and predict quite faithfully the dynamics of real bioprocesses exploiting microbial ecosystems. It is today an important tool for decision making in industrial world, such as for dimensioning bioreactors or designing efficient operating conditions [13, 20].

41

42

43

44

45

46

47

48

Several extensions of the original model of the chemostat, considering spatial heterogeneity, have been proposed to better cope reality (see for instance [19]). Lovitt and Wimpenny has proposed the "gradostat" experimental device as a collection of chemostats of same volume interconnected in series [22, 23], which has led to the so-called "gradostat model" representing in a more general framework a gradient of concentrations [37, 40]. The gradostat model has been further generalized as the "general gradostat model" representing more general interconnection graphs with tanks of different volumes [38, 39].

49

50

51

52

53

54

55

56

57

58

59

Efficient use of a chemostat in practice relies on the analysis of its performance. The performance is considered for different criteria studied in the literature [31], among which the most common are: the output substrate concentration, the residence time, the biogas flow rate and the biomass productivity. Particular interconnection structures have been investigated and compared for the properties in terms of input-output performances (see for instance [5, 7, 15, 28]). It has been notably shown that a series of reactors instead of a single perfectly mixed one can significantly improve the performances of the bioprocess (in terms of matter conversion) while preserving the same residence time, or equivalently that the same performance can be obtained with a smaller residence time considering several tanks in series instead of a single one [14, 17, 24, 25, 47].

60

61

62

63

64

65

66

67

68

69

72

On another hand, it is known that in real processes, various growth conditions can be met and that it could be difficult to setup exactly the same perfect conditions in different reactors. These conditions include toxicity levels of culture media, which means more concretely that the consideration of a bacterial mortality, although often neglected compared to the removal rate, might be non avoidable and could also be variable. To the best of our knowledge, the possible impacts of mortality in the design of series of chemostats has not been yet studied in the literature, which is the purpose of the present work. Its contributions also cover interests in theoretical ecology for a better grasp of the interplay between spatial heterogeneity and mortality in resource-consumers models. Indeed, considering different removal rates in the classical chemostat model or more general ones allows to consider additional mortality terms [21, 29, 34, 44]. However, these mathematical studies have mainly

73
74

75

concern analyses of equilibria and stability and not the performances of the
system in presence of mortality.

76

77

78

79

80

81

82

83

In view of providing clear messages to the practitioners, we investigate
how the operating diagram of a series of two interconnected chemostats in
series is modified when considering different or identical mortality rates in
both tanks. Operating diagrams have proven to be a good synthetic tool to
summarize the possible operating modes, emphasized in [26] for its importance
for bioreactors. Indeed, such diagrams are more and more often constructed
both in the biological literature [26, 36, 41, 45] and the mathematical literature
[1, 2, 4, 9-12, 18, 31-33, 35, 42, 43].

84

85

86

87

Then, we study the performances in terms of conversion ratio and byprod-
uct production (such as biogas). As we shall see, several aspects are not
intuitive, which show that the consideration of mortality can significantly
modify the favorable operating conditions.

88

89

Along the paper, we use the abbreviations LES for locally exponentially
stable and GAS for globally asymptotically stable in the positive orthant.

90

91

92

93

94

95

96

97

98

99

100

101

102

The paper is organized as follows. Section 2 includes the introduction
of the mathematical model corresponding to the serial configuration of two
chemostats with mortality rate. Afterwards, Section 3 focuses on the study of
performances of the serial configuration with respect of the output substrate
concentration. Then, Section 4 considers the performances of the serial con-
figuration with respect of the biogas production. Next, Section 5 is devoted to
illustrations and numerical simulations and a conclusion is given in Section 6.
Moreover, we set up the single chemostat with mortality in Appendix A, while
Appendix B is devoted to the existence and stability analysis of the steady
states of the serial chemostat and Appendix C to its operating diagram. These
results are extension of former results, in the case without mortality [7], but
that have required to revisit significantly the mathematical proofs. Finally,
Appendix D contains technical proofs.

103

## 2 Presentation of the model

104

105

106

107

108

109

110

We consider two serial interconnected chemostats where the total volume $V$ is
divided into $V_1 = rV$ and $V_2 = (1-r)V$, with $r \in (0,1)$, as shown in Fig. 1.
The substrate and the biomass concentrations in the tank $i$ are respectively
denoted $S_i$ and $x_i$, $i = 1,2$. The input substrate concentration in the first
chemostat is designated $S^{in}$, the flow rate is constant and is designated by $Q$.
The output substrate concentration is the concentration of substrate in the
second tank $S^{out} = S_2$.

111

The mathematical model is given by the following equations:
4 Performance study of two serial interconnected chemostats

<figure><img src="image_3.png" /><figcaption>Fig. 1 The serial configuration of two chemostats respectively of volumes rV and (1 − r ) V .</figcaption></figure>

$$
\begin{align*}
\dot{S}_1 &= \frac{D}{r}(S^{in} - S_1) - f(S_1)x_1 \\
\dot{x}_1 &= -\frac{D}{r}x_1 + f(S_1)x_1 - ax_1 \tag{1} \\
\dot{S}_2 &= \frac{D}{1-r}(S_1 - S_2) - f(S_2)x_2 \\
\dot{x}_2 &= \frac{D}{1-r}(x_1 - x_2) + f(S_2)x_2 - ax_2,
\end{align*}
$$

112 where $\dot{S}_i = \frac{dS_i}{dt}$, $\dot{x}_i = \frac{dx_i}{dt}$, $i = 1, 2$, $f$ is the growth function such that $f(S_i)$ is the growth function of the substrate in the tank $i = 1, 2$, $a$ is the mortality rate of the biomass and $D = Q/V$ is the dilution rate of the whole structure. 113 114 The dilution rate of the first tank is $Q/V_1 = D/r$. The dilution rate of the 115 second tank is $Q/V_2 = D/(1-r)$. 116

117 Note that these equations are not valid for r = 0 and r = 1, which corre- 118 spond to a single chemostat. For sake of completeness, the useful results on 119 the single chemostat are given in Appendix A. The considered growth function 120 satisfies the following properties. 121

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

We define

$$
m := \sup_{S>0} f(S), \quad (m \text{ may be } +\infty). \tag{2}
$$

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

$$
\lambda(D) := f^{-1}(D) \quad \text{when} \quad 0 \le D < m. \tag{3}
$$

The particular case without mortality of the biomass ( a = 0) is studied in [7]. The results on the existence and stability of steady states of system (1) are very similar to the case without mortality. The details are given in Appendix B. The system can have up to three steady states: 123 124 125 126

• The washout steady state E$_{0}$ = ( S $^{in}$, 0 , S $^{in}$, 0).

• The steady state E$_{1}$ = ( S $^{in}$, 0 , S$_{2}$ , x$_{2}$ ) of washout in the first chemostat but 128 129

• The steady state E$_{2}$ = ( S ∗ $_{1}$, x ∗ $_{1}$, S ∗ $_{2}$, x ∗ $_{2}$) of persistence of the species in both 130 131
As in the case without mortality, see Table C2 in the Appendix, for any operating condition ($S^{in}$, $D$), one and only one of the steady-states $E_0$, $E_1$ and $E_2$, is stable. It is then globally asymptotically stable (GAS).

The operating diagram of the system is described in Appendix C. The operating diagram has as coordinates the input substrate concentration $S^{in}$ and the dilution rate $D$, and shows how the solutions of the system behave for different values of these two parameters. The regions constituting the operating diagram correspond to different qualitative asymptotic behaviors. The operating diagram of system (1) is depicted in Fig. 2.

The aim of this work is to establish a comparison of the performance of the serial configuration with ones of the single chemostat. In the following, we compare both structures according to two different criteria; the output substrate concentration and the biogas flow rate.

## 3 Output substrate concentration

The output substrate concentration measures the biodegradation of the input substrate by the overall device. The reduction of the output substrate concentration is one of the main objectives of the biological wastewater treatment, and its minimization is often addressed in the literature, see for example [46]. We assume that the serial configuration is functioning at a stable steady state. The output substrate concentration at steady state depends on the parameters $D$, $S^{in}$ and $r$, and will be denoted $S_r^{out}(S^{in}, D)$.

**Proposition 1** Assume that Assumption 1 is satisfied. The output substrate concentration at steady state of system (1) is given by

$$S_r^{out}(S^{in}, D) = \begin{cases} S^{in} & \text{if } S^{in} \le \min\left(\lambda\left(\frac{D}{1-r} + a\right), \lambda\left(\frac{D}{r} + a\right)\right) \\ \bar{S}_2 & \text{if } \lambda\left(\frac{D}{1-r} + a\right) \le S^{in} \le \lambda\left(\frac{D}{r} + a\right) \\ S_2^* & \text{if } S^{in} > \lambda\left(\frac{D}{r} + a\right) \end{cases} \quad (4)$$

where $\bar{S}_2 = \lambda\left(\frac{D}{1-r} + a\right)$ and $S_2^*$ is the unique solution of equation $h(S_2) = f(S_2)$. In this equation, the function $h$ is defined by:

$$h(S_2) = \frac{D+(1-r)a}{1-r} \frac{S_1^*-S_2}{b-S_2}, \quad (5)$$
where $S_1^* = \lambda\left(\frac{D}{r} + a\right)$ and $b = \frac{D(S_1^{in}-S_1^*)}{D+ra} + S_1^*.$

*Proof* The output substrate concentration at steady state of system (1) is equal to $S^{in}$, if $E_0$ is the GAS steady state. It is equal to $\bar{S}_2$ if $E_1$ is the GAS steady state and to $S_2^*$ if $E_2$ is GAS. According to Theorem 3 in the Appendix, $E_0$ is GAS if and only if

$$D \ge \max(r, 1-r)(f(S^{in}) - a),$$
6 Performance study of two serial interconnected chemostats

162 which is equivalent to

$$
S^{in} \leq \min \left( \lambda \left( \frac{D}{1-r} + a \right), \lambda \left( \frac{D}{r} + a \right) \right).
$$

163 On the other hand, using Theorem 3, $\bar{S}_2$ depends on $D$ and $r$ and we have
164 $\bar{S}_2 = \lambda \left( \frac{D}{1-r} + a \right)$. $E_1$ is GAS if and only if

$$
165 r(f(S^{in}) - a) \le D \le (1-r)(f(S^{in}) - a),
$$

which is equivalent to

$$
\lambda \left( \frac{D}{1-r} + a \right) \le S^{in} \le \lambda \left( \frac{D}{r} + a \right).
$$

166 Finally, using Theorem 3, we know that $S_2^*$ depends on parameters $S^{in}, D, r$. It
167 is the unique solution of equation $h(S_2) = f(S_2)$, where $h$ is defined by (5). On the
168 other hand $E_2$ is GAS if and only if the condition $D < r(f(S^{in}) - a)$ is satisfied,
169 which is equivalent to the condition $S^{in} > \lambda \left(\frac{D}{r} + a\right)$. $\square$
170

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

$$
S_0^{\text{out}}(S^{\text{in}}, D) = S_1^{\text{out}}(S^{\text{in}}, D) = S^{\text{out}}(S^{\text{in}}, D). \quad (6)
$$

where $S^{\text{out}}(S^{\text{in}}, D)$, which is the output substrate concentration of the single chemostat, is given by

$$
S^{\text{out}}(S^{\text{in}}, D) = \begin{cases} S^{\text{in}} & \text{if } S^{\text{in}} \le \lambda(D+a), \\ \lambda(D+a) & \text{if } S^{\text{in}} > \lambda(D+a). \end{cases} \tag{7}
$$

174 For more information on $S^{\text{out}}(S^{\text{in}}, D)$, see Appendix A.

175 The proof of (6), comes from the following remarks. First, we have
176 $\bar{S}_2(D, 0) = \lambda(D + a)$ and second, according to Lemma 9 in the Appendix, we
177 can extend $S_2^*(S^{\text{in}}, D, r)$, by continuity, to $r = 1$, by

$$
S_2^*(S^{\text{in}}, D, 1) = \lambda(D + a).
$$

178 Our aim in this section is to compare $S_r^{\text{out}}$ defined by (4) and (6) and $S^{\text{out}}$
179 defined by (7).

3.1 The serial configuration can be more efficient than
the single chemostat

182 We fix $r$ and we describe the set of operating conditions ($S^{\text{in}}$, $D$) for which

$$
S_r^{\text{out}}(S^{\text{in}}, D) < S^{\text{out}}(S^{\text{in}}, D), \qquad (8)
$$

that is to say, the serial configuration with volumes $rV$ and $(1-r)V$, is more
efficient than the single chemostat of volume $V$. For $r \in (0, 1)$, let $g_r : [0, r(m-a)) \mapsto \mathbb{R}$ defined by

$$
g_r(D) := \lambda \left( \frac{D}{r} + a \right) + \frac{r(D+ar)}{(1-r)(D+a)} \left( \lambda \left( \frac{D}{r} + a \right) - \lambda(D+a) \right). \quad (9)
$$
186 **Lemma 1** For $r \in (0,1)$ we have $g_r(D) > \lambda(\frac{D}{r} + a)$.

187 **Proof** As $0 < r < 1$ and $\lambda$ is an increasing function then, we have $\lambda(D/r + a) > \lambda(D + a)$. Using (9), we have $g_r(D) > \lambda(D/r + a)$. $\square$

<figure><img src="image_4.png" /><figcaption>Fig. 2 The operating diagram of system (1) and the curve $\Gamma_r$ defined by (14) under which the serial configuration is more efficient than the single chemostat.</figcaption></figure>

189 **Theorem 1** Assume that Assumption 1 is satisfied. For any $r \in (0,1)$, we have

$$S_r^{\text{out}}(S^{\text{in}}, D) = S^{\text{out}}(S^{\text{in}}, D) \iff S^{\text{in}} = g_r(D).$$

Moreover,

$$S_r^{\text{out}}(S^{\text{in}}, D) < S^{\text{out}}(S^{\text{in}}, D) \iff S^{\text{in}} > g_r(D).$$

191 **Proof** Recall that $S_2^*(S^{\text{in}}, D, r)$ is the unique solution of equation $f(S_2) = h(S_2)$ with $h$ defined by (5). Let us first prove that

$$S_2^*(S^{\text{in}}, D, r) < \lambda(D + a) \iff S^{\text{in}} > g_r(D). \qquad (10)$$

193 Since $f$ is increasing, see Assumption 1, and $h$ is decreasing, see Lemma 8 in the Appendix, then the condition $S_2^*(S^{\text{in}}, D, r) < \lambda(D + a)$ is equivalent to the condition $h(\lambda(D + a)) < f(\lambda(D + a)) = D + a$. Using (5), a straightforward computation shows that the condition $h(\lambda(D + a)) < D + a$ is equivalent to $S^{\text{in}} > g_r(D)$, where $g_r$ is defined by (9). This proves (10).

198 Let us go now to the proof of the theorem. Assume that $S^{\text{in}} > g_r(D)$. Using Lemma 1, we have

$$S^{\text{in}} > \lambda(D/r + a) > \lambda(D + a).$$

200 Using (4) and (7), we have
8 Performance study of two serial interconnected chemostats

$$
\begin{align}
S_r^{\text{out}}(S^{\text{in}}, D) &= S_2^*(S^{\text{in}}, D, r), \nonumber \\
S^{\text{out}}(S^{\text{in}}, D) &= \lambda(D + a). \tag{11}
\end{align}
$$

From (10), we have $S_r^{\text{out}}(S^{\text{in}}, D) < S^{\text{out}}(S^{\text{in}}, D)$. Hence, we proved the following implication

$$
S^{\text{in}} > g_r(D) \implies S_r^{\text{out}}(S^{\text{in}}, D) < S^{\text{out}}(S^{\text{in}}, D). \quad (12)
$$

Assume now that $S^{\text{in}} \le g_r(D)$. When $r < 1/2$, three cases must be distinguished.
First, if

$$
\lambda(D + a) < \lambda\left(\frac{D}{r} + a\right) < S^{\text{in}} \le g_r(D),
$$

then, by (4) and (7), we obtain (11). Hence, using (10), we have $S_r^{\text{out}}(S^{\text{in}}, D) \ge S^{\text{out}}(S^{\text{in}}, D)$. Secondly, if

$$
\lambda(D + a) < \lambda\left(\frac{D}{1-r} + a\right) \le S^{\text{in}} \le \lambda\left(\frac{D}{r} + a\right),
$$

then, by (4) and (7), we have

$$
\begin{align*}
S_r^{\text{out}}(S^{\text{in}}, D) &= \lambda\left(\frac{D}{1-r} + a\right), \\
S^{\text{out}}(S^{\text{in}}, D) &= \lambda(D + a).
\end{align*}
$$

Therefore, we have $S_r^{\text{out}}(S^{\text{in}}, D) > S^{\text{out}}(S^{\text{in}}, D)$. Finally, if $S^{\text{in}} \le \lambda(D + a)$, then

$$
S_r^{\text{out}}(S^{\text{in}}, D) = S^{\text{out}}(S^{\text{in}}, D) = S^{\text{in}}.
$$

When $r \ge 1/2$, the proof is similar, excepted that we must distinguish only two
cases, $\lambda(D+a) < S^{\text{in}} \le \lambda(D/r+a)$ and $S^{\text{in}} \le \lambda(D+a)$. Hence, we have proved
the reciprocal implication of (12). This completes the proof of second equivalence in
the theorem.

The same calculations show the equivalence if inequalities are replaced by equal-
ities.
□

Theorem 1 asserts that the serial configuration is more efficient than the
single chemostat if and only if $S^{\text{in}} > g_r(D)$. Let us illustrate this result in the
operating diagram of system (1). Consider the curve of equation

$$
\Phi_r = \{(S^{\text{in}}, D) : S^{\text{in}} = \lambda(D/r + a)\}. \quad (13)
$$

According to the results given in Appendix C, the curves $\Phi_r$ and $\Phi_{1-r}$
defined by (13) separate the operating plane $(S^{\text{in}}, D)$ in four regions in which
the system has different asymptotic behaviour, see Table C2. To put it simply,
in the $I_0(r)$ region, $E_0$ is GAS, in $I_1(r)$, $E_1$ is GAS, and in $I_2(r) \cap I_3(r)$, $E_3$ is
GAS, see Fig. 2. This figure also shows the plot of the curve $\Gamma_r$, defined by

$$
\Gamma_r := \{(S^{\text{in}}, D) : S^{\text{in}} = g_r(D)\}. \quad (14)
$$

Using Lemma 1, we see that for all $r \in (0, 1)$, the curve $\Gamma_r$ is always at right of the curve $\Phi_r$. According to Theorem 1, the output substrate concentration of the serial configuration is smaller than the one of the single chemostat, if and only if $(S^{\text{in}}, D)$ is at right of the curve $\Gamma_r$ depicted in Fig. 2.
## 3.2 The output substrate concentration as a function of the volume fraction $r$

In this section we assume that $(S^{in}, D)$ is fixed and we look at the values of $r$ for which (8) holds. More precisely we are going to describe the function

$$r \mapsto S_r^{out}(S^{in}, D). \qquad (15)$$

**Proposition 2** Assume that Assumption 1 is satisfied. Let $D > 0$, $S^{in} > \lambda(a)$. We denote $r_0 = D/(f(S^{in}) - a)$.

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

2. If $\lambda(D + a) < S^{in} < \lambda(2D + a)$, then $1/2 < r_0 < 1$ and one has

$$S_r^{out}(S^{in}, D) = \begin{cases} \bar{S}_2 & \text{if } 0 \le r \le 1 - r_0 \\ S^{in} & \text{if } 1 - r_0 \le r \le r_0 \\ S_2^* & \text{if } r_0 \le r \le 1. \end{cases} \qquad (16)$$

3. If $\lambda(2D + a) \le S^{in}$, then $0 < r_0 \le 1/2$ and one has

$$S_r^{out}(S^{in}, D) = \begin{cases} \bar{S}_2 & \text{if } 0 \le r \le r_0 \\ S_2^* & \text{if } r_0 \le r \le 1. \end{cases} \qquad (17)$$

Here $\bar{S}_2 = \lambda(\frac{D}{1-r} + a)$ and $S_2^* = S_2^*(S^{in}, D, r)$ is the unique solution of equation $f(S_2) = h(S_2)$, where $h$ is defined by (5).

**Proof** If $S^{in} \le \lambda(D + a)$, then, for all $r \in (0, 1)$, one has

$$S^{in} \le \lambda(D + a) \le \min \left\{ \lambda \left( \frac{D}{1-r} + a \right), \lambda \left( \frac{D}{r} + a \right) \right\}.$$

Then, according to (4), one has $S_r^{out}(S^{in}, D) = S^{in}$. This proves item 1 of the proposition.

Let $r_0 = D/(f(S^{in}) - a)$, i.e. $S^{in} = \lambda(D/r_0 + a)$.

If $\lambda(D + a) < S^{in} < \lambda(2D + a)$, then $r_0 \in (1/2, 1)$, so that $1 - r_0 < r_0$. The interval $[0, 1]$ is subdivided into three sub-intervals. Firstly, if $0 \le r \le 1 - r_0 < r_0$, then $r < r_0 \le 1 - r$, so that

$$\lambda \left( \frac{D}{1-r} + a \right) \le S^{in} = \lambda \left( \frac{D}{r_0} + a \right) < \lambda \left( \frac{D}{r} + a \right).$$

Hence, according to (4), one has

$$S_r^{out}(S^{in}, D) = \lambda \left( \frac{D}{1-r} + a \right).$$

Secondly, if $1 - r_0 \le r \le r_0$, then $r_0 \ge \max\{r, 1 - r\}$, so that

$$S^{in} = \lambda \left( \frac{D}{r_0} + a \right) \le \min \left\{ \lambda \left( \frac{D}{1-r} + a \right), \lambda \left( \frac{D}{r} + a \right) \right\}.$$
10 Performance study of two serial interconnected chemostats

245 Hence, according to (4), one has $S_r^{out}(S^{in}, D) = S^{in}$. Finally, if $r_0 < r \le 1$, then one
246 has

247 Hence, according to (4), one has
$$S_r^{out}(S^{in}, D) = S_2^*(S^{in}, D, r).$$

248 This proves item 2 of the proposition.
249 If $\lambda(2D+a) \le S^{in}$, then $r_0 \in (0, 1/2]$. Therefore, $r_0 \le 1-r_0$. The proof of item 3
250 of the proposition is the same as the proof of item 2 excepted that now, the interval
251 $[0, 1]$ is subdivided now into two sub-intervals $[0, r_0]$ and $[r_0, 1]$, so that the interval
252 for which $S_r^{out}(S^{in}, D) = S^{in}$ is empty. $\square$

<figure><img src="image_5.png" /><figcaption>**Fig. 3** In each region $J_i$, $i = 0, \dots, 4$, the map $r \mapsto S_r^{out}(S^{in}, D)$ for fixed $(S^{in}, D)$ has a different behavior.</figcaption></figure>

253 We want to determine the values $r \in (0, 1)$ for which the condition (8) is satisfied. We need the following Assumption that is satisfied by any concave growth function but also by non concave growth functions, satisfying additional conditions, see Section 3.4.

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

259 Let $D < m - a$. Using $g_r(D) > \lambda(D/r + a)$, we have

$$ \lim_{r \to D/(m-a)} g_r(D) > \lim_{r \to D/(m-a)} \lambda(D/r + a) = +\infty. $$

260 On the other hand, using L'Hôpital's rule, we have
$$ \lim_{r \to 1} g_r(D) = g(D). \qquad (18) $$

where $g: [0, m-a) \to \mathbb{R}^+$ is defined by

$$ g(D) = \lambda(D + a) + D\lambda'(D + a). \qquad (19) $$

Therefore, from Assumption 2, the function $r \mapsto g_r(D)$ is decreasing from $(D/(m-a), 1)$ to $(g(D), +\infty)$. Hence, it admits an inverse function

$$ S^{in} \in (g(D), +\infty) \mapsto r_1(S^{in}, D) \in (D/(m-a), 1). $$

We use the notation $r_1(\cdot, D)$ to recall the dependence of the inverse function in $D$. For all $D \in (0, m-a)$, $r \in (D/(m-a), 1)$ and $S^{in} > g(D)$, we have

$$ r = r_1(S^{in}, D) \iff S^{in} = g_r(D), \qquad (20) $$

$$ r > r_1(S^{in}, D) \iff S^{in} > g_r(D). \qquad (21) $$

**Theorem 2** Assume that Assumptions 1 and 2 are satisfied. Let $g$ defined by (19).

* If $S^{in} \le g(D)$ then for any $r \in (0, 1)$, we have $S_r^{out}(S^{in}, D) > S^{out}(S^{in}, D)$. In addition, for $r = 0$ and $r = 1$ we have $S_r^{out}(S^{in}, D) = S^{out}(S^{in}, D)$.

* If $S^{in} > g(D)$ then $S_r^{out}(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 by (20). In addition, for $r = 0$, $r = r_1(S^{in}, D)$ and $r = 1$, we have $S_r^{out}(S^{in}, D) = S^{out}(S^{in}, D)$.

*Proof* The function $r \mapsto g_r(D)$ is decreasing and tends to $g(D)$ as $r$ tends to 1, as shown by (18). Thus, for all $r \in (0, 1)$, we have $g(D) < g_r(D)$. If $S^{in} \le g(D)$, then $S^{in} < g_r(D)$. According to Theorem 1, for all $r \in (0, 1)$, we have $S_r^{out}(S^{in}, D) > S^{out}(S^{in}, D)$.

Let $S^{in} > g(D)$. Let $r_1 = r_1(S^{in}, D)$. According to (21), for all $r > r_1$, we have $S^{in} > g_r(D)$. Thus, according to Theorem 1, we have $S_r^{out}(S^{in}, D) < S^{out}(S^{in}, D)$.

The equality $S_r^{out}(S^{in}, D) = S^{out}(S^{in}, D)$ is verified for the $r = 0$ and $r = 1$, see (6). In addition, we have $S^{in} = g_r(D)$, see (20). Hence, according to Theorem 1, we have $S_{r_1}^{out}(S^{in}, D) = S^{out}(S^{in}, D)$. $\square$

Let us now describe the subsets of the operational space $(S^{in}, D)$ for which the behaviour described in the three cases of Proposition 2 occurs. For a complete description we will also distinguish the sub-cases for which there exists $r_1 = r_1(S^{in}, D)$ such that, for $r_1 < r < 1$, (8) is satisfied, as shown in Theorem 2. Consider the curves $\Phi_1$ and $\Phi_{1/2}$, defined by (13), and the curve $\Gamma$ defined by

$$ \Gamma := \{(S^{in}, D) : S^{in} = g(D)\}, \qquad (22) $$

These three curves intersect at $(\lambda(a), 0)$ and, using the inequality $g(D) > \lambda(D + a)$, which is satisfied for all $D > 0$, one deduces that $\Gamma$ is at the right
12 Performance study of two serial interconnected chemostats

<figure><img src="image_6.png" /><figcaption>**Fig. 4** For $S^{\text{in}}$ and $D$ fixed, the output substrate concentration of the serial configuration, in red, compared to that of the single chemostat, in blue; $r_1(S^{\text{in}}, D)$ is defined by (20), $r_0 = D/(f(S^{\text{in}}) - a)$ and $J_1, J_2, J_3, J_4$ are depicted in Fig. 3.</figcaption></figure>

of $\Phi_1$. Therefore, the curves $\Phi_1, \Phi_{1/2}$ and $\Gamma$ separate the set of operating parameters $(S^{\text{in}}, D)$ into the following four subsets, see Fig. 3.

$$
\begin{align}
J_0 &= \{(S^{\text{in}}, D) : S^{\text{in}} \le \lambda(D+a)\}, \nonumber \\
J_1 &= \{(S^{\text{in}}, D) : \lambda(D+a) < S^{\text{in}} \le \min\{g(D), \lambda(2D+a)\}\}, \nonumber \\
J_2 &= \{(S^{\text{in}}, D) : g(D) < S^{\text{in}} < \lambda(2D+a)\}, \tag{23} \\
J_3 &= \{(S^{\text{in}}, D) : \max\{g(D), \lambda(2D+a)\} \le S^{\text{in}}\}, \nonumber \\
J_4 &= \{(S^{\text{in}}, D) : \lambda(2D+a) < S^{\text{in}} < g(D)\}. \nonumber
\end{align}
$$

Combining the results of Proposition 2 and Theorem 2, we find that the function $r \mapsto S_r^{\text{out}}(S^{\text{in}}, D)$ is as in Fig. 4. In the following we will comment on this figure.

• If $(S^{\text{in}}, D) \in J_1$, then when $S^{\text{in}} < \lambda(2D + a)$, $S_r^{\text{out}}(S^{\text{in}}, D)$ is given by (16) and when $S^{\text{in}} = \lambda(2D + a)$, $S_r^{\text{out}}(S^{\text{in}}, D)$ is given by (17). In addition, for all $r \in (0, 1)$, $S_r^{\text{out}}(S^{\text{in}}, D) > S_r^{\text{out}}(S^{\text{in}}, D)$. The equality is fulfilled for $r = 0$ and $r = 1$, see Fig. 4(a).

• If $(S^{\text{in}}, D) \in J_2$, then $S_r^{\text{out}}(S^{\text{in}}, D)$ is given by (16) and $S_r^{\text{out}}(S^{\text{in}}, D) < S_r^{\text{out}}(S^{\text{in}}, D)$ if and only if $r \in (r_1(S^{\text{in}}, D), 1)$, where $r_1(S^{\text{in}}, D)$ is defined
by (20). The equality is fulfilled for $r = 0$, $r = r_1(S^{in}, D)$ and $r = 1$, see Fig. 4(b).

* If $(S^{in}, D) \in J_3$ then $S_r^{out}(S^{in}, D)$ is given by (17) and $S_r^{out}(S^{in}, D) < S^{out}(S^{in}, D)$ if and only if $S^{in} > g(D)$ and $r \in (r_1(S^{in}, D), 1)$ where $r = r_1(S^{in}, D)$ is defined by (20). The equality is fulfilled for $r = 0$, $r = r_1(S^{in}, D)$ and $r = 1$, see Fig. 4(c).

* If $(S^{in}, D) \in J_4$ then $S_r^{out}(S^{in}, D)$ is given by (17) and for all $r \in (0, 1)$, $S_r^{out}(S^{in}, D) > S^{out}(S^{in}, D)$. The equality is fulfilled for $r = 0$ and $r = 1$, see Fig. 4(d).

Note that if $(S^{in}, D) \in J_0$, then case 1 of Proposition 2 occurs. One remarks that the lowest value of the red curve, corresponding to the lowest output substrate concentration of the serial configuration, is obtained for $(S^{in}, D) \in J_2 \cap J_3$ and $r > r_1(S^{in}, D)$. This lowest concentration is obtained with the best possible serial configuration.

Figures 2, 3 and 4 are made without graduations on the axes because they represent general situations where the growth function is only assumed to verify our hypotheses. It should be noticed that regions $J_0$, $J_1$ and $J_3$ always exist and are connected. However, regions the $J_2$ and $J_4$ do not always exist or are necessarily connected. This depends on the number of points of intersection between curves $\Phi_{1/2}$ and $\Gamma$. For a linear growth rate, $\Phi_{1/2} = \Gamma$ and hence, regions $J_2$ and $J_4$ do not exist, see Fig. 7(a). For a Monod growth function, curves $\Phi_{1/2}$ and $\Gamma$ intersect only at point $(\lambda(a), 0)$ and hence, region $J_3$ always exist and is connected but region $J_3$ does not exist, see Fig. 8(a). For a Hill growth function, curves $\Phi_{1/2}$ and $\Gamma$ always intersect at $(\lambda(a), 0)$ and also at a unique positive point, Lemma 6. Hence, regions $J_2$ and $J_4$ both exist and are connected, see Fig. 9(a,b,c).

## 3.3 The output substrate concentration as a function of the dilution rate

In this section we assume that $S^{in}$ and $r$ are fixed and we look at the values of the dilution rate $D$ for which (8) holds, i.e. the serial configuration, is more efficient than the single chemostat. More precisely we are going to describe the function

$$D \mapsto S_r^{out}(S^{in}, D). \qquad (24)$$

We want to determine the subset of values of $D$ for which the condition (8) is satisfied. We need the following Assumption that is satisfied by any concave growth function, but also by non concave growth functions, satisfying additional conditions, see Section 3.4.

**Assumption 3** For every $r \in (0, 1)$, the function $D \in [0, r(m-a)) \mapsto g_r(D) \in \mathbb{R}$ is increasing.
14 Performance study of two serial interconnected chemostats

Using $g_r(D) > \lambda(D/r + a)$, we have

$$ \lim_{D \to r(m-a)} g_r(D) > \lim_{D \to r(m-a)} \lambda(D/r + a) = +\infty. $$

From Assumption 3, the function $D \mapsto g_r(D)$ is increasing from $[0, r(m-a))$ to $[g_r(0) = \lambda(a), +\infty)$. Hence, its admits an inverse function

$$ S^{in} \in (\lambda(a), +\infty) \mapsto D_r(S^{in}) \in [0, r(m-a)). $$

For all $r \in (0, 1)$, $S^{in} \ge \lambda(a)$ and $D \in [0, r(m-a))$, we have

$$ D = D_r(S^{in}) \iff S^{in} = g_r(D), \qquad (25) $$

$$ D < D_r(S^{in}) \iff S^{in} > g_r(D). \qquad (26) $$

**Proposition 3** Assume that Assumptions 1 and 3 are satisfied. We have

$$ S_r^{out}(S^{in}, D) < S^{out}(S^{in}, D) \iff 0 < D < D_r(S^{in}), $$

where $D_r(S^{in})$ is defined by (25).

*Proof* Let $r \in (0, 1)$. According to (26), if $D < D_r(S^{in})$, then $S^{in} > g_r(D)$. Consequently, according to Theorem 1, we have $S_r^{out}(S^{in}, D) < S^{out}(S^{in}, D)$. □

339
340

## 3.4 How to check Assumptions 2 and 3

In this section we give sufficient conditions for Assumption 2 and 3 to be

satisfied. These conditions will be useful for the applications given in Section

5. For this purpose we consider the function $\gamma$ defined by

$$ \gamma(r, D) = g_r(D), \qquad (27) $$

defined on

$$ \mathrm{dom}(\gamma) = \{(r, D) : 0 < r < 1, 0 < D/r + a < m\}, $$

which consists simply in considering $g_r(D)$, given by (9), as a function of both

variables $r$ and $D$. If

$$ \frac{\partial \gamma}{\partial r}(r, D) < 0 \text{ for all } (r, D) \in \mathrm{dom}(\gamma), $$

then Assumption 2 is satisfied. Similarly, if

$$ \frac{\partial \gamma}{\partial D}(r, D) > 0 \text{ for all } (r, D) \in \mathrm{dom}(\gamma), $$

then Assumption 3 is satisfied. The following Lemmas give sufficient conditions,

for partial derivatives of $\gamma$ to have their signs as indicated above.
**Lemma 2** For $D \in (0, m-a)$, let $l_D$ be defined on $dom(l_D) = (D/(m-a), 1]$ by $l_D(r) = \lambda(D/r + a)$.

a) Assume that for $D \in (0, m-a)$ and $r \in dom(l_D)$ we have
$$l_D(1) > l_D(r) + (1-r)l'_D(r) \quad (28)$$

then, for all $(r, D) \in dom(\gamma)$, we have $\frac{\partial\gamma}{\partial r}(r, D) < 0$.

b) If, for $D \in (0, m-a)$, $l_D$ is strictly convex on $dom(l_D)$, then the condition (28) is satisfied.

c) If $f$ is twice derivable, then $l_D$ is twice derivable and the following conditions are equivalent

1. For $D \in (0, m-a)$ and $r \in dom(l_D)$, $l''_D(r) > 0$.

2. For $S > \lambda(a)$, $(f(S) - a)f''(S) < 2(f'(S))^2$.

**Proof** Notice first that $\gamma(r, D)$ can be written as follows

$$\gamma(r, D) = g_r(D) = \lambda(D + a) + \left(\frac{1}{1-r} - \frac{ra}{D+a}\right) \left(\lambda\left(\frac{D}{r} + a\right) - \lambda(D + a)\right). \quad (29)$$

Using the definition of $l_D$, $\gamma(r, D)$ is given then by

$$\gamma(r, D) = l_D(1) + \left(\frac{1}{1-r} - \frac{ra}{D+a}\right) (l_D(r) - l_D(1)).$$

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

$$\frac{\partial \gamma}{\partial r}(r, D) = \frac{a(1-2r)}{D+a} l'_D(r) + \left(\frac{1}{(1-r)^2} - \frac{a}{D+a}\right) (l_D(r) - l_D(1) + (1-r)l'_D(r)). \quad (30)$$

Notice that $\frac{1}{(1-r)^2} - \frac{a}{D+a} > 0$ for all $r \in (0, 1)$. From $l'_D(r) = -\frac{D}{r^2}\lambda'\left(\frac{D}{r} + a\right)$, it is deduced that $l'_D(r) < 0$. Therefore, if the condition (28) is satisfied, and, in addition $0 < r \le 1/2$, then, from (30), it is deduced that $\frac{\partial\gamma}{\partial r}(r, D) < 0$.

In the case $r \in (1/2, 1)$, we use the following expression of $\gamma(r, D)$ which is deduced from (29):

$$\gamma(r, D) = l_D(1) + B(r) \frac{l_D(r) - l_D(1)}{1-r},$$

where $B(r) = \frac{D+a-ar(1-r)}{D+a}$. Straightforward computation show that

$$\frac{\partial \gamma}{\partial r}(r, D) = \frac{D+ar(2-r)}{(D+a)(1-r)^2} (l_D(r) - l_D(1) + (1-r)C(r)l'_D(r)), \quad (31)$$

where $C(r) = \frac{D+a-ar(1-r)}{D+ar(2-r)}$. We have

$$C'(r) = \frac{a}{(D+ar(2-r))^2} (ar^2 + 2(a+2D)r - 3D - 2a).$$

Thus $C'(r) = 0$ for

$$r = r^* := \frac{1}{a} \left( \sqrt{3a^2 + 7aD + 4D^2} - a - 2D \right) \in (1/2, 1)$$
16 Performance study of two serial interconnected chemostats

and $(r - r^*)C'(r) > 0$ for $r \in (1/2, 1)$, $r \neq r^*$. Hence, from $C(1/2) = C(1) = 1$, we have $0 < C(r) < 1$ for all $r \in (1/2, 1)$. Now, if we assume that (28) is satisfied, for $1/2 < r < 1$ we have

$$l_D(1) > l_D(r) + (1-r)l'_D(r) > l_D(r) + (1-r)C(r)l'_D(r).$$

361 Hence, from (31), it is deduced that $\frac{\partial \gamma}{\partial r}(r, D) < 0$. This proves part a of the lemma.

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

$$l_D(s) > l_D(r) + (s-r)l'_D(r).$$

Taking $s = 1$ and $r \in dom(l_D)$ one obtains the condition (28). This proves part b of the lemma. Assume now that $f$, and hence $l_D$, are twice derivable. Using

$$\lambda'(D) = \frac{1}{f'(\lambda(D))}, \quad \lambda''(D) = - \frac{f''(\lambda(D))}{(f'(\lambda(D)))^3}, \qquad (32)$$

we can write

$$l''_D(r) = \frac{2D}{r^3} \lambda'\left(\frac{D}{r} + a\right) + \frac{D^2}{r^4} \lambda''\left(\frac{D}{r} + a\right) \\ = \frac{D\left(2\left(f'(\lambda(\frac{D}{r}+a))\right)^2 - \frac{D}{r}f''(\lambda(\frac{D}{r}+a))\right)}{r^3\left(f'(\lambda(\frac{D}{r}+a))\right)^3}.$$

362 Therefore, the condition 1 in item c in the lemma is equivalent to the following

363 condition: For all $D \in (0, m-a)$ and $r \in (D/(m-a), 1]$, we have

$$\frac{D}{r} f''\left(\lambda\left(\frac{D}{r} + a\right)\right) < 2f'\left(\lambda\left(\frac{D}{r} + a\right)\right)^2. \qquad (33)$$

364 Using the notation $S = \lambda(\frac{D}{r} + a)$, which is the same as $D/r = f(S) - a$, the

365 condition (33) is equivalent to : For all $S > 0$, $(f(S) - a)f''(S) < 2(f'(S))^2$, which

366 is the condition 2 in c in the lemma. $\square$

**Lemma 3** Assume that

$$f'\left(\lambda\left(\frac{D}{r} + a\right)\right) \leq \frac{1}{r}f'\left(\lambda(D + a)\right). \qquad (34)$$

368 Then, $\frac{\partial \gamma}{\partial D}(r, D) > 0$. Hence Assumption 3 is satisfied. If $f'$ is decreasing, then the

369 condition (34) is satisfied.

**Proof** From (29) we deduce that

$$\frac{\partial \gamma}{\partial D}(r, D) = \lambda'(D + a) + \frac{ra}{(D+a)^2} \left( \lambda\left(\frac{D}{r} + a\right) - \lambda(D + a) \right) \\ + \left( \frac{1}{1-r} - \frac{ra}{D+a} \right) \left( \frac{1}{r} \lambda'\left(\frac{D}{r} + a\right) - \lambda'(D + a) \right).$$

371 Notice that $\frac{1}{1-r} - \frac{ra}{D+a} > 0$, $\lambda'(D + a) > 0$ and $\lambda(\frac{D}{r} + a) > \lambda(D + a)$. Therefore

372 the condition

$$\frac{1}{r} \lambda' \left( \frac{D}{r} + a \right) - \lambda'(D + a) \geq 0$$

373 is sufficient to have $\frac{\partial \gamma}{\partial D}(r, D) > 0$. Using (32), this condition is equivalent to (34).

374 Note that if $f'$ is decreasing, then this condition is satisfied. Indeed, we have

$$f'\left(\lambda\left(\frac{D}{r} + a\right)\right) \leq f'\left(\lambda(D + a)\right) \leq \frac{1}{r}f'\left(\lambda(D + a)\right),$$

375 which is the condition (34). $\square$
Remark 1 Notice that:

i) The condition 2 in part c of Lemma 2 is equivalent to the condition

$$ \text{For all } S > \lambda(a), \frac{d^2}{dS^2} \left( \frac{1}{f(S)-a} \right) > 0. \qquad (35) $$

Therefore, if $f$ satisfies the condition (35), then it verifies Assumption 2.

ii) If the increasing growth function $f$ is twice derivable and satisfies $f''(S) \le 0$ for all $S > 0$, then the condition b in Lemma 2 and the condition (34) in Lemma 3 are satisfied. Thus, Assumptions 2 and 3 are satisfied and our results apply for any concave growth function.

iii) Assume that the increasing growth function $f$ is twice derivable and there exists $\hat{S} \in (0, +\infty)$ such that $f''$ is nonnegative on $(0, \hat{S})$ and nonpositive on $(\hat{S}, +\infty)$. If moreover the condition 2 in part c of Lemma 2 is verified for $a = 0$, then this condition is also verified for any $a > 0$ and $S \in (\lambda(a), \hat{S})$. Therefore, if $(1/f)'' > 0$ on $(0, \hat{S})$ then Assumption 2 is satisfied.

We will see in Section 5, how to use Remark 1 and Lemmas 2 and 3 to show that a linear growth function, a Monod function and a Hill function satisfy Assumptions 2 and 3.

## 4 Biogas flow rate

Microbial activity often produces by-products such as biogas, which can be a valuable source of energy in certain contexts. For instance, the anaerobic digestion of organic matter by microbial species produces methane and carbon dioxide. Valorizing biogas production while treating wastewater has received recently great attention, as a way of producing valuable energy and limiting the carbon footprint of the process [30].

We recall that the biogas flow rate is proportional to the microbial activity, as defined for instance in [3, 27]. We consider here the biogas flow rate as a function of the input substrate concentration $S^{in}$, the dilution rate $D$ and the parameter $r$.

For $r(f(S^{in}) - a) \le D < (1 - r)(f(S^{in}) - a)$, the biogas flow rate corresponding to the steady state $E_1$ is given by the expression

$$ G_1(S^{in}, D, r) = V_2 \bar{x}_2 f(\bar{S}_2), \qquad (36) $$

with $V_2 = (1 - r)V$, $\bar{x}_2$ and $\bar{S}_2$ defined in (B15).

For $D < r(f(S^{in}) - a)$, the biogas flow rate corresponding to the positive steady state $E_2$ is given by the expression

$$ G_2(S^{in}, D, r) = V_1 x_1^* f(S_1^*) + V_2 x_2^* f(S_2^*), \qquad (37) $$

with $V_1 = rV$, $V_2 = (1 - r)V$, $x_1^*$ and $S_1^*$ defined in (B17), $x_2^*$ defined by (B18) and $S_2^*$ the unique solution of $h(S_2) = f(S_2)$.
18 Performance study of two serial interconnected chemostats

**Proposition 4** 1. When $r(f(S^{in}) - a) \le D$ and $D < (1-r)(f(S^{in}) - a)$ then

$$G_1(S^{in}, D, r) = VD(S^{in} - \bar{S}_2). \qquad (38)$$

2. When $D < r(f(S^{in}) - a)$ then

$$G_2(S^{in}, D, r) = VD(S^{in} - S_2^*). \qquad (39)$$

*Proof* From system (B10), considering equation $\dot{S}_2 = 0$, one obtains $\bar{x}_2f(\bar{S}_2) = D(S^{in} - \bar{S}_2)/(1-r)$. Thus,

$$G_1(S^{in}, D, r) = V_1 \frac{D}{r} (S^{in} - \bar{S}_2) = VD(S^{in} - \bar{S}_2).$$

From system (B10), considering $\dot{S}_1 = 0$ and $\dot{S}_2 = 0$ gives respectively $x_1^*f(S_1^*) = D(S^{in} - S_1^*)/r$ and $x_2^*f(S_2^*) = D(S_1^* - S_2^*)/(1-r)$. Thus, one has

$$G_2(S^{in}, D, r) = V_1 \frac{D}{r} (S^{in} - S_1^*) + V_2 \frac{D}{1-r} (S_1^* - S_2^*) \\ = VD(S^{in} - S_2^*).$$

This ends the proof of the proposition.

□

410

411

412

Although $G_1(S^{in}, D, r)$ and $G_2(S^{in}, D, r)$, given by (36) and (37), respectively, are not defined for $r = 0$ or $r = 1$, the formulas (38) and (39) allow them to be extended to $r = 0$ and $r = 1$, as was done for $S_r^{out}$ in (6). We can write

$$G_1(S^{in}, D, 0) = G_2(S^{in}, D, 1) = G_{chem}(S^{in}, D),$$

where

$$G_{chem}(S^{in}, D) = VD(S^{in} - \lambda(D + a)), \qquad (40)$$

413

414

represents the biogas flow rate of the single chemostat when $0 < D < f(S^{in}) - a$. For more information on $G_{chem}(S^{in}, D)$, see (A7) in Appendix A.

<figure><img src="image_7.png" /><figcaption>**Fig. 5** For $r$ and $S^{in}$ fixed, the curves of the maps $D \mapsto G_1(S^{in}, D, r)$, in green, $D \mapsto G_2(S^{in}, D, r)$, in orange, and $D \mapsto G_{chem}(S^{in}, D)$, in black, where $G_1$, $G_2$ and $G_{chem}$ are given by (38), (39) and (40) respectively.</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>Fig. 6 The map $r \mapsto \bar{G}(r)$ with $\bar{G}$ defined by (42). (a) $f(S) = 4S$, $a = 0.6$ and $S^{in} = 1.5$. (b) $f(S) = 4S/(5+S)$, $a = 0.3$ and $S^{in} = 1.5$. (c) $f(S) = 4S^2/(25+S^2)$, $a = 0.3$ and $S^{in} = 10$.</figcaption></figure>

## 4.1 The serial configuration can be more efficient than the single chemostat

In this section, we prove that the biogas flow rate $G_1$ corresponding to the steady state $E_1$ is always smaller than the biogas flow rate of the single chemo- stat. However, the biogas flow rate $G_2$ corresponding to the steady state $E_2$ can be larger than the biogas flow rate of the single chemo- stat. More precisely, we have the following result.

**Proposition 5** Assume that Assumption 1 is satisfied. Let $r \in (0,1)$, $0 \le D < f(S^{in}) - a$ and $G_{chem}$ defined by (40).

1. If $r(f(S^{in}) - a) \le D$ and $D < (1-r)(f(S^{in}) - a)$, then $G_1(S^{in}, D, r) < G_{chem}(S^{in}, D)$, where $G_1$ is given by (38).

2. If $D < r(f(S^{in}) - a)$, then

$$G_2(S^{in}, D, r) > G_{chem}(S^{in}, D) \iff S^{in} > g_r(D),$$

where $G_2$ is given by (39) and $g_r$ is defined by (9).

* If, in addition, Assumption 2 is satisfied, and $S^{in} > g(D)$, then $G_2(S^{in}, D, r) > G_{chem}(S^{in}, D)$, if and only if $r > r_1(S^{in}, D)$, where $r_1(S^{in}, D)$ is defined by (20).

* If, in addition, Assumption 3 is satisfied, then $G_2(S^{in}, D, r) > G_{chem}(S^{in}, D)$, if and only if $D < D_r(S^{in})$, where $D_r(S^{in})$ is defined by (25).

**Proof 1.** Since $D/(1-r) > D$ and $\lambda$ is increasing, we have $\lambda(D/(1-r)+a) > \lambda(D+a)$. Then, using the formula for $G_1$ given in Proposition 4, this induces the inequality $G_1(S^{in}, D, r) < G_{chem}(S^{in}, D)$.
20 Performance study of two serial interconnected chemostats

2. According to Theorem 1, for any $r \in (0,1)$ and $D < r(f(S^{in}) - a)$ one has $S_2^*(S^{in}, D, r) < \lambda(D + a)$ if and only if $S^{in} > g_r(D)$. Consequently, using the formula for $G_2$ given in Proposition 4, one has $G_2(S^{in}, D, r) > G_{chem}(S^{in}, D)$ if and only if $S^{in} > g_r(D)$. If Assumption 2 is satisfied, then, using (21), we see that $G_2(S^{in}, D, r) > G_{chem}(S^{in}, D)$ if and only if $r > r_1(S^{in}, D)$. If Assumption 3 is satisfied, then, using (26), we see that $G_2(S^{in}, D, r) > G_{chem}(S^{in}, D)$ if and only if $D < D_r(S^{in})$.

This ends the proof of the proposition.

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

Let $S^{in}$ and $D$ be fixed. The graphs of the biogas flow rates functions

$$r \mapsto G_1(S^{in}, D, r), \text{ and } r \mapsto G_2(S^{in}, D, r),$$

are easily obtained from the graph of the output substrate concentration, $r \mapsto S_r^{out}(S^{in}, D)$, see Fig. 4. Indeed, the formulas given in Proposition 4 show that, whenever these functions are defined, we have

$$G_1(S^{in}, D, r) = VD(S^{in} - S_r^{out}(S^{in}, D)),$$

$$G_2(S^{in}, D, r) = VD(S^{in} - S_r^{out}(S^{in}, D)).$$

We will see in Section 5, some illustrative plots of the biogas flow rates $G_1$ and $G_2$ as functions of the parameter $r \in [0, 1]$, for linear growth, see Fig. 7, Monod growth, see Fig. 8 and Hill growth, see Fig. 9.

Let us illustrate the result of Proposition 5 by plotting the graphs of the biogas flow rates

$$D \mapsto G_1(S^{in}, D, r) \text{ and } D \mapsto G_2(S^{in}, D, r),$$

when $r$ and $S^{in}$ are fixed, see Fig. 5. This figure is made without graduations on the axes because it represents a general situation where the growth function is only assumed to verify our hypotheses. Indeed the behaviors of the functions, depicted in this figure, follow from our results and are not simply numerical illustrations.

Notice that for any $r \in (0,1)$, the graph of $G_1$ (plotted in green in the figure) is always below the graph $G_{chem}$ (plotted in black). This illustrates item 1 of Proposition 5. Assuming that Assumption 3 is satisfied, then for all $0 < D < D_r(S^{in})$, the graph of $G_2$ (plotted in orange) is above the graph of $G_{chem}$ (plotted in black). This illustrates item 2 of Proposition 5.

## 4.2 The maximal biogas of the serial configuration can exceed that of the single chemostat

In Figure 5(c) the plot shows that the maximum of $G_2$ (the red curve) is larger than the maximum of $G_{chem}$, as we want to emphasize that the following inequality is possible

$$\max_D G_2(S^{in}, D, r) > \max_D G_{chem}(S^{in}, D). \qquad (41)$$
Indeed we will show that there is a value $r^* \in (0, 1)$ such that this inequality is true for all $r \in (r^*, 1)$. The threshold $r^*$ obviously depends on $S^{in}$ and the rate of mortality $a$. It will be noted $r^*(S^{in}, a)$ when we want to highlight this dependence. This phenomenon never occurs in the case of no mortality, since we have $r^*(S^{in}, 0) = 1$. Indeed, in the case without mortality, we proved, see Proposition 6 of [7], that for all $S^{in} > 0$, and all $r \in (0, 1)$ we have

$$ \max_D G_2(S^{in}, D, r) < \max_D G_{chem}(S^{in}, D), $$

that is to say, the maximal biogas flow rate of the serial configuration never exceed the maximal biogas flow rate of the single chemostat.

Let us prove that, when $a > 0$, the inequality (41) is always true for $r$ sufficiently close to 1. Observe that for any fixed $S^{in} > \lambda(a)$ and $r \in (0, 1]$, the continuous function $D \mapsto G_2(S^{in}, D, r)$ is defined on the closed interval $[0, r(f(S^{in}) - a)]$. It is null at the extremities of this interval and positive on the open interval $(0, r(f(S^{in}) - a))$. Therefore, it reaches its maximum. For a given $S^{in} > \lambda(a)$, we then consider the function

$$ \bar{G}(r) := \max_{D \in [0, r(f(S^{in}) - a)]} G_2(S^{in}, D, r). \quad (42) $$

We want to ensure that this maximum is reached at a single value, denoted $\bar{D}(r)$. Note that $\bar{D}(1)$ represents the value, which we will assume to be unique, at which the function $D \mapsto G_{chem}(S^{in}, D)$ reaches its maximum. We need the following assumption.

**Assumption 4** The function $f$ is $C^2$ and increasing and, for $S^{in} > \lambda(a)$, there exists $\bar{D}(1) \in (0, f(S^{in}) - a)$ such that $D \mapsto G_{chem}(S^{in}, D)$ is

* strictly concave at $\bar{D}(1)$,

* increasing on $(0, \bar{D}(1))$,

* decreasing on $(\bar{D}(1), f(S^{in}) - a)$,

These conditions are related to the single chemostat model. They are verified for linear, Monod, or Hill growth functions, see Remark 3 in Appendix A.

If Assumption 4 is satisfied, then the maximum of the function $D \mapsto G_{chem}(S^{in}, D)$ is unique. The following lemma shows that the function $D \mapsto G_2(S^{in}, D, r)$ satisfies the same property for $r$ sufficiently close to 1.

**Lemma 4** Assume that Assumption 4 is satisfied, then for any $S^{in} > \lambda(a)$, there exists a neighborhood $V_1$ of 1, such that for any $r \in V_1 \cap \{r \le 1\}$, the maximum of the function $D \mapsto G_2(S^{in}, D, r)$ is unique. We denote it by $\bar{D}(r)$. Moreover, $\bar{D}$ is differentiable on $V_1 \cap \{r < 1\}$ with bounded derivative.
22 Performance study of two serial interconnected chemostats

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

**Proposition 6** Under Assumption 4, the function $\overline{G}$ admits left limits of its first and second derivatives at $r = 1$, which are

$$ \overline{G}'(1^{-}) = 0, \quad \overline{G}''(1^{-}) = \frac{2a\overline{D}(1)}{\overline{D}(1) + a}(S^{in} - \lambda(\overline{D}(1) + a)). \quad (43) $$

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

**Proposition 7** Under Assumption 4, there exists $r^*$ in $(0, 1)$ such that (41) is true for any $r \in (r^*, 1)$ and

$$ \max_D G_2(S^{in}, D, r^*) = \max_D G_{chem}(S^{in}, D). $$

Proof From Proposition 6, there exist $\varepsilon > 0$ such that for all $r \in (1 - \varepsilon, 1)$, we have $\overline{G}(r) > \overline{G}(1)$. Therefore, the subset $I$ of $(0, 1)$ defined by

$$ I = \{\rho \in (0, 1) : \forall r \in (\rho, 1), \overline{G}(r) > \overline{G}(1)\}, $$

494 is non empty. Let $r^*$ be the lower bound of $I$. We have $\overline{G}(r^*) = \overline{G}(1)$ and $\overline{G}(r) > \overline{G}(1)$

495 for $r \in (r^*, 1)$. Using (42), we deduce that the equality in the proposition is true and

496 (41) is true for any $r \in (r^*, 1)$. □

The function $r \mapsto \overline{G}(r)$ reaches its maximum at some $r^{max} \in (r^*, 1)$. Let $D^{max} = \overline{D}(r^{max})$ be the maximum of the function $D \mapsto G_2(S^{in}, D, r^{max})$. Therefore the maximal biogas flow rate of the serial chemostat is given by $G_2(S^{in}, D^{max}, r^{max})$. It satisfies

$$ G_2(S^{in}, D^{max}, r^{max}) > G_{chem}(S^{in}, \overline{D}(1)). $$

We have plotted the function $r \mapsto \overline{G}(r)$ for the linear, Monod, and Hill growth functions considered in Fig. 6. It is seen in this figure that the tangent at $r = 1$ is horizontal which corresponds to $\overline{G}'(1) = 0$. In addition, one remarks that $\overline{G}(r) > \overline{G}(1)$ for $r$ in some interval $(r^*, 1)$ and $\overline{G}(r^*) = \overline{G}(1)$. Thus, with presence of mortality rate, if practitioners are able to choose the dilution rate $D$, the good strategy consists in working with a serial configuration and choose $r$ in the interval $(r^*, 1)$. The serial configuration should be operated at $D = \overline{D}(r)$, where $\overline{D}(r)$ is defined in Lemma 4.

Remark 2 • If one is interested in increasing the flow of biogas, the best choice is $r = r^{max}$, $D = D^{max}$.

• If one is interested in reducing the dilution rate, the best choice is $r = r^*$ and $D = D^*$, where $D^* = \overline{D}(r^*)$.
Indeed, for the choice $r = r^*$ and $D = D^*$, we have

$$G_2(S^{in}, D^*, r^*) = G_{chem}(S^{in}, \bar{D}(1)),$$

but $D^*$ is expected to be significantly smaller than $\bar{D}(1)$, the dilution rate that maximises biogas for the simple chemostat. In fact, reducing $D$ means that the flow rate $Q$ has been reduced, and therefore energy has been saved to obtain the same result as with a simple chemostat

This result has an important message for practitioners: the serial configuration is worth considering when mortality is not negligible. To the best of our knowledge, this result is new in the literature. On the other hand, it is not intuitive. For more information on this issue, see Section 5.4. For biological comments on the heuristic underlying this non-intuitive behaviour, the reader is refered to [6].

# 5 Illustrations and numerical simulations

This section illustrates of results using three different growth functions. As concave functions, we choose the linear growth function and the Monod function. As a non concave function, we choose the Hill function.

<figure><img src="image_10.png" /><figcaption>Fig. 7 (a) The regions $J_0, J_1$ and $J_3$ of the operating plane with $f(S) = 4S$ and $a = 0.3$. The biogas flow rates corresponding to points $\eta_1 = (0.27, 0.6) \in J_1$ and $\eta_3 = (0.5, 0.6) \in J_3$ are depicted in panels (b) and (c) respectively. In these panels, the numbered curves (①) (in black), and ①, ② (in orange) are respectively defined by $y = G_{chem}(S^{in}, D)$, $y = G_1(S^{in}, D, r)$ and $y = G_2(S^{in}, D, r)$; $r_0(S^{in}, D) = D/(f(S^{in}) - a)$ and $r_1(S^{in}, D)$ is defined by (20). (b) $r_0 \approx 0.77$. (c) $r_0 \approx 0.35$ and $r_1 = 0.5$.</figcaption></figure>

## 5.1 Linear growth function

Let consider a linear function $f(S) = \alpha S$, $\alpha > 0$. As it is concave, according to item *ii* in Remark 1, the linear function verifies Assumptions 2 and 3. Therefore, our results apply for a linear function.

One has $\lambda(2D + a) = g(D) = (2D + a)/\alpha$ then, the curves $\Phi_{1/2}$, defined by (13), and $\Gamma$, defined by (22), are identical. Consequently, the operating plane
24 Performance study of two serial interconnected chemostats

<figure><img src="image_11.png" /><figcaption>**Fig. 8** (a) The regions $J_0, J_1, J_2$ and $J_3$ in the operating plane with $f(S) = 4S/(5+S)$ and $a = 0.3$. The biogas flow rates corresponding to points $\eta_1 = (3,0.7) \in J_1, \eta_2 = (3.45,0.7) \in J_2$ and $\eta_3 = (5,0.7) \in J_3$ are depicted in panels (b), (c) and (d) respectively. In these panels the curves are coloured and numbered as in Fig. 7, $r_0(S^{\text{in}}, D) = D/(f(S^{\text{in}})-a)$, and $r_1(S^{\text{in}}, D)$ is defined by (20) (b) $r_0 \approx 0.58$. (c) $r_0 \approx 0.53$ and $r_1 \approx 0.87$. (d) $r_0 \approx 0.41$ and $r_1 \approx 0.54$.</figcaption></figure>

$(S^{\text{in}}, D)$ is divided in three regions $J_i$, $i = 0, 1, 3$ defined in (23) that describe the behavior of the output substrate concentration and the biogas flow rate, see Figure 7(a).

Consider the operating points $\eta_1$ and $\eta_3$, fixed respectively in regions $J_1$ and $J_3$, as shown in Figure 7(a). The behavior of the biogas flow rate for these operating points is depicted in Figure 7(b,c). It should be noticed that for any other point $(S^{\text{in}}, D) \in J_1$, the curve representing the biogas flow rate with respect to $r$ should be similar to the curve shown in Figure 7(a), and corresponding to $(S^{\text{in}}, D) = \eta_1$. Similarly, for any other point $(S^{\text{in}}, D) \in J_3$, it should be similar to the curve shown in Figure 7(b), and corresponding to $(S^{\text{in}}, D) = \eta_3$.

In the linear case, the equation $S^{\text{in}} = g_r(D)$ is a second degree algebraic equation in $r$ that gives two solutions, one corresponds to $r_1(S^{\text{in}}, D)$ defined by (20) and the other one is not considered as it does not belong to $(0, 1)$.

Since the point $\eta_3 = (0.5, 0.6)$ satisfies the condition $S^{\text{in}} > g(D)$, as stated in item 2 of Proposition 5, the serial configuration has a higher biogas flow rate production than a single chemostat if and only if $r \in (r_1, 1)$, where $r_1(0.5, 0.6) \approx 0.5$, see Figure 7 (b).

## 5.2 Monod function

The Monod function is $f(S) = mS/(K+S)$. As it is concave, according to item ii in Remark 1, the Monod function verifies Assumptions 2 and 3. Therefore, our results apply for Monod function.

**Lemma 5** For any $D > 0$, the curve $\Gamma$, defined by (22), is at left of the curve $\Phi_{1/2}$, defined by (13).

*Proof* The curves $\Phi_{1/2}$ and $\Gamma$ are respectively defined by equations $S^{\text{in}} = \lambda(2D+a)$ and $S^{\text{in}} = g(D)$. Let the function $H: [0, (m-a)/2) \mapsto \mathbb{R}$ be defined by
$$H(D) = \lambda(2D + a) - g(D) = \frac{KmD^2}{(m-D-a)^2(m-a-2D)}.$$

Note that $H(0) = 0$ and, for any $D \in (0, (m-a)/2)$, one has $H(D) > 0$ i.e. $\lambda(2D + a) > g(D)$.

Hence, the curve $\Gamma$ is at left of the curve $\Phi_{1/2}$.

□

As a consequence of Lemma 5, the operating plane $(S^{in}, D)$ is divided in four regions $J_i$, $i = 0, 1, 2, 3$ defined in (23) that describe the behavior of the output substrate concentration and the biogas flow rate, see Fig. 8(a).

Consider the operating points $\eta_1, \eta_2$ and $\eta_3$, fixed respectively in regions $J_1, J_2$ and $J_3$, as shown in Fig. 8(a). The behavior of the biogas flow rate for these points is depicted in Fig. 8(b,c,d). It should be noticed that for any other point $(S^{in}, D) \in J_1$ (resp. $(S^{in}, D) \in J_2$ and $(S^{in}, D) \in J_3$), the curve representing the biogas flow rate with respect to $r$ should be similar to the curve shown in Fig. 8(b) (resp. 8(c) and 8(d)), and corresponding to $(S^{in}, D) = \eta_1$ (resp. $(S^{in}, D) = \eta_2$ and $(S^{in}, D) = \eta_3$).

In the Monod case, the equation $S^{in} = g_r(D)$ is a second degree algebraic equation in $r$ that gives two solutions, one corresponds to $r_1(S^{in}, D)$ defined by (20) and the other one is not considered as it does not belong to $(0, 1)$.

Since the point $\eta_2$ (resp. $\eta_3$) satisfies the condition $S^{in} > g(D)$, as stated in item 2 of Proposition 5, the serial configuration has a higher biogas flow rate production than a single chemostat if and only if $r \in (r_1, 1)$, with $r_1(3.45, 0.7) \approx 0.87$ in Fig. 8(c) and $r_1(5, 0.7) \approx 0.54$ in Fig. 8(c).

## 5.3 Hill function

The Hill function is $f(S) = mS^p/(K^p + S^p)$. Note that if $p = 1$ this function reduces to the Monod function. For $p > 1$ it is non-concave. We have

$$\lambda(a) = \left(\frac{a}{m-a}\right)^{1/p} K.$$

**Proposition 8** The Hill function satisfies the conditions (34) and (35). Therefore, according to item iii in Remark 1, it verifies Assumption 2 and according to Lemma 3, it satisfies Assumption 3.

**Proof** Let us first prove that the Hill function satisfies the condition (35). Straight-forward computation give

$$\frac{d^2}{dS^2} \left( \frac{1}{f(S)-a} \right) = mpK^p \frac{(p+1)(m-a)S^{2p-2} + (p-1)aK^p S^{p-2}}{((m-a)S^p - aK^p)^3}.$$

Therefore, $\frac{d^2}{dS^2} \left( \frac{1}{f(S)-a} \right) > 0$ for all $S > \lambda(a)$, that is to say, (35) is satisfied. This result can also be obtained without laborious calculations by using item iii of Remark 1. Let $\hat{S} \in (0, +\infty)$ be the inflexion point of the Hill function $f$. It is sufficient to show that $(1/f)'' > 0$ for all $S \in (0, \hat{S})$. One easily see that
26 Performance study of two serial interconnected chemostats

<figure><img src="image_12.png" /><figcaption>**Fig. 9** (a) The regions $J_0, J_1, J_2, J_3$ and $J_4$ in the operating plane with $f(S) = 4S^2/(25+S^2)$ and $a = 0.1$. The curves $\Gamma$ and $\Phi_{1/2}$ intersects for $D_1 = 0.69$ (see Lemma 6). (b,c) Zooms of (a) showing the region $J_4$. The biogas flow rates corresponding to points $\eta_1 = (7, 1.6) \in J_1$, $\eta_2 = (9, 1.6) \in J_2$, $\eta_3 = (12, 1.6) \in J_3$ and $\eta_4 = (2.33, 0.3) \in J_4$ are depicted in panels (d) to (g), respectively. In these panels curves are coloured and numbered as in Fig. 7, $r_0(S^{in}, D) = D/(f(S^{in}) - a)$, and $r_1(S^{in}, D)$ is defined by (20). (d) $r_0 \approx 0.63$. (e) $r_0 \approx 0.54$ and $r_1 \approx 0.81$. (f) $r_0 \approx 0.48$ and $r_1 \approx 0.61$. (g) $r_0 \approx 0.49$.</figcaption></figure>

$$ \left(\frac{1}{J}\right)''(S) = \frac{p(p+1)K^p}{mS^{p+2}} > 0, $$

for any $S > 0$. Consequently, for all $p > 1$, the Hill function verifies Assumption 2.

587

Let us now prove that the Hill function verifies the condition (34). Straightforward

588

computations give

589

$$ f'(\lambda(D+a)) = \frac{p}{Km}(D+a)^{\frac{p-1}{p}}(m-a-D)^{\frac{p+1}{p}}. \quad (44) $$

Therefore,

590
$$f'(\lambda \left(\frac{D}{r} + a\right)) = \frac{p}{Km} \left(\frac{D}{r} + a\right)^{\frac{p-1}{p}} \left(m - a - \frac{D}{r}\right)^{\frac{p+1}{p}}.$$

Since $p > 1$, $D + ra < D + a$ and

$$0 < rm - ra - D < m - a - D,$$

one has

$$(D + ra)^{\frac{p-1}{p}} < (D + a)^{\frac{p-1}{p}} \qquad (45)$$

$$(rm - ra - D)^{\frac{1}{p}} < (m - a - D)^{\frac{1}{p}}. \qquad (46)$$

From (45) one has

$$\left(\frac{D}{r} + a\right)^{\frac{p-1}{p}} = \left(\frac{1}{r}\right)^{\frac{p-1}{p}} (D + ra)^{\frac{p-1}{p}} < \left(\frac{1}{r}\right)^{\frac{p-1}{p}} (D + a)^{\frac{p-1}{p}}. \qquad (47)$$

On the other hand, we have

$$(m - a - \frac{D}{r})^{\frac{p+1}{p}} = \left(\frac{1}{r}\right)^{\frac{1}{p}} (m - a - \frac{D}{r}) A,$$

where $A = (rm - ra - D)^{\frac{1}{p}}$. From (46), and using

$$0 < m - a - D/r < m - a - D,$$

we then deduce

$$(m - a - \frac{D}{r})^{\frac{p+1}{p}} < \left(\frac{1}{r}\right)^{\frac{1}{p}} (m - a - D)^{\frac{p+1}{p}}. \qquad (48)$$

Therefore, using (44), (47) and (48) one obtains

$$f'\left(\lambda\left(\frac{D}{r} + a\right)\right) = \frac{p}{Km} \left(\frac{D}{r} + a\right)^{\frac{p-1}{p}} \left(m - a - \frac{D}{r}\right)^{\frac{p+1}{p}} \\ < \frac{p}{Km} \frac{1}{r} (D + a)^{\frac{p-1}{p}} (m - a - D)^{\frac{p+1}{p}} = \frac{1}{r} f'\left(\lambda(D + a)\right).$$

This ends the proof of (34). Consequently, according to Lemma 3, any Hill function

satisfies Assumption 3. □

Let now consider the case $p = 2$ of the Hill function: $f(S) = mS^2/(K^2+S^2)$.

**Lemma 6** Let $D_1 = (3m - 4a - \sqrt{m(5m - 4a)})/4$. If $0 < D < D_1$ then the curve

$\Phi_{1/2}$, defined by (13), at left of the curve $\Gamma$, defined by (22). In contrast, if $D_1 < D < (m - a)/2$ then the curve $\Phi_{1/2}$ is at right of the curve $\Gamma$.

**Proof** Let the function $H: [0, (m - a)/2) \to \mathbb{R}$ be defined by $H(D) := \lambda(2D + a) - g(D)$. We have

$$H(D) = K \left( \sqrt{\frac{2D+a}{m-a-D}} - \frac{(2D+a)(m-a-D) + (D+a)(m-a)}{2(m-a-D)^{3/2}\sqrt{D+a}} \right).$$

Straightforward computation shows that this function is positive if and only if the

polynomial
28 Performance study of two serial interconnected chemostats

<figure><img src="image_13.png" /><figcaption>Fig. 10 (a) The map $r \mapsto \bar{G}(r)$ defined by (42), with $f(S) = 4S/(5+S)$, $a = 0.1$ and $S^{in} = 1.5$, showing the values $r^*$ and $r^{max}$. (b) The corresponding maps $D \mapsto G_{chem}(S^{in}, D)$, in black, $D \mapsto G_2(S^{in}, D, r^*)$, in blue and $D \mapsto G_2(S^{in}, D, r^{max})$, in red, show the values $D^* < D^{max} < D_{chem}^{max}$. (c) The biogas flow rates for $a = 0, 0.1, 0.2, 0.3$ showing the effects of mortality.</figcaption></figure>

$$Q(D) := 4D^2 - 2(3m - 4a)D + 4a^2 - 5am + m^2$$

608 is negative. The solution of equation $Q(D) = 0$ are

$$D_1 = \frac{3m-4a-\sqrt{\Delta}}{4} \text{ and } D_2 = \frac{3m-4a+\sqrt{\Delta}}{4},$$

609 where $\Delta = 4m(5m - 4a) > 0$, as $a < m$. Notice that we have $0 < D_1 < (m - a)/2$ and $(m - a)/2 < D_2$. Thus, for any $D \in (D_1, (m - a)/2)$, we have $H(D) > 0$ and then the curve $\Phi_{1/2}$ at right of the curve $\Gamma$. $\square$

612 As a consequence of Lemma 6, the operating plane is divided in five regions $J_i$ $i = 0, 1, 2, 3, 4$ defined in (23), see Figure 9(a,b,c).

613 Consider the operating points $\eta_1, \eta_2, \eta_3$ and $\eta_4$ fixed respectively in regions $J_1, J_2, J_3$ and $J_4$, as shown in Figure 9(a,b,c). It should be noticed that for any other point $(S^{in}, D) \in J_1$ (resp. $(S^{in}, D) \in J_2, (S^{in}, D) \in J_3$ and $(S^{in}, D) \in J_4$), the curve representing the biogas flow rate with respect to $r$ should be similar to the curve shown in Fig. 9(a) (resp. (b), (c) and (d)), and corresponding to $(S^{in}, D) = \eta_1$ (resp. $(S^{in}, D) = \eta_2, (S^{in}, D) = \eta_3$ and $(S^{in}, D) = \eta_4)$. 616 615 614

620 619 621 Recall that $r_1(S^{in}, D)$ is defined by (20). It is obtained by solving numer- ically the equation $S^{in} = g_r(D)$. Since the point $\eta_2$ (resp. $\eta_3$) satisfies the condition $S^{in} > g(D)$, as stated in item 2 of Proposition 5, the serial configu- ration has a higher biogas flow rate production than a single chemostat if and only if $r \in (r_1, 1)$, with $r_1(9, 1.6) \approx 0.81$ in Fig. 9(e) and $r_1(12, 1.6) \approx 0.61$, in Fig. 9(f). 623 622 624 625 626

627 628 5.4 The serial configuration is worth considering when mortality is not negligible

629 In this section we numerically illustrate Remark 2. We fix $S^{in}$ and we adopt the following notations. 630

$$D_{chem}^{max} = \bar{D}(1), \quad G_{chem}^{max} = \bar{G}(1) = G_{chem}(S^{in}, D_{chem}^{max})$$
where $\bar{G}(r)$ is defined by (42) and $\bar{D}(r)$ is as in Lemma 4. Recall that $r^* \in (0, 1)$ satisfies

$$ \bar{G}(r^*) = \bar{G}(1) = G_{chem}^{max}, \qquad (49) $$

and $\bar{G}(r) > \bar{G}(1)$ for $r \in (r^*, 1)$, so that $\bar{G}(r)$ attains its maximum for $r = r^{max} \in (r^*, 1)$, see Fig. 10(a), obtained with a Monod function and $S^{in} = 1.5$. We adopt the following notations.

$$ \begin{aligned} D^{max} &= \bar{D}(r^{max}), & G^{max} &= G_2(S^{in}, D^{max}, r^{max}) \\ D^* &= \bar{D}(r^*), & G^* &= G_2(S^{in}, D^*, r^*) = G_{chem}^{max} \end{aligned} $$

<table><caption>**Table 1** Numerical values</caption><thead><tr><th></th><th>a = 0</th><th>a = 0.1</th><th>a = 0.2</th><th>a = 0.3</th></tr></thead><tbody><tr><td>D<sup>max</sup><sub>chem</sub></td><td>0.4918</td><td>0.4359</td><td>0.3806</td><td>0.3259</td></tr><tr><td>G* = G<sup>max</sup><sub>chem</sub></td><td>0.3930</td><td>0.3167</td><td>0.2478</td><td>0.1866</td></tr><tr><td>r*</td><td>1</td><td>0.839</td><td>0.717</td><td>0.631</td></tr><tr><td>D*</td><td>0.4918</td><td>0.3758</td><td>0.2969</td><td>0.2369</td></tr><tr><td>r<sup>max</sup></td><td>1</td><td>0.889</td><td>0.808</td><td>0.751</td></tr><tr><td>D<sup>max</sup></td><td>0.4918</td><td>0.3925</td><td>0.3190</td><td>0.2591</td></tr><tr><td>G<sup>max</sup></td><td>0.3930</td><td>0.3169</td><td>0.2490</td><td>0.1890</td></tr><tr><td>G<sup>max</sup> - G<sup>max</sup><sub>chem</sub></td><td>0</td><td>0.06%</td><td>0.5%</td><td>1.3%</td></tr><tr><td>D<sup>max</sup><sub>chem</sub> - D<sup>max</sup></td><td>0</td><td>10%</td><td>16.2%</td><td>20.5%</td></tr><tr><td>D<sup>max</sup><sub>chem</sub> - D*</td><td>0</td><td>13.6%</td><td>22%</td><td>27.3%</td></tr></tbody></table>

These notations are illustrated in Figs. 10(a,b). The zoom in Fig. 10(b) shows that $G^{max}$ exceeds $G^* = G_{chem}^{max}$ only slightly, but $D^*$ is significantly smaller than $D^{max}$, which is itself smaller than $D_{chem}^{max}$. We give in Table 1 the numerical values of $r^*$, $r^{max}$, $D^*$, $D^{max}$, $G^{max}$ and $G^* = G_{chem}^{max}$, for various values of the mortality rate $a$. The table also gives the relative gains

$$ \frac{G^{max} - G_{chem}^{max}}{G_{chem}^{max}}, \quad \frac{D_{chem}^{max} - D^{max}}{D_{chem}^{max}}, \quad \frac{D_{chem}^{max} - D^{*}}{D_{chem}^{max}}, $$

when replacing the single chemostat with the serial device using the ratios $r^*$ and $r^{max}$. The gain in biogas production is almost negligible, but the gain in bioreactor flow rate is significant.

The biogas flow rates $G_{chem}(S^{in}, D)$, $G_2(S^{in}, D, r^*)$ and $G_2(S^{in}, D, r^{max})$, for the various considered values of the mortality rate $a$, are depicted in Fig. 10(c), in black, blue and red, respectively. This figure shows that mortality is a real problem as it considerably reduces biogas production. Where mortality cannot be avoided or reduced, instead of using the single chemostat, by using a serial device, biogas production can be slightly improved while significantly reducing the bioreactor flow rate.
30 Performance study of two serial interconnected chemostats

651

# 6 Conclusion

652

653

654

655

656

657

658

659

660

661

662

663

664

665

666

667

668

669

670

671

672

673

674

675

In this work, an in-depth study is carried out on the mathematical model of two interconnected chemostats in serial with mortality. Equations contain a term representing the mortality rate of the species. Due to this added term characterizing the mathematical model, this paper is considered as an extension of the work done in [7], where the model does not consider the mortality rate. However, the mathematical analysis revealed that the proofs have had to be significantly revisited and reveal several new non intuitive differences compared to the case without mortality. Let us recall that without mortality, the dynamics admits a forward attractive invariant hyperplane related to the total mass conservation, which is no longer verified under mortality consideration. This at the core of the differences in the mathematical analysis. The study of the model is based on the analysis of the asymptotic behavior of its solutions, and is supported by an operating diagram which describes the number and stability of steady states. In a first step, we considered different mortality rates $a_1, a_2$ in each tank. Then, in view of comparing with the single configuration, we considered identical mortality rate $a = a_1 = a_2$. We analyzed the performances of the model at steady state for two different criteria: the output substrate concentration and the biogas flow rate (and compared them for the single chemostat with the same mortality rate $a$). Explicit expressions of criteria, depending on the dilution rate $D$ and the input substrate concentration $S^{in}$, are provided. These new results provide conditions that insure the existence of a serial configuration more efficient than a single chemostat, in the sense of minimizing the output substrate concentration or maximizing the biogas flow rate.

676

677

678

679

680

681

682

683

684

685

686

687

688

689

690

Along the paper, the similarities, specificities and differences of our model compared to the model without mortality (i.e. for $a = 0$) studied in [7] are highlighted. Among the differences that attract attention, on the one hand, we have the operating diagram with different mortality which presents many more cases than the diagram without mortality where it is reduced to only two cases. Thus, the presence of the four regions of stability on the same diagram is now possible. On the other hand, we have the biogas production of the serial device in its maximum state which can be significantly larger than the largest biogas production of the single chemostat. This never happens in the case without mortality. Finally, unlike the case without mortality, the biomass productivity and the biogas flow rate at steady state are not given by the same formulas. Therefore, if biomass productivity is taken into account as a performance criterion, the comparison between the serial chemostat and the single chemostat does not lead to the same conclusions. For more details on this issue the reader can refer to [8].

691

# Appendix A The single chemostat

692

693

In this section, we give a brief presentation of the mathematical model of the single chemostat with mortality rate. The mathematical equations are given by
<figure><img src="image_14.png" /><figcaption>Fig. A1 (a) The map $D \mapsto S^{\text{out}}(S^{\text{in}}, D)$ is increasing on $[0, \delta]$, where $\delta = f(S^{\text{in}}) - a$. (b) The map $D \mapsto x^{\text{out}}(S^{\text{in}}, D)$ with $f(S) = 4S/(5+S)$, $S^{\text{in}} = 10$ and $a = 0.6$. (c) The curve $\Gamma$ in the operating plane $(S^{\text{in}}, D)$ of the single chemostat.</figcaption></figure>

$$ \begin{aligned} \dot{S} &= D(S^{\text{in}} - S) - f(S)x, \\ \dot{x} &= -Dx + f(S)x - ax, \end{aligned} \qquad (A1) $$

where $S$ and $x$ denote respectively the substrate and the biomass concentration, $S^{\text{in}}$ the input substrate concentration, $a$ the mortality rate and $D = Q/V$ the dilution rate, with $Q$ the input flow rate and $V$ the volume of the tank. The specific growth rate $f$ of the microorganisms satisfies Assumption 1. It is well known (see [16, 39]) that, besides the washout $F_0 = (S^{\text{in}}, 0)$, this system can have a positive steady state

$$ F_1 = (S^*(D), x^*(S^{\text{in}}, D)), $$

where

$$ S^* = \lambda(D + a) \quad \text{and} \quad x^* = \frac{D}{D+a}(S^{\text{in}} - \lambda(D + a)). $$

See Fig. A1(a) for the plot of the function $D \mapsto S^*(D)$ and Fig. A1(b) for the plot of the function $D \mapsto x^*(S^{\text{in}}, D)$ for $0 \le D \le \delta$, where $\delta = f(S^{\text{in}}) - a$.

The washout steady state $F_0$ always exists. It is GAS if and only if $D \ge \delta$. It is LES if and only if $D > \delta$. The positive steady state $F_1$ exists if and only if $D < \delta$. It is GAS and LES whenever it exists. Therefore, the curve $\Phi$ defined by

$$ \Phi := \{(S^{\text{in}}, D) : D = f(S^{\text{in}}) - a\} \qquad (A2) $$

splits the set of operating parameters $(S^{\text{in}}, D)$ into two regions, denoted $I_0$ and $I_1$, as depicted in A1(c). These regions are defined by

$$ \begin{aligned} I_0 &:= \{(S^{\text{in}}, D) : D \ge f(S^{\text{in}}) - a\}, \\ I_1 &:= \{(S^{\text{in}}, D) : D < f(S^{\text{in}}) - a\}. \end{aligned} \qquad (A3) $$

The behavior of the system in each region is given in Table A1. Figure A1(c), together with A1 is called the operating diagram of the single chemostat.

The particularity of this operating diagram is that the curve limiting both regions $I_0$ and $I_1$ is translated from zero, unlike the case with mortality, as shown in Figure 2.5 of [16]. Thus, with presence of mortality rate, the region where the washout is GAS, is larger.
32 Performance study of two serial interconnected chemostats

<table><caption>**Table A1** Stability of steady states in the various regions of the operating diagram.</caption><thead><tr><th></th><th><i>I</i><sub>0</sub></th><th><i>I</i><sub>1</sub></th></tr></thead><tbody><tr><td><i>F</i><sub>0</sub></td><td>GAS</td><td>U</td></tr><tr><td><i>F</i><sub>1</sub></td><td></td><td>GAS</td></tr></tbody></table>

The output substrate concentration of the single chemostat, at its stable steady state is given by

$$S^{\text{out}}(S^{\text{in}}, D) = \begin{cases} S^{\text{in}} & \text{if } D \ge \delta \\ \lambda(D + a) & \text{if } D < \delta. \end{cases} \qquad (\text{A4})$$

Its output biomass concentration at steady state is then given by

$$x^{\text{out}}(S^{\text{in}}, D) = \frac{D}{D+a}(S^{\text{in}} - S^{\text{out}}(D, a)) \qquad (\text{A5})$$

For all $S^{\text{in}} > \lambda(a)$, one has

$$\frac{\partial S^{\text{out}}}{\partial D}(S^{\text{in}}, D) = \begin{cases} 0 & \text{if } D > \delta \\ \lambda'(D + a) & \text{if } D < \delta, \end{cases}$$

Thus, for any growth function satisfying Assumption 1 the function $D \mapsto S^{\text{out}}(S^{\text{in}}, D)$ is increasing on $[0, \delta]$, as shown in Figure A1(a). The function $D \mapsto x^{\text{out}}(S^{\text{in}}, D)$ is illustrated in Figure A1(b) for a Monod function.

The biogas flow rate of the single chemostat is defined, up to a multiplicative yield coefficient, by

$$G_{\text{chem}}(S^{\text{in}}, D) := V x^{\text{out}} f(S^{\text{out}}). \qquad (\text{A6})$$

Using the expressions (A4) and (A5) respectively of $S^{\text{out}}$ and $x^{\text{out}}$, the biogas flow rate of the single chemostat is given by:

$$G_{\text{chem}}(S^{\text{in}}, D) = \begin{cases} 0 & \text{if } D \ge \delta \\ VD(S^{\text{in}} - \lambda(D + a)) & \text{if } D < \delta. \end{cases} \qquad (\text{A7})$$

For a given $S^{\text{in}} > \lambda(a)$, the function $D \mapsto G_{\text{chem}}(S^{\text{in}}, D)$ is null for $D = 0$ or $D \ge \delta$, and is positive for $D \in (0, \delta)$. Therefore it admits a maximum in $(0, \delta)$, which is assumed to be unique. A characterization of the growth functions for which this uniqueness is satisfied can be found in [31].

**Proposition 9** Assume that for any $S^{\text{in}} > \lambda(a)$, the maximum of $D \mapsto G_{\text{chem}}(S^{\text{in}}, D)$ is unique, and define $\bar{D}(S^{\text{in}}) \in (0, \delta)$, such that

$$G_{\text{chem}}(S^{\text{in}}, \bar{D}(S^{\text{in}})) = \max_{D \ge 0} G_{\text{chem}}(S^{\text{in}}, D).$$

Then, the dilution rate $D = \bar{D}(S^{\text{in}})$ is the solution of the equation $S^{\text{in}} = g(D)$, where the function $g: [0, m - a) \mapsto \mathbb{R}$ is given by

$$g(D) := \lambda(D + a) + D\lambda'(D + a). \qquad (\text{A8})$$
735 Proof For any $S^{in} > \lambda(a)$ and $D \in (0, \delta)$, we have

$$ \frac{\partial G_{chem}}{\partial D}(S^{in}, D) = V(S^{in} - \lambda(D + a) - D\lambda'(D + a)) \quad (A9) $$

736 Therefore, $\frac{\partial G_{chem}}{\partial D}(S^{in}, D) = 0$ if and only if $S^{in} = g(D)$, where $g$ is defined by (A8). $\square$ 737

738 Notice that the function $g$ defined by (A8) is the same as the function 739 $g$, defined by (19), which was obtained as the limit, when $r$ tends to 1, to 740 the function $g_r$, defined by (9). Recall that $\Gamma$ is the curve of equation $S^{in} = 741 g(D)$, see (22). This curve is depicted in Fig. A1(c). It is the set of operating 742 conditions given the higher biogas of the single chemostat. More precisely, for 743 any $S^{in} > \lambda(a)$, the maximum $D = \bar{D}(S^{in})$ of the biogas satisfies the condition 744 $(S^{in}, D) \in \Gamma$. Therefore, a sufficient condition for the uniqueness of $\bar{D}(S^{in})$ is 745 that the mapping $g$ is increasing. If, in addition, $f$ is $C^2$, then, deriving (A9) 746 with respect of $D$, we have

$$ \frac{\partial^2 G_{chem}}{\partial D^2}(S^{in}, D) = -Vg'(D). $$

747 Hence, a sufficient condition for Assumption 4 to be satisfied is that $g'(D) > 748 0$ for $D \in [0, m - a)$. This last condition if satisfied whenever $f'' \le 0$ on 749 $(\lambda(a), +\infty)$, or $(\frac{1}{f-a})'' > 0$ on $(\lambda(a), +\infty)$, see Lemma 1 in [31]. Therefore we 750 can make the following remark.

751 Remark 3 Linear and Monod growth functions satisfy Assumption 4, since they 752 satisfy $f'' \le 0$ on $(0, +\infty)$. On the other hand the Hill function satisfy Assumption 753 4, since it satisfies $(\frac{1}{f-a})'' > 0$ on $(\lambda(a), +\infty)$, as shown in Proposition 8.

754 Appendix B The serial configuration

755 We consider a slight extension of system (1) with different mortality rates in 756 the two tanks. Indeed, we assume that the growth environment differs from 757 one tank to another one. This can lead to two different mortality rates in the 758 tanks. We denote by $a_1$ and $a_2$ the mortality rates. The mathematical model 759 is given by the following equations.

$$ \begin{aligned} \dot{S}_1 &= \frac{D}{r}(S^{in} - S_1) - f(S_1)x_1 \\ \dot{x}_1 &= -\frac{D}{r}x_1 + f(S_1)x_1 - a_1x_1 \\ \dot{S}_2 &= \frac{D}{1-r}(S_1 - S_2) - f(S_2)x_2 \\ \dot{x}_2 &= \frac{D}{1-r}(x_1 - x_2) + f(S_2)x_2 - a_2x_2. \end{aligned} \quad (B10) $$

760 The following result is classical in the mathematical theory of the chemo- 761 stat.
34 Performance study of two serial interconnected chemostats

762 Lemma 7 For any nonnegative initial condition, the solution of system (B10)

763 $(S_1(t), x_1(t), S_2(t), x_2(t))$ is nonnegative for any $t > 0$ and positively bounded.

764 Proof Since the vector field defined by (B10) is $C^1$, the uniqueness of the solution

765 to an initial value problem holds. From (B10) and using $f(0) = 0$, we have:

$$ \text{for } i = 1, 2, \quad S_i = 0 \implies \dot{S}_i > 0, $$

$$ x_1 = 0 \implies \dot{x}_1 = 0 $$

$$ x_1 \ge 0 \text{ and } x_2 = 0 \implies \dot{x}_2 \ge 0 $$

766 Therefore, for $i = 1, 2$, $S_i(t) \ge 0$ and $x_i(t) \ge 0$, for all $t \ge 0$, for which they are

767 defined, provided $S_i(0) \ge 0$ and $x_i(0) \ge 0$, for $i = 1, 2$, see Prop. B.7 in [39]. This

768 proves that the solutions of nonnegative initial conditions are always nonnegative.

769 Let $z_i = S_i + x_i$, $i = 1, 2$. From system (B10), we have

$$ \dot{z}_1 = \frac{D}{r}(S^{\text{in}} - z_1) - a_1 x_1, \quad \dot{z}_2 = \frac{D}{1-r}(z_1 - z_2) - a_2 x_2. $$

770 Consequently, we have the differential inequality

$$ \dot{z}_1 \le \frac{D}{r}(S^{\text{in}} - z_1), $$

771 It follows by comparison of solutions of ordinary differential equations (see for

772 instance [48]) that one has

$$ z_1(t) \le S^{\text{in}} + (z_1(0) - S^{\text{in}})e^{-\frac{D}{r}t} $$

773 Therefore, $z_1(t) \le Z_1$, where $Z_1 = \max(S^{\text{in}}, z_1(0))$. Then, we also have the

774 differential inequality

$$ \dot{z}_2 \le \frac{D}{1-r}(Z_1 - z_2). $$

775 It follows again by comparison of solutions of ordinary differential equations that one

776 has also

$$ z_2(t) \le Z_1 + (z_2(0) - Z_1)e^{-\frac{D}{1-r}t} $$

777 Therefore, $z_2(t) \le Z_2$, where $Z_2 = \max(Z_1, z_2(0))$. Hence, the solutions of (B10) are

778 positively bounded. Therefore, they are defined for all $t \ge 0$. $\square$

779 For the description of the steady states, we shall consider the following
780 function $h$ that will play a key role

$$ \begin{gathered} h(S_2) = \frac{D+(1-r)a_2}{1-r} \frac{S_2^*-S_2}{b-S_2}, \quad S_2 \in (0,b) \tag{B11} \\ \text{where } S_1^* = \lambda \left( \frac{D}{r} + a_1 \right), \quad b = \frac{D(S_1^{\text{in}}-S_1^*)}{D+ra_1} + S_1^*. \end{gathered} $$

781 This function satisfies the following property.

782 Lemma 8 Assume that $D/r + a_1 < f(S^{\text{in}})$. The function $h$ is decreasing from

783 $h(0) > 0$ to $h(S_1^*) = 0$, where $h(0)$ is given by

$$ h(0) = \frac{D+(1-r)a_2}{1-r} \frac{(D+ra_1)S_1^*}{DS^{\text{in}}+ra_1S_1^*} \tag{B12} $$
784 Proof From the condition $D/r + a_1 < f(S^{in})$ it is deduced that $S_1^* < S^{in}$. Note that

$$b = \frac{DS^{in} + ra_1S_1^*}{D + ra_1}.$$

785 Hence, b is a convex combination of $S^{in}$ and $S_1^*$, and we have $S_1^* < b < S^{in}$.

786 Therefore, the vertical asymptote $S_2 = b$ of $h$ is at right of $S_1^*$. The derivative of $h$ is

$$h'(S_2) = \frac{D+(1-r)a_2}{1-r} \frac{S_1^*-b}{(b-S_2)^2}.$$

787 Hence, we have $h'(S_2) < 0$ for all $S_2 < b$. Therefore, $h$ is defined on the interval

788 $(0, S_1^*)$ and is decreasing from $h(0)$, given by (B12) to $h(S_1^*) = 0$. □

789 Therefore, if $D/r + a_1 < f(S^{in})$, equation $f(S_2) = h(S_2)$ admits a unique

790 solution, denoted by $S_2^*(S^{in}, D, r)$, as shown in Fig. B2(a). This solution satisfy

791 the following property.

792 **Lemma 9** Considering $a_1 = a_2 = a$, for all $0 \le D < f(S^{in}) - a$, one has

$$\lim_{r \to 1} S_2^*(S^{in}, D, r) = \lambda(D + a).$$

793

794 Proof Let $0 \le D < f(S^{in}) - a$. Using (5), the condition $h(S_2) = f(S_2)$ is equivalent

795 to

$$(D + (1 - r)a)(S_1^* - S_2^*) = (1 - r) \left( \frac{D}{D + ra} (S^{in} - S_1^*) + S_1^* - S_2^* \right) f(S_2^*). \quad (B13)$$

796 For $r = 1$, we have $S_1^* = \lambda(D + a)$. As $\lim_{r \to 1} f(S_2^*) < +\infty$ then, (B13) gives

$$D(\lambda(D + a) - \lim_{r \to 1} S_2^*(S^{in}, D, r)) = 0.$$

797 Consequently, one has $\lim_{r \to 1} S_2^*(S^{in}, D, r)) = \lambda(D + a)$. □

<figure><img src="image_15.png" /><figcaption>**Fig. B2** (a) Existence and uniqueness of the solution $S_2^*$ of equation $f(S_2) = h(S_2)$. (b) Graphical illustration of Proposition 10: $S_2^*$ decreases when $S^{in}$ increases.</figcaption></figure>

798 The existence and stability of steady states of (B10) are given by the

799 following result.
36 Performance study of two serial interconnected chemostats

800 **Theorem 3** Assume that Assumption 1 is satisfied. The steady states of (B10) are:

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

$$D \geq \max\{r(f(S^{\text{in}}) - a_1), (1-r)(f(S^{\text{in}}) - a_2)\}. \quad (B14)$$

It is LES if and only if

$$D > \max\{r(f(S^{\text{in}}) - a_1), (1-r)(f(S^{\text{in}}) - a_2)\}.$$

• The steady state $E_1 = (S^{\text{in}}, 0, \bar{S}_2, \bar{x}_2)$ of washout in the first chemostat but not in the second one with

$$\bar{S}_2 = \lambda \left( \frac{D}{1-r} + a_2 \right), \quad \bar{x}_2 = \frac{D}{D+(1-r)a_2} (S^{\text{in}} - \bar{S}_2). \quad (B15)$$

It exists if and only if $D < (1-r)(f(S^{\text{in}}) - a_2)$. It is GAS if and only if

$$r(f(S^{\text{in}}) - a_1) \leq D \quad \text{and} \quad D < (1-r)(f(S^{\text{in}}) - a_2). \quad (B16)$$

It is LES if and only if

$$r(f(S^{\text{in}}) - a_1) < D < (1-r)(f(S^{\text{in}}) - a_2).$$

• The steady state $E_2 = (S_1^*, x_1^*, S_2^*, x_2^*)$ of persistence of the species in both chemostats with

$$S_1^* = \lambda \left( \frac{D}{r} + a_1 \right), \quad x_1^* = \frac{D}{D+ra_1} (S^{\text{in}} - S_1^*), \quad (B17)$$

$$x_2^* = \frac{D}{D+(1-r)a_2} \left( \frac{D}{D+ra_1} (S^{\text{in}} - S_1^*) + S_1^* - S_2^* \right) \quad (B18)$$

and $S_2^* = S_2^*(S^{\text{in}}, D, r)$ is the unique solution of the equation $h(S_2) = f(S_2)$ with $h$ defined by (B11). This steady state exists and is positive if and only if $D < r(f(S^{\text{in}}) - a_1)$. It is GAS and LES whenever it exists and is positive.

804 Proof The 4-dimensional system of ODEs (B10) has a cascade structure of two planar systems of ODEs, whose mathematical analysis is easy and well known in the 805 mathematical theory of the chemostat [16, 39]. Using this cascade structure, the 806 global behavior of the system is deduced from the global behaviour of planar systems 807 and Thieme's theory of asymptotically autonomous systems. 808

809 For the convenience of the reader the details of the proof are given in Appendix 810 D.1. □

811 **Proposition 10** The function $S^{\text{in}} \mapsto S_2^*(S^{\text{in}}, D, r)$ is decreasing.
## Appendix C Operating diagram

$$ \Phi_r^1 := \{(S^{in}, D) \in \mathbb{R}_+^2 : D = r(f(S^{in}) - a_1)\}, \quad (C19) $$

$$ \Phi_{1-r}^2 := \{(S^{in}, D) \in \mathbb{R}_+^2 : D = (1-r)(f(S^{in}) - a_2)\}. \quad (C20) $$

<footnote n="812">812 Proof D and r are fixed. Let $S^{in,1} > S^{in,2}$ and $h_i$ defined by (B11), with $S^{in} = S^{in,i}$,
813 $i = 1, 2$. Let $S_2^{*i}$, the solution of equation $f(S_2) = h_i(S_2)$, $i = 1, 2$. Using Lemma
814 8, $h_i$ is a decreasing hyperbola from $h_i(0)$ defined by (B12), with $S^{in} = S^{in,i}$, to
815 $h_i(S_1^*) = 0$. Since $h_1(0) < h_2(0)$, we have $h_1(S_2) < h_2(S_2)$ for all $S_2 \in (0, S_1^*)$.
816 Therefore, $S_2^{*1} < S_2^{*2}$, see Fig. B2(b). □</footnote>

<footnote n="817">817 This result means that the effluent steady state concentration of substrate
decreases when the influent concentration of substrate increases. This behavior
818 is very different from the single chemostat, where the effluent steady state
819 substrate concentration is independent of the influent substrate concentration.
820</footnote>

<footnote n="821">821 For the chemostat model, the operating diagram has as coordinates the input
822 substrate concentration $S^{in}$ and the dilution rate $D$, and shows how the solu-
tions of the system behave for different values of these two parameters. The
823 regions constituting the operating diagram correspond to different qualitative
824 asymptotic behaviors. Indeed, the main interest of an operating diagram is
825 to highlight the number and stability of the steady states for a given pair of
826 parameters ($S^{in}$, $D$). The input substrate concentration $S^{in}$ and the dilution
827 rate $D$ are the usual parameters manipulated by the experimenter of a chemo-
828 stat. Apart from these parameters, and the parameter $r$ that can be also chosen
829 by the experimenter but not easily changed as $S^{in}$ and $D$, all other parame-
ters have biological meaning and are fitted using experimental data from real
830 measurements of concentrations of micro-organisms and substrates. Therefore
831 the operating diagram is a bifurcation diagram, quite useful to understand the
832 possible behaviors of the solutions of the system from both the mathematical
833 and biological points of view.
834</footnote>

<footnote n="835">835 Here, we fix $r \in (0,1)$ and we depict in the plane ($S^{in}$, $D$) the regions
836 in which the solution of system (B10) globally converges towards one of the
837 steady state $E_0$, $E_1$ or $E_2$. From the results given in Theorem 3, it is seen that
838 these regions are delimited by the curves $\Phi_r^1$ and $\Phi_{1-r}^2$ defined by:</footnote>

<footnote n="841">841 When $a_1 = a_2 = 0$, as we have shown in [7], these curves meet only at one
point (the origin) and merge when $r = 1/2$. Therefore, in this case the curves
842 $\Phi_r^1$ and $\Phi_{1-r}^2$ separate the operating plane ($S^{in}$, $D$), in only three regions, see
843 [7, Figure 5]. This property continue to hold when $a_1 = a_2$, that is to say,
844 the curves intersect only at $(\lambda(a_1), 0)$ and merge when $r = 1/2$. In this case
845 the curves $\Phi_r^1$ and $\Phi_{1-r}^2$ separate the operating plane ($S^{in}$, $D$), in only three
846 regions, see Figure B3 (c) and (d). The novelty when $a_1$ and $a_2$ are different
847 and non null, is that the intersection of the curves $\Phi_r^1$ and $\Phi_{1-r}^2$ can lie outside
848 the $S^{in}$ axis. Therefore there can be four regions in the operating plane, as
849 depicted in Figure B3 (a) and (f). For the description of the intersection of
850</footnote>
38 Performance study of two serial interconnected chemostats

<figure><img src="image_16.png" /><figcaption>**Fig. B3** The operating diagram of (B10). The asymptotic behaviour in each region is depicted in Table C2.</figcaption></figure>

$$
\bar{r} := \frac{m-a_2}{2m-a_1-a_2}. \qquad (C21)
$$

<footnote n="851">851 the curves $\Phi_r^1$ and $\Phi_{1-r}^2$, we need some definitions and notations. Let $\bar{r} \in (0, 1)$
852 be defined by</footnote>

<footnote n="853">853 Note that if $a_1 < a_2$ then $\bar{r} < 1/2$, and if $a_1 > a_2$ then $\bar{r} > 1/2$. For $a_1 < a_2$
854 and $0 < r < \bar{r}$ (or $a_1 > a_2$ and $\bar{r} < r < 1$), we define the point $P = (S_P^{\text{in}}, D_P)$
855 of the operating plane by:</footnote>
<table><caption>**Table C2** The regions $I_k(r)$, $k = 0, 1, 2, 3$ of the operating diagram of (B10) and asymptotic behaviour in these various regions.</caption><thead><tr><th>Regions</th></tr></thead><tbody><tr><td>$I_0(r) = \{(S^{\text{in}}, D) : \max\{r(f(S^{\text{in}}) - a_1), (1-r)(f(S^{\text{in}}) - a_2)\} \le D\}$</td></tr><tr><td>$I_1(r) = \{(S^{\text{in}}, D) : r(f(S^{\text{in}}) - a_1) \le D \text{ and } D < (1-r)(f(S^{\text{in}}) - a_2)\}$</td></tr><tr><td>$I_2(r) = \{(S^{\text{in}}, D) : 0 < D < \min\{r(f(S^{\text{in}}) - a_1), (1-r)(f(S^{\text{in}}) - a_2)\}\}$</td></tr><tr><td>$I_3(r) = \{(S^{\text{in}}, D) : (1-r)(f(S^{\text{in}}) - a_2) \le D \text{ and } D < r(f(S^{\text{in}}) - a_1)\}$</td></tr></tbody></table>

<table><thead><tr><th></th><th>$I_0(r)$</th><th>$I_1(r)$</th><th>$I_2(r)$</th><th>$I_3(r)$</th></tr></thead><tbody><tr><th>$E_0$</th><td>GAS</td><td>U</td><td>U</td><td>U</td></tr><tr><th>$E_1$</th><td></td><td>GAS</td><td>U</td><td></td></tr><tr><th>$E_2$</th><td></td><td></td><td>GAS</td><td>GAS</td></tr></tbody></table>

$$S_P^{\text{in}} := \lambda \left( \frac{r a_1 (1-r) a_2}{2r-1} \right), \quad D_P := \frac{r (1-r) (a_2-a_1)}{1-2r}. \qquad (C22)$$

Note that $S_P^{\text{in}} > 0$ and $D_P > 0$. With these notations we can state the following result:

**Proposition 11** 1. If $a_1 < a_2$ then for all $r \in (0, \bar{r})$, the curves $\Phi_r^1$ and $\Phi_{1-r}^2$ intersect at the point $P$ and $\Phi_r^1$ is strictly below [resp. above] $\Phi_{1-r}^2$ for $S^{\text{in}} > S_P^{\text{in}}$ [resp. $S^{\text{in}} < S_P^{\text{in}}$], see Figure B3 (a). For all $r \in (\bar{r}, 1)$, $\Phi_r^1$ is strictly above $\Phi_{1-r}^2$, see Figure B3 (b).

2. If $a_1 > a_2$ then for all $r \in (\bar{r}, 1)$, the curves $\Phi_r^1$ and $\Phi_{1-r}^2$ intersect at the point $P$ and $\Phi_r^1$ is strictly above [resp. below] $\Phi_{1-r}^2$ for $S^{\text{in}} > S_P^{\text{in}}$ [resp. $S^{\text{in}} < S_P^{\text{in}}$], see Figure B3 (f). For all $r \in (0, \bar{r})$, $\Phi_r^1$ is below $\Phi_{1-r}^2$, see Figure B3 (e).

3. If $a_1 = a_2$ then, for $r = 1/2$, $\Phi_r^1 = \Phi_{1-r}^2$. Moreover, if $r < 1/2$ then $\Phi_r^1$ is strictly below $\Phi_{1-r}^2$, see Figure B3 (c) and, if $r > 1/2$ then $\Phi_r^1$ is strictly above $\Phi_{1-r}^2$, see Figure B3 (d).

*Proof* For $0 < r < 1$ and $S^{\text{in}} > \lambda(a_i)$ we define the function $\varphi_i$, $i = 1, 2$, by

$$\begin{aligned} \varphi_1(S^{\text{in}}, r) &= r(f(S^{\text{in}}) - a_1), \\ \varphi_2(S^{\text{in}}, r) &= (1-r)(f(S^{\text{in}}) - a_2). \end{aligned} \qquad (C23)$$

The curves $\Phi_r^1$ and $\Phi_{1-r}^2$, defined respectively by (C19) and (C20), intersect if and only if there exists $r \in (0, 1)$ and $S^{\text{in}} > \max(\lambda(a_1), \lambda(a_2))$ such that $\varphi_1(S^{\text{in}}, r) = \varphi_2(S^{\text{in}}, r)$, that is to say

$$f(S^{\text{in}}) = A(r), \quad \text{with} \quad A(r) := \frac{r a_1 (1-r) a_2}{2r-1}. \qquad (C24)$$

This equation has a solution $S^{\text{in}} > \max(\lambda(a_1), \lambda(a_2))$ if and only if

$$\max(a_1, a_2) < A(r) < m, \qquad (C25)$$
40 Performance study of two serial interconnected chemostats

874 where $m = \sup(f)$, as in (2). When these conditions are satisfied, the solution of (C24) is given by $S^{in} = \lambda(A(r))$, where $\lambda$ is the inverse function $pf$, i.e. the 875 break-even concentration defined by (3). Hence, $S^{in} = S_P^{in}$, given in (C22). The 876 corresponding intersection point of $\Phi_r^1$ and $\Phi_{1-r}^2$ is given by $D_P = r(f(S_P^{in}) - a_1)$, 877 which is the value given in (C22).

878 Let us determine now for which value of $r$, the conditions (C25) are satisfied. The 879 function $A$ is a homographic function. Its graphical representation is a hyperbola, 880 whose vertical asymptote is $r = 1/2$. Its derivative is given by

881 $A'(r) = \frac{a_2 - a_1}{(2r - 1)^2}$. (C26)
Note that $A(r) = m$ if and only if $r = \bar{r}$, where $\bar{r}$ is defined by (C21). Therefore 882 if $a_1 < a_2$ then, according to (C26), $A$ is increasing. Since $A(0) = a_2$, $A(\bar{r}) = m$, 883 and $\bar{r} < 1/2$, the condition (C25) is satisfied if and only if $0 < r < \bar{r}$. Similarly, if 884 $a_1 > a_2$, then, according to (C26), $A$ is decreasing. Since $A(1) = a_1$, $A(\bar{r}) = m$ and 885 $\bar{r} > 1/2$, the condition (C25) is satisfied if and only if $\bar{r} < r < 1$. Finally, if $a_1 = a_2$ 886 then $A(r) = a_1$ and the condition (C25) cannot be satisfied. 887

888 Suppose that $a_1 < a_2$. Note that for $0 < r < 1/2$, the condition $f(S^{in}) > A(r)$ 889 [resp. $f(S^{in}) < A(r)$] is equivalent to $\varphi_1(S^{in}, r) < \varphi_2(S^{in}, r)$ [resp. $\varphi_1(S^{in}, r) >$ 890 $\varphi_2(S^{in}, r)$]. Thus:

891 • If $r \in (0, \bar{r})$, then $f(S^{in}) < A(r)$ if and only if $S^{in} < S_P^{in}$, where $S_P^{in}$ is defined 892 by (C22). Hence, the curves $\Phi_r^1$ and $\Phi_{1-r}^2$ intersect at $P = (S_P^{in}, D_P)$ and 893 the curve $\Phi_r^1$ is strictly below [resp. above] the curve $\Phi_{1-r}^2$, for all $S^{in} > S_P^{in}$ 894 [resp. $S^{in} < S_P^{in}$].

895 • If $r \in [\bar{r}, 1/2)$ then $f(S^{in}) < A(r)$ for all $S^{in} > 0$, so that the curve $\Phi_r^1$ is 896 strictly above the curve $\Phi_{1-r}^2$.

897 • If $r \in [1/2, 1)$, then, using $r \ge 1 - r$ and $a_1 < a_2$, one has $\varphi_1(S^{in}, r) >$ 898 $\varphi_2(S^{in}, r)$. Therefore, the curve $\Phi_r^1$ is strictly above the curve $\Phi_{1-r}^2$.

899 If $a_1 > a_2$, the proof is similar to the case $a_1 < a_2$.

900 If $a_1 = a_2$ then $\varphi_1(S^{in}, r) = \varphi_2(S^{in}, r)$ is equivalent to $r(f(S^{in}) - a_1) = (1 - 901 r)(f(S^{in}) - a_1)$. Therefore, $r = 1 - r$, that is $r = 1/2$. In this case the curves $\Phi_r^1$ and 902 $\Phi_{1-r}^2$ merge. In addition, if $r < 1/2$ [resp. $r > 1/2$] then $r < 1 - r$ [resp. $r > 1 - r$] 903 and the curve $\Phi_r^1$ is strictly below [resp. above] the curve $\Phi_{1-r}^2$. This ends the proof 904 of the proposition. $\square$

905 For any $r \in (0, 1)$, the curves $\Phi_r^1$ and $\Phi_{1-r}^2$, defined by (C19) and (C20), 906 respectively split the plane $(S^{in}, D)$ in the regions denoted $I_0(r)$, $I_1(r)$, $I_2(r)$ 907 and $I_3(r)$ and defined in Table C2. These regions are depicted in Fig. B3 in 908 the cases $a_1 < a_2$, $a_1 = a_2$ and $a_1 > a_2$.

909 The behavior of the system in each region, when it is not empty, is given in 910 Table C2. Notice that $E_1$ exists in both regions $I_1(r)$ and $I_2(r)$, but is stable 911 only when $(S^{in}, D)$ is fixed in $I_1(r)$.

912 When $a_1 = a_2 = 0$ then $\lambda(a_1) = \lambda(a_2) = 0$ and the curves $\Phi_r^1$ and $\Phi_{1-r}^2$ of 913 the operating diagram start from the origin of the plane $(S^{in}, D)$ and merge 914 for $r = 1/2$. Therefore, the diagrams shown in panels (a), (b), (c), (d), (e) 915 and (f) of Fig. B3 are reduced to only two different cases characterized by 916 $0 < r < 1/2$ and $1/2 < r < 1$, as shown in Figure 5 of [7]. There is no changes
917

918

in the stability of the steady states and in the number of the regions depicted
in the operating diagram.

919

920

921

922

923

924

925

926

927

928

929

930

931

932

933

934

935

This result reveals an interplay between spatial heterogeneity (the ratio *r* of
volume distribution between tanks) and the mortality heterogeneity (difference
between *a*₁ and *a*₂). Indeed, panels (a) and (f) of Fig. B3 bring a particular
feature when mortality rates are different: domains *I*₁(*r*) and *I*₂(*r*) can appear
or disappear playing only with the spatial distribution *r*, a phenomenon which
does not happens when mortality is identical in each tank. This shows that
the existence of domains *I*₁(*r*) and *I*₂(*r*) is controlled by a relative toxicity in
the tanks, and not only the spatial distribution as it is the case for identical
mortality. This feature can have interest when practitioners can adjust pH or
other abiotic parameters having impacts on the mortality rate, independently
in each tank. Given operating parameters *S*<sup>*in*</sup>, *D* and *r*, panels (a) and (f) of
Fig. B3 show that it is theoretically possible to pass from domain *I*₂(*r*) to *I*₂(*r*)
when mortality parameter is diminished only in the second tank. In practice,
being in domain *I*₂(*r*) might be more desirable than *I*₂(*r*) with respect to some
dysfunctioning of the first tank that can drop suddenly its biomass to zero.
Indeed, in *I*₂(*r*), the second tank is no conducted to the wash-out differently
to the *I*₂(*r*) case.

936

937

938

939

When $a_1 = a_2 = a$, which is the case corresponding to the system (1)
considered in Section 2, only panels (c,d) of Fig. B3 are encountered, as shown
in Fig. 2. We describe hereafter the bifurcations that occur in this particular
case. The general case i.e. when $a_1 \neq a_2$ is similar.

940

941

942

943

944

945

946

Remark 4 Transcritical bifurcations occur in the limit cases $D = r(f(S^{\text{in}}) - a)$ and $D = (1-r)(f(S^{\text{in}}) - a)$, for system (1). If $0 < r < 1/2$ then, we have a transcritical bifurcation of $E_0$ and $E_1$ when $D = (1-r)(f(S^{\text{in}}) - a)$ and a transcritical bifurcation of $E_1$ and $E_2$ when $D = r(f(S^{\text{in}}) - a)$. If $1/2 < r < 1$ then, we have a transcritical bifurcation of $E_0$ and $E_1$ when $D = (1-r)(f(S^{\text{in}}) - a)$ and a transcritical bifurcation of $E_0$ and $E_2$ when $D = r(f(S^{\text{in}}) - a)$. If $r = 1/2$ and $D = (f(S^{\text{in}}) - a)/2$ then, we have transcritical bifurcations of $E_0$ and $E_1$, and $E_0$ and $E_2$, simultaneously.

947

# Appendix D Proofs

948

## D.1 Proof of Theorem 3

949

950

951

952

953

954

955

956

We begin by the existence of steady states. The steady states are the solutions
of the set of equations $\dot{S}_1 = 0$, $\dot{x}_1 = 0$, $\dot{S}_2 = 0$, $\dot{x}_2 = 0$. From equation $\dot{x}_1 = 0$,
it is deduced that $x_1 = 0$ or $f(S_1) = D/r + a_1$. Suppose first that $x_1 = 0$.
Then, from equation $\dot{S}_1 = 0$ it is deduced that $S_1 = S^{\text{in}}$ and from equation
$\dot{x}_2 = 0$ it is deduced that $x_2 = 0$ or $f(S_2) = D/(1-r) + a_2$. If $x_2 = 0$,
then from equation $S_2 = 0$ it is deduced that $S_2 = S^{\text{in}}$. Hence we obtain the
steady state $E_0 = (S^{\text{in}}, 0, S^{\text{in}}, 0)$, which always exist. On the other hand, if
$f(S_2) = D/(1-r) + a_2$, then $S_2 = \bar{S}_2$, defined in (B15). From equation $\dot{S}_2 = 0$,
it is deduced that $x_2 = \bar{x}_2$, defined in (B15). Hence we obtain the steady state

957
42 Performance study of two serial interconnected chemostats

958 $E_1 = (S^{\text{in}}, 0, \bar{S}_2, \bar{x}_2)$. This steady state exists if and only if $S^{\text{in}} > \bar{S}_2$, that is
959 $D < (1-r)(f(S^{\text{in}}) - a_2)$.

960 Suppose now that $f(S_1) = D/r + a_1$. Then $S_1 = S_1^*$, defined in (B17). From
961 equation $\dot{S}_1 = 0$, it is deduced that $x_1 = x_1^*$, defined in (B17). From equation
962 $\dot{S}_2 + \dot{x}_2 = 0$, it is deduced that

$$x_2 = \frac{D}{D + (1-r)a_2}(S_1^* + x_1^* - S_2). \qquad (\text{D27})$$

963 Replacing $x_2$ by this expression in the equation $\dot{S}_2 = 0$, it is deduced that
964 $f(S_2) = h(S_2)$, where $h$ is defined by (B11). Hence $S_2 = S_2^*$, which is the
965 unique solution of the equation $f(S_2) = h(S_2)$, as shown in Figure B2 (a).
966 Replacing $S_2$ by $S_2^*$ in (D27) gives $x_2 = x_2^*$, defined by (B18). Consequently,
967 we obtain the steady state $E_2 = (S_1^*, x_1^*, S_2^*, x_2^*)$. This steady state is positive
968 if and only if $S^{\text{in}} > S_1^*$, which is equivalent to $D < r(f(S^{\text{in}}) - a_1)$.

969 Let us now study the local stability. Since the system has a cascade struc-
970 the stability analysis reduces to the study of square 2×2 matrices. Indeed,
971 the Jacobian matrix associated to system (B10) is the lower triangular matrix
972 by blocs, $J = \begin{pmatrix} A & 0 \\ B & C \end{pmatrix}$ where $B$ is the diagonal matrix whose diagonal elements
973 are $D/(1-r)$, and $A$ and $C$ are given by:

$$A = \begin{pmatrix} -\frac{D}{r} - f'(S_1)x_1 & -f(S_1) \\ f'(S_1)x_1 & -\frac{D}{r} + f(S_1) - a_1 \end{pmatrix}, \\ C = \begin{pmatrix} -\frac{D}{1-r} - f'(S_2)x_2 & -f(S_2) \\ f'(S_2)x_2 & -\frac{D}{1-r} + f(S_2) - a_2 \end{pmatrix},$$

974 Hence, the eigenvalues of $J$ are the ones of $A$ and $C$.

975 For $E_0$, the eigenvalues are $-D/r$, $-D/r + f(S^{\text{in}}) - a_1$, $-D/(1-r)$ and
976 $-D/(1-r) + f(S^{\text{in}}) - a_2$. They are negative if and only if $D > r(f(S^{\text{in}}) - a_1)$
977 and $D > (1-r)(f(S^{\text{in}}) - a_2)$. Therefore, $E_0$ is LES if and only if the condition
978 in the theorem is satisfied.

979 For $E_1$, the eigenvalues of $A$ are $-D/r + f(S^{\text{in}}) - a_1$ and $-D/r$. The first
980 eigenvalue is negative if and only if $D > r(f(S^{\text{in}}) - a_1)$. On the other hand,
981 since the determinant of $C$ is positive, and its trace is negative, the eigenvalues
982 of $C$ are of negative real parts. Therefore, $E_1$ is LES if and only if the condition
983 in the theorem is satisfied.

984 For $E_2$, the determinant of $A$ is positive and its trace is negative. On the
985 other hand, using the notation $C_{E_2}$ for the matrix $C$ evaluated at $E_2$, we have

$$\det(C_{E_2}) = \left(-\frac{D}{1-r} - f'(S_2)x_2^*\right) \left(-\frac{D}{1-r} - a_2 + f(S_2^*)\right) + f(S_2^*)f'(S_2^*)x_2^*, \\ \operatorname{tr}(C_{E_2}) = -2\frac{D}{1-r} - a_2 - f'(S_2^*)x_2^* + f(S_2^*).$$

986 Note that $h(S_2) < D/(1-r) + a_2$ for all $S_2 \in (0, S_1^*)$. Therefore, from (B11),
987 we have $f(S_2^*) = h(S_2^*) < D/(1-r) + a_2$. Consequently, $\det(C_{E_2})$ and $\operatorname{tr}(C_{E_2})$
988 are respectively positive and negative. Therefore, $E_2$ is LES whenever it exists, that is $D < r(f(S^{in}) - a_1)$.

989 For the study of the global stability we use the cascade structure of the system (B10) and Thieme's Theorem (see Theorem A1.9 of [16]). In the rest of the proof, we denote by $(S_1(t), x_1(t), S_2(t), x_2(t))$ the solution of (B10) with the initial condition $(S_1^0, x_1^0, S_2^0, x_2^0)$. Then, $(S_1(t), x_1(t))$ is the solution of system

$$ 
\begin{aligned}
\dot{S}_1 &= \frac{D}{r}(S^{in} - S_1) - f(S_1)x_1 \\
\dot{x}_1 &= -\frac{D}{r}x_1 + f(S_1)x_1 - a_1x_1
\end{aligned}
\qquad (D28) $$

994 with initial condition $(S_1^0, x_1^0)$ and $(S_2(t), x_2(t))$ is the solution of the non-autonomous system of differential equations

$$ 
\begin{aligned}
\dot{S}_2 &= \frac{D}{1-r}(S_1(t) - S_2) - f(S_2)x_2 \\
\dot{x}_2 &= \frac{D}{1-r}(x_1(t) - x_2) + f(S_2)x_2 - a_2x_2
\end{aligned}
\qquad (D29) $$

996 with the initial condition $(S_2^0, x_2^0)$. The system (D28) is the classical model of a single chemostat. Its asymptotic behaviour is well known (see, for instance, Proposition 2.2 of [16]). This system admits the steady states:

$$ e_1^0 = (S^{in}, 0) \quad \text{and} \quad e_1^1 = (S_1^*, x_1^*) \qquad (D30) $$

999 where $S_1^*$ and $x_1^*$ are defined by (B17). Two cases must be distinguished.

1000 Firstly, if $\lambda(D/r + a_1) \ge S^{in}$, that is $D \ge r(f(S^{in}) - a_1)$ then, $e_1^0$, defined in (D30), is GAS for (D28) in the nonnegative quadrant. Hence, for any non-negative initial condition $(S_1^0, x_1^0)$,

$$ \lim_{t \to +\infty} (S_1(t), x_1(t)) = (S^{in}, 0). \qquad (D31) $$

1003 Therefore, the system (D29) is asymptotically autonomous with the limiting system

$$ 
\begin{aligned}
\dot{S}_2 &= \frac{D}{1-r}(S^{in} - S_2) - f(S_2)x_2 \\
\dot{x}_2 &= -\frac{D}{1-r}x_2 + f(S_2)x_2 - a_2x_2.
\end{aligned}
\qquad (D32) $$

1005 Recall that the solutions of (D29) are positively bounded. Therefore, we shall use Thieme's results which apply for bounded solutions.

1007 The system (D32) represents the classical model of a single chemostat. It admits the two steady states $e_0^2 = (S^{in}, 0)$ and $e_1^2 = (\bar{S}_2, \bar{x}_2)$, with $(\bar{S}_2, \bar{x}_2)$ defined by (B15). Two subcases must be distinguished.

1010 • If $\lambda(D/(1-r) + a_2) \ge S^{in}$, that is $D \ge (1-r)(f(S^{in}) - a_2)$ then, $e_0^2$ is GAS in the nonnegative quadrant. Using Thieme's Theorem, we deduce that for any nonnegative $(S_2^0, x_2^0)$, the solution $(S_2(t), x_2(t))$ of (D29) converges towards $e_0^2 = (S^{in}, 0)$. Using (D31) we deduce that, when $D \ge \max(r(f(S^{in}) - a_1), (1-r)(f(S^{in}) - a_2))$, the solution $(S_1(t), x_1(t), S_2(t), x_2(t))$ of (B10) converges towards $E_0 = (S^{in}, 0, S^{in}, 0)$, which proves (B14).
44 Performance study of two serial interconnected chemostats

• In contrast, if $\lambda(D/(1-r) + a_2) < S^{in}$, that is $D < (1-r)(f(S^{in}) - a_2)$ then, both steady states $e_0^0$ and $e_1^0$ exist and $e_1^0$ is GAS in the positive quadrant. Although system (D29) has the saddle point $e_0^0$, no polycycle can exist. Using Thieme's Theorem, for any positive $(S_2^0, x_2^0)$, the solution $(S_2(t), x_2(t))$ of (D29) converges towards $e_1^0 = (\bar{S}_2, \bar{x}_2)$. Using (D31) we deduce that, if $r(f(S^{in}) - a_1) \le D$ and $D < (1-r)(f(S^{in}) - a_2)$, then the solution $(S_1(t), x_1(t), S_2(t), x_2(t))$ of (B10) converges towards $E_1 = (S^{in}, 0, \bar{S}_2, \bar{x}_2)$, which proves (B16).

Secondly, if $\lambda(D/r + a_1) < S^{in}$, that is $D < r(f(S^{in}) - a_1)$ then, $e_1^0$, defined in (D30), is GAS for (D28) in the positive quadrant. Hence, for any positive initial condition $(S_1^0, x_1^0)$

$$ \lim_{t \to +\infty} (S_1(t), x_1(t)) = (S_1^*, x_1^*). \qquad (D33) $$

Therefore, the system (D29) is asymptotically autonomous with the limiting system

$$ \begin{aligned} \dot{S}_2 &= \frac{D}{1-r}(S_1^* - S_2) - f(S_2)x_2 \\ \dot{x}_2 &= \frac{D}{1-r}(x_1^* - x_2) + f(S_2)x_2 - a_2x_2. \end{aligned} \qquad (D34) $$

The system (D34) represents the classical model of a single chemostat with an input biomass. In this case, there is no washout and the system (D34) always admits one LES steady state $e_2 = (S_2^*, x_2^*)$ with positive biomass defined by (B18) and $S_2^*$ the unique solution of $h(S_2) = f(S_2)$.

Let us show that this steady state is GAS for (D34). Assume that $x_2 > 0$. Consider the change of variable $\xi = \ln(x_2)$. The system (D34) becomes as

$$ \begin{aligned} \dot{S}_2 &= \frac{D}{1-r}(S_1^* - S_2) - f(S_2)e^{\xi} \\ \dot{\xi} &= \frac{D}{1-r}(x_1^*e^{-\xi} - 1) + f(S_2) - a_2. \end{aligned} \qquad (D35) $$

The divergence of the vector field

$$ \psi(S_2, \xi) = \left[ \frac{D}{1-r}(S_1^* - S_2) - f(S_2)e^{\xi} \frac{D}{1-r}(x_1^*e^{-\xi} - 1) + f(S_2) - a_2 \right] $$

associated to (D35) is $\text{div}\psi(S_2, \xi) = -\frac{D}{1-r}(1 + x_1^*e^{\xi}) - f'(S_2)e^{\xi}$. It is negative. Thus, using Bendixon-Dulac criterion, system (D35) cannot have a periodic solution. Hence, system (D34) has no cycle in the positive quadrant. For any non negative initial condition $(S_2^0, x_2^0)$, the solution of (D34) is bounded. Hence, the $\omega$-limit set of $(S_2^0, x_2^0)$, denoted $\omega(S_2^0, x_2^0)$, is non-empty and included in the positive quadrant. If $e_2 \notin \omega(S_2^0, x_2^0)$ then, using Poincaré-Bendixon Theorem, $\omega(S_2^0, x_2^0)$ is a limit cycle, but the system does not present any, due to the divergence property. One then deduces $e_2 \in \omega(S_2^0, x_2^0)$ and, as $e_2$ is LES, then $\omega(S_2^0, x_2^0) = \{e_2\}$. Consequently, $e_2$ is GAS for (D34) in the positive quadrant. Using again Thieme's Theorem, for any positive $(S_2^0, x_2^0)$, the solution $(S_2(t), x_2(t))$ of (D29) converges towards $e_2 = (S_2^*, x_2^*)$. Using (D33) we deduce
## D.2 Proof of Lemma 4

$$G_2(S^{in}, D, r) < \max_{d \in (0, r\delta)} G_2(S^{in}, d, r), \quad (D, r) \in [0, \underline{D}] \times [\underline{r}, 1]. \quad (D36)$$

$$D_{\varepsilon} := \varepsilon \left( a + \max_{s \in [0, S^{in}]} f'(s)(S^{in} - s) \right) < \underline{D} \quad (D37)$$

$$\mathcal{D}_{\varepsilon} := \{ (D, r); \ D \in (D_{\varepsilon}, \delta), \ r \in \left( \max \left( \underline{r}, \frac{D}{\delta} \right), 1 + \varepsilon \right) \}.$$

$$\varphi(S_2, D, r) = (D + (1 - r)a)(\lambda(D/r + a) - S_2) - (1 - r)f(S_2)\left(\frac{DS^{in} + ra\lambda(D/r + a)}{D + ra} - S_2\right),$$

$$\varphi(S_2, D, r) = (1 - r) \left( \frac{DS^{in} + ra\lambda(D/r + a)}{D + ra} - S_2 \right) (h(S_2) - f(S_2))$$

$$\partial_{S_2} \varphi \Big|_{S_2=S_2^*} = (1-r) \left( \frac{DS^{in} + ra\lambda(D/r + a)}{D + ra} - S_2 \right) \left( \partial_{S_2} h - f' \right) \Big|_{S_2=S_2^*} < 0.$$

<footnote n="1047">1047 that, if $D < r(f(S^{in}) - a_1)$, then the solution $(S_1(t), x_1(t), S_2(t), x_2(t))$ of (B10) converges towards $E_2 = (S_1^*, x_1^*, S_2^*, x_2^*)$. This ends the proof of the theorem.</footnote>

<footnote n="1050">1050 Let us fix $S^{in}$ such that $\delta := f(S^{in}) - a > 0$. The proof consists in showing that the function $(D, r) \mapsto G_2(S^{in}, D, r)$ can be formally extended as a $C^2$ function for values of $r$ larger than 1 (although such values have no physical meaning). Recall first that for any $D \in (0, \delta)$, one has $G_2(S^{in}, D, 1) = G_{chem}(S^{in}, D)$. As $G_2(S^{in}, \bar{D}(1), 1) > 0$ and $G_2(S^{in}, 0, 1) = 0$, there exists by continuity of the function $G_2$, numbers $\underline{D} \in (0, \bar{D}(1))$, $\underline{r} \in (0, 1)$ such that</footnote>

<footnote n="1056">1056 Let $\varepsilon > 0$ be such that</footnote>

<footnote n="1057">1057 and consider the domain</footnote>

<footnote n="1058">1058 Note that for any $(D, r) \in \mathcal{D}_{\varepsilon}$, the number $\lambda(D/r + a) = f^{-1}(D/r + a)$ is well defined. Posit the function</footnote>

<footnote n="1060">1060 where $(S_2, D, r) \in (0, S^{in}) \times \mathcal{D}_{\varepsilon}$. As $f$ is $C^2$, $\varphi$ is $C^2$ on $(0, S^{in}) \times \mathcal{D}_{\varepsilon}$. For $r < 1$ and $(D, r) \in \mathcal{D}_{\varepsilon}$, one has</footnote>

<footnote n="1062">1062 where $h$ is the function defined in (5). According to Lemma 8, $h$ is positive decreasing on $(0, \lambda(D/r + a)$, and $h - f$ admits an unique zero $S_2^* = S_2^*(S^{in}, D, r)$ on $(0, \lambda(D/r + a)$. Then, one can write</footnote>

<footnote n="1065">1065 For $r \in [1, 1 + \varepsilon)$ and $(D, r) \in \mathcal{D}_{\varepsilon}$, on has</footnote>
46 Performance study of two serial interconnected chemostats

$$
\partial_{S_2} \varphi = -(D + (1-r)a) - (1-r)f'(S_2) \left( \frac{DS^{in} + ra\lambda(D/r + a)}{D + ra} - S_2 \right) \\
+ (1-r)f(S_2),
$$

which is negative for any $S_2 \in (0, S^{in})$ thanks to condition (D37). As $\varphi(0, D, r) > 0$ and $\varphi(S^{in}, D, r) < 0$, we deduce the existence of a unique $S_2^* = S_2^*(S^{in}, D, r)$ in $(0, S^{in})$ such that $\varphi(S_2^*, D, r) = 0$, which also verifies $\partial_{S_2}\varphi < 0$ at $S_2 = S_2^*$.

Then, by the Implicit Function Theorem, the function $(D, r) \mapsto S_2^*(S^{in}, D, r)$ is $C^2$ on $\mathcal{D}_\varepsilon$. Recall that for $r < 1$ and $D < r\delta$, on has the expression $G_2(S^{in}, D, r) = VD(S^{in} - S_2^*(S^{in}, D, r))$ (see Proposition 4). We extend now the function $(D, r) \mapsto G_2(S^{in}, D, r)$ with this last $C^2$ expression on $\mathcal{D}_\varepsilon$. As $G_2(S^{in}, D, 1) = G_{chem}(S^{in}, D)$ for any $D \in (0, \delta)$, one deduces, by continuity of the partial derivatives of $G_2$ with respect to $D$ and property (D36), the existence of $\mathcal{V}_D$, $\mathcal{V}_r$ as neighborhoods respectively of $\bar{D}(1)$ and $1$ with $\mathcal{V}_D \times \mathcal{V}_r \subset \mathcal{D}_\varepsilon$ such that for any $r \in \mathcal{V}_r$, the function $D \mapsto G_2(S^{in}, D, r)$ possesses the following properties

1. it is strictly concave on $\mathcal{V}_D$,

2. it is increasing on $(D_\varepsilon, \bar{D}(1)) \setminus \mathcal{V}_D$ and decreasing on $(\bar{D}(1), r\delta) \setminus \mathcal{V}_D$,

3. its maximum over $(0, r\delta)$ is not reached for $D \leq D_\varepsilon$.

We thus deduce that $D \mapsto G_2(S^{in}, D, r)$ admits a unique maximum $\bar{D}(r)$ on $(0, r\delta)$, for any $r \in \mathcal{V}_r$.

Finally, for any $r \in \mathcal{V}_r$, $\bar{D}(r)$ is characterized as the zero of the map $D \mapsto F(D, r)$ where $F$ is the $C^1$ function

$F(D, r) := \partial_D G_2(S^{in}, D, r)$

From property 1. above, one obtains

$$
\partial_D F(\bar{D}(r), r) = \partial_{\bar{D}D}^2 G_2(S^{in}, \bar{D}(r), r) < 0, \quad r \in \mathcal{V}_r
$$

and by the Implicit Function Theorem, there exists a neighborhood $\mathcal{V}_1 \subset \mathcal{V}_r$ of 1 such that $\bar{D}$ is $C^1$ on $\mathcal{V}_1$, which ends the proof of the lemma.

D.3 Proof of Proposition 6

$S^{in}$ being fixed, we shall drop the $S^{in}$ dependency in the expressions of $S_i^*$, $x_i^*$
$(i = 1, 2)$ and $G_2$. Thus, let us define

$$
\begin{align*}
G(D, r) &:= G_2(S^{in}, D, r), \\
F_i(D, r) &:= f(S_i^*(D, r))x_i^*(D, r), \quad i = 1, 2,
\end{align*}
$$

as functions of $D \ge 0$ and $r \in \mathcal{V}_1 \cap \{r < 1\}$. Remark from the expression of $F_1$,
that it is well defined as well as its partial derivatives at $r = 1$. In addition,
for the limiting case $r = 1$, using Lemma 9, for all $D \ge 0$, one has
$$ 
\begin{aligned}
S_2^*(D, 1) &= S_1^*(D, 1) = \lambda(D + a) \\
x_2^*(D, 1) &= x_1^*(D, 1) = \frac{D}{D+a}(S^{in} - \lambda(D + a)).
\end{aligned} 
\quad (D38) $$

Thus, for all $D \ge 0$, one has

$$ F_1(D, 1) = F_2(D, 1), \quad (D39) $$

and $F_2$ is also well defined for $r = 1$. Thus, according to (37), for all $D \ge 0$ and $r \in V_1 \cap \{r \le 1\}$, one has

$$ G(D, r) = rF_1(D, r) + (1-r)F_2(D, r), $$

and from Lemma 4, for $r \in V_1 \cap \{r < 1\}$, one has

$$ \overline{G}(r) = G(\overline{D}(r), r), \quad (D40) $$

with $\overline{G}$ defined by (42). For convenience, for a function $E$ of $(D, r)$ that is differentiable, we shall define the three following functions: $\overline{E}(r) := E(\overline{D}(r), r)$ and

$$ \partial_r E(r) := \frac{\partial E}{\partial r}(\overline{D}(r), r), \quad \partial_D E(r) := \frac{\partial E}{\partial D}(\overline{D}(r), r). $$

Therefore, the function $\overline{G}$ writes

$$ \overline{G}(r) = r\overline{F}_1(r) + (1-r)\overline{F}_2(r), \quad \text{for } r \in V_1 \cap \{r < 1\}. \quad (D41) $$

As the functions $F_i$ are differentiable and as $\overline{D}(r)$ is a maximizer of $D \mapsto rF_1(D, r)+(1-r)F_2(D, r)$ on the interior of the interval $[0, f(S^{in})-a]$, one has

$$ r\partial_D F_1(r) + (1-r)\partial_D F_2(r) = 0, \quad \text{for } r \in V_1 \cap \{r < 1\}, \quad (D42) $$

and $\partial_D F_1(1) = 0$. As $f$ is $C^2$ and $\overline{D}$ is assumed to be differentiable on $V_1 \cap \{r < 1\}$, $\overline{G}$ is differentiable and from (D41), for all $r \in V_1 \cap \{r < 1\}$, one has

$$ \overline{G}'(r) = \overline{F}_1(r) - \overline{F}_2(r) + r\partial_r F_1(r) + (1-r)\partial_r F_2(r) \\ + (r\partial_D F_1(r) + (1-r)\partial_D F_2(r))\overline{D}'(r), $$

and with (D42), for all $r \in V_1 \cap \{r < 1\}$, one has simply

$$ \overline{G}'(r) = \overline{F}_1(r) - \overline{F}_2(r) + r\partial_r F_1(r) + (1-r)\partial_r F_2(r). \quad (D43) $$

Let us now determine the limits of the terms of the right side of this last equality when $r$ tends to 1. Firstly, according to (D39), one has in particular

$$ \overline{F}_1(1) = \overline{F}_2(1). \quad (D44) $$

Secondly, remark that the dynamics of the first tank is parameterized by the single dilution rate $D_1 = D/r$, the other parameters being fixed (see the expression (B17)). The function $F_1$ takes then the form $F_1(D, r) = \tilde{F}_1(D/r)$ where $\tilde{F}_1$ is a smooth function. Therefore, one has

$$ \partial_D F_1(r) = -\frac{r}{\overline{D}(r)} \partial_r F_1(r). \quad (D45) $$
48 Performance study of two serial interconnected chemostats

1114 As $\partial_D F_1(1) = 0$ then one deduces

$$ \partial_r F_1(1) = 0. \qquad (D46) $$

1115 Finally, from $\dot{S}_2 = 0$, for all $r \in V_1 \cap \{r < 1\}$, one gets

$$ F_2(D, r) = \frac{D}{1-r} (S_1^*(D, r) - S_2^*(D, r)). \qquad (D47) $$

1116 Differentiating (D47) with respect to r gives

$$ \frac{\partial F_2}{\partial r}(D, r) = \frac{D}{1-r} \left( \frac{\partial S_1^*}{\partial r}(D, r) - \frac{\partial S_2^*}{\partial r}(D, r) \right) + \frac{D}{(1-r)^2} (S_1^*(D, r) - S_2^*(D, r)) $$

1117 which can be written equivalently as

$$ (1-r) \frac{\partial F_2}{\partial r}(D, r) = D \left( \frac{\partial S_1^*}{\partial r}(D, r) - \frac{\partial S_2^*}{\partial r}(D, r) \right) + F_2(D, r). $$

1118 Thus, for $D = \bar{D}(r)$, one has

$$ (1-r)\partial_r F_2(r) = \bar{D}(r)(\partial_r S_1^*(r) - \partial_r S_2^*(r)) + \bar{F}_2(r). $$

1119 Notice that for $D = \bar{D}(r)$, (D47) gives

$$ \bar{F}_2(r) = \frac{\bar{D}(r)}{1-r} (\bar{S}_1^*(r) - \bar{S}_2^*(r)), \quad \text{for all } r \in V_1 \cup \{r < 1\}. \qquad (D48) $$

1120 Using L'Hôpital's rule in (D48) when $r$ tends to 1, one gets

$$ \bar{F}_2(1) = \lim_{r \to 1^{-}} \frac{\frac{\partial'(r)(\bar{S}_1^*(r) - \bar{S}_2^*(r))}{\bar{D}(r)(\partial_r S_1^*(r) - \partial_r S_2^*(r))} - 1}{-1} $$

1121 and using (D38) and (D44), one obtains

$$ \bar{F}_1(1) = \lim_{r \to 1^{-}} -\bar{D}(r)(\partial_r S_1^*(r) - \partial_r S_2^*(r)). $$

1122 Consequently, one has

$$ \lim_{r \to 1^{-}} (1-r) \partial_r F_2(r) = 0. \qquad (D49) $$

1123 With (D44), (D46) and (D49), expression (D43) gives the existence of the limit of $\bar{G}'$ when $r$ tends to 1 with $r < 1$, which is

$$ \bar{G}'(1^{-}) = 0. \qquad (D50) $$

1125 Note that $\bar{G}''(1^{-})$ exists if and only if $\lim_{r \to 1^{-}} \frac{\bar{G}'(r) - \bar{G}'(1)}{r-1}$ exists. Using (D50) and (D43), one has

$$ \frac{\bar{G}'(r) - \bar{G}'(1^{-})}{r-1} = - \frac{\bar{G}'(r)}{1-r} = - \frac{\bar{F}_1(r) - \bar{F}_2(r) + r \partial_r F_1(r) + (1-r) \partial_r F_2(r)}{1-r} \qquad (D51) $$
1127 On the one hand, using L'Hôpital's rule, one has

$$ \lim_{r \to 1^{-}} \frac{\overline{F}_{1}(r) - \overline{F}_{2}(r)}{1 - r} = \lim_{r \to 1^{-}} \frac{\overline{F}'_{1}(r) - \overline{F}'_{2}(r)}{-1}. $$

1128 Recall that $\partial_r F_1(1) = 0$ and thus one has $\overline{F}'_1(1) = 0$. Consequently, one has

$$ \lim_{r \to 1^{-}} frac{\overline{F}_1(r) - \overline{F}_2(r)}{1-r} = \lim_{r \to 1^{-}} \overline{F}'_2(r) = \lim_{r \to 1^{-}} \partial_r F_2(r) + \partial_D F_2(r) \overline{D}'(r). \quad (D52) $$

1129 On the other hand, using (D42) and (D45), one has

$$ \frac{r}{1-r} \partial_r F_1(r) = \frac{\overline{D}(r)}{r} \partial_D F_2(r). \quad (D53) $$

1130 Thus, according to (D51), (D52) and (D53), one gets

$$ \lim_{r \to 1^{-}} \frac{\overline{G}'(r) - \overline{G}'(1^{-})}{r - 1} = \lim_{r \to 1^{-}} \quad (D54) $$

1131 Let us show now that the limit of $\partial_D F_2(r)$ is 0 when $r$ tends to 1. One has

$$ \frac{\partial F_2}{\partial D} = f'(S_2^*) \frac{\partial S_2^*}{\partial D} x_2^* + f(S_2^*) \frac{\partial x_2^*}{\partial D}. $$

1132 Let use the expression $G(D, r) = D(S^{in} - S_2^*(D, r))$ given by Proposition 4.

1133 As $\overline{D}(r)$ is a maximizer then one has

$$ \partial_D G(r) = S^{in} - \overline{S}_2^*(r) - \overline{D}(r) \partial_D S_2^*(r) = 0. $$

1134 Using (D38), one then deduces

$$ \partial_D S_2^*(1^{-}) = \frac{S^{in} - \lambda(\overline{D}(1) + a)}{\overline{D}(1)}. $$

1135 In addition, using expressions (B18) and (D38), one gets

$$ \partial_D x_2^*(1^{-}) = - \frac{\overline{D}(1)}{(\overline{D}(1) + a)^2} (S^{in} - \lambda(\overline{D}(1) + a)), $$

1136 and hence the limit of $\partial_D F_2$ when $r$ tends to 1 exists:

$$ \partial_D F_2(1^{-}) = \frac{S^{in} - \lambda(\overline{D}(1) + a)}{\overline{D}(1) + a} f'\left(\lambda(\overline{D}(1) + a)\right) A, $$

1137 where $A = S^{in} - \lambda(\overline{D}(1) + a) - \frac{\overline{D}(1)}{f'(\lambda(\overline{D}(1)+a))}$. Thus, one has

$$ \partial_D F_2(1^{-}) = \frac{S^{in} - \lambda(\overline{D}(1) + a)}{\overline{D}(1) + a} f'\left(\lambda(\overline{D}(1) + a)\right) (S^{in} - g(\overline{D}(1))), $$

1138 with $g$ defined by (A8). According to Proposition 9, one has $S^{in} - g(\overline{D}(1)) = 0$.

1139 Consequently, one has $\partial_D F_2(1^{-}) = 0$.
Performance study of two serial interconnected chemostats

1140 Finally, it remains to calculate the limit of $\partial_r F_2(r)$ when $r$ tends to 1. One
1141 has

$$ \frac{\partial F_2}{\partial r} = f'(S_2^*) \frac{\partial S_2^*}{\partial r} x_2^* + f(S_2^*) \frac{\partial x_2^*}{\partial r}. $$

1142 Let use again the expression $G(D, r) = D(S^{in} - S_{2}^{*}(D, r))$. According to (D41),
1143 one has

$$ \overline{G}'(r) = \partial_r G(r) + \partial_D G(r) \overline{D}'(r) $$

1144 where $\partial_D G(r) = 0$. According to (D50), we have $\partial_r G(1^{-}) = 0$, and thus
1145 $\partial_r S_{2}^{*}(1^{-}) = 0$. Using expression (B18), one gets

$$ \partial_r x_2^*(1^-) = -a \overline{D}(1) \frac{S^{in} - \lambda (\overline{D}(1) + a)}{(\overline{D}(1) + a)^2}, $$

1146 and then the limit of $\partial_r F_2$ when $r$ tends to 1 exists:

$$ \partial_r F_2(1^{-}) = -a \overline{D}(1) \frac{S^{in} - \lambda (\overline{D}(1) + a)}{\overline{D}(1) + a}. $$

1147 As $\overline{D}'$ is assumed to be bounded on $V_1 \cup \{r < 1\}$, we thus obtain from
1148 (D54) the existence of $\overline{G}''(1^{-})$ with

$$ \overline{G}''(1^{-}) = -2\partial_r F_2(1^{-}) $$

1149 which is given by expression (43).

Acknowledgements

1151 The authors thank Jérôme Harmand for valuable and fruitful comments.

1152 The authors thank the Euro-Mediterranean research network Treasure (http:
1153 //www.inra.fr/treasure).

Conflict of interest

1155 The authors have no relevant financial or non-financial interests to disclose.

Funding

1157 The authors declare that no funds, grants, or other support were received
1158 during the preparation of this manuscript.

Author Contributions

1160 All authors contributed to the study conception, methodology and mathe-

1161 matical analysis. The first draft of the manuscript was written by Manel

1162 Dali-Youcef. All authors commented on previous versions of the manuscript.

1163 All authors read and approved the final manuscript.
# Data Availability

Data sharing not applicable to this article as no datasets were generated or analysed during the current study.

## References

[1] N. Abdellatif, R. Fekih-Salem and T. Sari, Competition for a single resource and coexistence of several species in the chemostat, Math. Biosci. Eng., 13 (2016), 631–652.

[2] B. Bar and T. Sari, The operating diagram for a model of competition in a chemostat with an external lethal inhibitor, Discrete & Continuous Dyn. Syst. - B, 25 (2020), 2093–2120.

[3] G. Bastin and D. Dochain, On-line estimation and adaptive control of bioreactors: Elsevier, Amsterdam, 1991.

[4] A. Bornhöft, R. Hanke-Rauschenbach and K. Sundmacher: steady state analysis of the anaerobic digestion model no. 1 (adm1). Nonlinear Dynamics 73 (2013), 535–549.

[5] M. Crespo and A. Rapaport, About the chemostat model with a lateral diffusive compartment, Journal of Optimization, Theory and Applications, Vol. 185 (2020), 597–621.

[6] M. Dali-Youcef, J. Harmand, A. Rapaport, T. Sari. Some non-intuitive properties of serial chemostats with and without mortality. 2021. https://hal.archives-ouvertes.fr/03404740

[7] M. Dali-Youcef, A. Rapaport and T. Sari, Study of performance criteria of serial configuration of two chemostats, Math. Biosci. Eng., 17(6) (2020), 6278-6309.

[8] M. Dali-Youcef and T. Sari. The productivity of two serial chemostats (2021). https://hal.inrae.fr/hal-03445797

[9] Y. Daoud, N. Abdellatif, T. Sari and J. Harmand: Steady-state analysis of a syntrophic model: The effect of a new input substrate concentration. Math. Model. Nat. Phenom. 13 (2018), 31.

[10] M. Dellal, M. Lakrib and T. Sari, The operating diagram of a model of two competitors in a chemostat with an external inhibitor, Math. Biosci., 302 (2018), 27–45.
1196

1197

1198

1199

[11] R. Fekih-Salem, Y. Daoud, N. Abdellatif and T. Sari. A mathematical model of anaerobic digestion with syntrophic relationship, substrate inhibition and distinct removal rates. *SIAM Journal on Applied Dynamical Systems* 20 (2021), 621–1654.

1200

1201

1202

[12] R. Fekih-Salem, C. Lobry and T. Sari, A density-dependent model of competition for one resource in the chemostat, Math. Biosci., 286 (2017), 104–122.

1203

1204

[13] S. Fogler: Elements of Chemical Reaction Engineering, 4th edition. Prentice Hall, New-York (2008).

1205

1206

1207

[14] C. de Gooijer, W. Bakker, H. Beeftink and J. Tramper, Bioreactors in series: an overview of design procedures and practical applications. Enzyme Microb. Technol. 18 (1996), 202–219.

1208

1209

1210

[15] I. Haidar, A. Rapaport, A. and F. Gérard, Effects of spatial structure and diffusion on the performances of the chemostat. Mathematical Bioscience and Engineering. 8(4) (2011), 953–971.

1211

1212

1213

[16] J. Harmand, C. Lobry, A. Rapaport and T. Sari, The Chemostat: Mathematical Theory of Microorganism Cultures, John Wiley & Sons, Chemical Engineering Series, 2017.

1214

1215

1216

[17] J. Harmand, A. Rapaport and A. Trofino, Optimal design of two interconnected bioreactors–some new results. AIChE J. 49(6) (1999), 1433–1450.

1217

1218

1219

1220

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

1221

1222

1223

[19] C.M. Kung and B.C. Baltzis: The growth of pure and simple microbial competitors in a moving and distributed medium. Math. Biosci. 111 (1992), 295–313.

1224

1225

[20] O. Levenspiel, Chemical reaction engineering, 3rd edition. Wiley, New York (1999).

1226

1227

1228

[21] B. Li, Global asymptotic behavior of the chemostat : general response functions and differential removal rates. SIAM Journal on Applied Mathematics 59 (1998), 411–4.

1229

1230

1231

[22] R. W. Lovitt and J.W.T. Wimpenny, The gradostat: a tool for investigating microbial growth and interactions in solute gradients. Soc. Gen. Microbial Quart. 6 (1979), 80 .
1232

1233

1234

[23] R. W. Lovitt and J.W.T. Wimpenny, The gradostat: a bidirectional compound chemostat and its applications in microbiological research, *J. Gen. Microbiol.* 127 (1981), 261—268

1235

1236

1237

[24] K. Luyben and J. Tramper, Optimal design for continuously stirred tank reactors in series using Michaelis-Menten kinetics. *Biotechnol. Bioeng.* 24 (1982), 1217–1220.

1238

1239

[25] M. Nelson and H. Sidhu, Evaluating the performance of a cascade of two bioreactors. *Chem. Eng. Sci.* 61 (2006), 3159–3166.

1240

1241

[26] S. Pavlou, Computing operating diagrams of bioreactors, J. Biotechnol., 71 (1999), 7–16.

1242

1243

1244

[27] M. Polihronakis, L. Petrou and A. Deligiannis, Parameter adaptive control techniques for anaerobic digesters—real-life experiments, Elsevier, Computers & chemical engineering, 17(12) (1993), 1167-1179.

1245

1246

1247

[28] A. Rapaport, I. Haidar and J. Harmand, Global dynamics of the buffered chemostat for a general class of growth functions, J. Mathematical Biology, 71(1) (2015), 69–98.

1248

1249

1250

[29] A. Rapaport and J. Harmand, Biological control of the chemostat with nonmonotonic response and different removal rates. *Mathematical Biosciences and Engineering* 5, no. 3 (2008), 539–547.

1251

1252

1253

[30] T. Reh and J. Muller, CO2 abatement costs of greenhouse gas (GHG) mitigation by different biogas conversion pathways. *J. Environ. Manag.* 114, no. 15 (2013), 13–25.

1254

1255

[31] T. Sari. Best Operating Conditions for Biogas Production in Some Simple Anaerobic Digestion Models. *Processes* 2022, 10, 258.

1256

1257

[32] T. Sari and B. Benyahia. The operating diagram for a two-step anaerobic digestion model. *Nonlinear Dynamics* 2021, **105**, 2711–2737.

1258

1259

1260

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

1261

1262

[34] T. Sari and F. Mazenc, Global dynamics of the chemostat with different removal rates and variable yields. *Math Biosci Eng.* 8(3) (2011), 827–40.

1263

1264

[35] T. Sari and M.J. Wade, Generalised approach to modelling a three-tiered microbial food-web, *Math. Biosci.*, 291 (2017), 21–37.

1265

1266

[36] M. Sbarciog, M. Loccufier and E. Noldus, Determination of appropriate operating strategies for anaerobic digestion systems, *Biochem. Eng. J.*, 51
[37] H. Smith, The gradostat: A model of competition along a nutrient gradient. *Microbial Ecology*, 22(1) (1991), 207–26.

[38] H. Smith, B. Tang and P. Waltman: Competition in a n-vessel gradostat. *SIAM J. Appl. Math.* 91(5) (1991), 1451–1471.

[39] H. Smith and P. Waltman, The Theory of the Chemostat, Dynamics of Microbial Competition. Cambridge University Press, 1995.

[40] B. Tang, Mathematical investigations of growth of microorganisms in the gradostat, *J. Math. Biol.*, 23 (1986), 319–339.

[41] M.J. Wade, R.W. Pattinson, N.G. Parker and J. Dolfing, Emergent behaviour in a chlorophenol-mineralising three-tiered microbial ‘food web’, *J. Theor. Biol.*, 389 (2016), 171–186.

[42] M. Weedermann, G. Seo and G.S.K. Wolkowics: Mathematical model of anaerobic digestion in a chemostat: Effects of syntrophy and inhibition. *Journal of Biological Dynamics* 7 (2013), 59–85.

[43] M. Weedermann, G.S.K. Wolkowicz and J. Sasara: Optimal biogas production in a model for anaerobic digestion. *Nonlinear Dynamics* 81 (2015), 1097–1112.

[44] G.S.K. Wolkowicz, Z. Lu, *Global dynamics of a mathematical model of competition in the chemostat*: general response functions and differential death rates. *SIAM Journal on Applied Mathematics* 52 (1992), 222–23.

[45] A. Xu, J. Dolfing, T.P. Curtis, G. Montague and E. Martin, Maintenance affects the stability of a two-tiered microbial ‘food chain’?, *J. Theor. Biol.*, 276 (2011), 35–41.

[46] J. Zambrano and B. Carlsson, Optimizing zone volumes in bioreactors described by Monod and Contois growth kinetics, *Proceeding of the IWA World Water Congress & Exhibition*, (2014).

[47] J. Zambrano, B. Carlsson and S. Diehl, Optimal steady-state design of zone volumes of bioreactors with Monod growth kinetics. *Biochem. Eng. J.* 100 (2015), 59–66.

[48] W. Walter, Ordinary Differential Equations. Springer Graduate Texts in Mathematics, 182 (1998).