About the chemostat model with a lateral diffusive compartment

# About the chemostat model with a lateral diffusive compartment

María Crespo and Alain Rapaport
MISTEA, Univ. Montpellier, INRA, Montpellier SupAgro
2 pl. Viala 34060 Montpellier, France
E-mails: maria.crespo@umontpellier.fr,alain.rapaport@inra.fr
August 4, 2019
###### Abstract

We consider the classical chemostat model with an additional compartment connected by pure diffusion, and analyze its asymptotic properties. We investigate conditions under which this spatial structure is beneficial for species survival and yield conversion, compared to single chemostat. Moreover we look for the best structure (volume repartition and diffusion rate) which minimizes the volume required to attain a desired yield conversion. The analysis reveals that configurations with a single tank connected by diffusion to the input stream can be the most efficient.

Key-words. chemostat model, compartments, diffusion, stability, yield conversion, optimal design.

## 1 Introduction

The model of the chemostat has been developed as a mathematical representation of the chemostat apparatus invented in the fifities simultaneously by Monod [37] and Novick & Szilard [40], for studying the culture of micro-organisms, and is still today of primer importance (see for instance [24, 20, 60]). Its mathematical analysis has led to the so-called “theory of the chemostat” [25, 53, 19]. This model is widely used for industrial applications with continuously fed “bioreactors” for fermentations [43, 55] or waste-water treatments [14, 20, 60], but also in ecology for studying populations of micro-organisms (or plankton) in lakes, wetlands, rivers or aquaculture ecosystems [27, 1, 28, 29, 46, 2], The word “chemostat” is often used to describe continuous cultures of micro-organisms, even though it can be quite far from the the original experimental setup. The classical model of the chemostat assumes a perfectly mixed media which is generally verified for small volumes. For industrial bioreactors or natural lakes with large volumes, the validity of this assumption becomes questionable. This is why several extensions of this model with spatial considerations have been proposed and studied in the literature.

The classical approaches for modeling non ideally mixed chemostats or bioreactors rely on a continuous representation of the spatial dimension (with systems of p.d.e. as in [30, 8]) or on a finite number of interconnected compartments with different flow conditions (with systems of o.d.e. as in the “general gradostat” [51, 54, 16]). Although there are many works in the literature on p.d.e. models for unmixed bioreactors (where nutrient diffuses on the media from the boundary of the domain, see for instance [53, 10]), there are comparatively few works for bioreactors with an advection term that represents an incoming nutrient which is pushed inward (as in the chemostat apparatus). Moreover, most of the mathematical analysis available in the literature consider a spatial heterogeneity only in the axial dimension of the bioreactors (leading to 1d p.d.e. or compartments connected in series) as in tubular or “plug-flow” bioreactors [12, 7, 9, 61, 6] and (simple) gradostats [32, 33, 56, 52]. Surprisingly, configurations of tanks in parallel rather than in series have been much less investigated, apart simple considerations in chemical reaction engineering [31, 11].

In many cases, the axial direction appears to be the one that generates the larger heterogeneity between the input and output (when the main current lines are along this axis), especially for height and relatively thin tanks under significant flow rate. The modeling with compartments has become quite popular in the optimal design of bio-processes [59], as it allows to determine easily the optimal sizes of a given number of tanks in series for minimizing the residence time [34, 23, 4, 5, 22, 21, 18, 39, 42]. Several studies have shown the huge benefit that can be obtained from one to two tanks in series (and even more but marginally less with more than two). Such considerations are similar to patches models or islands models, commonly used in theoretical ecology [35, 17] (or lattice differential equations [48]). For instance, a recent investigation studies the influence of these structures on a consumer/resource model [15]. However, ecological consumer/resource models are similar to chemostat models apart the source terms that are modeled as constant intakes of nutrient, instead of dilution rates (or Robin boundary conditions) that are rather met in liquid media.

