The mathematical analysis of a syntrophic relationship between two microbial species in a chemostat

. A mathematical model involving a syntrophic relationship between two populations of bacteria in a continuous culture is proposed. A detailed qualitative analysis is carried out as well as the analysis of the local and global stability of the equilibria. We demonstrate, under general assump- tions of monotonicity which are relevant from an applied point of view, the asymptotic stability of the positive equilibrium point which corresponds to the coexistence of the two bacteria. A syntrophic relationship in the anaerobic digestion process is proposed as a real candidate for this model.


THE MATHEMATICAL ANALYSIS OF A SYNTROPHIC RELATIONSHIP BETWEEN TWO MICROBIAL SPECIES IN A CHEMOSTAT
Abstract. A mathematical model involving a syntrophic relationship between two populations of bacteria in a continuous culture is proposed. A detailed qualitative analysis is carried out as well as the analysis of the local and global stability of the equilibria. We demonstrate, under general assumptions of monotonicity which are relevant from an applied point of view, the asymptotic stability of the positive equilibrium point which corresponds to the coexistence of the two bacteria. A syntrophic relationship in the anaerobic digestion process is proposed as a real candidate for this model. 1. Introduction. A synthrophic relationship between two organisms refers to a situation where the species exhibit mutualism but where, in contrast to what happens in a purely symbiotic relationship, one of the species can grow without the other. Such a situation can be mathematically formalized as follows: assume that a first species denoted X 1 grows on a substrate S 1 forming an intermediate product S 2 . This intermediate product is required for its growth by a second species X 2 . Since the substrate needed for its growth by the second bacteria is the product of the first bioreaction, the second bacteria cannot grow if the first one is not present.
An important feature of our work is that it involves mutualism. Such an interaction is quite common in nature which explains why a number of models have already been proposed in the literature. Katsuyama et al. [11], describing pesticide degradation proposed a model involving two mutualistic species, while a more general case was considered by Kreikenbohm and Bohl [12]. Since mutualism generally involves species interacting through intermediate products, other studies have considered mutualistic relationships in food webs. For instance, Bratbak and Thingstad [3], and, more recently, Aota and Nakajima [1] have explored mutualism between phytoplankton and bacteria through carbon excretion by phytoplankton. A model studied by Freedman et al. [10] was proposed to explain the observed coexistence of such species. Another characteristic of the present work is that it involves several substrates. A considerable number of studies have appeared in the literature over the last years treating this subject. For a recent review of this field, the reader may refer for instance to Chase and Leibold [5] or Ballyk and Wolkowicz [2]. However, in previous studies the models have usually been very specific. In particular, the mathematical analysis of the models have been designed for specific growth rates that are explicitly given (in most cases as Monod functions).
To extend the study of mutualism to more general systems, we have recently considered more general assumptions, notably with respect to the growth rate functions considered in the models in using qualitative hypotheses, cf. [7] . Furthermore, wa have assumed that the species X 1 may be inhibited by the product S 2 that it has produced itself while the species X 2 was simply limited by S 2 . An example of such interaction is anaerobic digestion in which mutualistic relationships allow certain classes of bacteria to coexist. A mutualistic relation has also been considered in [6]. See [8] for another model of coexistence in the chemostat.
In this paper, following [9], we revisit the model proposed in [7] in incorporating two main changes which significantly extend the range of practical situations covered by the model. First, we assume that there is some S 2 in the influent or, to put it differently, the limiting substrate S 2 on which the species X 2 grows is not only produced by the species X 1 but is also available even if the species X 1 is not present. The second modification to the model is that the second species is assumed to be inhibited by an excess of S 1 , the limiting substrate on which the first species grows. To illustrate the usefulness of such extensions of the original model proposed by El Hajji et al. [7], the biological interpretation of these hypotheses within the context of the anaerobic process is given in the appendix.
The paper is organized as follows. In Section 2, we propose a modification of the model studied in [7] that involves four differential equations. In Section 3, the system is reduced to a planar system, for which equilibria are determined and their local stability properties established. In Section 4, when the system has a single positive equilibrium, the global asymptotic stability is demonstrated using Dulac's criterion (which rules out the possibility of the existence of periodic solutions for the reduced planar system), the Poincaré-Bendixon Theorem and the Butler-McGehee Lemma. Hence, in this case, for all positive initial conditions, the solutions converge to the positive equilibrium point which corresponds to the coexistence of the two bacterial species as observed in real processes. Application to the original model in [7] is given in Section 5. In this section, numerical simulations are presented in the case when the growth functions are of the Monod type. Concluding remarks are given in Section 6. An example of a syntrophic relationship is given in Appendix A as a candidate for this model. The mathematical proofs of the results are given in Appendix B.
2. Mathematical model. Let S 1 , X 1 , S 2 and X 2 denote, respectively, the concentrations of the substrate, the first bacteria, the intermediate product, and the second bacteria present in the reactor at time t. We neglect all species-specific death Author-produced version of the article published in Mathematical Biosciences and Engineering, 2012, 9(3), 627-645. The original publication is available at http://aimsciences. rates and take into account the dilution rate only. Hence our model is described by the following system of ordinary differential equations : where S in 1 > 0 denotes the input concentration of substrate, S in 2 > 0 the input concentration of the intermediate product and D > 0 the dilution rate.
Assume that the functional response of each species μ 1 , μ 2 : R 2 + → R + satisfies : Hypothesis A2 signifies that no growth can take place for species X 1 without the substrate S 1 and that the intermediate product S 2 is necessary for the growth of species X 2 . Hypothesis A3 means that the growth of species X 1 increases with the substrate S 1 but it is inhibited by the intermediate product S 2 that it produces. Hypothesis A4 means that the growth of species X 2 increases with intermediate product S 2 produced by species X 1 while it is inhibited by the substrate S 1 . Note that there is a syntrophic relationship between the two species. We first scale system (1) using the following change of variables and notations : Thus, the dimensionless equations obtained are : Where the functions f 1 , f 2 : R 2 + → R + are defined by Since the functions μ 1 and μ 2 satisfy hypostheses A1-A4, it follows that functions f 1 and f 2 satisfy H1: f 1 , f 2 : R 2 + → R + , of class C 1 , H2: For all (s 1 , s 2 ) ∈ R 2 + , f 1 (0, s 2 ) = 0 and f 2 (s 1 , 0) = 0. H3: For all (s 1 , s 2 ) ∈ R 2 + , ∂f 1 ∂s 1 (s 1 , s 2 ) > 0 and ∂f 1 ∂s 2 (s 1 , s 2 ) < 0.
is a positive invariant attractor of all solutions of system (2).
3. Restriction on the plane. The solutions of system (2) are exponentially convergent towards the set Ω and we are interested in the asymptotic behavior of these solutions. It is enough to restrict the study of the asymptotic behaviour of system (2) to Ω. In fact, thanks to Thieme's results [14], the asymptotic behaviour of the solutions of the restriction of (2) on Ω will be informative for the complete system (see Section 4). In this section, we study the following reduced system which is simply the projection on the plane (x 1 , x 2 ), of the restriction of system (2) on Ω. where Thus, for (3) the state-vector (x 1 , x 2 ) belongs to the following subset of the plane as illustrated in Fig. 1: Figure 1. The set S The point F 0 = (0, 0) is an equilibrium of (3). Besides this equilibrium point the system can have the following three types of equilibrium points: • Boundary equilibria Author-produced version of the article published in Mathematical Biosciences and Engineering, 2012, 9(3), 627-645. The original publication is available at http://aimsciences.org Doi : 10.3934/mbe.2012.9.627 We use the following notations The mapping Figure 2, right) such that: Figure 2. Existence and uniqueness ofx 1 . On the left, the case The mapping . We denote by D 4 ∈]D 1 , D 2 [ the sole real number (see Figure 3, right) such that: The nature of the trivial equilibrium point F 0 is given in the following lemma.
The conditions of existence of the boundary equilibria F 1 and F 2 , and their nature, are stated in the following lemmas.  Figure 3. Existence and uniqueness ofx 2 . On the left, the case Let us now discuss the conditions of existence of positive equilibria F * and their number.
which is a solution of (6) lying in S. this gives Assuming H3 is satisfied, this partial derivative is positive. Hence, equation Recall that x 1 =x 1 is the solution of (4) which, according to Lemma 2, exists and is unique, if and only if D < D 1 . This gives Hence the function F 1 is increasing. Since Φ 1 (s in 1 , 0) = 0, the graph Γ 1 of F 1 has no intersection with the right boundary of the domain S, defined by x 1 = s in 1 . This graph separates S in two regions denoted as the left and right sides of Γ 1 : see Figure  4. This also gives By assuming H4, this partial derivative is positive. Hence, equation x 2 =x 2 is the solution of (5) which, according to Lemma 3 exists and is unique, if and only if D < D 2 . This gives Hence the function F 2 is increasing. Since Φ 2 (x 1 , s in 2 + x 1 ) = 0, the graph Γ 2 of F 2 has no intersection with the top boundary of the domain S, defined by x 2 = s in 2 +x 1 . Thus the point at the very right of Γ 2 lies necessarily on the right boundary of S, defined by x 1 = s in 1 . Hence it lies on the right side of Γ 1 : see Figure 4. The graphs Γ 1 and Γ 2 can intersect or not (see Figures 4, 6 and 7). If they intersect at some point F * = (x * 1 , x * 2 ) then F * is a positive equilibrium. If the point A at the very left of Γ 2 lies on left side of Γ 1 then Γ 1 and Γ 2 intersect in at least one point F * = (x * 1 , x * 2 ). They can have multiple intersections. Generically, they have an odd number of intersections (see Figure 4, center). If the point A at the very left of Γ 2 lies on right side of Γ 1 then Γ 1 and Γ 2 can intersect or not. Generically, they have an even number of intersections (see Figure 4, right). The nature of a positive equilibrium F * is stated in the following lemmas.

