Effects of spatial structure and diffusion on the performances of the chemostat

Given hydric capacity and nutrient flow of a chemostat-like system, we analyse the influence of a spatial structure on the output concentrations at steady-state. Three configurations are compared: perfectly-mixed, serial and parallel with diffusion rate. We show the existence of a threshold on the input concentration of nutrient for which the benefits of the serial and parallel configurations over the perfectly-mixed one are reversed. In addition, we show that the dependency of the output concentrations on the diffusion rate can be non-monotonic, and give precise conditions for the diffusion effect to be advantageous. The study encompasses dead-zone models


Introduction
The chemostat is a popular apparatus, invented simultaneously by Monod [20] and Novick & Szilard [23], for the so-called continuous culture of micro-organisms.It has the advantage to study bacteria growth at steady state, in contrast to batch cultivation.In the classical experiments, the medium is assumed to be perfectly mixed, that justifies mathematical models described by systems of ordinary differential equations [29].The chemostat model is also used in ecology for studying populations of micro-organisms, such as lake plankton or wetlands ecosystems.In natural ecosytems, or in industrial applications that use large bioreactors, the assumption of perfectly mixed medium is questionable.This is why spatial considerations have been introduced in the classical model of the chemostat, such as the gradostat model [16] that is a series of interconnected chemostats (of identicalvolumes).Segregated habitats are also considered in lakes, where the bottom can be modeled as a dead zone and nutrient mixing between the two zones is achieved by diffusion rate [21].The consideration of dead zones is also often used in bioprocesses modelling [15,14,6,26,25,31,27].
Series of chemostats, instead of single chemostat, have shown to potentially improve the performances of bioprocesses, reducing the total residence time [13,17,10,11,12] or allowing species persistence [30,24].These properties have of course economical impacts for the biotechnological industry, and there is a significant literature on the design of series of reactors and comparison with plug-flow reactors (that can be seen as the limiting case of an arbitrary large number of tanks of arbitrary small volumes) [32,1,22,2,3,4,5].Sometimes a radial diffusion is also considered in plug-flow reactors [7], but surprisingly, configurations of tanks in parallel have been much less investigated [15].One can argue that knowing input rates and volumes of tanks in parallel, there dynamical characteristics can be studied separately, and there is no need of devoting a specific study for these configurations.This is no longer the case if one considers a passive communication between the tanks, through a membrane for instance.In saturated soils or wetlands, a spatial structure could be simply represented by separated domains with diffusive communication.This consideration is similar to patches models or islands models, commonly used in ecology [18,9], or lattice differential equations [28].For instance, a recent investigation studies the influence of such structures on a consumer/resource model [8].Consumer/resource models in ecology are similar to chemostat models, apart the source terms that are modeled as constant intakes of nutrient, instead of dilution rates that one rather met in liquid media.
In this paper, we propose to bring new insight on parallel configurations of chemostats with communication, in a spirit different than the one usually taken in bioprocesses design.One usually chooses a target for the output concentration of substrate, and looks for minimizing the total volume, or equivalently the residence time, among all the configurations that provide the same desired output at steady state.Here, we fix both the total hydric volume and the input flow and study the input-output map at steady-state, investigating the role of the spatial structure on the performances of the system.The performance is here measured by the level of substrate that is degraded by the system, and collected at the output.We draw precise comparisons between the three configurations: perfectly mixed, serial and parallel (with diffusion rate) with the same total hydric volume and flow rate.This set of configurations is far to be exhaustive, being limited to two compartments only, but it is a first attempt to grasp this input-output map of a structured chemostat, and study how a spatial structure can modify this map, and what are the key parameters.We believe that this study is of interest for the modelling of ecosystems such as saturated soils for which it is not easy to know the spatial structure, and where one has only access to input-output observations of the substrate degradation.
The paper is organized as follows.In Section 2, we present the three configurations under investigation and give the equations of the models.The main part of the paper is devoted to the analysis of the steady states, given in Section 3. The proofs of the global stability of the equilibriums are postponed to the Appendix, for lightening the presentation.Finally, discussion and numerical simulations are given in Section 4.

