Why is cyclic dominance so rare?
Abstract
Natural populations can contain multiple types of coexisting individuals. How does natural selection maintain such diversity within and across populations? A popular theoretical basis for the maintenance of diversity is cyclic dominance, illustrated by the rockpaperscissor game. However, it appears difficult to find cyclic dominance in nature. Why is this the case? Focusing on continuously produced novel mutations, we theoretically addressed the rareness of cyclic dominance. We developed a model of an evolving population and studied the formation of cyclic dominance. Our results showed that the chance for cyclic dominance to emerge is lower when the newly introduced type is similar to existing types compared to the introduction of an unrelated type. This suggests that cyclic dominance is more likely to evolve through the assembly of unrelated types whereas it rarely evolves within a community of similar types.
Introduction
Natural populations ranging from microbial communities to animal societies consist of many different individuals. Some individuals compete with each other to exploit a shared resource (Hardin, 1960; Connell, 1983), whereas others coexist (Morris et al., 2013). Interactions affect the death or reproduction of individuals and thus shape the composition of populations (Friedman and Gore, 2017). Different types of individuals are distinguishable at the interaction level and they have a complex interaction structure (Farahpour et al., 2018). Because interaction structures themselves can support the coexistence of multiple types, they have been extensively studied in ecology and evolution (Gross et al., 2009; Allesina and Levine, 2011). A particularly exciting type of interaction is cyclic dominance, in which each type dominates another one but is in turn dominated by a separate type, leading to a RockPaperScissors cycle (Maynard Smith, 1982; Hofbauer et al., 1998; Szab and Fth, 2007) as sketched in Figure 1A. None of the types fixates in the population, because each type is dominated by one type while it simultaneously dominates a third type. Thus, it has been argued that this type of interaction can support biodiversity (Reichenbach et al., 2007). Cyclic dominance has therefore attracted substantial attention and it has been extensively studied theoretically (Maynard Smith, 1982; Hofbauer et al., 1998; Frean and Abraham, 2001; Hauert et al., 2002; Reichenbach et al., 2007; Szab and Fth, 2007; Mathiesen et al., 2011; Jiang et al., 2011; Allesina and Levine, 2011; Szolnoki et al., 2014).
A famous example of this type of cyclic dominance in biology is toxin production in Escherichia coli (Kerr et al., 2002; Kirkup and Riley, 2004; Cascales et al., 2007). Toxinproducing (or colicinogenic) E. coli cells can purge cells that are sensitive to the toxin. However such toxin producers are dominated by resistant cells that do not produce the toxin. Once common, resistant cells are again dominated by sensitive cells, which avoid the costs of resistance. This leads to cyclic dominance in a RockPaperScissors manner, as shown in Figure 1B. Another example is the mating strategies of North American sideblotched lizards Uta stansburiana (Sinervo and Lively, 1996). The strategy of males guarding several females dominates the strategy of males guarding only a single female. However, sneaky strategies under which males secretly mate with guarded females can become dominant over the strategy of males guarding several females. Once such a sneaky strategy is common, the strategy of guarding only a single female can be successful again leading to cyclic dominance among the mating types as illustrated in Figure 1C. In addition to E. coli and sideblotched lizards, other examples of cyclic dominance have been described in ecology: Stylopoma spongites (Jackson and Buss, 1975), Drosophila melanogaster (Clark et al., 2000), European lizards Lacerta vivipara (Sinervo et al., 2007), and plant systems (Taylor, 1990; Lankau and Strauss, 2007; Cameron et al., 2009). These types of cyclic dominance arise because of competition, which can happen within and between species at the same trophical level. Mating or sperm competitions are the basis of cyclic dominance within species observed in U. stansburiana and D. melanogaster, whereas common resource competition can happen both between and within species.
However, traditional theoretical work assumes a set of predefined cyclic dominance types without asking how they developed or came together. In ecosystems, the introduction of a new species through migration can lead to such cyclic dominance. In this context, it is often termed intransitive competition (Allesina and Levine, 2011; Soliveres et al., 2015; Gallien et al., 2017). However, immigrating species can also disturb and destroy cyclic dominance. In evolving populations, new types can arise through mutation and recombination. In the same manner, mutation and recombination can lead to the formation of cyclic dominance but can also lead to types that do not fit into such types of dominance and break the cycle. A recent experimental study (Higgins et al., 2017) indicated that in the assembly of microbial ecosystems based on 20 bacterial strains found in a single grain of soil, only 3 of almost 1000 triplets exhibited cyclic dominance. Other soil bacterial species (Wright and Vetsigian, 2016; Friedman et al., 2017) also displayed a lack of cyclic dominance. This rareness is found in both soil bacteria and plant systems (Taylor, 1990). Why is it so difficult for cyclic dominance to emerge by assembly or evolution? In this study, we ask the following question theoretically: How frequent is cyclic dominance in situations in which new types constantly arise, providing an opportunity for new cycles but also breaking old cycles at the same time?
Forming a cyclic dominance from a single type is challenging because as soon as the second type arises, the dominance type will take over the whole population, driving extinction of the other type. Therefore, a third type must arise before the population loses either of the two previous types. Such a precise timing of the arrival of a new type is critical for developing cyclic dominance and it can occur when new types arise at a high frequency, either through high mutation rates, recombination, or immigration (Kotil and Vetsigian, 2018; Tarnita, 2018). This rapid evolution can be achieved through both high mutation rates per capita and large population sizes (Lanfear et al., 2014; Hague and Routman, 2016; Vahdati et al., 2017; Kotil and Vetsigian, 2018; Tarnita, 2018). Thus, we considered a model in which the population naturally evolves to a large population size (Park et al., 2019), which allows the development of cyclic dominances via an evolutionary process.
Once we introduced our model in more detail, we will show that the interaction of ecology and evolution leads to increasing population size. This increases the chances that cyclic dominance arises, but it remains rare. Next, we rationalize this finding: While the lifespans of cyclic and noncyclic dominance triplets are similar, it is more difficult to form cyclic dominance compared to noncyclic dominance. The underlying reason is that similarity between parental and offspring payoffs suppresses the formation of cyclic dominance. Finally, we discuss which genealogical structure can promote or suppress cyclic dominance.
Materials and methods
Interactions between individuals affect their death or birth. A traditional model for describing an interacting population is the generalized LotkaVolterra equation (Ginzburg et al., 1988; Bomze, 1995; Yoshida, 2003). In particular, some studies (Biancalani et al., 2015; Huang et al., 2015; Shtilerman et al., 2015; Barbier et al., 2018; Farahpour et al., 2018) assumed that the interaction determines the likelihood of death from a pairwise competition. These interaction parameters can be written as a form of a matrix, including selfinteraction. However, only a few studies considered novel mutations (Huang et al., 2012; Shtilerman et al., 2015; Farahpour et al., 2018). Drawing new interaction parameters for a new type and extending the interaction matrix, we considered such a novel mutation process. In addition, the size of the interaction matrix can be reduced when types go extinct. We traced an evolving population by dealing with this dynamically changing matrix.
We built the model based on individual reaction rules
where I and J are individuals of types i and j, respectively. We assumed that all types are in the same trophic level, and thus there is only competition and no predation. All types have the same background birth and death rates. Only competition makes a difference (Huang et al., 2015; Farahpour et al., 2018; Park et al., 2019). Because the population always collapses when $\lambda}_{b}\le {\lambda}_{d$, we only focused on $\lambda}_{b}>{\lambda}_{d$. For the sake of simplicity, we only considered wellmixed populations without any other highorder interactions.
Formulating the competition death rate $d}_{ij$ as a function of the payoff $A}_{ij$, we connected evolutionary game theory to the competitive LotkaVolterra type dynamics (Huang et al., 2012; Huang et al., 2015; Park and Traulsen, 2017; Park et al., 2019; Sidhom and Galla, 2020). Note that $A}_{ij$ is the payoff of an individual of type i from the interaction with an individual of type j. Because lower payoffs should increase the probability of death, we used an exponentially decaying function for the competition death rate as follows:
where $\alpha >0$ is the baseline death rate from competition, which ensures that the population remains bounded regardless of the value of the evolving payoffs $A}_{ij$. A larger payoff implies a lower death rate from competition. The overall competition death rate is always positive, such that we remain in the regime of the competitive LotkaVolterra equations.
For a large population size, the abundance $x}_{i$ of type i can be described using the competitive LotkaVolterra equation
where n is the number of types in the population, used as a diversity index, and T is the time. The stability of the population is determined by Equation (3). In parallel, the stability between only two types can be determined by the two associated equations in Equation (3), which are described by four payoff values of the two types.
Once a new mutant type arises during reproduction, new interactions occur. To describe these new interactions, we draw new payoff values from the parental payoff with Gaussian noise
where $\xi$ is a random variable sampled from a Gaussian distribution with zero mean and variance $\sigma}^{2$. This inheritance of payoffs with noise implies that the mutant type $i}^{\mathrm{\prime}$ slightly deviates from the parental type i. Here, we treat selfinteraction $A}_{ii$ and interaction with other types $A}_{ij(\ne i)$ in the same way. Because of new interactions, the population composition changes over time, as shown in Figure 2A. We let the population evolve from a single type with a randomly drawn initial payoff from the normal distribution with mean $\mathrm{ln}(1000)$ and standard deviation 1. As a natural consequence of evolving payoffs, the average population size also evolves. Because different types are fully described by the payoff matrix, we can trace the evolving population by tracking the payoff matrix, as shown in Figure 2B. We do not consider any tradeoff: having higher payoffs does not cost anything. Hence, the evolution tends to increase the payoffs constantly and thus drives the system into a regime where payoff differences become smaller.
To construct a pairwise interaction network, we used the stability between two types. Hence the term interaction refers to the pairwise relationship, considering the stability between two types. There are four possible scenarios for stability (see Appendix 1):
Dominance of type i:
represented by $\u25ef\phantom{\rule{negativethinmathspace}{0ex}}\phantom{\rule{negativethinmathspace}{0ex}}\phantom{\rule{negativethinmathspace}{0ex}}\phantom{\rule{negativethinmathspace}{0ex}}i\leftarrow \leftarrow \phantom{\rule{negativethinmathspace}{0ex}}\u25ef\phantom{\rule{negativethinmathspace}{0ex}}\phantom{\rule{negativethinmathspace}{0ex}}\phantom{\rule{negativethinmathspace}{0ex}}\phantom{\rule{negativethinmathspace}{0ex}}j\text{}$ for $A}_{ii}\phantom{\rule{negativethinmathspace}{0ex}}>\phantom{\rule{negativethinmathspace}{0ex}}{A}_{ji$ and $A}_{ij}\phantom{\rule{negativethinmathspace}{0ex}}>\phantom{\rule{negativethinmathspace}{0ex}}{A}_{jj$.
Dominance of type j:
represented by $\u25ef\phantom{\rule{negativethinmathspace}{0ex}}\phantom{\rule{negativethinmathspace}{0ex}}\phantom{\rule{negativethinmathspace}{0ex}}\phantom{\rule{negativethinmathspace}{0ex}}i\to \to \phantom{\rule{negativethinmathspace}{0ex}}\u25ef\phantom{\rule{negativethinmathspace}{0ex}}\phantom{\rule{negativethinmathspace}{0ex}}\phantom{\rule{negativethinmathspace}{0ex}}\phantom{\rule{negativethinmathspace}{0ex}}j\text{}$ for $A}_{ii}\phantom{\rule{negativethinmathspace}{0ex}}<\phantom{\rule{negativethinmathspace}{0ex}}{A}_{ji$ and $A}_{ij}\phantom{\rule{negativethinmathspace}{0ex}}<\phantom{\rule{negativethinmathspace}{0ex}}{A}_{jj$.
Bistability:
represented by $\u25ef\phantom{\rule{negativethinmathspace}{0ex}}\phantom{\rule{negativethinmathspace}{0ex}}\phantom{\rule{negativethinmathspace}{0ex}}\phantom{\rule{negativethinmathspace}{0ex}}i\leftarrow \to \phantom{\rule{negativethinmathspace}{0ex}}\u25ef\phantom{\rule{negativethinmathspace}{0ex}}\phantom{\rule{negativethinmathspace}{0ex}}\phantom{\rule{negativethinmathspace}{0ex}}\phantom{\rule{negativethinmathspace}{0ex}}j\text{}$ for $A}_{ii}\phantom{\rule{negativethinmathspace}{0ex}}>\phantom{\rule{negativethinmathspace}{0ex}}{A}_{ji$ and $A}_{ij}\phantom{\rule{negativethinmathspace}{0ex}}<\phantom{\rule{negativethinmathspace}{0ex}}{A}_{jj$.
Coexistence:
represented by $\u25ef\phantom{\rule{negativethinmathspace}{0ex}}\phantom{\rule{negativethinmathspace}{0ex}}\phantom{\rule{negativethinmathspace}{0ex}}\phantom{\rule{negativethinmathspace}{0ex}}i\to \leftarrow \phantom{\rule{negativethinmathspace}{0ex}}\u25ef\phantom{\rule{negativethinmathspace}{0ex}}\phantom{\rule{negativethinmathspace}{0ex}}\phantom{\rule{negativethinmathspace}{0ex}}\phantom{\rule{negativethinmathspace}{0ex}}j\text{}$ for $A}_{ii}\phantom{\rule{negativethinmathspace}{0ex}}<\phantom{\rule{negativethinmathspace}{0ex}}{A}_{ji$ and $A}_{ij}\phantom{\rule{negativethinmathspace}{0ex}}>\phantom{\rule{negativethinmathspace}{0ex}}{A}_{jj$.
Constructing the interaction network, we can examine the formation and the collapse of cyclic dominance, as shown in Figure 2C. If the links are drawn from a random matrix, each stability scenario described above occurs with the same probability. Thus, we find a proportion of 0.50 dominance links and a proportion of 0.25 proportions for bistability and coexistence links, respectively. Because the networks can contain three different link types (dominance, bistability, and coexistence), both cyclic dominance and other types of triplets can be found. However, in the main text, we only focused on cyclic and noncyclic dominance triplets which are composed of only dominance because the proportions of each triplet strongly depend on the proportions of link types. Hence at a given link composition (three dominance links), we investigate how often we can observe cyclic dominance compared to noncyclic one.
Results
Evolution leads to increasing population size
The population dynamics described in Equation (1) appears simple, but its tracing is complicated because of the novel mutations. Due to the emergence of a new mutant and its consequences, the payoff matrix dynamically changes. As large payoffs lower competition death rates, types with higher payoffs are more likely to survive. Therefore, payoffs evolve to larger values, which increases the population size (Park et al., 2019). Since there is no tradeoff on the payoffs, the average payoff increases monotonically. This makes types become more similar, enhancing diversity (Scheffer and van Nes, 2006). However, the population size saturates at a certain level because of the baseline death rate $\alpha$ corresponding to resource limitation and enters a steady state (see Appendix 2).
For small $\alpha$ values (rich environments) in particular, the population size N at the steady state becomes large, containing many different types (see Figure 3A). This evolution toward a large population induces rapid mutation. Once the populations size becomes large, new mutant types are generated faster than in smaller populations given a fixed mutation rate per individual. In this rapid mutation regime, a new mutant can arise before the population equilibrates, thereby establishing a cyclic dominance from a timely emerged mutant (Kotil and Vetsigian, 2018). Thus, cyclic dominance can be established when the populations enters the rapid mutation regime, as shown in Figure 3B. The proportions of triplets were averaged over all surviving realizations. In principle, we can observe a triplet from $t=2$ even though there were fewer than three types on average. Because of the smaller average diversity $\u27e8n\u27e9$, there were large fluctuations in measuring the proportions of triplets in the early regime ($t\lesssim 100$). However, the measurement became more accurate as diversity increased.
Cyclic dominance triplets are rare
The proportions of cyclic and noncyclic dominance triplets increase in the early dynamics and quickly saturate. Whereas population dynamics illustrates the formation of both cyclic and noncyclic dominance, cyclic dominance is much rarer than noncyclic dominance. To quantify this rareness of cyclic dominance, we measured the fraction $\chi$, which is defined by the fraction of cyclic dominance triplets among all dominance triplets (cyclic and noncyclic) (see Figure 3C). In steady state, the fraction yielded $\chi \approx 0.033$, indicating that one cyclic dominance triplet can be found among 30 dominancecomposed triplets. If the pairwise relationships are random (all four possible links appear in the same probability 1/4, called a random network), then the fraction $\chi$ of cyclic dominance should be 0.25 because there are only two configurations of cyclic dominance triplets, whereas six configurations produce noncyclic dominance triplets. Hence cyclic dominance that developed from our population dynamics is much rarer than expected from a random choice of interactions. To elucidate why it is the case, we checked how long each triplet can be sustained in the population and how often they emerge. First, we focused on the lifetime of cyclic and noncyclic dominance and then moved on the formation of each triplet. When it comes to the chance to emerge cyclic and noncyclic dominance, we argued that the genealogy structure shaped by ecoevolutionary dynamics will enhance or suppress the formation of cyclic dominance.
The lifespans of cyclic and noncyclic dominance triplets are similar
The rareness of cyclic dominance triplets may be caused by their shorter lifespan compared with that of noncyclic dominance triplets. Thus, we investigated the lifespan of triplets first to understand the rareness of cyclic dominance. Once triplets arise in populations, we can identify them, and trace how long they persist. Lifespan distributions in the steady state of both cyclic and noncyclic dominance triplets decayed algebraically. We plotted the complementary cumulative distribution functions (CCDFs), clearly revealing a power law decay, as shown in Figure 4A. Surprisingly, there was no difference in the lifespan of both triplets. Both cyclic and noncyclic dominance triplets were destroyed in five mutation events on average. The noncyclic dominance triplet has a higher chance of persisting longer, although the difference is small. In addition, the median is the same for both distributions because almost all probabilities are concentrated on short lifespans. In conclusion, lifespan does not explain why cyclic dominance is rarer than noncyclic dominance. Hence, the lower chance for cyclic dominance to emerge is the reason.
The condition to form cyclic dominance is more strict than that for noncyclic dominance
Why is it more difficult for cyclic dominance to emerge than for noncyclic dominance? One factor is that the conditions needed for an interaction matrix to provide cyclic dominance are more restrictive than those for noncyclic dominance. For a matrix to reveal cyclic dominance, it is necessary that in each of the three columns, the three payoffs $A}_{ij$, $A}_{jj$, and $A}_{kj$ are ordered ($A}_{ij}<{A}_{jj}<{A}_{kj$ or $A}_{ij}>{A}_{jj}>{A}_{kj$). Conversely, the formation of noncyclic dominance requires this condition to be satisfied only in a single column, whereas the other two columns should satisfy a less restrictive condition ($A}_{ij}<{A}_{jj$ and $A}_{kj}<{A}_{jj$ or vice versa). For example, in random payoff matrices where all payoffs are randomly drawn from the standard normal distribution, the fraction of cyclic dominance is $1/13\approx 0.077$ (as shown in Appendix 4), which is smaller than the value of $2/8=0.25$ expected in a random network of directed links. This is because in the matrix approach, unlike a random network, links in a triplet are interdependent; selfinteraction payoffs contribute to the character of several links at once.
Similarity between parental and offspring payoffs suppresses the formation of cyclic dominance
The fraction $\chi \approx 0.077$ in the random matrix is still larger than that obtained from our population dynamics $\chi \approx 0.033$, implying there are other factors suppressing the development of cyclic dominance. A key reason is the correlation between payoffs. In our model, the elements of the payoff matrix are not fully independent because of heredity. Offsprings payoffs are derived from their parents payoffs. We found that this correlation between payoffs plays an important role in suppressing the formation of cyclic dominance. For example, let us imagine two different uncorrelated preexisting types represented by a 2 × 2 random matrix. All elements are drawn from the standard normal distribution. If a third type emerges from a preexisting one by a mutation, the average fraction of cyclic dominance becomes $\chi \approx 0.02$ (the standard normal distribution is used for new payoffs and we average across 5 10^{5} samples), which is lower than that of the assembly of thee uncorrelated types. This is because the payoffs of a new offspring are similar to its parental payoffs: the offspring is likely to have the same relationships with other types in a population as the parental type and any type dominated by the parent will also be dominated by the offspring. The triplets including these offspring and parental types are more likely to form noncyclic dominance than a cyclic one. Hence the correlation between payoffs affects the fractions by which cyclic dominance emerges compared with noncyclic dominance.
To check the effect of the correlation on emerging triplets, we measured the similarity between types as a proxy of the correlation between payoffs. We defined the trait vector $\overrightarrow{T}}_{l$ of the type l using the row capturing with its own payoff and the column of the others payoff against it in the payoff matrix, ${\overrightarrow{T}}_{l}=({A}_{l,i},{A}_{l,j},{A}_{l,k},{A}_{i,l},{A}_{j,l},{A}_{k,l})$, similarly to Farahpour et al., 2018. Because the average payoff increases over time, we shifted all elements in those vectors by a constant to ensure that the average of all values is zero. Then, using normalized vectors after shifting we calculated the scalar product for all pairs of trait vectors as a similarity measure. Larger values indicate that the two types have more similar payoff values. Each triplet has three trait vectors and thus has three similarity measures. Taking the mean and standard deviation of these three similarities, we found that in cyclic dominance the three types tend to be less similar compared to noncyclic dominance, as shown in Figure 4B. The inheritance plays a key role in the emergence of cyclic and noncyclic triplets, giving rise to a correlation between payoffs. Because the payoff correlation is determined by a genealogy structure, in the following we investigated which of these structures have higher or lower chances to promote the emergence of cyclic dominance.
Genealogical structure can promote or suppress cyclic dominance
Between the last common ancestor and the present types, there are intermediary types accumulating mutations between them. Genealogies tell us who is whose parent, tracing back to the common ancestor of the observed types. From the genealogy, we can infer how many mutations were accumulated by each type and the time at which they diverged. If two types have only accumulated a few mutations from the most recent common ancestor, their payoffs are likely to be similar. Hence, the genealogy structure shapes the correlation between payoffs and affects the value of fraction $\chi$. To study the role of genealogies, we first characterized genealogies for three types based on accumulated mutational distances, see Figure 5AB. We found that for a triplet of types, the distribution of payoff elements could be characterized by five parameters of the genealogy, which all measure the number of mutations between types w_{1}, w_{2}, x, y, and z (see Appendix 5).
We analyzed those mutational distances for cyclic and noncyclic triplets found in the simulation at the steadystate. From the analysis, we inferred that the crucial parameter influencing the fraction $\chi$ of cyclic dominance is the fraction $F}_{l$ of mutations accumulated after each type evolves independently from others (the time $t}_{d$ in Figure 5B).
Payoff correlations must be developed before lineages become independent at $t}_{d$, depending on the mutational history (w_{1} and x). With more mutations being accumulated after the types divergence, the correlation between types becomes weaker, increasing the chance of the cyclic dominance to emerge, see Figure 5C. On the other hand, with only few mutations after $t}_{d$ the strong correlation between parental and offsprings payoffs remains, which mostly induces the emergence of noncyclic dominance. This means that the mutations accumulated after all three types diverge are important for decoupling their payoff correlations, increasing the chance to form the cyclic dominance. This finding agrees well with numerically investigated genealogies that gave the minimal and maximal values of the fraction $\chi$, called minimizer and maximizer genealogies, respectively. For minimizer genealogies, the majority of mutations were accumulated before $t}_{d$ (see Appendix 6). Conversely, for maximizer genealogies, the majority of mutations were accumulated after all three lineages diverged from each other.
In addition, numerically calculated minimizer genealogies almost completely suppressed the emergence of cyclic dominance ${\chi}_{\text{min}}<0.001$. For the maximizer, the fraction could be as high as ${\chi}_{\text{max}}=1/6\approx 0.167$. Both the fraction of cyclic dominance arising from the random matrix $\chi \approx 0.077$ and that found in our simulations $\chi \approx 0.033$ fell between the minimal and the maximal values possible with the genealogy structure. The fractions $\chi$ from population dynamics are closer to the minimal value, as shown in Figure 6. Therefore, we can infer that in the genealogies occurred in population dynamics, the similarity of the new types to their parental type prevents the emergence of cyclic dominance.
Discussion
Cyclic dominance is extremely interesting from a conceptual and theoretical perspective and it has thus been analyzed in great detail in mathematical biology (Hofbauer et al., 1998; Hofbauer and Sigmund, 2003; Szab and Fth, 2007). However, the theoretical literature typically refers to only a handful of examples in nature. Moreover, recent experiments have revealed that it is difficult for cyclic dominance to emerge in microbial populations (Wright and Vetsigian, 2016; Friedman et al., 2017; Higgins et al., 2017). Why is the establishment of cyclic dominance so difficult? To address this question, we used an evolutionary process with evolving interactions for the formation of such cyclic dominance instead of following the more conventional approach of using a predefined set of interactions. For example, Kotil and Vetsigian, 2018 observed the formation of cyclic dominance in the fast evolution regime with adaptation, but the involved traits were predefined. However, we argue that it is difficult for cyclic dominance to emerge even in the presence of rapid evolution. In addition, our investigation shows no correlation between diversity and the probability to find cyclic dominance. The diversity observed in our model does not originate from the cyclic dominance structure, but from the generation on a time scale that is fast compared to equilibration of the system. Our results indicate cyclic dominance in general could be a mechanism to support diversity, but it is most probably not essential in a situation where it can emerge, as such situations are characterized by a high diversity in the first place. Even if we rescale the payoffs to prevent going towards neutrality which increases diversity, cyclic dominance remains rare as long as the cyclic dominance is not a main driver for the diversity. This result is robust even when we consider the population frequency weights, see Appendix 7. In our model, the noncyclic dominance can reach a majority of the population, while the cyclic dominance typically remains bounded to a much lower frequency.
We also examined the circumstances under which cyclic dominance can appear more frequently. While the probability of assembling such an interaction structure by chance in a payoff matrix with uncorrelated random entries is small, the probability to evolve such an interaction structure is even smaller. The inheritance of interactions from parent to offspring is a key mechanism shaping the correlations between payoffs and determines the formation of cyclic dominance. From our approach, we found that the introduction of an uncorrelated type is crucial for the formation of cyclic dominance triplets. Because the migration of new species can be interpreted as such an introduction, our results suggest that cyclic dominance might be more frequent on an interspecies basis than on an intraspecies basis. As widespread intransitive competition is found in ecological systems (Soliveres et al., 2015; Gallien et al., 2017), our manuscript nicely supports the basic idea that assembly of unrelated types is more likely to lead to cyclic triplets than evolution, in which emerging types are closely related.
Our approach, which reduces the complexity from continuous values to a categorical classification, may help to bridge the model dynamics and experimental data more easily. Experimental work has provided data regarding both the constituents of a microbial community but also the interactions between them. However for large communities, parameterizing all interactions in the model numerically makes it difficult to identify the fundamental factors shaping the dynamics. Reducing the complexity may permit study of the large scales of experimental data connecting the underlying model dynamics and large datasets.
An important limitation of our work is the assumption of global interactions. In our model, all individuals can interact with each other, ignoring the spatial population structure. A spatial model could localize the interactions and lead to the more frequent formation of cyclic dominance. Such a localization can foster cyclic dominance for a predefined cyclic set (Durrett and Levin, 1997; Durrett and Levin, 1998; Frean and Abraham, 2001; Reichenbach et al., 2007; Szab and Fth, 2007; Mathiesen et al., 2011; Jiang et al., 2011; Mitarai et al., 2012; Szolnoki et al., 2014; Kelsic et al., 2015; Sneppen, 2017). However, before moving into spatial models it appears essential to investigate this issue in the absence of all potentially confounding factors. Such models appear necessary for explaining why cyclic dominance within one species is not found often in nature, and they may open a new direction for the extensive theoretical work on this topic.
Appendix 1
Pairwise relationship based on stability
For large population sizes, the change of abundances of types can be described by deterministic equations. Linear stability analysis reveals which types will persist in equilibrium (Strogatz, 2000). Using this stability, we determine the relationship between two types i and j. In this section, we perform linear stability analysis first and then determine the pairwise relationships based on that stability.
To determine the stability between two types, type 1 and type 2, we focus on the two equations
where $\lambda ={\lambda}_{b}{\lambda}_{d}$. Since the death rates d_{ij} are all positive, the population would always go extinct for $\lambda <0$. Thus, we only consider a positive $\lambda$. There are four fixed points $({x}_{1}^{\ast},{x}_{2}^{\ast})$: One is the extinction of both types $(0,0)$, two are indicating each singletype population $(\lambda /{d}_{11},0)$ and $(0,\lambda /{d}_{22})$, and the fourth one is the coexistence point $\left(\frac{({d}_{22}{d}_{12})\lambda}{{d}_{11}{d}_{22}{d}_{12}{d}_{21}},\frac{({d}_{11}{d}_{21})\lambda}{{d}_{11}{d}_{22}{d}_{12}{d}_{21}}\right)$. To check the stability of each fixed point, we calculate the Jacobian matrix $\mathbf{J}$,
and then calculate the eigenvalues for the fixed point.
The extinction point (0,0) is always unstable, because both eigenvalues are identical to $\lambda >0$. On the other hand, the three other fixed points have one negative eigenvalue $\lambda$, indicating that they are saddle or stable fixed points; If the second eigenvalue is positive, the fixed point becomes a saddle and is unstable. If the second eigenvalue is negative, the fixed point is stable. Hence the fixed point $(\lambda /{d}_{11},0)$ becomes stable when $d}_{21}>{d}_{11$, because the eigenvalue is $({d}_{11}{d}_{21})\lambda /{d}_{11}$. In the same manner, the condition $d}_{12}>{d}_{22$ can be obtained for a stable fixed point $(\lambda /{d}_{11},0)$. Lastly, the coexistence point has an eigenvalue $\frac{({d}_{11}{d}_{21})({d}_{22}{d}_{12})}{{d}_{11}{d}_{22}{d}_{12}{d}_{21}}\lambda$. For meanginful values of abundances, stable coexistence points have to satisfy the conditions ${x}_{1}^{\ast}>0$ and ${x}_{2}^{\ast}>0$. Together with these restrictions, we can find the coexistence fixed point is stable only if both $d}_{11}>{d}_{21$ and $d}_{22}>{d}_{12$ are satisfied.
We classify the pairwise relationship based on these stabilities. Dominance relationships are given if only a singletype fixed point is stable while all other fixed points are unstable. When both singletype fixed points are stable, a bistability relationship is drawn. A coexistence relationship is achieved when only the coexistence 1xed point is stable. We summarized the stabilities of 1xed points at a given condition and named the pairwise relationship in Appendix 1—table 1.
Appendix 2
Population size and diversity at steadystate for various baseline death rates
In this section, we show how the population size N and diversity n evolve in time t with various baseline death rates $\alpha$. Since larger payoffs imply lower death rates, individuals with larger payoffs tend to survive. As a result of evolution, the overall death rates decrease, which enlarges the population size. This can be interpreted that species improve their efficiency to consume a resource. However, if resources are limited, the population size is confined. We implement this resource limit by introducing the baseline death rate $\alpha$. The death rates from the competition cannot be lower than the baseline death rate, and thus the population size saturates at a certain level at the end, $\lambda /\alpha$. For various $\alpha$, we run 5000 independent stochastic simulations and measure the average population size and diversity in time t, see Appendix 2—figure 1A and B. The average runs over the surviving samples at a given time $t=10000$, and we denote the averaged quantity O as $\u27e8O\u27e9$. As we expected, the population sizes evolve to $\lambda /\alpha$, and thus larger $\alpha$ leads to smaller population sizes.
Large population size reduces the time until new types occur in the population, and at the same time it takes longer to equilibrate. Thus, a large population size increases the number of different types in the population via two processes: (i) Even when the most competitive type is expected to fixate in the population, new types which are less fit can constantly emerge due to mutation. (ii) A second mechanism is ecological tunneling (Kotil and Vetsigian, 2018): Even though the population dynamics yields the extinction of certain types, it can be rescued by the emergence of a new type which supports the coexistence of types predicted to be extinct. From those effects, larger populations have higher diversity Appendix 2—figure 1B. We compare the average population size and diversity in the stationary regime as well, see Appendix 2—figure 1C.
We also take a closer look at the average payoff values. Because of selection, the average payoff increases in time, but the increment becomes smaller as the average increases, see Appendix 2—figure 1D. After a transient time, the average payoffs increases logarithmically. It means that the natural selection induced by the payoff difference becomes weaker and is almost neural at the end. Hence, in the long run, the diversity mainly originates from the mutationselection balance.
Appendix 3
Link proportions
From the payoff matrix, we can determine the pairwise relationship between types, constructing a network. We represent the relationship between two types i and j with three different link types, dominance, bistability, and coexistence:
As a basic element of a network, we investigate the frequencies of link types first.
As a reference, we consider a random payoff matrix model wherein all payoffs are randomly drawn from the standard normal distribution. In this case, the probability that one payoff is larger than another is 0.5 because all payoffs are independently sampled from the same distribution. Thus each condition in Appendix 1—table 1 happens with the same probability, 0.25. Thus, the probability to observe a dominance link is 0.5 (since there are two directions), and the others are 0.25. We use these values as a reference to compare the link proportions obtained from population dynamics.
With population dynamics, the behavior of link frequencies differs for different $\alpha$. For small $\alpha$, the ensembleaveraged proportions are stable in time, while with large $\alpha$ values the average proportions fluctuate a lot, see Appendix 3——figure 1A and B. For large $\alpha$, the number of types rarely exceeds one. Hence, the link proportions fluctuate a lot, but the average values are well predicted from a random payoff matrix because the links appear by chance. On the other hand, for small $\alpha$, many links can be formed due to the presence of multiple types. In this case, the average link proportions are stable in the stationary regime, and coexistence links are favored compared to the random matrix case. We also plot the mean of ensembleaveraged link proportions in the stationary regimes $9500\le t\le 10000$ for various $\alpha$ in Appendix 3——figure 1A. As we can see, the link proportions get closer to that of the random matrix as increasing the baseline death rates $\alpha$.
Appendix 4
Triplet proportions
Let us now consider triplets that are constructed by three links. With four kinds of link relationships, in total we can find 64 possible triplets. However, if we take into account symmetries, this reduces to 16 different triplet structure, see Appendix 4—figure 1. We call the number of triplets, which have the same triplet structure degeneracy. For example, with type indices i, j, and k, the cyclic dominance triplet structure can be found in two different triplets: $i\to \to j\to \to k\to \to i$ and $i\leftarrow \leftarrow j\leftarrow \leftarrow k\leftarrow \leftarrow i$. In contrast, there is no degeneracy for the triplet types marked by 1 and 4 in Appendix 4—figure 1A because changing the indices of types does not give any difference. Different triplet structures have different degeneracy because they have different mirror and rotation symmetries, summarized in Appendix 4—figure 1. Since our interest is the triplet structure, we investigate proportions of 16 triplet structures considering those degeneracies.
Again we use a random payoff matrix to attain the reference triplet proportions first. For three types i, j, and k, we can construct a 3 3 random matrix,
All elements are independently drawn from the standard normal distribution $\mathcal{\mathcal{N}}(x0,1)$. Here, $\mathcal{\mathcal{N}}(xm,{\sigma}^{2})$ indicates a Gaussian distribution with mean m and variance $\sigma}^{2$. Any triplet is fully defined by three relationships, and each pairwise game determines one relationship. Mathematically, this corresponds to three $2\times 2$ submatrices of $\mathbf{A}$. However, these matrices are not independent, as it would require $3\times 2\times 2=12$ independent parameters while the payoff matrix has only $3\times 3=9$ entries.
To determine the type of triplet structure, we focus on one type of individual j. For this type j, we can characterize invasion possibilities against the two other types, giving four possible situations,
$i\phantom{\rule{1em}{0ex}}\to j\to \phantom{\rule{1em}{0ex}}k$, when $A}_{ij}<{A}_{jj}<{A}_{kj$
$i\phantom{\rule{1em}{0ex}}\leftarrow j\leftarrow \phantom{\rule{1em}{0ex}}k$, when $A}_{ij}>{A}_{jj}>{A}_{kj$
$i\phantom{\rule{1em}{0ex}}\to j\leftarrow \phantom{\rule{1em}{0ex}}k$, when $A}_{jj}>{A}_{ij$ and $A}_{jj}>{A}_{kj$
$i\phantom{\rule{1em}{0ex}}\leftarrow j\to \phantom{\rule{1em}{0ex}}k$, when $A}_{jj}<{A}_{ij$ and $A}_{jj}<{A}_{kj$
where $i\phantom{\rule{1em}{0ex}}\to j$ means that j is stable with respect to invasion of i from rare (either j dominates i, $i\to \to j$ or they are bistable $i\leftarrow \to j$, hence the notation). Mathematically, this corresponds to three $3\times 1$ columns of $\mathbf{A}$. Since all columns are independent, the probability of finding a certain triplet is the product of three probabilities to find a certain set of stubs for each type.
Thus we calculate the probability to find a certain set of subs for a type. The probability $2\times 2$ that three random variables a, b, and c satisfy the condition $a<b<c$ is given by convolution,
where $\mathrm{\Phi}(x)$ is the probit function, $\mathrm{\Phi}(x)={\int}_{\mathrm{\infty}}^{x}\mathcal{\mathcal{N}}(\theta 0,1)d\theta$. Hence, the probability that a type has in and outstubs is 1/6. On the other hand, the probability $P(c<a\text{and}cb)$ is
In the same way, $P(c>a\phantom{\rule{thinmathspace}{0ex}}\mathrm{a}\mathrm{n}\mathrm{d}\text{}\mathrm{c}\mathrm{b})$ is 1/3. Hence, the probabilities to observe each set of stubs are
$i\phantom{\rule{1em}{0ex}}\to j\to \phantom{\rule{1em}{0ex}}k$ appears with probability 1/6
$i\phantom{\rule{1em}{0ex}}\leftarrow j\leftarrow \phantom{\rule{1em}{0ex}}k$ appears with probability 1/6
$i\phantom{\rule{1em}{0ex}}\to j\leftarrow \phantom{\rule{1em}{0ex}}k$ appears with probability 1/3
$i\phantom{\rule{1em}{0ex}}\leftarrow j\to \phantom{\rule{1em}{0ex}}k$ appears with probability 1/3
With those probabilities and the degeneracies, we can calculate the triplet proportions expected in a random payoff matrix. The results are summarized in Appendix 4—table 1, where we find the cyclic dominance (Triad ID 15) is the rarest one. At the same time, the noncyclic dominance (Triad ID 16) is one of the most abundant triplets. From the proportions we can also obtain the fraction of cyclic dominance among cyclic and noncyclic dominance triplets as $\chi =1/13$, indicating 12 noncyclic dominances can occur while only a single cyclic dominance appears.
Now we move to triplet proportions in population dynamics, see Appendix 4—figure 1B. Since the triplet proportions strongly depend on link proportions, the results for $\alpha =5\cdot {10}^{6}$ and $\alpha =5\cdot {10}^{4}$ are very different. For example, coexistence links are norm for $\alpha =5\cdot {10}^{6}$ and thus the triplet with Triad ID 4 is more abundant than in other cases. Thus, it is hard to directly compare proportions of triplets, which have different link composition. Instead of that, we compare the triplets, which have the same link composition. This comparison allows us to find which structure is more abundant, eliminating the effect of link proportions. There are four such sets of triplets, see Appendix 4—figure 1B. One of the sets consists of cyclic and noncyclic dominance triplets, which are composed of three dominance links. If we look at the fraction $\chi$ of cyclic dominance, we can find that noncyclic dominance is more suppressed in population dynamics than in the random matrix. For the other three sets of triplets, we can also find that the types which have bistability usually dominate another type while the types with coexistence are dominated by others. However, this tendency becomes weaker in population dynamics compared to the random payoff matrix analysis.
Appendix 5
Properties of payoff matrices for three types and a given genealogy
The emergence of a new mutant type l induces a new row and a new column in the payoff matrix. If we denote the parental type of type l as $r(l)$, the new payoffs are written by
where $\xi ({\sigma}^{2})$ are random values sampled from the Gaussian distribution with zero mean and variation $\sigma}^{2$, $\mathcal{\mathcal{N}}(x0,{\sigma}^{2})$. Due to this inheritance, the genealogical structure which tells us who is whose parent type shapes the payoffs. In this section, we show how payoff elements are determined by a given genealogy.
Genealogies can be schematically depicted with time and mutational distance axes, see Appendix 5—figure 1A. For three focal types i, j, and k, a number of possible genealogies exists. We begin with a scenario shown in Appendix 5—figure 1B, where the lineage leading to the type i diverges from the last common ancestor o first and a different lineage leading to the two types j and k diverges from o later. We index the last common ancestor of k and j in their lineage as n. Also, we give an index m for the progenitor of type i at the moment when lineages leading to types j and k diverged.
We begin with a payoff matrix of types i, j, and k and trace it back to the moment of the last common ancestor type o of all three types. All nine payoffs between types i, j, and k are derived from the single payoff value $A}_{oo$.
For a pair of types j and k, all four payoffs $A}_{jj$, $A}_{jk$, $A}_{kj$, and $A}_{kk$ can be traced back to the payoff $A}_{nn$. Let us imagine that y and z mutations happen leading the emergence of type j and k from the type n respectively. If y mutation events occur first before z mutation events, we can write the value of $A}_{jk$ as
Since the random variable $\xi$ follows the normal distribution, we can simply write
Thus, multiple mutation events can be written as a single mutation event with larger variance. For any other order of the mutation events, the final expression is the same even though all intermediate terms will be different.
In the same way, we can proceed with the other three payoffs. Consequently, all four payoffs can be written as
Selfinteractions, $A}_{jj$ and $A}_{kk$, do not change when the other type accumulates a mutation, so their mutational distances from $A}_{nn$ are smaller.
Next, we consider type i which is diverged from the type o but is a different lineage from j and k. Then, the payoff $A}_{ii$ can be traced back to $A}_{mm$ as
where w_{2} is the number of mutations accumulated in the lineage of type i from type m. The payoffs between i and two other types can be traced back to payoffs between their progenitors $A}_{mn$ and $A}_{nm$
Hence, all nine payoffs describing interactions between types i, j, and k can be traced back to four payoffs, $A}_{mm$, $A}_{mn$, $A}_{nm$, and $A}_{nn$, characterizing interactions between m and n. These four payoffs, in turn can be traced back to $A}_{oo$ in the same way,
where x is the number of mutations accumulated from the common ancestor o to the type n, and w_{1} is the number of mutations accumulated from the type o to the type m. In summary, we show how payoffs can be calculated from their common ancestors payoffs, see Appendix 5—figure 2.
The same calculations apply to any other genealogy, for example ones presented in Appendix 5—figure 1CE, if we label the types in a specific way. The labels j and k should indicate the pair of types, which diverged the last, that is, the last common ancestor of j and k is the most recent one among all three pairwise the last common ancestors. By exclusion, the rest one becomes type i. Type n is the last common ancestor of j and k (for genealogy on Appendix 5—figure 1E and F, n is the same as o). Finally, type m is the progenitor of i existed at the moment of divergence of j and k (for genealogy on Appendix 5—figure 1D, m is the same as o).
Appendix 6
Minimizer and maximizer genealogies
In this section, we numerically identify those genealogies which produce the largest or smallest proportions of certain triplets. First of all, we note that the fractions of triplets do not depend on scales of the set of mutational distances $D}_{m}=\{{w}_{1},{w}_{2},x,y,z\$ without population dynamics; rescaling the distances results in the same triplet proportions. Therefore, without loss of generality, we can assume that all five mutational distances sum up to one and are thus located on the simplex ${w}_{1}+{w}_{2}+x+y+z=1$.
We search for extreme genealogies by numerical optimization, implemented by a hill climbing algorithm. The optimization starts from a random set of D_{m} and identifies the local optimum. Since our method can find only local optima, we use multiple initial values to find optimal set. After we identified the mutational distance sets which give local optima, we cluster them using principle component analysis (PCA). We denote the Xth class of genealogies that maximize and minimize the proportion of triplet tri by ${D}^{tri}\text{}X$ and ${D}_{tri}\text{}X$, respectively. For different scenarios, we find different numbers of distinct classes of genealogies.
Maximization of the number of cyclic dominances
First, we ask which genealogies maximize the number of cyclic dominance triplets. Performing PCA analysis, we found five different classes of such genealogies, see Appendix 6—table 1. We denote cyclic dominance triplets as T_{15} and noncyclic dominance ones as T_{16}. Classes ${D}^{15}\text{}1$ and ${D}^{15}\text{}2$ are effectively the same genealogy, where all three types diverge from each other very early and then one of these types accumulates all the mutations (the most isolated type i in ${D}^{15}\text{}1$ or one of the more recently diverged types k in ${D}^{15}\text{}2$, see the end of Appendix 5 for the details of types labelling). Classes ${D}^{15}\text{}3$ and ${D}^{15}\text{}5$ also are also similar each other. All three types diverge from each other as early as possible and then two of them equally share the subsequent mutations (j and k in ${D}^{15}\text{}3$, or i and k in ${D}^{15}\text{}5$). In the class ${D}^{15}\text{}4$, all three types diverge from each other as early as possible. In this case, each of three types have the similar number of accumulated mutations (33%).
Minimization of the number of noncyclic dominances
Next, we ask which genealogies instead minimize the number of noncyclic dominance triplets, with the expectation that they will be similar to the ones above. Performing PCA analysis, we again find five different classes, see Appendix 6—table 2. The obtained genealogies minimizing T_{15} have the same characteristics of mutational distances with the genealogies maximizing T_{16}, see section F. The genealogy class ${D}_{16}\text{}1$ is the same with ${D}^{15}\text{}1$. The class ${D}_{16}\text{}2$ is the same as ${D}^{15}\text{}2$, and ${D}_{16}\text{}3,4,5$ are the same as ${D}^{15}\text{}4$. Classes ${D}_{16}\text{}3$, ${D}_{16}\text{}4$, and ${D}_{16}\text{}5$ differ only in their values of w_{1} and x. Thus these can be considered as a fine structure of a single class. Classes equivalent to ${D}^{15}\text{}3$ and ${D}^{15}\text{}5$ are not found in D_{16}.
Minimization of the number of cyclic dominances
Alternatively, we can also ask for which genealogies it is hardest to obtain cyclic dominances. Performing PCA analysis, we find three different classes, see Appendix 6—table 3. In the class ${D}_{15}\text{}1$, most of mutations are accumulated in the lineage of the type i before types j and k are diverged. The divergence of types j and k tends to be the last mutational event in this genealogy. In the class ${D}_{15}\text{}2$, most mutations are accumulated in the lineage of the common progenitor of types j and k before their divergence. The divergence of types j and k tends to be the last mutational event in this genealogy. In the class ${D}_{15}\text{}3$, mutations are equally shared between the lineage of type i and the common progenitors of types j and k before their divergence. Similar to above two cases, the divergence of types j and k tends to be the last mutational event in this genealogy.
Maximization of the number of noncyclic dominances
Now we ask for which genealogies it is easiest to obtain noncyclic dominances. Performing PCA analysis, we find again three different classes of geneaolgies, see Appendix 6—table 4. Again, the obtained three classes are the same as minimizing cyclic dominances T_{15}: The class ${D}^{16}\text{}1$ is the same as ${D}_{15}\text{}1$, and ${D}^{16}\text{}2$ is the same as ${D}_{15}\text{}2$, and ${D}^{16}\text{}3$ is the same as ${D}_{15}\text{}3$.
Summary for optimization
From all optimization results we find two extreme genealogies for maximizing and minimizing the fraction $\chi$ of cyclic dominances. In the first class, most of mutations occur after all three lineages separate and these mutations are accumulated in a single lineage, see Appendix 6—table 5. We call them maximizer genealogies. These genealogies promote cyclic dominances and suppress noncyclic dominances. In the second class, the most of mutations occur before types j and k diverge. These genealogies suppress cyclic dominances, while promoting noncyclic dominances. We call them maximizer genealogies.
Analytics for maximizer genealogies
We calculate the proportions of cyclic and noncyclic dominance triplets at maximizer genealogies. In that case, the cyclic dominance triplets occur
$i\to j\to k\to i$ with probability $1/96$
$i\to k\to j\to i$ with probability $1/96$
Altogether, this results in $p({T}_{15})=1/48\approx 0.0208$. The noncyclic dominance triplet T_{16} has degeneracy six, and they occur
$j\to k\to i$ with probability $1/96$
$i\to k\to j$ with probability $1/96$
$k\to j\to i$ with probability $1/48$
$k\to i\to j$ with probability $1/48$
$j\to i\to k$ with probability $1/48$
$i\to j\to k$ with probability $1/48$
Altogether, this results in $p({T}_{16})=5/48\approx 0.104$, yielding $\chi =1/6\approx 0.167$. The results agree well with numerical results in Appendix 6—table 5.
Minimal model
Besides genealogies, there is also another way to shape the payoff correlation. Mutational distances are controlled by not only how many steps proceed but also how big the jump is related variance of noise. In a minimal model, we can directly control the closeness of newborn type by using a different variance $\sigma}^{2$ of Gaussian noise. We generate the first mutant type M_{1} from the original type o with the normal distribution, and we use variance $\sigma}_{1}^{2$ when we generate the type M_{1}. Then, we assume that M_{1} is the parental type of type M_{2}, and the variance $\sigma}_{2}^{2$ is used in this case. This procedure is schematically drawn in Appendix 6—figure 1A. Note that only the relative mutational distances are important, and thus we introduce the fraction $F}_{l$ that measures the distances before and after the divergence of the third type. When $F}_{l$ is small, the close type to the resident type arises which is corresponding to the minimizer genealogy. On the other hand, if we use large $F}_{l$, it will reproduce the maximizer by introducing uncorrelated types in the population. The minimum and maximum average fractions $\chi$ are obtained by changing $F}_{l$, yielding almost zero and 0.165 respectively, see Appendix 6—figure 1B. These results agree well with the results of the genealogy approach.
Appendix 7
Correlation between diversity and triplet fraction
We investigated the correlation between diversity and the probability of having cyclic and noncyclic dominance triplets in the system at the steadystate. We measured the Pearson correlation coefficient $r}_{xy$ between two variables x and y defined as
where $\overline{x}$ and $\overline{y}$ indicate the averages of the variables x and y. The index i runs over all samples (in our case, 1944 surviving realizations are used). The correlation coefficient is bound from 1 to 1, representing fully anticorrelated to fullycorrelated variables.
We used three different diversity indices: (1) richness n as the number of different types, (2) Shannon diversity index $H=\sum _{i=1}^{n}{f}_{i}\mathrm{ln}{f}_{i}$, and (3) Simpson diversity index $S={\left(\sum _{i=1}^{n}{f}_{i}^{2}\right)}^{1}$. Note that a population frequency of type i is denoted by $f}_{i$. Correlations were measured between each diversity index and the probabilities to having cyclic or noncyclic dominance triplets, $P}_{cyc$ and $P}_{ncyc$, see Appendix 7—figure 1. All correlation coefficients are summarized in Appendix 7—table 1. The results show weak anticorrelation between diversity and cyclic dominance and stronger anticorrelation between diversity and noncyclic dominance. This tendency does not change even when population frequencies are considered. Finally we measured the correlation coefficient between the fraction of cyclic dominance $\chi$ and diversity, and surprisingly the coefficient showed almost zero. This implies that there is no significant correlation between cyclic dominance and diversity.
Anticorrelation between the fraction of noncyclic dominance $P}_{ncyc$ and diversity could be understood from the population frequencies of noncyclic dominance $f}_{ncyc$. In the steadystate, noncyclic dominance usually occupies the majority of the population, taking around 0.83 of the population on average, see Appendix 7—figure 2AB. It seems that the noncyclic dominance emerges from the majority, and thus their fraction could be high when the winner takes over almost the whole population, leading to low diversity. On the other hand, cyclic dominance usually takes only a minority, showing a peak near ${f}_{cyc}\approx 0$ in Appendix 7—figure 2C.
Data availability
Our simulation code is available at https://github.com/ParkHyeJin/CyclicDominance (copy archived at https://github.com/elifesciencespublications/CyclicDominance).
References

LotkaVolterra equation and replicator dynamics: new issues in classificationBiological Cybernetics 72:447–453.https://doi.org/10.1007/BF00201420

Colicin biologyMicrobiology and Molecular Biology Reviews 71:158–229.https://doi.org/10.1128/MMBR.0003606

On the prevalence and relative importance of interspecific competition: evidence from field experimentsThe American Naturalist 122:661–696.https://doi.org/10.1086/284165

Allelopathy in spatially distributed populationsJournal of Theoretical Biology 185:165–171.https://doi.org/10.1006/jtbi.1996.0292

Spatial aspects of interspecific competitionTheoretical Population Biology 53:30–43.https://doi.org/10.1006/tpbi.1997.1338

Rockscissorspaper and the survival of the weakestProceedings of the Royal Society of London. Series B: Biological Sciences 268:1323–1327.https://doi.org/10.1098/rspb.2001.1670

Community structure follows simple assembly rules in microbial microcosmsNature Ecology & Evolution 1:109.https://doi.org/10.1038/s415590170109

Ecological systems biology: the dynamics of interacting populationsCurrent Opinion in Systems Biology 1:114–121.https://doi.org/10.1016/j.coisb.2016.12.001

The effects of intransitive competition on coexistenceEcology Letters 20:791–800.https://doi.org/10.1111/ele.12775

Evolution of community structure: competitionJournal of Theoretical Biology 133:513–523.https://doi.org/10.1016/S00225193(88)803382

The competitive exclusion principleScience 131:1292–1297.https://doi.org/10.1126/science.131.3409.1292

Replicator dynamics for optional public good gamesJournal of Theoretical Biology 218:187–194.https://doi.org/10.1006/jtbi.2002.3067

BookEvolutionary Games and Population DynamicsCambridge, United Kingdom: Cambridge University Press.https://doi.org/10.1017/CBO9781139173179

Evolutionary game dynamicsBulletin of the American Mathematical Society 40:479–520.https://doi.org/10.1090/S0273097903009881

Emergence of evolutionarily stable communities through ecoevolutionary tunnellingNature Ecology & Evolution 2:1644–1653.https://doi.org/10.1038/s4155901806557

Population size and the rate of evolutionTrends in Ecology & Evolution 29:33–41.https://doi.org/10.1016/j.tree.2013.09.009

Ecosystems with mutually exclusive interactions selforganize to a state of high diversityPhysical Review Letters 107:188101.https://doi.org/10.1103/PhysRevLett.107.188101

BookEvolution and the Theory of GamesCambridge: Cambridge University Press.https://doi.org/10.1017/CBO9780511806292

Emergence of diversity in a model ecosystemPhysical Review E 86:011929.https://doi.org/10.1103/PhysRevE.86.011929

Microbial syntrophy: interaction for the common goodFEMS Microbiology Reviews 37:384–406.https://doi.org/10.1111/15746976.12019

Emergence of structured communities through evolutionary dynamicsJournal of Theoretical Biology 383:138–144.https://doi.org/10.1016/j.jtbi.2015.07.020

Models of densitydependent genic selection and a new rockpaperscissors social systemThe American Naturalist 170:663–680.https://doi.org/10.1086/522092

Models of life: epigenetics, diversity and cyclesReports on Progress in Physics 80:042601.https://doi.org/10.1088/13616633/aa5aeb

BookNonlinear Dynamics and Chaos: With Applications to Physics, Biology, Chemistry, and Engineering (Studies in Nonlinearity)Westview Press.

Evolutionary games on graphsPhysics Reports 446:97–216.https://doi.org/10.1016/j.physrep.2007.04.004

Cyclic dominance in evolutionary games: a reviewJournal of the Royal Society Interface 11:20140735.https://doi.org/10.1098/rsif.2014.0735

Fast evolution unlocks forbidden communitiesNature Ecology & Evolution 2:1525–1526.https://doi.org/10.1038/s415590180688y

AlleleFrequency change in a ClassStructured populationThe American Naturalist 135:95–106.https://doi.org/10.1086/285034

Effect of population size and mutation rate on the evolution of RNA sequences on an adaptive landscape determined by RNA foldingInternational Journal of Biological Sciences 13:1138–1151.https://doi.org/10.7150/ijbs.19436

Inhibitory interactions promote frequent bistability among competing BacteriaNature Communications 7:11274.https://doi.org/10.1038/ncomms11274

Evolutionary dynamics of species diversity in an interaction web systemEcological Modelling 163:131–143.https://doi.org/10.1016/S03043800(02)004179
Article and author information
Author details
Funding
JRG Program at the APCTP through Science and Technology Promotion Fund and Lottery Fund of the Korean Government and Korean Local Government Gyeongsangbukdo Province and Pohang City
 Hye Jin Park
The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.
Acknowledgements
HJP, YP, and AT thank the Max Planck Society for generous funding. HJP was also supported by (1) an appointment to the JRG Program at the APCTP through the Science and Technology Promotion Fund and Lottery Fund of the Korean Government and by (2) the Korean Local Governments  Gyeongsangbukdo Province and Pohang city.
Version history
 Received: April 14, 2020
 Accepted: August 1, 2020
 Version of Record published: September 4, 2020 (version 1)
Copyright
© 2020, Park et al.
This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.
Metrics

 2,155
 Page views

 212
 Downloads

 14
 Citations
Article citation count generated by polling the highest count across the following sources: Crossref, PubMed Central, Scopus.
Download links
Downloads (link to download the article as PDF)
Open citations (links to open the citations from this article in various online reference manager services)
Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)
Further reading

 Ecology
 Evolutionary Biology
