Performance Study of Two Serial Interconnected Chemostats with Mortality

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.


Introduction
The mathematical model of the chemostat has received a great attention in the literature for many years [see for instance Harmand et al. (2017) 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 (Fogler 2008;Levenspiel 1999).
Several extensions of the original model of the chemostat, considering spatial heterogeneity, have been proposed to better cope reality (see for instance Kung and Baltzis (1992)). Lovitt and Wimpenny have proposed the "gradostat" experimental device as a collection of chemostats of same volume interconnected in series Wimpenny 1979, 1981), which has led to the so-called gradostat model representing in a more general framework a gradient of concentrations (Smith 1991;Tang 1986). The gradostat model has been further generalized as the "general gradostat model" representing more general interconnection graphs with tanks of different volumes Smith and Waltman 1995).
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 (Sari 2022), 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 inputoutput performances [see for instance Crespo and Rapaport (2020); Dali-Youcef et al. (2020); Haidar et al. (2011); Rapaport et al. (2015)]. 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 (de Gooijer et al. 1996;Harmand et al. 1999;Luyben and Tramper 1982;Nelson and Sidhu 2006;Zambrano et al. 2015).
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 (Li 1998;Rapaport and Harmand 2008;Sari and Mazenc 2011;Wolkowicz and Lu 1992). However, these mathematical studies have mainly concern analyses of equilibria and stability and not the performances of the system in presence of mortality.
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 Pavlou (1999) for its importance for bioreactors. Indeed, such diagrams are more and more often constructed both in the biological literature (Pavlou 1999;Sbarciog et al. 2010;Wade et al. 2016;Xu et al. 2011) and the mathematical literature (Abdellatif et al. 2016; Bar and Sari 2020;Bornhöft et al. 2013;Daoud et al. 2018;Dellal et al. 2018;Fekih-Salem et al. 2021Khedim et al. 2018;Sari 2022;Sari and Benyahia 2021;Sari and Harmand 2016;Sari and Wade 2017;Weedermann et al. 2013Weedermann et al. , 2015. Then, we study the performances in terms of conversion ratio and byproduct 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.
Along the paper, we use the abbreviations LES for locally exponentially stable and GAS for globally asymptotically stable in the positive orthant.
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. Afterward, Sect. 3 focuses on the study of performances of the serial configuration with respect of the output substrate concentration. Then, Sect. 4 considers the performances of the serial configuration with respect of the biogas production. Next, Sect. 5 is devoted to illustrations and numerical simulations and a conclusion is given in Sect. 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 (Dali-Youcef et al. 2020), but that have required to revisit significantly the mathematical proofs. Finally, Appendix D contains technical proofs.

Presentation of the Model
We consider two serial interconnected chemostats where the total volume V is divided into V 1 = r V and V 2 = (1 − r )V , with r ∈ (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 .
The mathematical model is given by the following equations: whereṠ i = dS i dt ,ẋ i = 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. The dilution rate of the first tank is Q/V 1 = D/r . The dilution rate of the second tank is Q/V 2 = D/(1−r ).
Note that these equations are not valid for r = 0 and r = 1, which correspond to a single chemostat. For sake of completeness, the useful results on the single chemostat are given in Appendix A. The considered growth function satisfies the following properties.
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), (m may be + ∞). (2) As f is increasing then the break-even concentration is defined by The particular case without mortality of the biomass (a = 0) is studied in Dali-Youcef et al. (2020). 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: • 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 not in the second one. • The steady state E 2 = (S * 1 , x * 1 , S * 2 , x * 2 ) of persistence of the species in both chemostats.
As in the case without mortality, see Table 3 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 The operating diagram of system (1) and the curve r defined by (14) under which the serial configuration is more efficient than the single chemostat (Color figure online) 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.

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, Zambrano and Carlsson (2014). 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 out r (S in , D). Proposition 1 Assume that Assumption 1 is satisfied. The output substrate concentration at steady state of system (1) is given by where S 2 = λ D 1−r + a and S * 2 is the unique solution of equation h(S 2 ) = f (S 2 ).

In this equation, the function h is defined by
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 S 2 if E 1 is the GAS steady state and to S * 2 if E 2 is GAS. According to Theorem 3 in Appendix, E 0 is GAS if and only if which is equivalent to On the other hand, using Theorem 3, S 2 depends on D and r and we have which is equivalent to Finally, using Theorem 3, we know that S * 2 depends on parameters S in , D, r . It is the unique solution of equation h(S 2 ) = f (S 2 ), where h is defined by (5). On the other hand E 2 is GAS if and only if the condition D < r ( f (S in ) − a) is satisfied, which is equivalent to the condition S in > λ D r + a .
Although S out r (S in , D) is defined only for 0 < r < 1, we can extend it, by continuity, for r = 0 and r = 1 by where S out (S in , D), which is the output substrate concentration of the single chemostat, is given by For more information on S out (S in , D), see Appendix A.
The proof of (6), comes from the following remarks. First, we have S 2 (D, 0) = λ(D + a) and second, according to Lemma 9 in the Appendix, we can extend S * 2 (S in , D, r ), by continuity, to r = 1, by Our aim in this section is to compare S out r defined by (4) and (6) and S out defined by (7).

The Serial Configuration can be More Efficient Than the Single Chemostat
We fix r and we describe the set of operating conditions (S in , D) for which that is to say, the serial configuration with volumes r V and (1 − r )V , is more efficient than the single chemostat of volume V . For r ∈ (0, 1), let g r : Lemma 1 For r ∈ (0, 1) we have g r (D) > λ D r + a .
Theorem 1 Assume that Assumption 1 is satisfied. For any r ∈ (0, 1), we have Moreover, Since f is increasing, see Assumption 1, and h is decreasing, see Lemma 8 in Appendix, then the condition S * 2 (S in , D, r ) < λ(D + a) is equivalent to the condition h(λ(D + a)) < f (λ(D + a)) = D + a. Using (5), a straightforward computation shows that the condition h(λ(D + a)) < D + a is equivalent to S in > g r (D), where g r is defined by (9). This proves (10).
Let us go now to the proof of the theorem. Assume that S in > g r (D). Using Lemma 1, we have Using (4) and (7), we have From (10), we have S out r (S in , D) < S out (S in , D). Hence, we proved the following implication Assume now that S in ≤ g r (D). When r < 1/2, three cases must be distinguished. First, if then, by (4) and (7), we obtain (11). Hence, using (10), we have S out then, by (4) and (7), we have When r ≥ 1/2, the proof is similar, excepted that we must distinguish only two cases, λ(D + a) < S in ≤ λ(D/r + a) and S in ≤ λ(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 equalities.
Theorem 1 asserts that the serial configuration is more efficient than the single chemostat if and only if S in > g r (D). Let us illustrate this result in the operating diagram of system (1). Consider the curve of equation According to the results given in Appendix C, the curves r and 1−r defined by (13) separate the operating plane (S in , D) in four regions in which the system has different asymptotic behavior, see Table 3. 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 ) ∪ I 3 (r ), E 3 is GAS, see Fig. 2. This figure also shows the plot of the curve r , defined by Using Lemma 1, we see that for all r ∈ (0, 1), the curve r is always at right of the curve 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 in , D) is at right of the curve r depicted in Fig. 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 Proposition 2 Assume that Assumption 1 is satisfied. Let D > 0, S in > λ(a). We 1. If S in ≤ λ(D + a), then for any r ∈ [0, 1], one has S out r (S in , D) = S out (S in , D) = S in . 2. If λ(D + a) < S in < λ(2D + a), then 1/2 < r 0 < 1 and one has 3. If λ(2D + a) ≤ S in , then 0 < r 0 ≤ 1/2 and one has Here Proof If S in ≤ λ(D + a), then, for all r ∈ (0, 1), one has Then, according to (4), one has S out r (S in , D) = S in . This proves item 1 of the proposition.
We want to determine the values r ∈ (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 Sect. 3.4. On the other hand, using L'Hôpital's rule, we have
We use the notation r 1 (·, D) to recall the dependence of the inverse function in D. For all D ∈ (0, m − a), r ∈ (D/(m − a), 1) and S in > g(D), we have Theorem 2 Assume that Assumptions 1 and 2 are satisfied. Let g defined by (19).
• If S in ≤ g(D) then for any r ∈ (0, 1), we have S out r (S in , D) > S out (S in , D). In addition, for r = 0 and r = 1 we have S out (20). In addition, for r = 0, r = r 1 (S in , D) and r = 1, Proof The function r → g r (D) is decreasing and tends to g(D) as r tends to 1, as shown by (18). Thus, for all r ∈ (0, 1), we have g(D) < g r (D). If S in ≤ g(D), then S in < g r (D). According to Theorem 1, for all r ∈ (0, 1), we have S out r (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 out r (S in , D) < S out (S in , D). The equality S out r (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 1 (D), see (20). Hence, according to Theorem 1, we have S out Let us now describe the subsets of the operational space (S in , D) for which the behavior 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 1 and 1/2 , defined by (13), and the curve defined by These three curves intersect at (λ(a), 0) and, using the inequality g(D) > λ(D + a), which is satisfied for all D > 0, one deduces that is at the right of 1 . Therefore, the curves 1 , 1/2 and separate the set of operating parameters (S in , D) into the Combining the results of Proposition 2 and Theorem 2, we find that the function r → S out r (S in , D) is as in Fig. 4. In the following we will comment on this figure. Note that if (S in , D) ∈ 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) ∈ J 2 ∩ 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 1/2 and . For a linear growth rate, 1/2 = and hence, regions J 2 and J 4 do not exist, see Fig. 7a. For a Monod growth function, curves 1/2 and intersect only at point (λ(a), 0) and hence, region J 2 always exist and is connected but region J 4 does not exist, see Fig. 8a. For a Hill growth function, curves 1/2 and always intersect at (λ(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. 9a-c.

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 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 Sect. 3.4.
Hence, its admits an inverse function Proposition 3 Assume that Assumptions 1 and 3 are satisfied. We have

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 Sect. 5. For this purpose we consider the function γ defined by defined on which consists simply in considering g r (D), given by (9), as a function of both variables r and D. If then Assumption 3 is satisfied. The following Lemmas give sufficient conditions, for partial derivatives of γ to have their signs as indicated above. (28) is satisfied. c) If f is twice derivable, then l D is twice derivable and the following conditions are equivalent Proof Notice first that γ (r , D) can be written as follows Using the definition of l D , γ (r , D) is given then by The partial derivative, with respect to r of γ , is given then by (28) is satisfied, and, in addition 0 < r ≤ 1/2, then, from (30), it is deduced that ∂γ ∂r (r , D) < 0. In the case r ∈ (1/2, 1), we use the following expression of γ (r , D) which is deduced from (29): where Hence, from (31), it is deduced that ∂γ ∂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 Taking s = 1 and r ∈ 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 we can write Therefore, the condition 1 in item c in the lemma is equivalent to the following condition: Using the notation S = λ D r + a , which is the same as is the condition 2 in c in the lemma.
Lemma 3 Assume that (32), this condition is equivalent to (34). Note that if f is decreasing, then this condition is satisfied. Indeed, we have which is the condition (34).

Remark 1 Notice that:
i) The condition 2 in part c of Lemma 2 is equivalent to the condition 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) ≤ 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 existŝ S ∈ (0, +∞) such that f is non-negative on (0,Ŝ) and non-positive on (Ŝ, +∞).
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 ∈ (λ(a),Ŝ). Therefore, if (1/ f ) > 0 on (0,Ŝ), then Assumption 2 is satisfied.
We will see in Sect. 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. 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 (Reh and Muller 2013).
We recall that the biogas flow rate is proportional to the microbial activity, as defined for instance in Bastin and Dochain (1991), Polihronakis et al. (1993). 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 , the biogas flow rate corresponding to the steady state E 1 is given by the expression , the biogas flow rate corresponding to the positive steady state E 2 is given by the expression Proof From system (B10), considering equationṠ 2 = 0, one obtains From system (B10), consideringṠ 1 = 0 andṠ 2 = 0 gives, respectively, This ends the proof of the proposition.
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 out r in (6). We can write represents the biogas flow rate of the single chemostat when 0

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 chemostat. 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 chemostat. More precisely, we have the following result.
where G 2 is given by (39) and g r is defined by (9).
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). 2. According to Theorem 1, for any r ∈ (0, 1) and . Consequently, using the formula for G 2 given in Proposition 4, one has (D). If Assumption 2 is satisfied, then, using (21), we see that This ends the proof of the proposition.
Let S in and D be fixed. The graphs of the biogas flow rates functions For r and S in fixed, the curves of the maps where G 1 , G 2 and G chem are given by (38), (39) and (40), respectively (Color figure online) are easily obtained from the graph of the output substrate concentration, r → S out r (S in , D), see Fig. 4. Indeed, the formulas given in Proposition 4 show that whenever these functions are defined, we have We will see, in Sect. 5, some illustrative plots of the biogas flow rates G 1 and G 2 as functions of the parameter r ∈ [0, 1], for linear growth, see Let us illustrate the result of Proposition 5 by plotting the graphs of the biogas flow rates when r and S in are fixed, see Fig. 5. This figure is made without graduations on the axes because its 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 ∈ (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.

The Maximal Biogas of the Serial Configuration Can Exceed That of the Single Chemostat
In Fig. 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 Indeed we will show that there is a value r * ∈ (0, 1) such that this inequality is true for all r ∈ (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 Dali-Youcef et al. (2020), that for all S in > 0, and all r ∈ (0, 1) we have 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 > λ(a) and r ∈ (0, 1], the continuous function 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 > λ(a), we then consider the function We want to ensure that this maximum is reached at a single value, denoted D(r ). Note that D(1) represents the value, which we will assume to be unique, at which the function D → 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 > λ(a), there exists 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 → G chem (S in , D) is unique. The following lemma shows that the function D → 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 > λ(a), there exists a neighborhood V 1 of 1, such that for any r ∈ V 1 ∩{r ≤ 1}, the maximum of the function Proof The proof is given in in Appendix D.2.
Proposition 6 Under Assumption 4, the function G admits left limits of its first and second derivatives at r = 1, which are Proof The proof is given in Appendix D.3.
The function r → G(r ) reaches its maximum at some r max ∈ (r * , 1). Let D max = D(r max ) be the maximum of the function D → 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 We have plotted the function r → 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 G (1) = 0. In addition, one remarks that G(r ) > G(1) for r in some interval (r * , 1) and G(r * ) = 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 = D(r ), where 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 * = D(r * ).
Indeed, for the choice r = r * and D = D * , we have but D * is expected to be significantly smaller than D(1), the dilution rate that maximizes biogas for the single 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 single 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 Sect. 5.4. For biological comments on the heuristic underlying this non-intuitive behavior, the reader is referred to Dali-Youcef et al. (2021).

Illustrations and Numerical Simulations
This section illustrates our results using three different growth functions. As concave functions, we choose the linear growth function and the Monod function. As a nonconcave function, we choose the Hill function.

Linear Growth Function
Let consider a linear function f (S) = αS, α > 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.
Consider the operating points η 1 and η 3 , fixed, respectively, in regions J 1 and J 3 , as shown in Fig. 7a. The behavior of the biogas flow rate for these operating points is depicted in Fig. 7b In the linear 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 η 3 = (0.5, 0.6) 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 ∈ (r 1 , 1), where r 1 (0.5, 0.6) ≈ 0.5, see Fig. 7b.

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.
Hence, the curve is at left of the curve 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)  Consider the operating points η 1 , η 2 and η 3 , fixed, respectively, in regions J 1 , J 2 and J 3 , as shown in Fig. 8a. The behavior of the biogas flow rate for these points is depicted in Fig. 8b-d. It should be noticed that for any other point (S in , D) ∈ J 1 (resp. (S in , D) ∈ J 2 and (S in , D) ∈ J 3 ), the curve representing the biogas flow rate with respect to r should be similar to the curve shown in Fig. 8b (resp. 8c, d), and corresponding to (S in , D) = η 1 (resp. (S in , D) = η 2 and (S in , D) = η 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).

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 (34) and (35). Therefore, according to item iii in Remark 1, it verifies Assumption 2 and according to Lemma 3, it satisfies Assumption 3.

Proposition 8 The Hill function satisfies the conditions
Proof Let us first prove that the Hill function satisfies the condition (35). Straightforward computation gives . Therefore, d 2 dS 2 1 f (S)−a > 0 for all S > λ(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Ŝ ∈ (0, +∞) be the inflexion point of the Hill function f . It is sufficient to show that (1/ f ) > 0 for all S ∈ (0,Ŝ). One easily see that for any S > 0. Consequently, for all p > 1, the Hill function verifies Assumption 2.
Let us now prove that the Hill function verifies the condition (34). Straightforward computations give Therefore, Since p > 1, D + ra < D + a and From (45) On the other hand, we have Therefore, using (44), (47) and (48) one obtains D + a)).
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 ). 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 Fig. 9a-c.
Recall that r 1 (S in , D) is defined by (20). It is obtained by solving numerically the equation S in = g r (D). Since the point η 2 (resp. η 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 ∈ (r 1 , 1), with r 1 (9, 1.6) ≈ 0.81 in Fig. 9e and r 1 (12, 1.6

The Serial Configuration is Worth Considering When Mortality is not Negligible
In this section we numerically illustrate Remark 2. We fix S in , and we adopt the following notations.  where G(r ) is defined by (42) and D(r ) is as in Lemma 4. Recall that r * ∈ (0, 1) satisfies and G(r ) > G(1) for r ∈ (r * , 1), so that G(r ) attains its maximum for r = r max ∈ (r * , 1), see Fig. 10a, obtained with a Monod function and S in = 1.5. We adopt the following notations.
These notations are illustrated in Fig. 10a,b. The zoom in Fig. 10b shows that G max exceeds G * = G max chem only slightly, but D * is significantly smaller than D max , which is itself smaller than D max chem . We give in Table 1 the numerical values of r * , r max , D * , D max , G max and G * = G max chem , for various values of the mortality rate a. The table also gives the relative gains 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. 10c, 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.

Conclusion
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 Dali-Youcef et al. (2020), 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 is 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.
Along the paper, the similarities, specificities and differences of our model compared to the model without mortality (i.e., for a = 0) studied in Dali-Youcef et al. (2020) 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 Dali-Youcef and Sari (2021).
Acknowledgements The authors thank Jérôme Harmand for valuable and fruitful comments. The authors thank the Euro-Mediterranean research network Treasure (http://www.inra.fr/treasure).

Author Contributions
All authors contributed to the study conception, methodology and mathematical analysis. The first draft of the manuscript was written by Manel Dali-Youcef. All authors commented on previous versions of the manuscript. All authors read and approved the final manuscript.
Funding The authors declare that no funds, grants, or other support was received during the preparation of this manuscript.
Data Availability Data is sharing not applicable to this article as no datasets were generated or analyzed during the current study.

Conflict of interest
The authors have no relevant financial or non-financial interests to disclose.

Appendix A The Single Chemostat
In this section, we give a brief presentation of the mathematical model of the single chemostat with mortality rate.
The mathematical equations are given bẏ where S and x denote, respectively, the substrate and the biomass concentration, S 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 Harmand et al. (2017); Smith and Waltman (1995)] that, besides the washout F 0 = (S in , 0), this system can have a positive steady state See Fig. 11a for the plot of the function D → S * (D) and Fig. 11(b) for the plot of the The washout steady state F 0 always exists. It is GAS if and only if D ≥ δ. It is LES if and only if D > δ. The positive steady state F 1 exists if and only if D < δ. It is GAS and LES whenever it exists. Therefore, the curve defined by splits the set of operating parameters (S in , D) into two regions, denoted I 0 and I 1 , as depicted in 11c. These regions are defined by The behavior of the system in each region is given in Table 2. Figure 11c, together with 2 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 Harmand et al. (2017). Thus, with presence of mortality rate, the region where the washout is GAS, is larger.
The output substrate concentration of the single chemostat, at its stable steady state is given by Its output biomass concentration at steady state is then given by Thus, for any growth function satisfying Assumption 1 the function D → S out (S in , D) is increasing on [0, δ], as shown in Fig. 11a. The function D → x out (S in , D) is illustrated in Fig. 11b for a Monod function.
The biogas flow rate of the single chemostat is defined, up to a multiplicative yield coefficient, by Using the expressions (A4) and (A5), respectively, of S out and x out , the biogas flow rate of the single chemostat is given by: For a given S in > λ(a), the function D → G chem (S in , D) is null for D = 0 or D ≥ δ, and is positive for D ∈ (0, δ). Therefore it admits a maximum in (0, δ), which is assumed to be unique. A characterization of the growth functions for which this uniqueness is satisfied can be found in Sari (2022). Notice that the function g defined by (A8) is the same as the function g, defined by (19), which was obtained as the limit, when r tends to 1, to the function g r , defined by (9). Recall that is the curve of equation S in = g(D), see (22). This curve is depicted in Fig. 11c. It is the set of operating conditions given the higher biogas of the single chemostat. More precisely, for any S in > λ(a), the maximum D = D(S in ) of the biogas satisfies the condition (S in , D) ∈ .

Proposition 9 Assume that for any S in > λ(a), the maximum of D
Therefore, a sufficient condition for the uniqueness of D(S in ) is that the mapping g is increasing. If, in addition, f is C 2 , then, deriving (A9) with respect of D, we have Hence, a sufficient condition for Assumption 4 to be satisfied is that g (D) > 0 for D ∈ [0, m − a). This last condition if satisfied whenever f ≤ 0 on (λ(a), +∞), or 1 f −a > 0 on (λ(a), +∞), see Lemma 1 in Sari (2022). Therefore we can make the following remark.

Appendix B The Serial Configuration
We consider a slight extension of system (1) with different mortality rates in the two tanks. Indeed, we assume that the growth environment differs from one tank to another one. This can lead to two different mortality rates in the tanks. We denote by a 1 and a 2 the mortality rates. The mathematical model is given by the following equations.
The following result is classical in the mathematical theory of the chemostat.
Proof Since the vector field defined by (B10) is C 1 , the uniqueness of the solution to an initial value problem holds. From (B10) and using f (0) = 0, we have: Therefore, for i = 1, 2, S i (t) ≥ 0 and x i (t) ≥ 0, for all t ≥ 0, for which they are defined, provided S i (0) ≥ 0 and x i (0) ≥ 0, for i = 1, 2, see Prop. B.7 in Smith and Waltman (1995). This proves that the solutions of non-negative initial conditions are always non-negative. Let z i = S i + x i , i = 1, 2. From system (B10), we havė Consequently, we have the differential inequalitẏ It follows by comparison of solutions of ordinary differential equations (see, for instance, Walter (1998)) that one has Therefore, z 1 (t) ≤ Z 1 , where Z 1 = max(S in , z 1 (0)). Then, we also have the differential inequalityż It follows again by comparison of solutions of ordinary differential equations that one has also z 2 (0)). Hence, the solutions of (B10) are positively bounded. Therefore, they are defined for all t ≥ 0.
For the description of the steady states, we shall consider the following function h that will play a key role This function satisfies the following property.

Lemma 8 Assume that D/r +a 1 < f (S in ). The function h is decreasing from h
Hence, b is a convex combination of S in and S * 1 , and we have S * 1 < b < S in . Therefore, the vertical asymptote Hence, we have h (S 2 ) < 0 for all S 2 < b. Therefore, h is defined on the interval (0, S * 1 ) and is decreasing from h(0), given by (B12) to h(S * 1 ) = 0.
admits a unique solution, denoted by S * 2 (S in , D, r ), as shown in Fig. 12a. This solution satisfies the following property.
The existence and stability of steady states of (B10) are given by the following result.
Theorem 3 Assume that Assumption 1 is satisfied. The steady states of (B10) are: • The washout steady state E 0 = (S in , 0, S in , 0) which always exists. It is GAS if and only if

It is LES if and only if
• The steady state E 1 = (S in , 0, S 2 , x 2 ) of washout in the first chemostat but not in the second one with

It is LES if and only if
• The steady state E 2 = (S * 1 , x * 1 , S * 2 , x * 2 ) of persistence of the species in both chemostats with and S * 2 = S * 2 (S 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 in ) − a 1 ). It is GAS and LES whenever it exists and is positive.
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 mathematical theory of the chemostat (Harmand et al. 2017;Smith and Waltman 1995). Using this cascade structure, the global behavior of the system is deduced from the global behavior of planar systems and Thieme's theory of asymptotically autonomous systems, all the solutions being bouded.
For the convenience of the reader the details of the proof are given in Appendix D.1.