The models
The flow rate is labeled Q and V is the total capacity of the system.The three simple patterns we analyze are depicted on Figure 1: We recall the dynamical equations of resource (nutrient) and biomass concentrations, respectively denoted by S i and X i in a compartment i of volume V i fed from a compartment i − with a flow rate Q i and connected by diffusion rate d to a compartment i d (see Figure 2).
For sake of simplicity of the analytical analysis, we assume that the growth function µ(•) is a linear function of the resource concentration: In Section 4, we shall consider Monod growth function and show that the qualitative results of our study are not changed.The yield coefficient y of the bio-conversion is kept equal to one (this is always possible by choosing the unit measuring the biomass).It is convenient to write dimensionless concentrations: for each concentration C i in the compartment i (C i can denote S i or X i ), we define We shall also consider that the time t is measured in units such that Q = V .Finally, we assume that the input concentration S in is large enough to avoid the (trivial) wash-out equilibrium to be the only steady-state in each compartment.
3 Steady-state analysis of the three configurations

Configuration with one compartment
The dynamical equations of the configuration with a single compartment are The non-trivial equilibrium is (1, s in − 1) under the condition s in > 1.Then, one has Remark.This is a well known property from the theory of the chemostat that the output concentration at steady state is independent of the input concentration, provided this latter to be large enough (i.e.s in ≥ 1).

