Rock-Paper-Scissors in the Chemostat

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.


Introduction
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. [1] 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 [2], expanded on the results of [1].
Another interesting study (in silico study) is by Karolyi et al. [3], 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 [4], 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. [3], and Kerr et al. [1], 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 [4].
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, [5]. Many variations of the basic model model have been studied. For example Li [6], 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 [7], 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. [8], 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 [9], 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

Formation of the Continuous Model
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 [5], and then the competition models and, specially, the model when one competitor produces a toxin studied by Hsu and Waltman [16]. We restate the basic competition equations (scaled) for two competitors in the chemostat as stated in Smith and Waltman [5] (note that '=d/dt): The competition models (scaled) when one competitor produces a toxin is studied in by Hsu and Waltman [16]: In (1) and (2), m 1 represents the maximal growth rates, a 1 are the Michaelis-Menten constants, and k represents the fraction of nutrient availability to the x 2 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, -γPx 1 (is the proportionality constant). In both equations, S(t) denotes the concentration of the nutrient in the chemostat. In the first equation, x 1 (t) and x 2 (t) denote the concentration of the competitors. In the second equation, x 1 (t) denotes the density of the toxin sensitive organism and x 2 (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 [16] 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 x 1 (t), x 2 (t), and x 3 (t): the x 1 -organism produces a toxin against the x 2 -organism (P 1 (t)), which is represented by k 1 as a fraction of its nutrient consumption; the x 2 -organism produces a toxin against the x 3 -organism (P 2 (t)), which is represented by k 2 as a fraction of its nutrient consumption; and the x 3 -organism produces a toxin against the x 1 -organism (P 3 (t)), which is represented by k 3 as a fraction of its nutrient consumption for a true rock-paper-scissors relationship. For simplicity, we rest assume that k i (but k 1 . k 2 , and k 3 may have different values) is constant, but then generalize the approach followed by Braselton and Waltman [21], when k is not constant.
As stated above, let S(t) denote the concentration of the nutrient at time t, x 1 (t) the concentration of the microorganism susceptible to the toxin secreted by the organism with concentration x 3 (t), x 2 (t) the concentration of the microorganism susceptible to the toxin secreted by the organism with concentration x 1 (t), x 3 (t) the concentration of the microorganism susceptible to the toxin secreted by the organism with concentration x 2 (t). The concentrations of the toxin producing organisms are given by P 1 (t) (for x 1 ), P 2 (t) (for x 2 ) and P 3 (t) (for x 3 ). The underlying assumption is that the chemostat is well stirred so the nutrient is equally available to all competitors. The model takes the form (1 ) 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 ηi are 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, [4]) 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 k i represents 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 [9]. Also, note that our assumption that the chemostat is well stirred differs significantly from the assumptions of Kerr et al. [1] and Karolyi et al. [3] described previously. The interaction between the allelopathic agent and the sensitive microorganism have been taken to be of mass action form, γ i P j x k . 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.
To linearize at the equilibrium points (rest points), we let ( ) and then compute the Jacobian for system (4): Where, The boundary rest point that we require in which follows.

Boundary Rest Points
Other boundary rest points exist when x 1 =0, x 2 =0, or x 3 =0. Observe that if system (4) has boundary rest points, the resulting system is equivalent to that studied by Hsu and Waltman [16]. To illustrate why this is true, we assume that one of the population densities is 0. For convenience we choose x 1 =0. Then, system (4) becomes For system (7) to be biological meaningful, we must have that k 3 =0, which reduces the system to which is equivalent to the system (2) studied by Hsu and Waltman [16]. The same argument applies if x 2 =0 or x 3 =0.

Interior Rest Points
When interior rest points exist of system (4) exist, they take the form EA=(SA; where SA is a positive solution of The coefficients c i are given by ) )))), ))), Note that for interior rest points (9) to exist, we must have

An Example
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.  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, E A , is unstable. Two of the boundary rest points are locally stable.

Choosing
To see that E A exists so that coexistence is possible, we illustrate the coexistence in Figure 1. However, it is important to remember that E A is unstable. Changing the initial conditions slightly leads to very di_erent results. (We choose initial conditions near E A .) 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 E A , we adjust the parameter values that we started with initially in (12). First, we vary k 1 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 k 2 and k 3 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 k i 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 [4], except that in our model all competitors are producing colicins.    The Case for Non-Constant k The model (4) does not appear to exhibit the dynamics and coexistence observed by Hsu and Roeger [9], or Cameron et al. [8], over short periods of time. We were not able to produce cyclical coexistence as observed by Karolyi et al. [3], 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. [1] and Karolyi et al. [3]. 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. [22], 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 [23], 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 k i is not constant but rather k i= k i (x 1 ; x 2 ; x 3 ). The problem is magnified in what follows so that little rigorous analytical results can be obtained for a general k i (x 1 ; x 2 ; x 3 ) and so it is necessary to consider special cases. We note first that there are two undesirable cases. If k i =0, no agent is produced and since y is assumed to be the weaker competitor, it becomes extinct. Similarly, if k i= 1, all uptake is devoted to toxin production and none to growth, so again the organism cannot survive. k i (x 1 ; x 2 ; x 3 ) 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 k i (x 1 ; x 2 ; x 3 ). First, we consider  For example, looking at the k 1 (x 1 ; x 2 ) equation in (13), if x 2 is large, x 1 devotes more of its resources to producing the toxin. On the other hand, if x 2 is small or nonexistent, x 1 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 k 2 and k 3 equations. The fraction of resources devoted to toxin production depends on the organism producing the toxin and the organism affected by the toxin.
With k i =k i (x 1 ; x 2 ; x 3 ) given by (13), system (4) As with system (4), if S(0), x 1 (0), x 2 (0), x 3 (0), P 1 (0), P 2 (0), and P 3 (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.

An Example
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.
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.

My Enemy's Enemy Is My Friend
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 k i , 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 k i 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 x 1 , might produce its toxin against the x 2 organism at a rate inversely proportional to the density of the x 3 organism because the x 3 organism produces a toxin against the x 1 organism. Thus, when the x 1 organism decreases its toxin production against the x 2 organism, the x 2 organism is allowed to devote more of its resources to producing its toxin against the x3 organism, which then benefits the x 1 organism by lowering the x 3 population and less toxin is produced against the x 1 organism. On the other hand, when the x 2 population (density) is high but the x 3 population (density) is low, x 1 devotes more of its resources    Table 3. to producing its toxin against the x 2 organism. We then continue a similar process for the other organisms. We try to choose simple choices of k i (x 1 ; x 2 ; x 3 ) that take this into consideration. In our discussion, we choose x With k i given by (16), system (4) becomes As with systems (4) and (14), if S(0), x 1 (0), x 2 (0), x 3 (0), P 1 (0), P 2 (0), and P 3 (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.

An Example
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.
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 [4]. 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. [3]. 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   The maximum value of the real part of the eigenvalues of the Jacobian evaluated at each rest point listed in Table 5. 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 x 1 organism produces a toxin against the x 2 organism (rock crushes scissors), the x 2 organism produces a toxin against the x 3 organism (scissors cuts paper), and the x 3 organism produces a toxin against the x 1 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, [4].
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, x 1 ) incorporates quorum sensing and depends on the population (density) of the organisms sensitive to the given organisms toxin (x 2 ) as well as the density of the organism producing a toxin against the given organism (x 3 ).
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 < k i < 1.

Computational notes
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 jbraselton@georgiasouthern.edu 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=k i , 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=k i is constant and then two cases when k=k i 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, [1], 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 [9], 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.