Strong cooperation or tragedy of the commons in the chemostat

1.   Introduction

Many organisms exhibit cooperative behavior that benefits the whole group. Examples of such behavior in microbial populations are secreted products such as extracellular enzymes that digest nutrients, siderophores that acquire iron, or exopolysaccharides that support biofilm growth ^[16]. Such products are known as public goods. While common goods production is advantageous, and often critical to the growth of the population, it is also costly. Cheating individuals that choose not to cooperate could benefit from the fruits of labor of individuals that do, and would have a competitive advantage over the cooperators. One would expect that cheaters should be able to successfully invade a population of cooperators. This situation reflects a Prisoner's dilemma in evolutionary game theory, where cheating is the evolutionary stable strategy ^[12]. However, when cheaters start to dominate cooperators, the common good would become scarce, ultimately leading to the demise of cooperators and cheaters alike. This idea is known as the Tragedy of the Commons ^[8,6]. It has since then spurred a body of research that continues to grow up to this day, devoted to identifying mechanisms that explain the evolution of cooperation.

Although there has been extensive experimental work on the topic of the Tragedy of the Commons, perhaps somewhat surprisingly, there exist very few mathematical models based on first principles that express population growth, public good production and cheating behavior, for which the tragedy of the commons can be verified by means of a mathematical proof. In ^[11] we proposed a minimal chemostat model that incorporates cooperative and cheating behavior in microbial populations, and we established mathematically that the tragedy of the commons does indeed occur. The main goal of this paper is to generalize the chemostat model of ^[11] to investigate the fate of the population when the cooperator has the ability to evolve its ecological characteristics, specifically its per capita uptake rate function and yield constant. Although we find that the tragedy continues to hold in many cases, we also find scenarios where the cooperator not only survives, it even outcompetes the cheater by driving it to extinction. We will show that his happens when the break-even concentration for the growth nutrient of the cooperator evolves to become lower than that of the cheater. Interestingly, under no circumstances can the cooperator and cheater coexist in this generalized chemostat model: either the tragedy of the commons occurs, or the cooperator outcompetes the cheater, hereby favoring a very strong form of cooperation.

We propose the generalized chemostat model and prove our main result in Section 2. Section 3 contains a discussion of our main result, and conclusions are found in Section 4.

2.   Model and main results

We consider a general chemostat model with positive dilution rate $D$ and positive input nutrient concentration $S^0$. Nutrient, cooperator and cheater concentration are denoted by $S$, $X_1$ and $X_2$ respectively. The cooperator produces an enzyme which has concentration $E$, and this enzyme is used to convert the nutrient into a processed nutrient with concentration $P$. Once processed, nutrient is available for growth of the cooperator and cheater, but the cooperator also diverts a fraction towards production of the enzyme; the cheater does not produce the enzyme. The enzymatic reaction converting nutrient into processed nutrient is given by:

The rate of this reaction is proportional to $E$ and to a possibly nonlinear $C^2$ function $G(S)$ which is assumed to be zero at zero, strictly increasing ($G'(S)>0$ for all $S>0$) and concave ($G''(S)\leq 0$ for all $S$). The prototypical example is a mass action rate, corresponding to $G(S) = kS$ for some positive rate constant $k$. The per capita uptake rates of cooperator and cheater are given by $F_1(P)/\gamma_1$ and $F_2(P)/\gamma_2$, where the $\gamma_i$ are positive yield constants, and the functions $F_i(P)$ are assumed to be $C^1$, zero at $P = 0$, and strictly increasing ($F_i'(P)>0$ for all $P>0$). In applications, one often picks Michaelis-Menten functions for the $F_i(P)$ which have the form $mP/(a+P)$ where $m$ and $a$ are positive constants, but Hill functions of the form $mP^n/(a^n+P^n)$ where $n$ is a positive integer, and $m$ and $a$ also are positive constants, are allowed here as well.

As mentioned, the cooperator only uses a fraction of the available processed nutrient towards growth, and we denote this fraction by a positive constant $q<1$. The remaining fraction $1-q$ is used to produce the enzyme. The conversion efficiency for the conversion of $P$ into $E$ by the cooperator, is a positive constant $\eta$.

These considerations lead to the following model:

defined on the forward invariant set $\mathbb{R}^5_+$.