Serial connection of two compartments
The dynamical equations of the model with two compartments in series (see Figure 1), assuming r to be different to 0 and 1, are Proposition 1.When s in > 1/r, there exists an unique equilibrium (s of (1) on the positive orthant.One has necessarily s ⋆ 1 = 1/r and s ⋆ 2 < min(1/r, 1/(1 − r)).Furthermore, one has Proof.One can readily check that there exists a non-trivial equilibrium (1/r, s in − 1/r) for the first compartment exactly when s in > 1/r.Furthermore, this equilibrium is unique.Then, any equilibrium for the overall system (1) has to be (s ⋆ 2 , s in − s ⋆ 2 ) for the second compartment, with s ⋆ 2 solution of the equation with s ⋆ 2 < 1/r.One can easily verify that there exists a unique s ⋆ 2 solution of (2) on (0, 1/r).Graphically, s ⋆ 2 is the abscissa of the intersection of the graphs (see Figure 3) of the polynomial function and the affine function Remark that s in > 1/r implies the inequality φ 2 with the value obtained in the configuration of one compartment:

Parallel interconnection of two compartments
The dynamical equations of the model with two compartments in parallel and diffusion (see Figure 1), assuming r to be different to 0 and 1, are the following where the output concentration s out is given by The wash-out in both compartments corresponds to the trivial equilibrium (s in , 0, s in , 0), that leads to the trivial steady-state s ⋆ out = s in .For convenience, we posit and assume, without any loss of generality that one has α 2 ≥ α 1 (if it is not the case one can just exchange indexes 1 and 2).
When d = 0 (no diffusion), the equilibrium of the system can be determined independently in the two compartments as simple chemostats.In this case, there is an unique globally stable equilibrium (s in the non-negative orthant, where s * i = min(α i , s in ) (i = 1, 2).When d > 0, we define the functions Proposition 2. When s in > 1 and d > 0, there exists a unique equilibrium (s is the unique solution of the system on the domain (0, s in ) × (0, s in ), with when α 2 > α 1 .
Proof.At equilibrium, one has which amounts to write, from equations ( 3) and deduces the property Consequently, an equilibrium in the positive orthant has to fulfill 3) at equilibrium, one obtains the equations which amounts to write that (s ⋆ 1 , s ⋆ 2 ) is solution of the system (4) (see Figure 4) or equivalently s ⋆ 1 is a zero of the function g(•).When (4).When α 2 > α 1 , one has necessarily α 1 < 1 and the condition s in > 1 implies g(α 1 ) < 0. We distinguish now two cases: In both cases, one deduces by the Mean Value Theorem the existence of s Rolle and Mean Value Theorems allow to conclude the existence of s ⋆ 1 ∈ (α 1 , s in ) such that g(s ⋆ 1 ) = 0.
In any case, we obtain the existence of (s ⋆ 1 , s ⋆ 2 ) solution of (4) with Thus, the inequalities (5) are fulfilled.
Finally, notice that functions φ 1 (•), φ 2 (•) are both strictly concave, and steady states (s ⋆ 1 , s ⋆ 2 ) are intersections of G 1 , the graph of the function φ 1 (•), and G 2 the symmetric of the graph of φ 2 (•) with respect to the first diagonal.Consequently, if (s ⋆ 1 , s ⋆ 2 ) is a steady state different from (s in , s in ), G 1 and G 2 are respectively above and below the line segment (s ⋆ 1 , s ⋆ 2 ) − (s in , s in ).We conclude that there exists at most one non-trivial equilibrium.
Proof.When α 1 = α 2 , one has s ⋆ 1 = s ⋆ 2 = 1 and one can easily check When α 2 > α 1 , one has g(α 1 ) < 0 and we recall from the proof of former Proposition that s ⋆ 1 is the unique zero of g(•) on (α 1 , min(α 2 , s in )).We conclude that g is non decreasing at s ⋆ 1 .Notice that φ 1 and φ 2 are concave functions and that ) cannot be equal to zero, and consequently one has g ′ (s ⋆ 1 ) > 0.
The global stability of the non-trivial equilibrium is proved in the Appendix.
Proposition 2 defines properly the map d → s ⋆ out = αs ⋆ 1 + (1 − α)s ⋆ 2 for the unique non-trivial steady-state, that we aim at studying as a function of d.Accordingly to Proposition 2, s ⋆ out is equal to one for any value of the parameter d in the non-generic case α 2 = α 1 .We shall focus on the case α 2 = α 1 (and without loss of generality we shall consider α 2 > α 1 ).We start by the two extreme situations: no diffusion and infinite diffusion.
Lemma 1.For the non trivial equilibrium, one has Proof.Under the assumptions s in > 1 and α 2 ≥ α 1 , we distinguish two cases when Then, one can write (recall that assuming α 2 ≥ α 1 imposes to have α < 1, and s 0 in is well defined).Notice that the number s 0 in is necessarily larger than one because α 1 ≤ 1, and one has also Consequently one concludes that s ⋆ out ≥ 1 exactly when s in ≥ s 0 in .Finally, remark that one has Lemma 2. For s in > 1, the non trivial equilibrium fulfill Proof.For any d > 0, Proposition guarantees the existence of a unique non trivial equilibrium (s ⋆ 1 , s ⋆ 2 ) ∈ (0, s in ) × (0, s in ) that is solution of (6).When d is arbitrary large, one obtains from (6) From equations ( 6), one deduces also the following equality valid for any d that can rewritten, taking into account the equality rα Consequently, one has If α 2 < s in , the property s ⋆ 1 < α 2 valid for any d > 0 implies that s ⋆ 1 cannot converges to s in .If α 2 ≥ s in and lim s ⋆ 1 = lim s ⋆ 2 = s in , there exists d such that rs ⋆ 1 + (1 − r)s ⋆ 2 > (s in + 1)/2.Then, one has that contradicts Corollary 2. Finally, one has lim s ⋆ 1 = lim s ⋆ 2 = 1 and consequently lim s ⋆ out = 1.
We present now our main result concerning properties of the map d → s ⋆ out (d) defined at the non-trivial steady-state.
-When s in ≥ 2, the map d → s ⋆ out (d) the non trivial equilibrium) is decreasing and s ⋆ out (d) > 1 for any d ≥ 0.
-When s in < 2, the map d → s ⋆ out (d) (for the non trivial equilibrium) admits a minimum in d ⋆ < +∞, that is strictly less than one.Furthermore, one has Proof.Let differentiate with respect to d the equations ( 6) at steady state: that can rewritten as follows 1 1 Remark that one has From the Corollary 2, one has det(Γ) < 0 and one deduces that the derivatives Notice from inequalities (5) that we obtain B > 0 and deduce ∂ d s ⋆ 1 > 0 for any d.With Lemma 2 we conclude that s ⋆ 1 (d) < 1 for any d.From equations (7), we can write When s in ≥ 2, one has A > 0 and then ∂ d s ⋆ 2 < 0. With Lemma 2 we conclude that s ⋆ 2 (d) > 1 for any d.Then, one obtain the inequality σ < (s in − 2)(α 1 − α 2 ) ≤ 0 which proves with Lemma 2 that s ⋆ out is a decreasing function of d that converges to one.
When s in < 2, we write As s ⋆ 1 and s ⋆ 2 tend to one when d takes arbitrary large values, we conclude that there exists d < +∞ such that σ > 0 for any d > d and consequently s ⋆ out is smaller than one and increasing for d > d.We conclude that the map d → s ⋆ out (d) admits a minimum, say at d ⋆ < +∞, that is strictly less than one.
, for which we conclude Remark that this case is feasible because of the inequality 2α 1 α 2 < min(2, α 2 )(α 1 + α 2 ).We conclude that for s in larger than this last value, d ⋆ is necessarily strictly positive.
Remark.The particular case α = 0 corresponds to a configuration of a perfectly mixed tank of volume (1 − r)V connected to a dead-zone of volume rV .This is a way to approximate a non well-mixed tank or segregated bioreactors of total volume V , estimating the fraction of the volume occupied by the highly agitated area.