Proposition 10
The function S in → S * 2 (S in , D, r ) is decreasing.

Proof
for all S 2 ∈ (0, S * 1 ). Therefore, S * 1 2 < S * 2 2 , see Fig. 12b. This result means that the effluent steady-state concentration of substrate decreases when the influent concentration of substrate increases. This behavior is quite different from the single chemostat, where the effluent steady-state substrate concentration is independent of the influent substrate concentration.

Appendix C Operating Diagram
For the chemostat model, 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. Indeed, the main interest of an operating diagram is to highlight the number and stability of the steady states for a given pair of parameters (S in , D). The input substrate concentration S in and the dilution rate D are the usual parameters manipulated by the experimenter of a chemostat. Apart from these parameters, and the parameter r that can be also chosen by the experimenter but not easily changed as S in and D, all other parameters have biological meaning and are fitted using experimental data from real measurements of concentrations of micro-organisms and substrates. Therefore the operating diagram is a bifurcation diagram, quite useful to understand the possible behaviors of the solutions of the system from both the mathematical and biological points of view.
Here, we fix r ∈ (0, 1) and we depict in the plane (S in , D) the regions in which the solution of system (B10) globally converges toward one of the steady state E 0 , E 1 or E 2 . From the results given in Theorem 3, it is seen that these regions are delimited by the curves 1 r and 2 1−r defined by: (e) (f) Fig. 13 The operating diagram of (B10). The asymptotic behavior in each region is depicted in Table 3 (color figure online) When a 1 = a 2 = 0, as we have shown in Dali-Youcef et al. (2020), these curves croos only at one point (the origin) and merge when r = 1/2. Therefore, in this case the curves 1 r and 2 1−r separate the operating plane (S in , D), in only three regions, see (Dali-Youcef et al. 2020, Figure 5). This property continues to hold when a 1 = a 2 , that is to say, the curves intersect only at (λ(a 1 ), 0) and merge when r = 1/2. In this case the curves 1 r and 2 1−r separate the operating plane (S in , D), in only three regions, see Fig. 13c, d. The novelty when a 1 and a 2 are different and non-null, is that the intersection of the curves 1 r and 2 1−r can lie outside the S in axis. Therefore there can be four regions in the operating plane, as depicted in Fig. 13a, f. For the description of the intersection of the curves 1 r and 2 1−r , we need some definitions and notations. Let r ∈ (0, 1) be defined by Note that if a 1 < a 2 , then r < 1/2, and if a 1 > a 2 , then r > 1/2. For a 1 < a 2 and 0 < r < r (or a 1 > a 2 and r < r < 1), we define the point P = S in P , D P of the operating plane by: Note that S in P > 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 ∈ (0, r ), the curves 1 r and 2 1−r intersect at the point P and 1 r is strictly below [resp. above] 2 1−r for S in > S in P [resp. S in < S in P ], see Fig. 13a. For all r ∈ (r , 1), 1 r is strictly above 2 1−r , see Fig. 13b. 2. If a 1 > a 2 then for all r ∈ (r , 1), the curves 1 r and 2 1−r intersect at the point P and 1 r is strictly above [resp. below] 2 1−r for S in > S in P [resp. S in < S in P ], see Fig. 13f. For all r ∈ (0, r ), 1 r is below 2 1−r , see Fig. 13e. 3. If a 1 = a 2 then, for r = 1/2, 1 r = 2 1−r . Moreover, if r < 1/2 then 1 r is strictly below 2 1−r , see Fig. 13c and, if r > 1/2, then 1 r is strictly above 2 1−r , see Fig. 13d.
Proof For 0 < r < 1 and S in > λ(a i ) we define the function ϕ i , i = 1, 2, by The curves 1 r and 2 1−r , defined, respectively, by (C19) and (C20), intersect if and only if there exists r ∈ (0, 1) and S in > max (λ(a 1 ), λ(a 2 )) such that ϕ 1 (S in , r ) = ϕ 2 (S in , r ), that is to say This equation has a solution S in > max (λ(a 1 ), λ(a 2 )) if and only if where m = sup( f ), as in (2). When these conditions are satisfied, the solution of (C24) is given by S in = λ(A(r )), where λ is the inverse function of f , i.e., the break-even concentration defined by (3). Hence, S in = S in P , given in (C22). The corresponding intersection point of 1 r and 2 1−r is given by D P = r f (S in P ) − a 1 , which is the value given in (C22).
Let us determine now for which value of r , the conditions (C25) are satisfied. The function A is a homographic function. Its graphical representation is a hyperbola, whose vertical asymptote is r = 1/2. Its derivative is given by Note that A(r ) = m if and only if r = r , where r is defined by (C21). Therefore if a 1 < a 2 , then, according to (C26), A is increasing. Since A(0) = a 2 , A(r ) = m, and r < 1/2, the condition (C25) is satisfied if and only if 0 < r < r . Similarly, if a 1 > a 2 , then, according to (C26), A is decreasing. Since A(1) = a 1 , A(r ) = m and r > 1/2, the condition (C25) is satisfied if and only if r < r < 1. Finally, if a 1 = a 2 then A(r ) = a 1 and the condition (C25) cannot be satisfied. Suppose that a 1 < a 2 . Note that for 0 where S in P is defined by (C22). Hence, the curves 1 r and 2 1−r intersect at P = (S in P , D P ) and the curve 1 r is strictly below [resp. above] the curve 2 1−r , for all S in > S in P [resp. S in < S in P ]. • If r ∈ [r , 1/2), then f (S in ) < A(r ) for all S in > 0, so that the curve 1 r is strictly above the curve 2 1−r . • If r ∈ [1/2, 1), then, using r ≥ 1−r and a 1 < a 2 , one has ϕ 1 (S in , r ) > ϕ 2 (S in , r ).
Therefore, the curve 1 r is strictly above the curve 2 1−r . If a 1 > a 2 , the proof is similar to the case a 1 < a 2 . If a 1 = a 2 , then ϕ 1 (S in , r ) = ϕ 2 (S in , r ) is equivalent to r ( f (S in ) − a 1 ) = (1 − r )( f (S in ) − a 1 ). Therefore, r = 1 − r , that is r = 1/2. In this case the curves 1 r and 2 1−r merge. In addition, if r < 1/2 [resp. r > 1/2], then r < 1 − r [resp. r > 1 − r ] and the curve 1 r is strictly below [resp. above] the curve 2 1−r . This ends the proof of the proposition.
For any r ∈ (0, 1), the curves 1 r and 2 1−r , defined by (C19) and (C20), respectively, split the plane (S in , D) in the regions denoted I 0 (r ), I 1 (r ), I 2 (r ) and I 3 (r ) defined in Table 3. These regions are depicted in Fig. 13 in the cases a 1 < a 2 , a 1 = a 2 and a 1 > a 2 .
The behavior of the system in each region, when it is not empty, is given in Table 3. Notice that E 1 exists in both regions I 1 (r ) and I 2 (r ), but is stable only when (S in , D) is in I 1 (r ).
When a 1 = a 2 = 0, λ(a 1 ) = λ(a 2 ) = 0 and the curves 1 r and , 2 1−r of the operating diagram start from the origin of the plane (S in , D) and merge for r = 1/2. Therefore, the diagrams shown in panels (a), (b), (c), (d), (e) and (f) of Fig. 13 are reduced to only two different cases characterized by 0 < r < 1/2 and 1/2 < r < 1, Table 3 The regions I k (r ), k = 0, 1, 2, 3 of the operating diagram of (B10) and asymptotic behavior in these various regions as shown in Figure 5 of Dali-Youcef et al. (2020). There is no change in the stability of the steady states and in the number of the regions depicted in the operating diagram. This result reveals an interplay between spatial heterogeneity (the ratio r of volume distribution between tanks) and the mortality heterogeneity (difference between a 1 and a 2 ). Indeed, panels (a) and (f) of Fig. 13 bring a particular feature when mortality rates are different: domains I 1 (r ) and I 3 (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 1 (r ) and I 3 (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 in , D and r , panels (a) and (f) of Fig. 13 show that it is theoretically possible to pass from domain I 3 (r ) to I 2 (r ) when mortality parameter is diminished only in the second tank. In practice, being in domain I 2 (r ) might be more desirable than I 3 (r ) with respect to some dysfunctioning of the first tank that can drop suddenly its biomass to zero. Indeed, in I 2 (r ), the second tank is no conducted to the wash-out differently to the I 3 (r ) case.
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. 13 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 = a 2 is similar.

Remark 4 Transcritical bifurcations occur in the limit cases
, 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 in )−a) and a transcritical bifurcation of E 1 and E 2 when D = r ( f (S 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 in )−a) and a transcritical bifurcation of E 0 and E 2 when D = r ( f (S in ) − a). If r = 1/2 and D = ( f (S in ) − a)/2 then, we have transcritical bifurcations of E 0 and E 1 , and E 0 and E 2 , simultaneously.
of C is positive, and its trace is negative, the eigenvalues of C have negative real parts. Therefore, E 1 is LES if and only if the condition in the theorem is satisfied.
For E 2 , the determinant of A is positive and its trace is negative. On the other hand, using the notation C E 2 for the matrix C evaluated at E 2 , we have Note that h(S 2 ) < D/(1 − r ) + a 2 for all S 2 ∈ (0, S * 1 ). Therefore, from (B11), we have f (S * 2 ) = h(S * 2 ) < D/(1 − r ) + a 2 . Consequently, det(C E 2 ) and tr(C E 2 ) are, respectively, positive and negative. Therefore, E 2 is LES whenever it exists, that is 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 Harmand et al. (2017)). 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 0 1 , x 0 1 , S 0 2 , x 0 2 ). Then, (S 1 (t), x 1 (t)) is the solution of systeṁ with initial condition (S 0 1 , x 0 1 ) and (S 2 (t), x 2 (t)) is the solution of the non-autonomous system of differential equationṡ with the initial condition (S 0 2 , x 0 2 ). The system (D28) is the classical model of a single chemostat. Its asymptotic behavior is well known (see, for instance, Proposition 2.2 of Harmand et al. (2017)). This system admits the steady states: where S * 1 and x * 1 are defined by (B17). Two cases must be distinguished. Firstly, if λ (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 non-negative quadrant. Hence, for any non-negative initial condition (S 0 1 , x 0 1 ), Therefore, the system (D29) is asymptotically autonomous with the limiting systeṁ Recall that the solutions of (D29) are positively bounded. Therefore, we shall use Thieme's results which apply for bounded solutions. The system (D32) represents the classical model of a single chemostat. It admits the two steady states e 2 0 = (S in , 0) and e 2 1 = (S 2 , x 2 ), with (S 2 , x 2 ) defined by (B15). Two subcases must be distinguished.
, then e 1 1 , defined in (D30), is GAS for (D28) in the positive quadrant. Hence, for any positive initial condition Therefore, the system (D29) is asymptotically autonomous with the limiting systeṁ 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 ξ = ln(x 2 ). The system (D34) becomes aṡ The divergence of the vector field where h is the function defined in (5). According to Lemma 8, h is positive decreasing on (0, λ(D/r +a), and h− f admits an unique zero S 2 = S 2 (S in , D, r ) on (0, λ(D/r + a). Then, one can write For r ∈ [1, 1 + ε) and (D, r ) ∈ D ε , on has which is negative for any S 2 ∈ (0, S in ) thanks to condition (D37). As ϕ(0, D, r ) > 0 and ϕ(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 ϕ(S 2 , D, r ) = 0, which also verifies ∂ S 2 ϕ < 0 at S 2 = S 2 . Then, by the implicit function theorem, the function (D, r ) → S 2 (S in , D, r ) is C 2 on D ε . Recall that for r < 1 and D < r δ, on has the expression G 2 (S in , D, r ) = V D(S in − S 2 (S in , D, r )) (see Proposition 4). We extend now the function (D, r ) → G 2 (S in , D, r ) with this last C 2 expression on D ε . As G 2 (S in , D, 1) = G chem (S in , D) for any D ∈ (0, δ), one deduces, by continuity of the partial derivatives of G 2 with respect to D and property (D36), the existence of V D , V r as neighborhoods, respectively, of D(1) and 1 with V D × V r ⊂ D ε such that for any r ∈ V r , the function D → G 2 (S in , D, r ) possesses the following properties 1. it is strictly concave on V D , 2. it is increasing on (D ε , D(1)) \ V D and decreasing on (D(1), r δ) \ V D , 3. its maximum over (0, r δ) is not reached for D ≤ D ε .
We thus deduce that D → G 2 (S in , D, r ) admits a unique maximum D(r ) on (0, r δ), for any r ∈ V r .
Finally, for any r ∈ V r , D(r ) is characterized as the zero of the map D → F(D, r ) where F is the C 1 function F(D, r ) := ∂ D G 2 (S in , D, r ) From property 1. above, one obtains ∂ D F(D(r ), r ) = ∂ 2 D D G 2 (S in , D(r ), r ) < 0, r ∈ V r and by the implicit function theorem, there exists a neighborhood V 1 ⊂ V r of 1 such that D is C 1 on V 1 , which ends the proof of the lemma.
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 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 ) =F 1 (D/r ) whereF 1 is a smooth function. Therefore, one has As ∂ D F 1 (1) = 0 then one deduces ∂ r F 1 (1) = 0.
With (D44), (D46) and (D49), expression (D43) gives the existence of the limit of G when r tends to 1 with r < 1, which is On the one hand, using L'Hôpital's rule, one has Recall that ∂ r F 1 (1) = 0 and thus one has F 1 (1) = 0. Consequently, one has On the other hand, using (D42) and (D45), one has r 1 − r ∂ r F 1 (r ) = D(r ) r ∂ D F 2 (r ).
Thus, according to (D51), (D52) and (D53), one gets Let us show now that the limit of ∂ D F 2 (r ) is 0 when r tends to 1. One has Let us use the expression G(D, r ) = D(S in − S * 2 (D, r )) given by Proposition 4. As D(r ) is a maximizer, then one has ∂ D G(r ) = S in − S * 2 (r ) − D(r )∂ D S * 2 (r ) = 0.