By scaling $s = S, p = P, e = E/(\eta \gamma_1), x_1 = X_1/\gamma_1, x_2 = X_2/\gamma_2$ and defining $g(s) = \eta \gamma_1G(S), f_i(p) = F_i(P)$ for $i = 1, 2$, and $d = D$ and $s^0 = S^0$, we obtain the scaled model:

Defining two new variables

which satisfy the following equations:

hence $m(t)\to s^0$ and $v(t)\to 0$ as $t\to +\infty$. To understand the behavior of $(2.6)-(2.10)$, we therefore first investigate the following limiting system:

which is defined on the forward invariant state space $\{x_1\geq 0, x_2\geq 0, p\geq 0, p+x_1/q+x_2 \leq s^0\}$. It will be more convenient to analyze this system with one more change of variable:

instead of using the variable $p$:

which is a system with forward invariant state space $\Omega = \{x_1\geq 0, x_2\geq 0, x_1/q+x_2\leq w\leq s^0\}$. The Jacobian matrix of this system is:

where we have suppressed the arguments of the functions $g$, $f_i$ and their derivatives to avoid a cumbersome notation. Our assumptions about these functions imply that the Jacobian matrix has the sign structure displayed below, where $*$ means that the sign is not fixed, but depends on the state where the Jacobian is evaluated, and where $+$ and $-$ indicate a non-negative and non-positive value for every state in the state space $\Omega$ respectively:

The key observation is that this sign structure implies that system $(2.11)-(2.13)$ is a 3-dimensional competitive system, see ^[13] for more on this particular class of systems. This means that backward-time solutions of this system remain ordered with respect to the partial order generated by the cone $K = \{(w, x_1, x_2)\, |\, w\leq 0, x_1\geq 0, x_2\geq 0\}$. The relevance of this fact is that a Poincaré-Bendixson theory is available for this class of systems which makes a global analysis feasible in ways comparable to planar systems. Note also that the sign structure of the Jacobian implies that the subsystem $(2.11)-(2.12)$ on the part of the boundary of $\Omega$ where $x_2 = 0$ is a cooperative system in the usual sense.

We begin our investigation of the dynamics of system $(2.11)-(2.13)$ by recalling from ^[11] the dynamics on the part of the boundary of $\Omega$ where $x_2 = 0$, which is easily seen to be a forward invariant set. Let $p_1$ denote the solution to the equation $qf_1(p) = d$. Since $f_1(p)$ is increasing in $p$, this solution is unique whenever it exists. When it does not exist, which happens if and only if $\lim_{p\to \infty} qf_1(p)\leq d$, we define $p_1 = +\infty$. For example, when $f_1(p) = m_1p/(a_1+p)$ (Monod function with $a_1$ and $m_1$ positive), and assuming that $qm_1>d$, we can easily compute $p_1$:

The quantity $p_1$ represents the break-even concentration of the processed nutrient for the cooperator. We introduce the assumption that:

which says that the break-even concentration of the cooperator should be less than the input nutrient concentration. If ${\bf H}$ does not hold, then $x_1(t)\to 0$, implying the same for $w(t)$ and $x_2(t)$; this trivial extinction case is not of interest. Next define the function $h:[0, s^0)\to \mathbb{R}_+$, by

We showed in ^[11] that $h$ is zero at zero, that it is strictly increasing and strictly convex (i.e. $h''(w)>0$) with a vertical asymptote at $w = s^0$. Consider the following equation:

Then strict convexity of the function on the left, and linearity of the function on the right, implies that this equation has at most $2$ solutions in the interval $(0, s^0)$, and will have exactly $2$ solutions for all sufficiently small values of $p_1$. In fact, we showed in ^[11] that:

Lemma 2.1. Let ${\bf H}$ hold. If equation $(2.15)$ has two solutions $w_1<w_2$, then system $(2.11)-(2.13)$, restricted to the invariant set $\{x_2 = 0\}$ has 3 steady states, $e_0 = (0, 0)$, $e_1 = (w_1, h(w_1))$ and $e_2 = (w_2, h(w_2))$. The steady states $e_0$ and $e_2$ are locally asymptotically stable, and $e_1$ is a saddle with one-dimensional stable manifold $W_s$, and one-dimensional unstable manifold $W_u$. The stable manifold $W_s$ intersects the boundary of $\Omega\cap \{x_2 = 0\}$ in two points, one on the boundary $x_1 = qw$, the other on the boundary $w = s^0$, forming a separatrix: Initial conditions below $W_s$ give rise to solutions converging to $e_0$, whereas initial conditions above $W_s$ give rise to solutions converging to $e_2$, yielding a bistable system.

