Abstract
Evolution is simultaneously driven by a number of processes such as mutation, competition and random sampling. Understanding which of these processes is dominating the collective evolutionary dynamics in dependence on system properties is a fundamental aim of theoretical research. Recent works quantitatively studied coevolutionary dynamics of competing species with a focus on linearly frequencydependent interactions, derived from a gametheoretic viewpoint. However, several aspects of evolutionary dynamics, e.g. limited resources, may induce effectively nonlinear frequency dependencies. Here we study the impact of nonlinear frequency dependence on evolutionary dynamics in a model class that covers linear frequency dependence as a special case. We focus on the simplest nontrivial setting of two genotypes and analyse the coaction of nonlinear frequency dependence with asymmetric mutation rates. We find that their coaction may induce novel metastable states as well as stochastic switching dynamics between them. Our results reveal how the different mechanisms of mutation, selection and genetic drift contribute to the dynamics and the emergence of metastable states, suggesting that multistability is a generic feature in systems with frequencydependent fitness.
1. Introduction
Selection, random genetic drift and mutations are the processes underlying Darwinian evolution. For a long time, population geneticists have analysed the dynamics in the simplest setting consisting of two genotypes evolving under these processes [1]. In those studies, a genotype represents an individual's genetic makeup, completely determining all relevant properties of the individual. A key concept is the socalled fitness of a genotype that represents the selection pressure for the individuals. The fitness defines the expected number of offspring an individual will produce. Thus, selection acts on fitness differences preferring individuals with higher fitness over individuals with lower fitness. Usually it is assumed that individuals have fixed fitnesses defined by their genotype alone [1–3]. Yet, experimental studies have revealed that many natural systems exhibit frequencydependent selection [4–6], which means that an individual's fitness not only depends on its genotype, but also on its interactions with other individuals and hence on the frequency of the different genotypes in the population. Although such frequencydependent selection had already been studied early by Crow & Kimura [3], only recently has it received more attention [7–14]. In these theoretical and computational studies, individuals' interactions are represented by interaction matrices from game theory. This leads to a frequency dependence where the fitness depends directly on the interaction parameters in a linear way. However, fitness may depend on many diverse factors such as cooperation (i.e. individuals acting together to increase their fitness [12,13]) and resource competition, so that certain systems may exhibit frequencydependent fitness that is nonlinear. For example, in experiments, certain hermaphrodites exhibit such nonlinear fitnessdependence [6]. To the best of our knowledge, the impact of such nonlinear dependencies on coevolutionary dynamics has not been investigated theoretically.
In this article, we show that nonlinear frequency dependence [6,15] may induce new stable configurations of the evolutionary dynamics. Furthermore, we study the impact of asymmetric mutation probabilities on the dynamics [16,17], which was also neglected in most models until now [7,9,13,14]. As in previous works on coevolutionary dynamics, we base our work on the Moran process in a nonspatial environment, which is a wellestablished model to study evolutionary dynamics and was already used in many applications [18]. The Moran process is a stochastic birth–death process that keeps the population size constant [19]. Therefore, in a twogenotype model, the system becomes effectively onedimensional, so that the dynamics may be described by a onedimensional Markov chain with transition rates defined by the Moran process. We derive the stationary probability distribution of the system dynamics via its Fokker–Planck equation [20]. Sharp maxima of the distribution reveal metastable points of the dynamics and a multitude of such maxima lead to stochastic switching dynamics between multiple stable points.
The article is structured as follows. In §2, we introduce the model details, and in §3, we derive the Fokker–Planck equation describing the probabilistic dynamics of the population. Using this equation, we derive the stationary probability distribution that describes the longtime behaviour of the system. In §4, we analyse this probability distribution, which yields information about the impact of nonlinear frequencydependent selection and of different mutation probabilities on the coevolutionary dynamics. In §5, we give a summary and discuss our results.
2. The model
Consider a population of N individuals evolving in a nonspatial environment, meaning that the population is wellmixed so that each individual interacts with all other individuals at all times. In this population, the individuals may assume one out of two genotypes A and B. The population sizes and () evolve according to the timecontinuous Moran process described in the following (cf. [19]). The number of individuals of genotype A completely determines the state of the system as . At all times, the interactions of the individuals determine the actual (frequencydependent) fitness, so that an individual's fitness of genotype A or B is defined by a fitness function or respectively. The fitness functions and may be any functions with the only condition that for all as negative fitness is not defined. At rate , an individual of type A produces an identical offspring that may mutate to genotype B with probability . This applies analogously to genotype B. Then one individual of the population is chosen randomly to die, so that the population size N stays constant and the variables change maximally by 1. This is the socalled Moran process, which was originally introduced for a population of two genotypes with fixed fitnesses and no mutations occurring (cf. [19]). However, this process is easily generalizable to more genotypes and frequencydependent selection in the way described already (cf. [11,21]).
Note that in our model, the rate of reproduction of, e.g. type A is directly determined by the term as, e.g. in [22]. In other models, the fitness is first normalized so that the rate of reproduction is given by [13], where 2.1is the population's mean fitness. While in both of these models the events occur with the same probability, the times between the events differ by a common factor determined by the mean fitness (2.1). Thus, these models exhibit a quantitatively different time course, but the same event sequences.
Until now, usually linear functions have been considered for . However, in many applications, nonlinear functions seem more appropriate [6,15]. For example, cooperation effects in game theory induce functions linearly increasing in , so that [11,13]. This would mean that fitness increases infinitely with the population size of genotype A. Yet, in any habitat there is only a limited amount of resources available for living. Therefore, fitness should decline when the population of one genotype becomes too large as all of the individuals will compete for the same resources. We conclude that a fitness function including these two effects has to contain a nonlinear factor, e.g. 2.2where and . Here we choose a quadratic nonlinear factor as a simple example. Naturally, all other nonlinear functions could be applicable as well.
Note that linearly increasing fitness functions such as and lead to quadratic relations for the mean fitness (see equation (2.1)) of the population, which becomes 2.3The quadratic factor originates from the fact that the linear fitness function in (2.1) is multiplied by the number of individuals, leading to a quadratic dependence on , e.g. . This factor enters the replicator equation in evolutionary game theory, thus leading to a quadratic fitness dependence in standard gametheoretic problems [11,22]. However, this does not imply nonlinear interactions and thus fitness effects such as presented in equation (2.2). Therefore, the fitness functions that can be generated by evolutionary game theory are a special case of the functions that occur in the theory presented here. Fitness functions as for example presented in equation (2.2) will—as we show in this article—lead to dynamics with more complex stability properties.
Usually, when analysing the dynamics of such a model system as described earlier, it was assumed either that no mutations occur at all () [10,11,19] or that the mutation probabilities are equal ) [13,14,21]. However, usually mutation probabilities can be highly diverse [16,17]. This has been considered in some examples [3], but until now the effect of different mutation probabilities in the introduced model system has not been studied systematically. Here, we explicitly study the effects caused by asymmetric mutation rates (see figure 2).
3. Analysis
In the following, we assume the fitness of the individuals to be given by and 3.1 with the interaction functions and . The model system is effectively onedimensional with the variable completely describing the system state [12]. At each event k may change by at most , depending only on the actual state of the system, so that the dynamics are Markovian. Let the transition occur with the rate . Then this rate is determined by the fitnesses and mutation probabilities in the following way.