Recent mathematical studies have revealed that considering heterogeneity in directions transverse to the axial one (with 2d or 3d p.d.e models or compartments models with non serial interconnections) could have a significant impacts on the performances of the bioreactor and the input-output behavior [3]. From an operational view point, it is often reported that “dead zones” are observed in bioreactors and that the effective volumes of the tanks have to be corrected in the models to provide accurate predictions [31, 26, 13, 45, 44, 58, 47]. 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 [38]. In a similar way, stagnant zones are well-known to occur in porous media such as soils, at various extents depending on soil structure. The effect of these dead zones on reactive and conservative mass transport, and thus in turn on the biogeochemical cycles of elements, can also be significant [57, 49]. However there have been relatively few analysis of models with explicit “dead-zones”. The wording “dead-zone” might be slightly miss-leading as it can make believe that a part of the tank where there is no advection stream from the input flow has no biological activity. But this does not necessarily mean that these “dead-zones” are entirely disconnected from other parts of the reactor. It is likely to be influenced by diffusion rather than convection. This is why we prefer to qualify these zones as “lateral-diffusive compartments”.

The aim of the present work is to analyze the chemostat model with two compartments (or two tanks), one of them being connected by “lateral-diffusion”, and to investigate conditions under which having this compartment could be beneficial for the yield conversion compared to the chemostat model with a single compartment of the same total volume. Although this structure falls into a particular case of the general gradostat, the existing results in the literature (see, e.g., [56, 53]) are too general to give accurate conditions for the existence of a positive equilibrium and its stability, and do not compare the performances of each configuration. The structure considered here can be also seen as a limiting case of the pattern “chemostats in parallel with diffusion connection” studied in [16], with only one vessel receiving an input flow rate. Nevertheless, this later reference imposes some restrictions such as linear reaction between species and removal rate large enough to avoid washout with a single tank. In the present work, we conduct a deeper model analysis and investigate the impact of lateral diffusion from two view points:

1. From an ecological perspective, we study the effect of the diffusion on a given volume repartition in terms of resource conversion. In particular, we aim at characterizing situations for which having a structure with lateral diffusion is better than having a single perfectly mixed volume.

2. From an engineering perspective, we look for the best volume repartition which, for a given diffusion rate, minimizes the total volume required to attain a desired resource conversion. This allows us to revisit the optimal design problem with such configurations, that was previously tackled but considering tanks connected in series (see, e.g., [4, 22, 18]). Additionally, we aim at determining the diffusion rate parameter that gives the best volume reduction.

The article is organized as follows: in Section 2 we introduce the model describing the dynamics in the chemostat composed of two compartments, one of them being connected by diffusion, and give results for the nonnegativity and boundedness of its solutions. In Section 3 we determine the steady states and analyze its asymptotic stability. Section 4 investigates the resource conversion rate and characterizes when this structure is better than the single chemostat (i.e. a single perflectly mixed volume). Section 5 is dedicated to optimal design questions, firstly when the diffusion rate is fixed and then when it can be tuned. Finally, Section 6 discusses and interprets the results.

## 2 Modeling and preliminaries results

We consider configurations of one tank of volume interconnected by Fickian diffusion with a tank of volume , as depicted on Figure 1.

We denote by , the concentrations of substrates and biomass in tank and write the equations of the chemostat model for such interconnections:

 (1)

where we have assumed, without any loss of generality, that the yield conversion factor of substrate into biomass is equal to . The parameters and denote the flow rate and substrate concentration of the input stream, while the parameter is the diffusion coefficient between the two tanks (that we assume to be identical for the substrate and the micro-organisms). The specific growth rate function of the micro-organisms is denoted and fulfills the classical following assumption.

###### Hypothesis 1.

The growth function is increasing concave function with .

A typical such instance of function is given by the Monod law (see, e.g., [19, 53]):

 μ(s)=μmaxsK+s,

where is the maximum specific growth rate and is the half-saturation constant.