In the unpredictable Anthropocene, a particularly pressing open question is how certain species invade urban environments. Sexbiased dispersal and learning arguably influence movement ecology, but their joint influence remains unexplored empirically, and might vary by space and time. We assayed reinforcement learning in wildcaught, temporarily captive core, middle, or edgerange greattailed grackles—a bird species undergoing urbantracking rapid range expansion, led by dispersing males. We show, across populations, both sexes initially perform similarly when learning stimulusreward pairings, but, when reward contingencies reverse, male—versus female—grackles finish ‘relearning’ faster, making fewer choiceoption switches. How do male grackles do this? Bayesian cognitive modelling revealed male grackles’ choice behaviour is governed more strongly by the ‘weight’ of relative differences in recent foraging payoffs—i.e., they show more pronounced risksensitive learning. Confirming this mechanism, agentbased forward simulations of reinforcement learning—where we simulate ‘birds’ based on empirical estimates of our grackles’ reinforcement learning—replicate our sexdifference behavioural data. Finally, evolutionary modelling revealed natural selection should favour risksensitive learning in hypothesised urbanlike environments: stable but stochastic settings. Together, these results imply risksensitive learning is a winning strategy for urbaninvasion leaders, underscoring the potential for life history and cognition to shape invasion success in humanmodified environments.

 Ecology
The understanding of ecoevolutionary dynamics, and in particular the mechanism of coexistence of species, is still fragmentary and in need of test bench model systems. To this aim we developed a variant of SELEX in vitro selection to study the evolution of a population of ∼10^{15} singlestrand DNA oligonucleotide ‘individuals’. We begin with a seed of random sequences which we select via affinity capture from ∼10^{12} DNA oligomers of fixed sequence (‘resources’) over which they compete. At each cycle (‘generation’), the ecosystem is replenished via PCR amplification of survivors. Massive parallel sequencing indicates that across generations the variety of sequences (‘species’) drastically decreases, while some of them become populous and dominate the ecosystem. The simplicity of our approach, in which survival is granted by hybridization, enables a quantitative investigation of fitness through a statistical analysis of binding energies. We find that the strength of individual resource binding dominates the selection in the first generations, while inter and intraindividual interactions become important in later stages, in parallel with the emergence of prototypical forms of mutualism and parasitism.