Department of Mathematical Sciences, P. O. Box 8093, Georgia Southern University, Statesboro, Georgia 30460-8093, USA
Received date: May 17, 2013; Accepted date: June 11, 2013; Published date: June 17, 2013
Citation: Braselton JP, Abell ML, Braselton LM (2013) Rock-Paper-Scissors in the Chemostat. J Comput Sci Syst Biol 6:118-131. doi:10.4172/jcsb.1000109
Copyright: © 2013 Braselton JP, et al. This is an open-access article distributed under the terms of the Creative Commons Attribution License,which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.
Visit for more related articles at Journal of Computer Science & Systems Biology
Rock-Paper-Scissors is a game played by two players to deter-mine a single winner. Biological relationships of Rock-Paper-Scissors are documented. In this paper, we form a continuous model of Rock-Papers-Scissors in the chemostat that coincides with the biology of such relationships. The basic models that we develop coincide with the observed phenomena. Be-cause the model involves a system of seven nonlinear differential equations, global results are difficult to obtain. We present several numerical studies that are the result of a substantial number of numerical trials to illustrate the various possibilities that might occur in the context of the problem discussed here.
Chemostat; Competition; Rock-paper-scissors; Inhibitor; Quorum sensing; Colicin; Antibiotic
Traditionally, rock-paper-scissors is a children’s game. We describe how to play the game with your hand. Rock is represented by a fist, scissors by two fingers extended to represent a pair of scissors, and paper by a flattened hand. At the count of three, the two players extend their hands into one of the forms described. Rock breaks scissors so defeats scissors, scissors cut paper so scissors defeat paper, and paper covers rock so paper defeats rock. Viewing the relationship between the rock, paper, and scissors, each has an advantage over the other. Although rock- paper-scissors may be a relatively easy way to settle simple disputes between two friends, "rock-paper-scissors" occurs in interesting and meaningful biological contexts as well.
Kerr et al.  describes a relationship between three populations of Escherichia coli that produce toxins against each other resulting in a biological rock-paper- scissors relationship. Their mathematical simulations (in silico study) using a chemostat model showed that coexistence of all three populations could occur if the chemostat was not well stirred so that interactions between the competitors were local. On the other hand, if the chemostat was well stirred so that interactions between the competitors were global, one would win and eliminate both of the other competitors. They carefully examined their study in the lab (using in vitro techniques) and generally saw that one population dominated the mixture after one week or so. In their article, Roelke and Eldrige , expanded on the results of .
Another interesting study (in silico study) is by Karolyi et al. , who examine a similar problem and demonstrate that the competition is affected by the distribution of the competitors. They further show that by chaotic mixing of the nutrient (with density S(t) in which that follows) that there are cyclic oscillations in the population densities of the three competitors.
Of particular biological interest Kirkup and Riley , perform an experiment (in vivo) using mice whose intestinal bacteria are in a rockpaper- scissors relationship. In their experiment, they examine colicin (antibiotic) producing strains of Escherichia coli against its close relatives in a rock-paper-scissors relationship. From the results of their experiment, they conclude that a diverse family of antibiotics [colicins] serves to promote microbial diversity in one of the dominant niches exploited by enteric bacteria, the mammalian colon.
In the models studied by Karolyi et al. , and Kerr et al. , the “rock-paper-scissors” relationship is formed by a killer (toxin producing or more specifically, colicin producing) organism, toxin sensitive organism, and toxin resistant organism. The colicin producing strains kill the sensitive strains that outcompete the resistant strains which outcompete the colicin producing strains. Unlike these previously studied models, in the model we develop and study in this paper we give each competitor a specific method of attack by giving it the ability to produce a toxin against its enemy. We expand on these previous results by considering the case when each organism produces a toxin (or colicin) to which one of its competitors is sensitive and the other resistant, which is more similar to the experiment performed by Kirkup and Riley .
By viewing the rock-paper-scissors relationship as a competition problem, the natural question to ask is who wins?
Competing species problems have been extensively studied theoretically and experimentally. We do not list an exhaustive reference of such studies but rather mention a few that seem particularly interesting in the context of the model we develop here. In a general sense, competing species in the chemostat are studied in Smith and Waltman, . Many variations of the basic model model have been studied. For example Li , studies a competition model with three competing species competing for three nutrients. His study is based on the experimental results seen by Husiman and Weissing , which show that cyclic competition, such as in a rock-paper-scissors relationship, can result in periodic oscillations (or limit cycles). More recently, Cameron et al. , examine a rock-paper-scissors relationship involving a parasitic plant and experimentally show that coexistence of the three competitors occurs and confirm the conclusion with biological observations. More recently, in Hsu and Roeger , use the May-Leonard competition models to study more complex competition relationships in the chemostat.
The production of anti-competitor toxins is of interest when the weaker competitor can devote some of its resources to the production of a toxin (or inhibitor) against its competitors at some expense to its own growth. In Husiman and Weissing , this was modeled as a constant proportion of the resources. Biologically, this makes sense. For example, it is documented that mixing strains of Escherichia coli, some of which produce bacteriocins, which are poisons that have a negative a effect on competing organisms, the end result can lead to co-existence of all competitors [10-13]. Inhibitors (including those added to the environment as well as those produced by the organisms) in the chemostat have been studied [14-20].
In this paper we consider a basic, resource-based model of competition of which the chemostat is the standard example. Such models have applications in ecology to model a simple lake, a simple digestive system, and in biotechnology to model the commercial bioreactor. Experimental verification of the match between theory and experiment can be found in Hanson and Hubbel  study. For a general discussion of competition models in the chemostat see Frederickson and Stephanopoulos , or Smith .
To formulate a continuous model of a biological rock-paper-scissors relationship in the chemostat, we start with the basic models that are summarized by Smith and Waltman , and then the competition models and, specially, the model when one competitor produces a toxin studied by Hsu and Waltman . We restate the basic competition equations (scaled) for two competitors in the chemostat as stated in Smith and Waltman  (note that ‘=d/dt):
The competition models (scaled) when one competitor produces a toxin is studied in by Hsu and Waltman :
In (1) and (2), m1 represents the maximal growth rates, a1 are the Michaelis- Menten constants, and k represents the fraction of nutrient availability to the x2 organism allocated to the production of the toxin. We assume that that interaction between the toxin and the aected micro organism is of mass-action form, -γPx1 (is the proportionality constant). In both equations, S(t) denotes the concentration of the nutrient in the chemostat. In the first equation, x1(t) and x2(t) denote the concentration of the competitors. In the second equation, x1(t) denotes the density of the toxin sensitive organism and x2(t) denotes the density of the toxin producing organism. P(t) denotes the density of the concentration of the toxin in the chemostat.
To generalize the Hsu and Waltman study in  to model a “Rock-Paper-Scissors” relationship, the number of equations increases but then we also must indicate how each organism affects the other. In the model we develop, we assume that we are dealing with organisms with densities x1(t), x2(t), and x3(t): the x1- organism produces a toxin against the x2-organism (P1(t)), which is represented by k1 as a fraction of its nutrient consumption; the x2-organism produces a toxin against the x3-organism (P2(t)), which is represented by k2 as a fraction of its nutrient consumption; and the x3-organism produces a toxin against the x1-organism (P3(t)), which is represented by k3 as a fraction of its nutrient consumption for a true rock-paper-scissors relationship. For simplicity, we rest assume that ki (but k1, k2, and k3 may have different values) is constant, but then generalize the approach followed by Braselton and Waltman , when k is not constant.
As stated above, let S(t) denote the concentration of the nutrient at time t, x1(t) the concentration of the microorganism susceptible to the toxin secreted by the organism with concentration x3(t), x2(t) the concentration of the microorganism susceptible to the toxin secreted by the organism with concentration x1(t), x3(t) the concentration of the microorganism susceptible to the toxin secreted by the organism with concentration x2(t). The concentrations of the toxin producing organisms are given by P1(t) (for x1), P2(t) (for x2) and P3(t) (for x3). The underlying assumption is that the chemostat is well stirred so the nutrient is equally available to all competitors. The model takes the form
In (3), S(0) is the input concentration of the nutrient, D is the washout rate, mi the maximal growth rates, ai are the Michaelis-Menten constants, and ηiare the yield constants. S(0) and D are usually controlled by the experimenter although in real life situations such as studying how Escherichia coli interact in a (human or mouse, ) digestive system might not be controllable by the individual involved. The equations (3) are usually called the Monod Model or the model with Michaelis- Menten dynamics. For the purposes here, the constant kirepresents the fraction of potential growth devoted to producing the toxin. If one had chosen to not start with a chemostat model, an alternative choice might have been starting with the May-Leonard competition equations as in Hsu and Roeger . Also, note that our assumption that the chemostat is well stirred differs significantly from the assumptions of Kerr et al.  and Karolyi et al.  described previously. The interaction between the allelopathic agent and the sensitive microorganism have been taken to be of mass action form, γiPjxk. This is common in modeling when an interaction depends on the concentrations.
For chemostat problems, we scale to dimensionless variables. To do so, we first assume η=η1=η2=η3 and then we let
After scaling, dropping the bars, and replacing τ with t system (3) simplifies to
In this reduced form, the system is dimensionless and we have fewer variables to work with so is preferred to the original system.
Observe that if S(0), x1(0), x2(0), x3(0), P1(0), P2(0), and P3(0) are non negative, then the solutions are also non negative.
To see this, let Σ=S + x1 + x2 + x3 + P1 + P2 + P3. Then,
so lim supt→oe Because each component is non-negative, it follows that system (4) is dissipative.
To linearize at the equilibrium points (rest points), we let and then compute the Jacobian for system (4):
The boundary rest point
E0=(S0; x10; x20; x30; P10; P20; P30)=(1; 0; 0; 0; 0; 0; 0)
always exists. Evaluated at the Jacobian we have which has eigenvalues along the diagonal: -1 (with multiplicity 4) and (1-k1)f1(1)-1, (1 - k2) f2(1) - 1, and (1 - k3)f3(1) - 1. E0 will be unstable when at least one of the last three eigenvalues is positive:
that we require in which follows.
Boundary Rest Points
Other boundary rest points exist when x1=0, x2=0, or x3=0. Observe that if system (4) has boundary rest points, the resulting system is equivalent to that studied by Hsu and Waltman . To illustrate why this is true, we assume that one of the population densities is 0. For convenience we choose x1=0. Then, system (4) becomes
For system (7) to be biological meaningful, we must have that k3=0, which reduces the system to
which is equivalent to the system (2) studied by Hsu and Waltman . The same argument applies if x2=0 or x3=0.
Interior Rest Points
When interior rest points exist of system (4) exist, they take the form EA=(SA; x1A; x2A; x3A; P1A; P2A; P3A) where
where SA is a positive solution of
The coefficients c1 are given by
Note that for interior rest points (9) to exist, we must have
If this occurs when this occurs when
To investigate the stability of any interior rest points, we chose parameter values that resulted in an interior rest point. The following example is typical of what we observed.
results in the equilibrium solutions shown in Table 1. Observe that these parameter values result in one interior rest point, EA.
Table 1: Equilibrium solutions using the parameter values .
Next, we investigate the stability of each rest point by computing the eigenvalues of the Jacobian, J, at each rest point. In Table 2, we list the maximum of the real part of the eigenvalues of the Jacobian, J, evaluated at each rest point. Observe that the interior rest point, EA, is unstable. Two of the boundary rest points are locally stable.
|Label||Maximum of Real Part of Eigenvalues of Jacobian||Classi cation|
Table 2: The maximum value of the real part of the eigenvalues of the Jacobian evaluated at each rest point listed in Table 1.
To see that EA exists so that coexistence is possible, we illustrate the coexistence in Figure 1. However, it is important to remember that EA is unstable. Changing the initial conditions slightly leads to very di_erent results. (We choose initial conditions near EA.) We illustrate a few of the possible situations in Figure 2. In the plots shown in the figures observe that depending on the initial conditions; generally, the results indicate that one or two species will dominate after long periods of time.
To investigate the stability of EA, we adjust the parameter values that we started with initially in (12). First, we vary k1 and leave the remaining parameter values the same as in (12). When an interior rest point exists, we compute the eigenvalues of the Jacobian evaluated at the rest point, then compute the maximum value of the real part of any eigenvalue, and plot the result. See the black plot in Figure 1 of the Jacobian evaluated at each rest point listed in Table 1. We then repeated the procedure for k2 and k3 in Figure 3 (a), a1, a2, and a3 in Figure 3 (b), 1, 2, and 3 in Figure 3 (c), and m1, m2, and m3 in Figure 3 (d). In all cases, observe that the maximum value of the real part of any of the eigenvalues is always positive: the interior rest point is unstable.
In fact, using constant ki mutual exclusion always occurred except at the exact densities of the unstable interior rest point as shown in Figure 1. This makes sense and agrees with other models: the coexistence of colicin producer and a sensitive strain in a well-mixed culture is not possible , except that in our model all competitors are producing colicins.
The model (4) does not appear to exhibit the dynamics and coexistence observed by Hsu and Roeger , or Cameron et al. , over short periods of time. We were not able to produce cyclical coexistence as observed by Karolyi et al. , probably because of our assumption that the chemostat is well stirred, while they allowed the nutrient density to fluctuate randomly and our underlying assumptions regard ing the relationships between the competing organisms are very different from the assumptions made by Kerr et al.  and Karolyi et al. . On the other hand, the model (4) indicates that over long periods of time, biological systems do stabilize unless external influences cause them to become unstable. Intuitively this makes sense for several reasons. For example, one might argue that an unstable biological relationship would propagate and then influence other biological systems, leading to a very unpredictable biological world. Nevertheless, in the last two papers mentioned, interaction between the organisms was considered or observed. In particular, this would imply that the organisms have some ability to sense those around them. Through mechanisms known as quorum sensing bacteria are able to control the expression of their genes in response to density of other bacteria in their environment. Quorum sensing mechanisms have been demonstrated to play a role in the control of gene expression associated with diverse activities like bioluminescence, in massing to form biofilms, and the expression of the gene code for characters responsible for the pathogenicity of these organisms. For example, Sandoz et al. , describe a situation in which mutations of the bacteria Pseudomonas aeruginosa use quorum sensing to cheat and coexist with non-mutated strains of the bacteria. In a review on quorum sensing by Bassler , provides a large annotated bibliography and reading list on this fascinating subject. In her review she states that recent studies show that quorum sensing modulates both intra- and inter-species cell-cell communication. In discussing a particular case she notes that the capacity to respond to both intra- and inter-species signals could allow V. harveyi to know not only its own cell density, but also the relative frequency of bacteria of its type in mixed populations.
Thus, if a bacterium has the ability to sense the current state of its habitat and the presence of other bacteria, we believe it is reasonable to conclude that ki is not constant but rather ki=ki(x1; x2; x3). The problem is magnified in what follows so that little rigorous analytical results can be obtained for a general ki(x1; x2; x3) and so it is necessary to consider special cases. We note first that there are two undesirable cases. If ki=0, no agent is produced and since y is assumed to be the weaker competitor, it becomes extinct. Similarly, if ki=1, all uptake is devoted to toxin production and none to growth, so again the organism cannot survive. ki(x1; x2; x3) must fall between these extremes of 0 (no toxin produced) and 1 (the growth rate of the organism is then 0 and the organism faces extinction). We consider two special cases that we feel may be extremes for the toxin producing functions ki(x1; x2; x3). First, we consider
For example, looking at the k1(x1; x2) equation in (13), if x1 is large, x1 devotes more of its resources to producing the toxin. On the other hand, if x2 is small or nonexistent, x1 devotes less effort to toxin production. This strategy has the advantage that if there is no competitor, no energy is wasted on toxin production. The interpretation follows the same pattern in the k2 and k3 equations. The fraction of resources devoted to toxin production depends on the organism producing the toxin and the organism affected by the toxin.
With ki=ki(x1; x2; x3) given by (13), system (4) becomes
As with system (4), if S(0), x1(0), x2(0), x3(0), P1(0), P2(0), and P3(0) are non-negative, the solutions are also non-negative because (5) is the same as in the previous model. The choice of non constant k results in significantly more complicated algebra and calculus computations. Therefore, we do not explicitly state formulas for the rest points for system (17) because they are symbolically complicated. Similarly, the Jacobian, J, is also symbolically complicated so is not shown here for length considerations.
For a basic model of quorum sensing and using equations (13), we choose parameter values that are nearly identical to those used in the first example to observe differences between the two models.
Results in the equilibrium solutions shown in Table 3. We then compute the real part of the eigenvalues of the Jacobian at each rest point and list the maximum value in Table 4. Observe that these parameter values result in one unstable interior rest point, EA.
Table 3: Equilibrium solutions using the parameter values .
|Label||Maximum of Real Part of Eigenvalues of Jacobian||Classi cation|
Table 4: The maximum value of the real part of the eigenvalues of the Jacobian evaluated at each rest point listed in Table 3.
To generate Figure 4, we used initial conditions near EA as we did in the previous example. Of particular interest, note that even though we are using a non constant k, the results are quite similar as to the results obtained using constant k. In particular, notice that mutual exclusion tends to occur. Also, as with constant k the outcome depends on the parameter values and initial conditions.
In equations (13), the organism’s toxin production depend on the density of the organism and the density of the organism affected by its toxin (inhibitor). However, in the rock-paper-scissors relationship, studied here a given organism is affected by both the organism affected by its toxin as well as by the organism that produces a toxin against the organism. For our second choices of ki, we consider a natural strategy to explore is my enemy’s enemy is my friend, which is not considered in equations (13) in the sense that in those choices of ki depend only on the density of the organism and its ability to detect the organism that it produces its toxin against. However, a given organism might want to detect the densities of both competitors. To take this into consideration in the competition model, the toxin producing organism, say x1, might produce its toxin against the x2 organism at a rate inversely proportional to the density of the x3 organism because the x3 organism produces a toxin against the x1 organism. Thus, when the x1 organism decreases its toxin production against the x2 organism, the x2 organism is allowed to devote more of its resources to producing its toxin against the x3 organism, which then benefits the x1 organism by lowering the x3 population and less toxin is produced against the x1 organism. On the other hand, when the x2 population (density) is high but the x3 population (density) is low, x1 devotes more of its resources to producing its toxin against the x2 organism. We then continue a similar process for the other organisms. We try to choose simple choices of ki(x1; x2; x3) that take this into consideration. In our discussion, we choose
With ki given by (16), system (4) becomes
As with systems (4) and (14), if S(0), x1(0), x2(0), x3(0), P1(0), P2(0), and P3(0) are non-negative, the solutions are also non-negative because (5) is the same as in the previous models. Because each component is non-negative, it follows that system (17) is dissipative. The choice of non constant k results in significantly more complicated algebra and calculus computations. Therefore, we do not explicitly state formulas for the rest points for system (17) because they are symbolically complicated. Similarly, the Jacobian, J, is also symbolically complicated so is not shown here for length considerations.
For our second model of quorum sensing and using equations (16), we choose parameter values that are nearly identical to those used in the first two examples to observe the similarities and differences between the three models.
Results in the equilibrium solutions shown in Table 5. Observe that these parameter values result in two interior rest points, EA and EB. As in the previous cases, we then evaluate the Jacobian at each rest point, compute its eigenvalues, compute the real part of each eigenvalue and look at the maximum to determine the stability of each rest point in Table 6. Of particular interest, notice that the interior rest point EB is locally stable.
Table 5: Equilibrium solutions using the parameter values .
|Label||Maximum of Real Part of Eigenvalues of Jacobian||Classi cation|
Table 6: The maximum value of the real part of the eigenvalues of the Jacobian evaluated at each rest point listed in Table 5.
In Figure 5, we choose the initial conditions first to be at EA to illustrate the unstable coexistence. Changing the initial conditions slightly results in competitive exclusion, as indicated in Figure 5 (a).
Of more interest is when we graph solutions with initial conditions near EB in Figure 6. EB is locally stable so many initial conditions result in coexistence as shown Figures 6 (a) and (c), which reaffirms the conclusions reached by Kirkup and Riley . On the other hand, it is possible to change the conditions so that competitive exclusion occurs as shown in Figures 6 (b) and (d).
Our simulations did not result in cyclic oscillations of the density of competitors observed by Karolyi et al. . This may be because of our assumption that the chemostat is well stirred, while they allowed the nutrient density to fluctuate randomly in their study. Another explanation is that their rock-paper-scissors relationship is different from ours. In their study, they observed cyclic competition between toxin producing, sensitive, and resistant organisms. On the other hand, our model specifically gives each organism the ability to produce a toxin (or colicin in the context of bacteria) against one of its competitors that is resistant to a different competitor. That is, the x1 organism produces a toxin against the x2 organism (rock crushes scissors), the x2 organism produces a toxin against the x3 organism (scissors cuts paper), and the x3 organism produces a toxin against the x1 organism (paper covers rock). Further, model (14) indicates that over long periods of time, biological systems do stabilize unless external influences causes them to become unstable. We also mathematically see that in a biological rock-paper-scissors relationship coexistence can occur as has been observed experimentally in studies such as the one done by Kirkup and Riley, .
Comparing the similarities and differences between the three models, we believe they indicate why coexistence occurs in a biological rockpaper- scissors relationship. A given organism’s toxin production (for example, x1) incorporates quorum sensing and depends on the population (density) of the organisms sensitive to the given organisms toxin (x2) as well as the density of the organism producing a toxin against the given organism (x3).
Of course, one can conceive of many other strategies that could be investigated. For example, equations (13) and (16) are specific cases of a more general strategy that takes the form.
Incorporating (19) into system (4) and examining various combinations and interpretations of the parameter values might result in an interesting study. We used Mathematica to perform our simulations and the notebooks are available from the authors. Our approach is quite general so other choices of the coefficient functions could be investigated with relative small changes in the Mathematica notebooks used in the calculations and plots here and are available from the authors as described next. Finally, we note that the positive constant in the denominator keeps the function differentiable at the origin. A smaller constant mimics stricter ratio dependence. In both cases k(0; 0; 0)=0: We also chose _ so that 0 < ki < 1.
For length considerations, we have not listed all parameter values used in our calculations, all initial conditions, or the eigenvalues of all the equilibrium solutions in our stability calculations. Generally, the initial conditions we used were near the interior rest points, EA or EB. The Mathematical notebooks that the authors used to carry out all of the calculations as well as to generate the figures are available from the authors by sending a request to Jim Braselton at [email protected]
In this paper we have numerically analyzed a seven-dimensional nonlinear system of differential equations that models a biological rock-paper-scissors relationship. In previous models with fewer competitors or different assumptions, a reduction process was able to be carried out to reduce the dimension of the system. In the cases studied here, reducing the dimension of the system (seven equations) is not possible because of the interpretation of the agent, k=ki, that is allowed to produce a toxin (inhibitor) against its competitors.
We have considered several extreme cases of the cost of the metabolic load of producing a toxin against a competitor. We examined the case when k=ki is constant and then two cases when k=ki was not constant. To see the coexistence that is biologically documented, we used a non constant k, which indicates that when coexistence occurs in such situations, quorum sensing might be involved in the biological relationship between the organisms.
Several of the studies mentioned here have incorporated multiple competitors competing for multiple nutrients. One could incorporate those considerations into our examples by including multiple Si (nutrient availability) for the competing species.
Other interesting steps might be to investigate the results seen by Kerr et al. , where they observed multiple mutations after several generations. Although incorporating genetics into the equations developed here might only result in numerical results, such results could have meaningful biological significance. Another interesting study would be to adjust the May-Leonard competition equations as in Hsu and Roeger , repeat the analysis, and analyze any unexpected results.
We think that studying the situation when there are multiple nutrients that the competitors are competing for would be particularly interesting, especially if the availability of the nutrients was dependent on seasons or possibly even random.
The late Professor Paul Waltman motivated the first author to study this model in 2002. Therefore, the first author dedicates this paper to him.