Figure 1 in the S3 Appendix in ^[11] illustrates the geometry of the nullclines. Using the fact that this planar system is cooperative and irreducible, one may conclude that $W_s$ is an unordered set, meaning that no two of its points are related by the usual component-wise ordering while $W_u$ is a totally ordered curve consisting of a heteroclinic orbit connecting $e_1$ to $e_0$ and a heteroclinic orbit connecting $e_1$ to $e_2$.

The dynamics of system $(2.11)-(2.13)$ restricted to the part of the boundary of $\Omega$ where $x_1 = 0$ (also a forward invariant set) is trivial: $w(t)\to 0$ as $t\to \infty$, and hence $x_2(t)\to 0$ as well, since the argument of the function $f_2$ in the $x_2$-equation becomes arbitrarily small for all sufficiently large times. The remaining parts of the boundary of $\Omega$ are not forward invariant. In fact, solutions starting there, enter ${\rm int}(\Omega)$ instantaneously. This feature will be used later.

We are now ready to investigate the global behavior of system $(2.11)-(2.13)$ in $\Omega$. Before stating the precise result, we also define $p_2$, the break even concentration of the processed nutrient for the cheater, as the unique solution $p$ to the equation $f_2(p) = d$ (as before, we define $p_2 = +\infty$ if there is no solution).

Theorem 2.2. Assume that the conditions of Lemma $2.1$ hold. Then system (2.11)-(2.13) has 3 steady states $e_0 = (0, 0, 0)$, $e_1 = (w_1, h(w_1), 0)$ and $e_2 = (w_2, h(w_2), 0)$. Moreover,

1. Strong cooperation or tragedy: If $p_1<p_2$, then $e_0$ and $e_2$ are locally asymptotically stable, and $e_1$ is a saddle point with 2-dimensional stable, and 1-dimensional unstable manifold. All solutions converge to one of these 3 steady states. In particular, $x_2(t)\to 0$ for every solution.

2. Tragedy: If $p_2<p_1$, then $e_0$ is locally asymptotically stable, $e_1$ is a saddle point with 1-dimensional stable and 2-dimensional unstable manifold, and $e_2$ is a saddle point with 2-dimensional stable manifold contained in the boundary of $\Omega$ where $x_2 = 0$, and 1-dimensional unstable manifold. All solutions with $x_2(0)>0$ converge to $e_0$.

Proof. Using Lemma $2.1$, it is easy to verify that there are exactly 3 steady states $e_0$, $e_1$ and $e_2$ which are all contained in the part of the boundary of $\Omega$ where $x_2 = 0$. In particular, there are no steady states where both cooperator and cheater coexist at positive levels, a feature that will reveal its importance later on.

1. Suppose that $p_1<p_2$. Then the eigenvalues of the Jacobian evaluated at $e_i = (w^i, x_1^i, 0)$, $i = 0, 1, 2$ are given by the eigenvalues of the matrix

and the number

Evaluating these at $e_0$, we obtain the eigenvalue $-d$ with multiplicity 3. It was shown in ^[11] that $J_s^1$ has one positive and one negative eigenvalue; moreover, $\lambda_3^1 = f_2(p_1)-d<f_2(p_2)-d = 0$ because $p_1<p_2$. Similarly, we showed in ^[11] that $J_s^2$ has two eigenvalues with negative real part; moreover, $\lambda_3^2 = f_2(p_1)-d<f_2(p_2)-d = 0$. These calculations prove the statements regarding the nature of the steady states and the dimension of the stable and unstable manifold of $e_1$.