— Each individual in the population receives offspring with rate with . This means that genotype A receives one offspring with rate .

— Such an offspring increases the population A with probability , and population B in case it mutates with probability .

— One individuum is chosen uniformly to die belonging to genotype B with probability .

— Equivalently, the number of individuals of genotype A can increase if genotype B receives one offspring, which mutates to genotype A with probability , and one individuum of B is chosen to die.
Taken these processes together, the transition rate for is 3.2Decreasing k to is equivalent to increasing the number of individuals of genotype B, so that the transition rate 3.3is obtained analogously. The master equation of the Markov chain is then given by 3.4where we define . As also , the system has reflective boundaries at and . If (), then () is an absorbing state of the dynamics.
It is possible to derive the exact solution to such an equation, similar to the derivations in [23], which we do in appendix B. However, the obtained solutions are not very explicit and their implications are hard to grasp. It is thus more useful to advance to the Fokker–Planck equation [20] describing the system in the limit of large N. As previous works have shown, this approximation already leads to very good results for moderate population sizes of (cf. [13]). We use the transformation 3.5and 3.6yielding the Fokker–Planck equation. The scaling factors are derived in appendix A, where the derivation of the Fokker–Planck equation in the scaling limit is shown in detail (equation (A 5)). This leads to 3.7with the normalization condition 3.8for all . Here and are the rescaled interaction functions defined in (3.6), is the rescaled mean mutation rate and is the rescaled mutation rate difference. These functions have to be rescaled to obtain a nondegenerate limit, the socalled weak selection regime.
Notice, that the drift term 3.9contains three different effects. The fitness difference at each point x causes a selective drift, the mean mutation rate causes a drift directed to and the mutation rate difference causes a onedirectional drift towards one of the boundaries x = 0 or x = 1. The diffusion term 3.10reflects the undirected genetic drift.
If fitness differences are on a scale of , then mutation, selection and genetic drift all act on the same scale [13,24]. Then the dynamics are characterized by an interplay of these different effects. On the other hand, in the strong selection regime—where fitness differences are on a scale —genetic drift becomes negligible and the dynamics in the bulk of the system become deterministic for large N [25]. However, as the factor vanishes for and the effects of mutations play an important role. Clearly, without mutations, x = 0 and x = 1 are absorbing states, but even in the strong selection regime, they become nonabsorbing for any positive mutation rate. This is reflected in the form of the drift term , as only the factors reflecting mutations remain nonzero in the limits and .
For a onedimensional Fokker–Planck equation (3.7), the stationary distribution is [20] 3.11with the potential 3.12which is defined up to a constant irrelevant to our calculations as it is cancelled by the normalization 3.13Thus, we obtain where we used 3.14The stationary solution becomes 3.15where the constant C defined in (3.13) has to be calculated numerically for given parameters.
This stationary distribution contains contributions from the earliermentioned four different effects.

