Bistability in a system of two species interacting through mutualism as well as competition: Chemostat vs. Lotka-Volterra equations

We theoretically study the dynamics of two interacting microbial species in the chemostat. These species are competitors for a common resource, as well as mutualists due to cross-feeding. In line with previous studies (Assaneo, et al., 2013; Holland, et al., 2010; Iwata, et al., 2011), we demonstrate that this system has a rich repertoire of dynamical behavior, including bistability. Standard Lotka-Volterra equations are not capable to describe this particular system, as these account for only one type of interaction (mutualistic or competitive). We show here that the different steady state solutions can be well captured by an extended Lotka-Volterra model, which better describe the density-dependent interaction (mutualism at low density and competition at high density). This two-variable model provides a more intuitive description of the dynamical behavior than the chemostat equations.

(3) These obey the differential equations dPi dt = Φ(S i − P i ) (i = 0, 1, 2) for which the solutions are given by exponentially decreasing functions P i (t) =S i + (P i (0) −S i )e −Φt . After a transient time, when t 1 Φ , the equations reduce to P i (t) =S i . Hence, we obtain the relations given in Eq.(4).
We can use these equations to eliminate S 0 , S 1 and S 2 from the chemostat equations (Eq.(1)), and obtain a 2 dimensional system that only depends on the species concentrations (Eq.(5)). 1/6 This system has a different transient behavior than the original chemostat system, but it converges to the same steady-state attractor as the original system. Indeed, the variables P i are decreasing exponentially with time constant 1 Φ , implying the same qualitative dynamics in both systems for t 1 Φ . However, as the nutrient concentrations cannot be negative this gives some constraints for the concentrations in order to be realistic. When imposing S i > 0 in Eq.4, we obtain the restrictions Eq.(6): Using (4) to eliminate the nutrients S i (i = 0, 1, 2) in the growth rates (2), it is clear that the growth rates are highly nonlinear functions and difficult to interpret. Looking at the behavior of the species in extreme situations, for example when the nutrient concentrations are very low or very high, gives important insights about the dynamics of this system.
The Monod terms S/(K + S) increase from zero to one with S. For high nutrient concentrations, S K, this term approximates one, whereas for low concentrations, S K, it behaves as a linear function S/K. Both growth rates in the chemostat system (2) depend on two nutrients, hence we can consider four different limiting We show that in the first three cases we obtain equations that are normal growth equations, similar to LV systems. The last situation is the most interesting case to obtain a system with the same dynamics as the chemostat system.
This is the simplest case. We obtain exponential growth equations Eq.(7) for both species since there are enough nutrients available. They do not depend on each other in a competitive nor mutualistic manner. As long as the growth rate is larger than the flow rate the populations increase exponentially. This ends when the consumption of nutrients has increased so much that the condition of high nutrient concentrations is no longer valid.
Due to the relations Eq.(4) the growth rates reduce to linear functions of the species concentrations: This is a competitive Lotka-Volterra system. The growth is limited by the compound S 0 in this case, so that the competition for S 0 dominates the system. Only one species survives due to competitive exclusion.
The growth rates are linear functions as well. However, here the growth is limited by the cross-feeding nutrients S 1 and S 2 , so that we obtain mutualistic Lotka-Volterra equations: Here, the concentrations keep increasing until the condition K 0i S 0 is no longer valid and the competition for S 0 becomes significant.
The growth is limited by S 0 , as well as by S 1 and S 2 . Therefore both competition and mutualism play a significant role. The growth rates are now given by: with r i = µ i /(K 0i K ji ). As we are interested in explaining bistability in the chemostat system, where competition as well as mutualism are important, we use these equations to simplify the system.
Theoretically, the simplification is valid if the conditions K 0i S 0 & K ji S j are fulfilled. However, if this is not the case, it is possible to find a parameter set with higher values of K for which the system has the same steady states. To understand this, we consider the steady state equations dXi dt = 0 and dSi dt = 0 of Eq.(1). It is possible to use a different parameter set to obtain the same solutions for the steady state concentrations of X i and S i . Such a change will make the dynamics faster or slower, but as the system converges to the same fixed points, it is the same qualitative behavior.
Therefore, as long as Eq.(11) still yields the solutions for S 0 and S j , we can change the parameters µ, Φ and K such that we find values for K for which K 0i S 0 and K ji S j . We show this with an example in Fig.1. In Fig.1(a) we show the dynamics simulated with parameter set as defined in the main text, with µ = 1600 1600 , which leads to a bistable situation. As K 0i = 200 the conditions K 0i S 0 are not satisfied. In Fig.1(b) the values of K are increased, so that K = 100K, µ is then changed according to Eq.(12), which follows from Eq.(11). Therefore the steady state concentrations remain the same, as can be seen in Fig.1(b). The species X 1 and X 2 either coexist (state 1) or die out (state 2). The end concentration of S 0 is S 0 =S 0 = 50 or S 0 = 26. Therefore the condition S 0 K 0i is not satisfied. The end concentrations of S 1 , equal to those of S 2 , are very small so that S i K ij is satisfied. (b) If we increase the values K such that K = 100K and we change µ so that Eq.(12) is satisfied, then the qualitative dynamics is not changed. The steady state concentrations are the same as in (a), but the transient dynamics is (moderately) different. The condition S 0 K 0i is well satisfied for this parameter set, so that the approximation to Eq.(13) is valid.
By substituting the growth rates Eq.(10) in Eq.(5), we obtain the following system (with i, j = 1, 2 and i = j): The relationships between these new parameters and those of the chemostat are 4/6 given by: The dynamics of system Eq.(13) can be visualized in the phase plane. As an example we did this for the situation where the behavior is bistable (Fig.2). The nullclines of the system are hyperbolic functions (as well as the x-and y-axis). The realistic region of the phase plane is limited by the restrictions given by Eq.(6), hence the only physical regime is situated inside the red lines in Fig.2. To further simplify the equations, we examine the effect of the different nonlinear terms in the growth: c ij X 2 j , f ii X 2 i and f ij X i X j . All these nonlinear terms become particularly relevant when the concentrations are large. In the realistic regime, this means that the growth is negative for large concentrations, which is caused by the competition. This is modeled via the term c ij X 2 j , which is therefore necessary for the dynamics. The term f ii X 2 i causes the divergence in the growth when the branches of the hyperbola in the unrealistic regime are crossed, however they do not affect the qualitative behavior in the realistic regime.
The last term f ij X i X j has the effect of increased competition or self-inhibition at high concentrations, where the growth already is negative. Therefore, the qualitative dynamics in the physical regime is not affected by elimination of the two terms f ii X 2 i and f ij X i X j , as the competition at high concentrations is incorporated in the term c ij X 2 j . The minimal model we use to describe the dynamics in the physical regime is therefore given by Eq.(15), which has parabolic nullclines. In this article, we refer to these equations as the extended Lotka-Volterra model. dX 1 dt = r 1 (a 1 − b 11 X 1 + b 12 X 2 − c 1 X 2 2 )X 1 − dX 1 dX 2 dt = r 2 (a 2 + b 21 X 1 − b 22 X 2 − c 2 X 2 1 )X 2 − dX 2