Numerical computation and discussion
Propositions 1 and 3 reveal the existence of a threshold on the value of the input concentration s in (equal to 2 for our choice of the parameters units) that reverses the performances of the serial and parallel configurations in terms of s ⋆ out , compared to the single tank case (for which s ⋆ out = 1): -for s in > 2, there exist serial configurations such that s ⋆ out < 1 for r large enough (i.e. the first tank has to be large enough), but any parallel configuration produces s ⋆ out > 1, -for s in < 2, there exists parallel configurations such that s ⋆ out < 1, while any serial configuration has s ⋆ out > 1.There exists another threshold s 0 in ∈ (1, 2) such that configurations with s ⋆ out < 1 require to have d large enough when s in > s 0 in (cf Lemma 1).Furthermore, the best performance of the parallel configuration is obtained -for arbitrary large values of d when s in > 2 , -for a finite positive d ⋆ when s in ∈ (s in , 2) (where the expression of s in is given in Proposition 3).For the parallel interconnection, we depict on Figure 6 the two kind of configurations that occur, depending on whether the number s in is larger than one or not.The analytic analysis of Section 3 has been conducted under the assumption of the linearity of the function µ(•).It is often in microbiology that the growth rate µ(•) presents a concavity, as described by the usual Monod (or Michaelis-Menten) function.We have computed numerically the same curves s ⋆ out (•) than Figures 5 and 6, considering the Monod function µ(S) = 6S 5 + S instead of the linear function (see Figure 7).This function has been chosen to fulfill s ⋆ out = 1 for the single tank configuration, guaranteeing the same steady state than the linear growth for this configuration.

The values of the parameters are given on the table below
On Figures 8 and 9, we observe that the concavity of the growth function does not change qualitatively the theoretical results and the existence of threshold for s in that favourites one of the configuration.
We notice on all the figures that the yield is better for the Monod function in the parallel configuration and worst for the serial one.This implies that the threshold on s in , that was determined to be equal to 2 for the linear case is higher when the growth function is concave.
Remarks.The serial configuration for the limiting value r = 1 is equivalent to a single tank.This explains why all the curves on Figures 5 and 8   For the parallel configuration with α = 0.1 and r = 0.9 one has α 2 = 9.This implies that for the limiting value d = 0 the only equilibrium in the second tank is the wash-out when s in < 9.This is not the case for the first tank but the flow rate αQ being small, the output s ⋆ out remains closed to s in in any case, as one can see on Figures 6 and 9 for small values of the parameter d.

Conclusion
Given a flow rate and the total volume of a chemostat system, this study shows the existence of a threshold on the value of the input concentration s in such that above and below this threshold, serial and the parallel configurations are respectively the best ones with respect to the criterion of minimizing the output concentration s ⋆ out at steady state.For the parallel scheme, the best performances are obtained for a precise value of the diffusion parameter that is proved to be positive when s in is not too small.This study concerns also dead-zone configurations, as particular cases of the parallel configurations.
Whatever are the data of the problem, there always exists a configuration that is better than a single perfectly mixed tank.We have shown that the non-trivial steady states are unique and globally exponentially stable under the assumption of a linear increasing growth rate.
Finally, this study reveals the role of the structure of the space on the performances of simple ecosystems or bioprocesses.The possibly non-monotonic influence of the diffusion parameter on the output steady state is not intuitive, and leave further investigations open for understanding or taking benefit of this property for natural ecosystems (such as saturated soils or wetlands) as well as for bioprocesses (such as waste-water treatments).This result can be also of interest for reverse engineering when deciding which among serial or parallel configurations is better fit for the modeling of chemostat-like ecosystems, providing that one has an estimation of the hydric capacity of the system.
Acknowledgments.This work has been achieved within the VITELBIO (VIRtual TELluric BIOreactors) program, sponsored by INRA and INRIA.The authors are grateful to this support.The work is also part of the PhD thesis of the first author.
6 Appendix: global exponential stability of the non-trivial equilibrium First, one can easily check that the domain D = IR 4 + is invariant by the dynamics ( 1) and ( 3).We consider the 2-dimensional vector z of variables z i = s in − x i − s i (i = 1, 2) whose dynamics are respectively for the serial and parallel configurations Notice that matrices A s and A p are Hurwitz : So z converges exponentially toward 0 for both systems, which implies that dynamics (1) and ( 3) are dissipative, in the sense that any solution of (1) or (3) in D converge exponentially to the compact set K = {(s 1 , x 1 , s 2 , x 2 ) ∈ D s.t.x 1 + s 1 = s in and x 2 + s 2 = s in }.
We recall a result from [19,Theorem 1.8] that shall be useful in the following.