###### Lemma 1.

The non-negative orthant is invariant by dynamics (1) and any solution in is bounded.

###### Proof.

Define for each tank and consider the dynamics (1) in coordinates:

 ⎧⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎩˙z1=−QV1z1−dV1(z1−z2),˙s1=−μ(s1)(sin−s1−z1)+QV1(sin−s1)+dV1(s2−s1),˙z2=−QV2(z2−z1),˙s2=−μ(s2)(sin−s2−z2)+dV2(s1−s2). (2)

This system has a cascade structure with a first independent sub-system linear in

 ˙z=⎡⎢ ⎢ ⎢⎣−Q+dV1dV1dV2−dV2⎤⎥ ⎥ ⎥⎦Az, (3)

where one has

 tr(A)=−Q+dV1−dV2<0anddet% (A)=QdV1V2>0.

Therefore the matrix is Hurwitz and any solution of (3) converges exponentially to . Then, the solution can be written as the solution of the non autonomous dynamics

 ˙s=F(t,s)=⎡⎢ ⎢ ⎢⎣(QV1−μ(s1))(sin−s1)+dV1(s2−s1)+μ(s1)z1(t)−μ(s2)(sin−s2)+dV2(s1−s2)+μ(s2)z2(t)⎤⎥ ⎥ ⎥⎦. (4)

Notice that, for any , one has

 ∂F1(t,s)∂s2=dV1>0and∂F2(t,s)∂s1=dV2>0,