We are left to show that all solutions converge to one of the 3 steady states. This follows from Lemma $2.1$ for solutions on the boundary where $x_2 = 0$, and the same is true for solutions on the boundary where $x_1 = 0$, as explained earlier. So we consider a solution $y(t) = (w(t), x_1(t), x_2(t))$ with $x_1(0)>0$ and $x_2(0)>0$. Our goal is to show that the omega limit set of this solution, $\omega(y(0))$, is a singleton consisting of one of the 3 steady states. If either $e_0$ or $e_2$ belong to $\omega(y(0))$, then $\omega(y(0))$ is indeed the singleton $\{e_0\}$ or $\{e_2\}$ respectively, because both are locally asymptotically stable. So we assume that $\omega(y(0))$ does not contain $e_0$, nor $e_2$. Then either $e_1\in \omega(y(0))$ or $e_1\notin \omega(y(0))$. If $e_1\in \omega(y(0))$, then there are two possibilities: $y(0)\in W_s(e_1)$, the stable manifold of $e_1$, and then $\omega(y(0)) = \{e_1\}$; or $y(0)\notin W_s(e_1)$, but then since $e_1\in \omega(y(0))$, the Butler-McGehee Lemma (p.12 in ^[14]) implies that $\omega(y(0))$ must intersect the 1-dimensional unstable manifold $W_u(e_1)$ of $e_1$ in the $x_2 = 0$ face at a point distinct from $e_1$. But as noted following Lemma 1, $W_u(e_1)\setminus \{e_1\}$ consists of two monotone orbits, one connecting to $e_0$, the other to $e_2$. It follows that $\omega(y(0))$ contains either $e_0$ or $e_1$, a contradiction to earlier arguments. Thus, we are left to consider the case that $e_1\notin \omega(y(0))$. Then $\omega(y(0))$ does not contain any of the steady states, and it follows from Hirsch's Theorem (Theorem 4.1 in chapter 3 of ^[13]) that $\omega(y(0))$ must be a periodic orbit $O$. Moreover, it is easy to see that any periodic orbit of system $(2.11)-(2.13)$ must belong to ${\rm int}(\Omega)$: the invariant sets on the boundary where $x_1 = 0$, and where $x_2 = 0$ do not contain periodic orbits, see Lemma $2.1$ and the discussion following it; the remaining parts of the boundary are instantaneous repellors. We aim to force a contradiction to the existence of a periodic orbit by using Proposition 4.3 in chapter 3 of ^[13] but first we must prepare its application by extending system $(2.11)-(2.13)$ to the extended state space

by extending the functions $f_i$, $i = 1, 2$ and their domains, to domains $\mathbb{R}$ and extended $C^1$ functions such that $f'_{i, e}(p)>0$ for all $p$ in $\mathbb{R}$. This implies that the extended system on $\Omega_e$ is still a 3-dimensional competitive system because the sign structure of the Jacobian matrix remains unchanged by construction. Moreover, it is not difficult to check that this extension does not introduce new steady states in ${\rm int}(\Omega_e)$ because the $x_1$ and $x_2$ steady state equations would imply that $w-x_1/q-x_2$ must equal both $p_1$ and $p_2$, contradicting our assumption that $p_1<p_2$.

As noted earlier, the time-reversed extending system $(2.11)-(2.13)$ preserves the partial order generated by the cone $K = \{(w, x_1, x_2)\, |\, w\leq 0, x_1\geq 0, x_2\geq 0\}$. We write $y\ge_K 0$ for $y\in K$, $y\le_K y'$ if $y'-y\ge_K 0$, and $y\ll_K y'$ if $y'-y\in {\rm int}(K)$. The point $A' = (s^0, 0, 0)\in \Omega$ is easily seen to satisfy $A'\ll_K o$ for every $o\in O$ so the same holds for point $A = (s^0-\epsilon, \epsilon, \epsilon)\in {\rm int}(\Omega)$ for small $\epsilon>0$. Similarly, $B = (\epsilon, \epsilon^{-1}, \epsilon^{-1})$ belongs to the interior of $\Omega_e$ and satisfies $o\ll_K B$ for all $o\in O$ for all small $\epsilon>0$. Fix suitably small $\epsilon>0$. Let $[A, B]_K$ denote the box-set $\{y = (w, x_1, x_2)\in \mathbb{R}^3: A\le_K y\le_K B\} = [\epsilon, s^0-\epsilon]\times[\epsilon, \epsilon^{-1}]^2$.

Consequently, we have