— Genetic drift is reflected by , which diverges for at the boundaries.

— The mean mutation rate causes the balancing term .

— The asymmetry in the mutation probabilities causes the term , which diverges at one boundary and vanishes at the other.

— All frequencydependent selection effects are contained in the exponential factor , which can take various shapes.
4. Multiple stable points
What are the possible shapes for the stationary distribution in the form given by equation (3.15)? Until now, the term representing the asymmetric mutation probabilities has not been considered to our knowledge, and the interaction functions and in the selection term were considered to be at most linear in x [13]. We should therefore be interested in the effects of nonlinear interaction functions and asymmetric mutation rates.
Let us first analyse the dynamics for nonlinear interaction functions [6,15] describing the effects of cooperation and limited resources as described by equation (2.2), so that and 4.1Already with such interaction functions induce dynamics stochastically switching between three metastable points. An example is shown in figure 1, where the theoretically calculated stationary distribution from equation (3.15) is shown together with data from simulations with a population of N = 1000 individuals which is enough to obtain almost perfect fitting (cf. also figure 4). The fitness functions (figure 1b) both show at first an increase on increasing the number of individuals of genotype A or B from 0 owing to cooperation effects and then a strong decrease owing to resource competition. As the resulting fitness functions are asymmetric, also the stationary distribution is asymmetric. There is a maximum at x = 0 owing to genetic drift. Selection drives the dynamics towards a metastable state at , because for genotype A is fitter than B thus increasing in frequency and for genotype A is less fit than B, thus decreasing in frequency (cf. figure 1b). The maximum of the stationary distribution is thus exactly at the point where . For again genotype A is fitter than B, and thus the dynamics are driven towards x = 1 by selection as well as genetic drift. The mutational force induced by increases the height of the maximum at , driving the system away from the maxima at x = 0 and x = 1. So, in this example, genetic drift, mutation and selection all significantly influence the dynamics.
Let us now study the influence of the mutation rates in detail. Interestingly, for asymmetric mutation rates ) the factor always diverges in the interval [0, 1]. For , it diverges for ; otherwise for . This can cause the emergence of a maximum of the stationary distribution at x = 0 (or, respectively, x = 1). Figure 2 shows an example where owing to asymmetric mutation rates a maximum occurs at x = 1, when for there is an absolute minimum at x = 1 as this stable state minimizes the population's mean fitness. This means that the system dynamics are mutation dominated near x = 1. Furthermore, the maximum located at x = 0 for is shifted for to values , causing a minimum of the stationary distribution at x = 0. Thus, in this example, both selection and the asymmetry in mutation probabilities are the driving forces of the system's dynamics.
For low mutation rates, the population has a high tendency to become dominated by one of the two genotypes for long times, which is illustrated in figure 3a. Thus, the dynamics stay close to the edges x = 0 or x = 1 waiting for one of the few mutations to occur. For low overall mutation rates, differences in the mutation rates have only weak effects. High mutation rates cause the population to be drawn towards a mixture of both genotypes (). Even more, for high mutation rates the differences in the mutation rates can have an important influence on the dynamics. For instance, figure 3b illustrates a shift of an existing stable point, which completely vanishes for high . We conclude that asymmetric mutation rates can cause the emergence of new stable states (as in figure 2), the disappearance of existing ones (as in figure 3), and shift existing stable states to new positions.
As we saw earlier, the fit of the theoretically calculated stationary distribution to simulation data is almost perfect for weak selection, i.e. the maximal fitness difference satisfies (cf. [13]). To quantify the quality of the fit, we measured the stationary distribution with the same parameters as mentioned earlier for different population sizes N. This is shown in figure 4, demonstrating that for increasing N the data get ever closer to the theoretically obtained distribution. For this, we define the empirical distribution 4.2where is the discrete process defined in equation (3.4). is a time large enough for the process to reach stationarity and depends on the system size, as does the measurement time which is specified in figure 4. Using this definition, we quantify the quality of the fit using the distance measure 4.3which takes the mean distance of the measured empirical distribution from the theoretical distribution for every point except k = 0 and k = N where the density diverges. Here, the theoretical distribution is calculated by integrating the theoretical density over bins of the size 1/N around the points k/N. Figure 4b shows that the mean distance decays with increasing N.
The decay of this quantity is dominated by the slow convergence at the domain boundaries k = 0 and k = N. Therefore, we additionally defined the maximum distance measure 4.4which gives the maximum distance between the measured probability distribution from the theoretical distribution for all points except k = 0 and k = N due to the same argument as above. Figure 4b shows that the maximum distance decreases with increasing N, however slower than the mean distance owing to the divergence of at the boundaries. Altogether we conclude that the theoretical curve fits the data very well already for . However, near the domain boundaries, special care has to be taken, as divergences of the theoretical curve lead to larger deviations.
Furthermore, as the Fokker–Planck approximation only holds for weak selection, we study the quality of the solution obtained through the Fokker–Planck equation in dependence of the selection strength. For this, we introduce a scaling factor to the interaction functions, so that they become 4.5We find that the stationary solution of the Fokker–Planck equation well approximates the solution of the Master equation 4.6for selection strengths up to of the order of 1 (cf. figure 5). We derive the solution (4.6) of the Master equation in appendix B. Note that for the derivation of this solution, no approximation is necessary and it hence fits data from simulations perfectly up to an error due to finite sampling in simulations. As figure 5 illustrates, for strong selection , the solution of the Fokker–Planck equation does not fit the exact solution of the Master equation perfectly, but still catches the overall trend of the dynamics. The error, quantified by the mean distance measure analogous to (4.3), increases with as a power law both in comparison with the solution from the Master equation and with data from simulations (cf. figure 5d). All in all, for weak selection the Fokker–Planck approximation works well, while for strong selection, it does not perfectly predict the stationary distribution, but still reflects the overall trend of the dynamics. If this approximation is not satisfying, the direct solution (4.6) of the Master equation yields the exact distribution fitting the data for all selection strengths.
We have shown earlier that the dynamics of a twogenotype system can exhibit multiple stable points induced by nonlinear selection. Even more, there is no theoretical limit to the number of stable points in the system. Assume for example and , so that mutation exactly balances genetic drift. Then, the stationary distribution (3.15) has a maximum at each point, where has a maximum. Theoretically, and can be any function with an arbitrary amount of extreme points in [0, 1], so that there is no limit to the stationary distribution's number of maxima, if selection dominates the dynamics. However, for finite N the number of possible maxima is naturally limited by N/2. Figure 6 shows an example, where we used periodic interaction functions 4.7 Although this is not a realistic interaction function in most applications, it demonstrates what is theoretically possible in the introduced system.
5. Summary
In this article, we analysed a twogenotype system in a very general setting with (possibly) asymmetric mutation probabilities and nonlinear fitness functions [6,15] in finite populations. The underlying Moran process is a wellestablished model [1,12,13] to gain an understanding of the interplay of selection, mutation and genetic drift in evolutionary dynamics. However, the Moran process is studied mostly with symmetric mutation probabilities and at most linear interaction functions. We reasoned that neither need mutation probabilities be symmetric—as experiments have shown that mutation probabilities are often asymmetric [16]—nor can all interaction effects be described by linear interaction functions. For example, cooperation in game theory leads to an interaction function increasing linearly in the frequency of the cooperating genotype. Yet, in many applications, also cooperators in the end compete for the same type of resource that is limited. Therefore, owing to limited resources, a population being too large cannot be sustained, leading to a decrease in the fitness. There is no linear function that can reflect both of these effects at the same time.
We derived the Fokker–Planck equation describing the dynamics of the number of individuals k of genotype A in the limit of large population sizes N. We quantified the quality of the Fokker–Planck approach for an example (cf. figure 4), where the difference of simulation data and theoretical solution became almost not detectable for population sizes larger than . Actually, if the system exhibits absorbing states, then the Fokker–Planck method does not work to study the corresponding quasistationary distributions. Instead, Wentzel–Kramers–Brillouin methods are more appropriate to describe the system dynamics as, for example, in [10], where fixation resulting from large fluctuations was studied. In our model system, no such absorbing states exist, as long as the mutation rates are positive and therefore the Fokker–Planck equation is appropriate to describe the system dynamics.
We identified the individual effects of selection, mean mutation rate and mutation difference as well as genetic drift and derived the stationary probability distribution as determined by the Fokker–Planck equation. Analysing the distribution, we found that asymmetries in the mutation probabilities may not only induce the shifting of existing stable points of the dynamics to new positions, but also lead to the emergence of new stable points. Thus, a genotype that has a selective disadvantage can anyway have a stable dynamical state where its individuals dominate the population owing to a higher mutational stability (figure 2). Further, we found that dynamic fitness leads to multiple stable points of the dynamics induced by selection and also genetic drift. We showed an example (figure 1) where three stable points exist, two caused by selection and one by genetic drift.
We conclude that frequencydependent fitness together with asymmetric mutation rates induce complex evolutionary dynamics—in particular, if the interactions imply nonlinear fitness functions. Theoretically, there is no limit for the number of stable points that the dynamics can exhibit (figure 6). All in all, we interpret our results such that in real biological systems, multiple metastable equilibria may exist whenever species interact in a way complex enough to imply a fitness that nonlinearly depends on frequency. As a consequence, one species may exhibit a certain frequency for a long time before a sudden shift occurs and then a new frequency prevails. Such a change may thus occur even in the absence of changes of the environment; it may be induced as well by a stochastic switching from one metastable state to another, owing to complex interspecies interactions.
The Moran process is a standard tool to gain theoretical insights into experimental data [18]. Of course, we do not propose here that any experimental setup may be exactly described by the Moran process. However, we think that it should be feasible to develop an experiment where two different mutants evolve with asymmetric mutation rates. To find an experimental setup where the two genotypes also exhibit nonlinear fitness could however prove more difficult. Rather, our study is a theoretical study indicating that nonlinear fitness may be the cause of multiple stable states when observed in experimental data.
For further studies on systems with more genotypes, it should be useful to combine our considerations presented here with the work of Traulsen et al. [13], where an analysis of systems with more than two genotypes was carried out. Extending those results, it may be possible to gain a better understanding of the effects of nonlinear interactions for many different genotypes. Also, it may be interesting to study the effects of changing interactions, where the interactions change according to the system dynamics [26]. Thus, our study might serve as a promising starting point to investigate how nonlinear frequency dependencies impact evolutionary dynamics in complex environments.
Acknowledgements
We thank Steven Strogatz for fruitful discussions during project initiation and Stefan Eule for helpful technical discussions. Stefan Grosskinsky acknowledges support by EPSRC, grant no. EP/E501311/1.
Appendix A. The Fokker–Planck equation
The master equation (3.4) may be transformed to a Fokker–Planck equation in the limit of large N [20]. We introduce the transformation A1together with the rescaled functions A2We fix the scaling functions F(N), G(N) and H(N) such that in the limit , all terms in equation (3.4) remain finite so that mutation, selection and genetic drift all act on the same scale. We further define A3and substituting all this into the master equation (3.4), we obtain
A4
We choose F(N)=1/(N + 1) and G(N) = H(N) = N so that in the limit the terms stay finite. Further, we introduce the mean mutation rate and the mutation rate difference . To not overload the notation, we drop the time argument s of in the following calculation. This leads to
In the limit , the different terms with N^{2} in front become secondorder derivatives with respect to x, while the other terms become firstorder derivatives. The terms that have a 1/N factor vanish in the limit and thus the earliermentioned equation becomes A5which is the Focker–Planckequation of the system.
Appendix B. Stationary solution of the master equation
We directly derive the stationary solution of the master equation (3.4) using the detailed balance equation B1which applies to any chain with only nearest neighbour transitions [20] (cf. also [23]). Thus, by using the rewritten balance equation B2iteratively, we obtain B3Finally, we may use the normalization condition B4of the stationary distribution to eliminate the factor . We then obtain the exact stationary solution of the master equation B5which can be evaluated numerically for any transition rates and . For more details on the exact solution of the Master equation in the Moran process, see, for example, the work by Claussen & Traulsen [23].
 Received June 11, 2012.
 Accepted July 13, 2012.
 This journal is © 2012 The Royal Society