and so the dynamics (4) is cooperative (see, e.g., [50]).
Define , for which it follows that for any . Proposition 2.1 in [50] allows to state that any solution of (4) with () satisfies ( for any , where is solution of the dynamics

 ˙ˇs=ˇF(t,ˇs)=[ˇF1(t,ˇs)F2(t,ˇs)],ˇs(0)=0

As one has for any , the solution is identically null and one obtains that () stays non-negative for any positive .

Similarly, can be written as a solution of a non-autonomous cooperative dynamics

 ˙x=H(t,x)=⎡⎢ ⎢ ⎢ ⎢⎣(μ(s1(t))−Q+dV1)x1+dV1x2dV2x1+(μ(s2(t))−dV2)x2⎤⎥ ⎥ ⎥ ⎥⎦

with , which allows to conclude that () stays non-negative for any positive .

Finally, the convergence of to provides the boundedness of the solutions , . ∎

Let us discuss the modeling of the two limiting cases that are not covered by system (1):

• . Physically, this corresponds to a single tank (of volume ) connected by diffusion to the input pipe with flow rate (see Figure 2).

There is no biological activity in the pipe but simply a dilution given by the mass balance at the connection point:

 {Q(sin−sout)=d(sout−s2)−Qxout=d(xout−x2)⇒sout=Qsin+ds2Q+d,xout=dx2Q+d (5)

Then the dynamics in the tank with instead of is given by the equations

This is equivalent to have a single tank of volume with input flow rate but with an output given by . Equivalently, one can consider the system (1) as a slow-fast dynamics when the parameter is small. Then the fast dynamics is

 {ϵ˙s1=−ϵμ(s1)x1+Q(sin−s1)+d(s2−s1)ϵ˙x1=ϵμ(s1)x1−Qx1+d(x2−x1)

and the slow manifold is given exactly by the system of equations (5) with .

• . The dynamics of is the one of the single chemostat model with volume

 ⎧⎪ ⎪ ⎪⎨⎪ ⎪ ⎪⎩˙s1=−μ(s1)x1+QV1(sin−s1)˙x1=μ(s1)x1−QV1x1 (6)

This is also equivalent to having no diffusion () between the tanks.

## 3 Study of the equilibria

For the analysis of the steady sates, it is convenient to introduce the function

 β(s)=μ(s)(sin−s) (7)

that verifies the following property.

###### Lemma 2.

Under Hypothesis 1, the function is strictly concave on . Thus, one can define the unique value

 ^s=argmaxs∈(0,sin)β(s). (8)
###### Proof.

One has and , which is negative for any . ∎

One can check that the washout state is always a steady-state of the dynamics (1). In the next Propositions, we characterize the existence of other equilibrium and their global stability.

###### Proposition 1.

The washout equilibrium is the unique steady state of (1) exactly when satisfies the condition

 μ(sin)≤QV1 and P(μ(sin))≥0 (9)

where is defined as

 P(X)=V1V2X2−(dV1+(Q+d)V2)X+dQ

When the condition (9) is not fulfilled, there exists an unique positive steady state of (1) distinct from .

###### Proof.

From the two last equations of (1), one has at steady-state, and from the two first ones . The values , at steady state are solutions of the system of two equations

 0 = (QV1−μ(s1))(sin−s1)+dV1(s2−s1) (10) 0 = −μ(s2)(sin−s2)+dV2(s1−s2) (11)

and , at steady state are uniquely defined from each solution of (10)-(11).

Clearly is a solution of (10)-(11). We look for (positive) solutions different to .

Posit

 λ1(sin):=max{s1∈[0,sin]|μ(s1)≤QV1}

From equations (10)-(11), a solution different to has to verify and then from equation (10), one has also . Define then the functions:

 ϕ1(s1) := s1−Q−V1μ(s1)d(sin−s1)=s1−Qd(sin−s1)+V1dβ(s1), ϕ2(s2) := s2+V2μ(s2)d(sin−s2)=s2+V2dβ(s2).

It is straightforward to see that any solution of (10)-(11) fulfills and . One has

 ϕ′1(s1)=1+V1dμ′(s1)(sin−s1)+Q−V1μ(s1)d.

Therefore is increasing on , with and . Thus, is invertible on with

 ϕ−11(0)∈(0,λ1(sin)).

From Lemma 2, it follows that and are strictly concave functions on . Consider then the function

 γ(s2)=ϕ2(s2)−ϕ−11(s2)s2∈[0,sin],

which is also strictly concave on . Then, a solution can be written as a solution of

 γ(s2)=0,s1=ϕ2(s2) with s2∈[0,λ1(sin)].

Notice that one has , and as is strictly concave, it cannot have more than two zeros. Therefore there is at most one solution different to . Furthermore, one has . Now, distinguish two different cases:

When (or equivalently ), one has

 γ(λ1(sin))=QV2dV1(sin−λ1(sin))>0.

By using the Mean Value Theorem, one concludes that there exists such that .

When (that is when ), the function takes positive value on the interval if and only if ( being strictly concave on ), or equivalently when the condition

 ϕ′2(sin)<1ϕ′1(sin)

is fulfilled. Notice that one has because . So the condition can be also written as

 ϕ′1(sin)ϕ′2(sin)<1.

From the expressions of and , one can write this condition as

 (d+Q−V1μ(sin))(d−V2μ(sin))d2<1,

and easily check that this amounts to require to satisfy .

We conclude that there exists a positive steady state if and only if or and that this steady state (when it exists) is unique. ∎

###### Proposition 2.

When the washout equilibrium is the unique steady state, it is globally asymptotically stable on .

When the positive steady state exists, for any initial condition except on a set of null measure, the solution of (1) converges asymptotically to , which is moreover locally exponentially stable.

###### Proof.

As shown in the proof of Lemma 1, the dynamics (1) has a cascade structure in coordinates, where converges exponentially to . Then converges to the set and is solution of the non-autonomous system (4) in the plane, which is asymptotically autonomous with limiting dynamics

 (12)

We study now the asymptotic behavior of dynamics (12) in the domain . Denote by the subset of its boundary. Accordingly to Proposition 1, this dynamics has (which is the projection of on the -plane) as an equilibrium, and at most another equilibrium (which is the projection of on the -plane) that has to belong to . On this last subset, we consider the variables , whose dynamics are given by

 ˙σ=~Fa(σ)=⎡⎢ ⎢ ⎢⎣−QV1+μ(sin−eσ1)−dV1(1−eσ2−σ1)μ(sin−eσ2)−dV2(1−eσ1−σ2)⎤⎥ ⎥ ⎥⎦ (13)

One obtains

 div(~Fa)=−μ′(sin−eσ1)eσ1−μ′(sin−eσ2)eσ2−dV1eσ2−σ1−dV2eσ1−σ2<0

From Dulac criteria and Poincaré-Bendixon theorem (see for instance [41]), we conclude that bounded trajectories of (13) cannot have limit cycle or closed path and necessarily converge to an equilibrium point. Consequently, any trajectory of (12) in either converges to the equilibrium (if it exists) or approaches the boundary . But

 if there exists t such that si(t)=sin and sj(t)

and so the only possibility for approaching is to converge to . This shows that the only non-empty closed connected invariant chain recurrent subsets of are the isolated points and .

We use the theory of asymptotically autonomous systems (see Theorem 1.8 in [36]) to deduce that any trajectory of system (1) in converges asymptotically to or .

Due to the cascade structure of the dynamics (1) that is made explicit in the proof of Lemma 1, the Jacobian matrix in the coordinates is

where the matrix defined in (3) is Hurwitz. Thus, the eigenvalues of are the ones of the matrix plus the ones of . One has

 tr(Ja(s))=−d(ϕ′1(s1)V1+ϕ′2(s2)V2)anddet(Ja(s))=d2V1V2(ϕ′1(s1)ϕ′2(s2)−1).

Accordingly to Proposition 1, the equilibrium exists when or .

When , one has or equivalently . Then is a saddle point (with a stable manifold of dimension one) and we deduce that almost any trajectory of (1) converges to .

When , notice that the equilibrium is not necessarily hyperbolic (as one can have which implies then ) and we cannot conclude its stability properties directly. We setup a proof by contradiction, inspired by [16]. Consider a solution of (12) with initial condition different to that converges asymptotically to , and define the function

 v(t)=min(~x(t),x1(t)) with ~x=V1x1+V2x2V1+V2,

which verifies for any and has to converge to when tends to . The condition implies the existence of positive numbers and such that for any . If with , one has and

 ˙x1≥(μ(s1(t))−QV1)x1≥ηx1≥0

If with , one has and

 ˙~x=V1V1+V2(μ(s1(t))−QV1)x1+V2V1+V2μ(s2(t))x2≥V1V1+V2ηx1≥0

The positive function is thus non decreasing for , in contradiction with its convergence to . We conclude that any solution of (1) with initial condition different to converges to .

Finally, when the equilibrium exists, the analysis conducted in the proof of Proposition 1 allows us to deduce the inequalities and , which in turn imply and so . Then, one has and i.e. is Hurwitz, which proves that the attractive equilibrium is also locally exponentially stable. ∎

## 4 Influence of lateral diffusion on the performances

Here, we investigate conditions under which having a second compartment (as proposed in Section 2) is beneficial for the yield conversion compared to the chemostat model with a single compartment of the same total volume. To this aim, we fix the hydric volumes and , the input flow , and analyze the output map at steady state, as function of the diffusion parameter . The benefits of the structured chemostat in terms of resource conversion are discussed in Section 6.

Proposition 1 defines properly the map for the unique non-trivial steady-state of system (1), that we study here as a function of . We start by deducing the range of existence of this steady-state.

###### Proposition 3.

Let and . It follows that:

1. If , then the non-trivial equilibrium exists when .

2. If , then the non-trivial equilibrium exists when .

3. If , then the non-trivial equilibrium exists when .

###### Proof.

When (that is, when the volume is detached) the classical equilibria analysis of the single chemostat model with volume (see, e.g., [53]) assures that the positive equilibrium exists when , which corresponds to the case (iii) on the proposition statement.
When , we prove cases (i)-(iii) by taking into account that they correspond to three different scenarios where condition (9) is not fulfilled. For ease of reasoning, we rewrite as

 P(μ(sin))=V2μ(sin)(V1μ(sin)−Q)A+d(Q−(V1+V2)μ(sin))B. (14)

(i). In this case, the non-trivial equilibrium exists when . It is straightforward to see that , and so exists when .
(ii). In this case, the non-trivial equilibrium exists when . It is straightforward to see that , and so exists for all .
(iii). In this case, the non-trivial equilibrium exists for all values of , and so exists for all . ∎

We now study the two extreme situations: no diffusion and infinite diffusion.

###### Lemma 3.

It follows that

1. When , the non trivial equilibrium of system (1) fulfills

 s⋆1(0)=s⋆,01,

where is the non-trivial steady state of the single chemostat model with volume .
In other case .

2. When , the non trivial equilibrium of system (1) fulfills

 limd→+∞s⋆1(d)=s⋆,∞1,

where is the non-trivial steady state of the single chemostat model with volume .

###### Proof.

(i). This result is a direct consequence of the classical equilibria analysis of the single chemostat model with volume (see, e.g., [53]), which assures that exists when .

(ii). For any , Proposition 1 guarantees the existence of a unique non trivial equilibrium that is solution of

 ⎧⎪ ⎪⎨⎪ ⎪⎩d(s⋆2−s⋆1)=(V1μ(s⋆1)−Q)(sin−s⋆1),d(s⋆1−s⋆2)=V2μ(s⋆2)(sin−s⋆2). (15)

When is arbitrary large, one obtains

 limd→+∞s⋆1−s⋆2=0.

From equations (15), one can also deduce the following equality (valid for any )

 (V1μ(s⋆1)−Q)(sin−s⋆1)=−V2μ(s⋆2)(sin−s⋆2). (16)

Consequently, one has

 limd→+∞s⋆1(d)=limd→+∞s⋆2(d)=sin or limd→+∞s⋆1(d)=limd→+∞s⋆2(d)=s⋆,∞1as V=V1+V2,

where the classical equilibria analysis of the single chemostat model with volume assures that exists when . But Proposition 3 shows that, under the assumptions of the lemma, cannot converge to . ∎

We now present our main result concerning properties of the map , defined at the non-trivial steady state.

###### Proposition 4.

Let be defined in (8) and . It follows that:

1. If , then the map admits a minimum in that is strictly less than .

2. If and , then the map admits a minimum in , that is strictly less than .

3. If and , then the map is decreasing and for any .

###### Proof.

If one differentiates system (15) with respect to , it follows that

 (s⋆2−s⋆1)+d(∂ds⋆2−∂ds⋆1)=∂ds⋆1(Q+V1μ′(s⋆1)(sin−s⋆1)−V1μ(s⋆1))A,
 (s⋆1−s⋆2)+d(∂ds⋆1−∂ds⋆2)=∂ds⋆2(V2μ′(s⋆2)(sin−s⋆2)−V2μ(s⋆2))B,

which can be rewritten as

Remark that

 A+d=dϕ′1(s⋆1),B+d=dϕ′2(s⋆2),det(Γ)=d2(1−ϕ′1(s⋆1)ϕ′2(s⋆2)).

As seen in the proof of Proposition 2, one has that and so the derivatives and can be defined as

Firstly, we prove that by showing that .
From Proposition 1, one has that the positive steady-state fulfills

 0

Since is concave (equivalently, is decreasing) on , one has that . Thus, we prove that by showing :

• If , then and .

• If , then and .

Therefore one has i.e. is an increasing map.

Now, notice that and its sign depends on the relative position of with respect to parameter . The cases considered on the proposition statement are treated separately.

(i) Since is increasing, , and , by using the Mean Value Theorem it follows that there exists a unique value (denoted by ) such that , with for and for . Consequently,