Then Proposition 4.3 in chapter 3 of ^[13] implies that $[A, B]_K$ contains a steady state. However, this contradicts the fact, established above, that the extended system has no steady states in ${\rm int}(\Omega_e)$. This completes the proof in case $p_1<p_2$.

2. Suppose that $p_1>p_2$. Similar calculations of the Jacobian matrices as in the previous case show the nature of the 3 steady states and their stable and unstable manifold. The main difference is that here $\lambda_3^i = f_2(p_1)-d>f_2(p_2)-d = 0$ for $i = 1, 2$ because now $p_1>p_2$.

We are left to show that all solutions with $x_2(0)>0$, converge to $e_0$. Pick such a solution $y(t) = (w(t), x_2(t), x_2(t))$. Then $y(0)$ does not belong to the stable manifolds of $e_1$ or $e_2$ (as these are contained in the part of the boundary of $\Omega$ where $x_2 = 0$). We claim that neither $e_1$, nor $e_2$ can belong to $\omega(y(0))$. Indeed, if $e_1\in \omega(y(0))$ then the Butler-McGehee Lemma implies that $\omega(y(0))$ intersects $W_s(e_1)\setminus\{e_1\}$, but then backward time invariance of omega limit sets would imply that $\omega(y(0))$ contains points outside of $\Omega$ which is impossible. A similar argument rules out that $e_2$ belongs to $\omega(y(0))$. Therefore, either $e_0\in \omega(y(0))$, in which case also $\omega(y(0)) = \{e_0\}$ because $e_0$ is asymptotically stable, establishing the desired result; or $e_0\notin \omega(y(0))$, but then again by Hirsch's Theorem, $\omega(y(0))$ must be a periodic orbit. We can now argue as in the proof of the previous case, by extending the system without introducing any steady states in the interior of its state sp ace, and showing that this leads to a contradiction as above.

Having proved Theorem $2.2$, we can now invoke standard results from the theory of asymptotically autonomous systems, see e.g. Appendix F in ^[14] (note that the reduced system has no cycle of steady states), to characterize the global dynamics of the scaled system $(2.6)-(2.10)$, and then also of the unscaled system $(2.1)-(2.5)$ in terms of the break-even concentrations $P_1$ of the cooperator, and $P_2$ of the cheater, which are defined as the unique solutions of the equations $qF_1(P) = D$ and $F_2(P) = D$ respectively (and defined as $+\infty$ if no solution exists). Of course, $P_i = p_i$ since $F_i = f_i$. The Main Result of this paper is:

Theorem 2.3. Assume that $P_1<S^0$, and that equation (2.15) has two distinct solutions $w_1$ and $w_2$ with $P_1<w_1<w_2<S^0$. Then system (2.1)-(2.5) has exactly 3 steady states $SS_0 = (S^0, 0, 0, 0, 0)$, $SS_1 = (S^0-w_1, P_1, \eta (1-q) \gamma_1 (w_1-P_1), q \gamma_1(w_1-P_1), 0)$ and $SS_2 = (S^0-w_2, P_1, \eta (1-q)\gamma_2(w_2-P_1), q\gamma_2(w_2-P_1), 0)$. Moreover,

1. Strong cooperation or tragedy: If $P_1<P_2$, then $SS_0$ and $SS_2$ are locally asymptotically stable, and $SS_1$ is a saddle point with 4-dimensional stable, and 1-dimensional unstable manifold. All solutions converge to one of these 3 steady states. In particular, $X_2(t)\to 0$ as $t\to +\infty$ for every solution. The system is bistable: With the exception of solutions starting on the stable manifold of $SS_1$, all solutions converge to either $SS_2$ (strong cooperation), or to $SS_0$ (tragedy).

2. Tragedy: If $P_2<P_1$, then $SS_0$ is locally asymptotically stable, $SS_1$ is a saddle point with 3-dimensional stable and 2-dimensional unstable manifold, and $SS_2$ is a saddle point with 4-dimensional stable manifold contained in the boundary of $\mathbb{R}^5_+$ where $X_2 = 0$, and 1-dimensional unstable manifold. All solutions with $X_2(0)>0$ converge to $SS_0$ (tragedy).