Lemma 4. If an equilibrium
It is a saddle point if the opposite inequality is satisfied. The number of equilibria of (3) and their nature are summarized in the following theorem.

Proposition 2.
There are no periodic orbits nor polycycles inside S.
Using the Butler MacGehee Theorem and Poincaré-Bendixon theory, we obtain the asymptotic behaviour of the reduced system (3). (3) has at most one positive equilibrium F * , then for every initial condition in S, the trajectories of system (3) converge asymptotically to :

Theorem 2. Assuming that system
Using Thieme's results [14] on asymptotically autonomous systems, we can deduce that the asymptotic behaviour of the solution of the complete system (2) is the same as the asymptotic behaviour described for the reduced system (3).

A particular case.
In this section we consider the case where s in 2 = 0 and f 2 depends only on s 2 . We obtain the model considered in [7]. The assumptions about f 1 are the same as our assumptions.
Assumption H4 about f 2 becomes H4: For all s 2 ∈ R + , f 2 (s 2 ) > 0. The reduced system is Thus, for (10) the state-vector (x 1 , x 2 ) belongs to (see Figure 5) (7) becomes (8) gives the sole real number D 3 ∈]0, D 1 [ such that (see Figure 2, right): . Hence, the graphs Γ 1 and Γ 2 can intersect at most at one point: see Figure 5. The following proposition which gives the number and nature of equilibria of (10) is a consequence of Theorems 1 and 2.

Proposition 3.
1. If D < D 3 , then (10) has three equilibria: F 0 , which is an unstable node, F 1 , which is a saddle point and F * , which is a globally stable node. 2. If D 3 < D < D 1 , (10) has two equilibria: F 0 , which is a saddle point, and F 1 , which is a globally stable node. 3. If D 1 < D, then (10) has one equilibrium, F 0 , which is a globally stable node.

Growth functions of generalized Monod type.
In this section we consider growth functions f 1 and f 2 of the following form Such functions are simply the product of a Monod function in s 1 by a decreasing function of s 2 . Such functions are currently used in biotechnology when the growth of a functional species is limited by one substrate while inhibited by another. Such situations are common in water treatment technology such as the denitrification (limited by the nitrate and inhibited by the dissolved oxygen) or in anoxic or anaerobic hydrolysis (limited by the slowly-biodegradable substrates while inhibited by an excess of oxygen), processes which can be modeled this way (cf. [15]). One can easily check that (12) satisfy Assumptions H1 to H4. Straightforward calculations give Hence equation F 1 (x 1 ) = F 2 (x 1 ), giving the abscissa of positive equilibria, is an algebraic equation of degree 2. Thus, it cannot have more than two solutions. Therefore, the situation shown in the center of Figure 4 with three positive equilibria, is excluded. However, the situation depicted on the right in Figure 4, showing two positive equilibria, can occur. For instance, consider the following values of the parameters: Then D 1 = 6/5, D 3 = 8/9, D 2 = 3/5. There is another bifurcation value, D = 1, which correspond to the case when the graphs Γ 1 and Γ 2 are tangent: see Figure 7. For this example, five cases can occur (see Figure 6): (3) where f 1 and f 2 are given by (12) with parameters (13).
1. If D < 3/5, then the system has four equilibria: F 0 which is an unstable node, F 1 and F 2 , which are saddle points, and F * , which is a stable node. 2. If 3/5 < D < 8/9, then the system has three equilibria: F 0 and F 1 , which are saddle points and F * , which is a stable node. 3. If 8/9 < D < 1, then the system has four equilibria: F 0 and F * 1 , which are saddle points and F 1 and F * 2 , which are stable nodes. 4. If 1 < D < 6/5, then the system has two equilibria: F 0 , which is a saddle point, and F 1 which is a stable node. 5. If D > 6/5, then the system has one equilibrium: F 0 , which is a stable node. Figure 7. The non-hyperbolic cases. When D = 6/5, F 0 and F 2 coalesce. When D = 8/9, F * 1 and F 1 coalesce (saddle node bifurcation). When D = 1, F * 1 and F * 2 coalesce (saddle node bifurcation). When D = 6/5, F 0 and F 2 coalesce. (3) where f 1 and f 2 are given by (12) with parameters (13). According to Proposition 4, in the case when 8/9 < D < 1, a bistability phenomenon occurs. According to the initial condition, both species can coexist at equilibrium F * 2 , or species x 2 goes to extinction at equilibrium F 1 . This phenomenon is illustrated numerically, with D = 0.95, in Figure 8. In the center, the phase portrait. On the right, the isoclines. (3) where f 1 and f 2 are given by (12) with parameters now given by

Consider again system
The bifurcation values are D 1 = 4/3 and D 2 = 21/16. If D > max(D 1 , D 2 ), for instance for D = 3/2, one obtains a bistability phenomenon corresponding to case (3) of Theorem 1, with two positive equilibria. Depending on the initial conditions, both species can coexist at equilibrium F * 2 , or both species go to extinction at equilibrium F 0 . This phenomenon is illustrated numerically in Figure 9. the input. For one of the populations, one resource is needed for growth and the other is inhibitory; for the other population, in contrast, the roles of the resources are reversed. One of the populations produces as a by-product the resource that is inhibitory to itself but needed for growth by the other population. The analysis of this model predicts that, under general and natural assumptions of monotonicity on the functional responses, the stable asymptotic coexistence of the two bacteria is possible. Extending the model studied in [7] by considering that there may have some S 2 in the influent and using a more general class of kinetics functions, we show that the qualitative behavior of the system can be significantly modified. In particular, if S in 2 > 0, the situation where X 1 is washed out but where X 2 > 0 is possible. In addition, positing the possible inhibition of the growth of X 2 by S 1 yields bistability. This work offers the perspective of providing an "elegant" formulation for synthesizing all results concerning the analysis of two-steps models of bioprocesses currently used for monitoring and control purposes.
Appendix A. The anaerobic digestion process : An example of a synthrophic relationship. "Methane fermentation" or "anaerobic digestion" is a process that converts organic matter into a gaseous mixture composed mainly of methane and carbon dioxide (CH 4 and CO 2 ) through the action of a complex bacterial ecosystem (cf. Fig.10). It is often used for the treatment of concentrated wastewaters or to convert the excess sludge produced in wastewater treatment plants into more stable products. There is also considerable interest in digesters fed with plant biomass, since the methane produced can be used profitably as a source of energy. It is usually considered that a number of metabolic groups of bacteria are involved sequentially in the fermentation process.
One specific characteristic of the anaerobic process is that within the involved groups of bacteria, there exist populations exhibiting obligatory mutualistic relationships. Such a syntrophic relationship is necessary for the biological reactions to be thermodynamically possible. In the first stages of the reactions (which together are called "acidogenesis"), some hydrogen is produced. However, in El Hajji et al. [7], this production of hydrogen at this stage was ignored (compare Fig.10 with Fig.1 of [7]). This hypothesis constitutes the first novelty with respect to [7]. It is to be noticed that an excess of hydrogen in the medium inhibits the growth of another bacterial group called "acetogenic bacteria". Their association with H 2consuming bacteria is thus necessary for the second stage of the reaction to occur. Such a syntrophic relationship has been pointed out in a number of experimental works (cf. for instance the seminal work by [4]). Let us consider the subsystem of the anaerobic system where the VFA (for Volatile Fatty Acids) are transformed into H 2 , CH 4 and CO2. The corresponding biological reactions can be formalized as a first bacterial consortium X 1 (the acetogens) transforming S 1 (the VFA) into S 2 (the hydrogen) and acetate (cf. Fig.10). Then, a second species X 2 (the hydrogenotrophic-methanogenic bacteria) grows on S 2 . In practice, acetogens are inhibited by an excess of hydrogen and methanogens by an excess of VFA. Thus, it is further assumed that X 1 is inhibited by S 2 and X 2 by S 1 . The last inhibiting relationship constitutes the second innovation with respect to [7]. This situation is precisely that dealt with by the model (1).  Next, it must be proved that the solution is bounded. Let z 1 = s 1 + x 1 , theṅ z 1 = −D(z 1 − s in 1 ) from which it follows that : Thus s 1 (t) and x 1 (t) are positively bounded. Let z 2 = s 2 + x 2 − x 1 , thenż 2 = −D(z 2 − s in 2 ) from which it can be deduced that: Thus s 2 (t) and x 2 (t) are positively bounded. Hence, the solution is defined for all positive t. From (15) and (16), it can be deduced that the set Ω is an invariant set which is an attractor.
Proof of Lemma 1. The Jacobian matrix J of (3), at point (x 1 , x 2 ), is given by: where the functions are evaluated at s in The Jacobian matrix at F 0 is given by: Proof of Lemma 2. An equilibrium F 1 = (x 1 , 0) exists if and only if x 1 =x 1 ∈ ]0, s in 1 [ is a solution of (4). Let ψ 1 (x 1 ) = Φ 1 (x 1 , 0). Then By assuming H3, ψ 1 (x 1 ) < 0. Since ψ 1 (0) = D 1 , and ψ 1 (s in 1 ) = 0, equation (4) admits a solution in the interval ]0, s in 1 [ if and only if D < D 1 . If this condition is satisfied then (4) admits one single solution since the function ψ 1 (.) is decreasing (see Figure 2). The Jacobian matrix at F 1 is given by: where the functions are evaluated at (s in 1 −x 1 , s in 2 +x 1 ). The eigenvalues are Thus F 1 is a saddle point if Φ 2 (x 1 , 0) > D. If D 1 < D 2 , this condition is satisfied for all D < D 1 . If D 2 < D 1 , it is satisfied for all 0 < D < D 3 (see Figure 2). F 1 is a stable node if D 3 < D < D 1 and D 2 < D 1 .
Given assumptions H3 and H4, the product of the partial derivatives is negative. Therefore, the determinant is positive if F 1 (x * 1 ) > F 2 (x * 1 ) and negative otherwise.
) and a saddle point otherwise.
Proof of Theorem 1. We give the details of the proof in the case where The other cases can be proved similarly. The washout equilibrium point F 0 always exists. According to Lemma 1, it is an unstable node if D < min(D 1 , D 2 ) and a stable node for D > max(D 1 , D 2 ).
Consider first the case D < min(D 1 , D 2 ). Lemmas 2 and 3 show that the system admits two boundary equilibria F 1 and F 2 , which are saddle points. Since Γ 2 starts on the left side of Γ 1 and ends on its right, the graphs Γ 1 and Γ 2 have at least one intersection. In the generic case, where all intersections are transverse, the graphs can have an odd number of intersections (see Figure 4, center). According Author-produced version of the article published in Mathematical Biosciences and Engineering, 2012, 9(3), 627-645. The original publication is available at http://aimsciences.org Doi : 10.3934/mbe.2012. 9.627 to Lemma 4, these intersections are alternatively saddle points and stable nodes. The one at the very left of these positive equilibria is a stable node since at this point the condition F 1 (x * 1 ) > F 2 (x * 1 ) is satisfied. Consider now the case D > max(D 1 , D 2 ). Lemmas 2 and 3 show that the system admits no boundary equilibria F 1 and F 2 . Since Γ 2 starts and ends on the right side of Γ 1 , the graphs Γ 1 and Γ 2 can either intersect or not. In the generic case, where all intersections are transverse, the graphs can have an even number of intersections (see Figure 4, right). According to Lemma 4, these intersection are alternatively saddle points and stable nodes. The one at the very left of these positive equilibria is a saddle point since at this point the condition F 1 (x * 1 ) < F 2 (x * 1 ) is satisfied.
From the Dulac criterion [13], it can be deduced that the system (17) has no periodic trajectory. Hence (3) has no periodic orbit in S.
Proof of Theorem 2. We give the details of the proof in the case when D < min(D 1 , D 2 ).
The other cases can be proved similarly. Let x 1 (0) > 0, x 2 (0) > 0 and ω the ω-limit set of (x 1 (0), x 2 (0)). ω is an invariant compact set and ω ⊂S. Assume that ω contains a point M on the x 1 x 2 axis : • M cannot be F 0 because F 0 is an unstable node and cannot be a part of the ω-limit set of (x 1 (0), x 2 (0)), . ω is not reduced to F 1 (respectively to F 2 ).