The serial configuration
Proposition 5.Under the condition s in > 1/r, any trajectory of (1) with initial condition in D such that (s 1 (0), x 1 (0)) = (s in , 0) converges exponentially to the unique non-trivial steady-state (s 2 ) given by Proposition 1.
Proof.Dynamics (1) has a cascade structure.It is straightforward to check that the solutions of the (s 1 , x 1 ) sub-system converges asymptotically towards the non-trivial equilibrium (1/r, s in − 1/r) from any initial condition away from the wash-out equilibrium (s in , 0).From the convergence of z 2 toward 0, we deduce that the s 2 variable has to converge to the bounded interval [0, s in ] and that its dynamics can be written as a scalar non autonomous differential equation: This last dynamics has the property to be asymptotically autonomous with the limiting differential equation: Statement of Proposition 1 implies that this last scalar dynamics has a unique equilibrium s ⋆ 2 that belongs to [0, s in ].Furthermore, one has f (0) > 0 and f (s in ) < 0. Consequently any solution of ( 9) in [0, s in ] converges asymptotically to s ⋆ 2 .Then applying Theorem 4, we conclude that any bounded solution of (8) converges to s ⋆ 2 .Finally any solutions of the (s 2 , x 2 ) sub-system converges asymptotically to (s ⋆ 2 , s in − s ⋆ 2 ).
The Jacobian matrix of dynamics (1) at the non-trivial equilibrium (s 2 ) that are both negative numbers, accordingly to Proposition 1.The exponential stability of the non-trivial equilibrium is thus proved.

The parallel configuration
Proposition 6.When s in > 1 and d > 0, any trajectory of (3) with initial condition in D such that x 1 (0) > 0 and x 2 (0) > 0 converges exponentially to the unique non-trivial steady-state (s given by Proposition 2.
Proof.Considering the time vector z(•), the (s 1 , s 2 ) sub-system of dynamics (3) can be written as solution of a non-autonomous planar dynamics ṡ1 = s 1 (z We know that z converges to 0 and consequently the vector S of variables s From Poincaré-Bendixon theorem and Dulac criterion, we conclude that bounded trajectories of (12) cannot have limit cycle or closed path and necessarily converge to an equilibrium point.Consequently, any trajectory of (11) in S either converges to the rest point S ⋆ = (s ⋆ 1 , s ⋆ 2 ) or approaches the boundary B. Notice that one has s i = s in , s j < s in ⇒ ṡi < 0 (i = j) So the only possibility for approaching B is to converge to the other rest point S 0 = (s in , s in ).This shows that the only non-empty, closed, connected, invariant and chain recurrent subsets of S are the singletons {S ⋆ } and {S 0 }.
Applying Theorem 4 we conclude that any trajectory of (10), issued from initial condition of dynamics (3) in D, converges asymptotically to S ⋆ or S 0 .Consider now any initial condition with x 1 (0) > 0 and x 2 (0) > 0. We show that the solution (s 1 (•), s 2 (•)) of (3) cannot converge to S 0 .If it is the case, there exists T < +∞ such that one has s 1 (t) > α 1 and rs 1 (t) + (1 − r)s 2 (t) > 1 for any t ≥ T under the assumption s in > 1.Let us consider the function V (x 1 , x 2 ) = min(rx 1 + (1 − r)x 2 , x 1 ) (see Figure 10) and v(t) = V (x 1 (t), x 2 (t)) that is positive and tends to 0 when t tends to +∞ We conclude that the function t → v(t) is non-decreasing for t ≥ T and consequently cannot converge to zero, thus a contradiction.
The Jacobian matrix of dynamics (3) at the non-trivial equilibrium (s ⋆ 1 , x ⋆ 1 , s ⋆ 2 , x ⋆ 2 ) is of the following form in (z 1 , z 2 , s 1 , s 2 ) coordinates Recall that A p is Hurwitz.One has

Figure 1 :
Figure 1: The set of configurations under investigation.

For
the serial configuration, the graph of the function s ⋆ out is plotted as function of r ∈ [1/s in , 1] on Figure5for different values of the input concentration s in .

Figure 5 :
Figure 5: Comparison of s ⋆ out for the serial configuration.

Figure 6 :
Figure 6: Comparison of s ⋆ out for the parallel configuration (s in > 1 on the left and s in < 1 on the right) .

Figure 9 :
Figure 9: Comparison of s ⋆ out for the parallel configuration (s in > 1 on the left and s in < 1 on the right) .

Figure 10 :
Figure 10: Iso-value of the V function.
coincide for this value of the parameter r.