Figure 1 illustrates the case that $P_1<P_2$ with $G(S) = kS$ and Monod uptake functions $F_i(P) = \frac{m_iP}{a_i+P}$ where $m_1 = 6, a_1 = 0.025, m_2 = 5, a_2 = 0.05$ and where $q = 0.8$ and $k = 20, D = 1, S^0 = 1, \gamma_i = 1, \eta = 1$. Initial conditions for the Tragedy outcome are $S = 0.3, P = 0.0, E = 0.01, X_1 = 0.10, X_2 = 0.55$. Initial conditions for the Strong Cooperation outcome are $S = 0.48, P = 0.0, E = 0.08, X_1 = 0.34, X_2 = 0.1$. Parameters and initial data are chosen merely for illustrating the two possible outcomes; they do not have biological significance.

3.   Conclusions

Theorem $2.3$ represents a generalization of Theorem 1 in ^[11] in the following sense: In ^[11] we considered the case where $\gamma_1 = \gamma_2$ and $F_1(P) = F_2(P)$, which means that the cooperator and cheater are indistinguishable as far as their per capita uptake rates $F_i(P)/\gamma_i$, and their yield constants $\gamma_i$ are concerned. In this special case $P_2<P_1$ always holds, and therefore the tragedy occurs (2nd case in Theorem $2.3$), and thus we recover the result of Theorem 1 in ^[11] as a special case.

The main purpose of this paper was to investigate whether a tragedy can be avoided if the cooperator evolves by either changing its yield constant $\gamma_1$, and/or its growth function $F_1(P)$, relative to the cheater. Interestingly, only changing its yield $\gamma_1$, but keeping $F_1(P) = F_2(P)$, is insufficient: In this case, the inequality $P_2<P_1$, remains valid and a tragedy cannot be avoided. This is supported by experimental observations with different metabolic strategies. In direct competitions, low-yield but high-flux fermentation out-competes high-yield but low-flux respiration ^[9].

Therefore, the cooperator must modify $F_1(P)$ to have a chance at survival. The key to cooperator survival lies in the reversal of the inequality $P_2<P_1$ into $P_1<P_2$. In other words, the cooperator must evolve in such a way that it has a lower break-even concentration than the cheater, thereby becoming better adapted to the environment of the chemostat. We can gain some insight into how the cooperator can achieve this by considering the case that $F_1(P)$ is a Monod function. Then formula $(2.14)$ (and recalling that $p_1 = P_1$) shows that $P_1$ can be decreased by either increasing $m_1$ or by decreasing $a_1$; that is, by either increasing its maximal uptake rate through an increase in $m_1$, or by decreasing its half-saturation constant $a_1$, the cooperator can evolve to become the superior competitor that drives the cheater to extinction.

This has been observed during experimental evolution of cooperating microbial populations, where adaptations that improve nutrient uptake stabilize cooperative behavior ^[15,1]. More broadly, any principle or mechanism that disproportionally increases the benefit of cooperation in cooperating individuals, including kin selection or spatial structuring, will favor its evolution and maintenance ^[17].

Assuming that $P_1<P_2$ holds, then depending on the initial condition of the system, Theorem $2.3$ shows that essentially there are two possible scenarios: Either the tragedy continues to hold, or, more strikingly, the cooperator outcompetes the cheater and drives the cheater to extinction. Under no circumstances can there be a coexistence of cooperators and cheaters. There is some experimental evidence that time and/or space heterogeneities promote cooperation ^[2,7,3]. In our model time heterogeneities can be introduced by letting the dilution rate and input nutrient concentration fluctuate explicitly with respect to time, by replacing $D$ by $D(t)$, and $S^0$ by $S^0(t)$ in our model. However, as we have shown in ^[11], such fluctuations still give rise to a Tragedy of the Commons, when $F_1(P) = F_2(P)$ and $\gamma_1 = \gamma_2$, confirming that time heterogeneities cannot induce coexistence of cooperator and cheater. We have not investigated scenarios with spatial heterogeneities in the environment, but in principle one could extend our model to reflect the gradostat (when space is discrete), or the unstirred chemostat (when space is continuous), see ^[14] for more on models of this kind.

Acknowledgments

Patrick De Leenheer is supported in part by NSF-DMS-1411853. Martin Schuster is supported by NSF-MCB-1158553. Hal Smith is supported by Simons Foundation Grant 355819.

Conflict of interest

All authors declare no conflicts of interest in this paper.