document: selection , random genetic drift and mutations are the processes underlying darwinian evolution . for a long time 
   population geneticists have analyzed the dynamics in the simplest setting consisting of two genotypes evolving under these processes @xcite . in those studies 
   , a genotype represents an individual s genetic makeup , completely determining all relevant properties of the individual . 
   a key concept is the so - called fitness of a genotype which 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 @xcite . yet 
   , experimental studies have revealed that many natural systems exhibit frequency - dependent selection @xcite , 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 frequency - dependent selection had already been studied early by crow and kimura @xcite , only recently has it received more attention @xcite . 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 @xcite ) and resource competition , so that certain systems may exhibit frequency - dependent fitness that is nonlinear . 
   for example , in experiments certain hermaphrodites exhibit such nonlinear fitness - dependence @xcite . 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 @xcite may induce new stable configurations of the evolutionary dynamics . 
   furthermore , we study the impact of asymmetric mutation probabilities on the dynamics @xcite , which was also neglected in most models until now @xcite . as in previous works on coevolutionary dynamics we base our work on the moran process in a non - spatial environment which is a well established model to study evolutionary dynamics and 
   was already used in many applications @xcite . 
   the moran process is a stochastic birth - death process which keeps the population size constant @xcite . 
   therefore , in a two - genotype model the system becomes effectively one - dimensional , so that the dynamics may be described by a one - dimensional 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 @xcite . 
   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 section 
   ii we introduce the model details and in section iii 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 long - time behavior of the system . in section 
   iv we analyze this probability distribution , which yields information about the impact of nonlinear frequency - dependent selection and of different mutation probabilities on the coevolutionary dynamics . in section 
   v we give a summary and discuss our results . 
   consider a population of @xmath0 individuals evolving in a non - spatial environment , meaning that the population is well - mixed so that each individual interacts with all other individuals at all times . in this population 
   the individuals may assume one out of two genotypes @xmath1 and @xmath2 . 
   the population sizes @xmath3 and @xmath4 ( @xmath5 ) evolve according to the time - continuous moran process described in the following , cf . 
   the number of individuals @xmath3 of genotype @xmath1 completely determines the state of the system as @xmath6 . at all times 
   the interactions of the individuals determine the actual ( frequency - dependent ) fitness , so that an individual s fitness of genotype @xmath1 or @xmath2 is defined by a fitness function @xmath7 or @xmath8 respectively . 
   the fitness functions @xmath7 and @xmath8 may be any functions with the only condition that @xmath9 for all @xmath10 $ ] as negative fitness is not defined . at rate 
   @xmath11 an individual of type @xmath1 produces an identical offspring which may mutate to genotype @xmath2 with probability @xmath12 . 
   this applies analogously to genotype @xmath2 . 
   then one individual of the population is chosen randomly to die , so that the population size @xmath0 stays constant and the variables @xmath13 change maximally by @xmath14 . 
   this is the so - called moran process which was originally introduced for a population of two genotypes with fixed fitnesses and no mutations occurring , cf . 
   however , this process is easily generalizable to more genotypes and frequency - dependent selection in the way described above , cf . 
   @xcite .    note that in our model the rate of reproduction of e.g. type @xmath1 is directly determined by the term @xmath15 as e.g. in @xcite . in other models 
   the fitness is first normalized so that the rate of reproduction is given by @xmath16 @xcite , where @xmath17 \label{eq : meanfitness}\ ] ] is 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 ( [ eq : meanfitness ] ) . 
   thus , these models exhibit a quantitatively different time course , but the same event sequences .    until now , 
   usually linear functions have been considered for @xmath18 . however , in many applications nonlinear functions seem more appropriate @xcite . 
   for example , cooperation effects in game theory induce functions linearly increasing in @xmath3 , so that @xmath19 @xcite . 
   this would mean that fitness increases infinitely with the population size of genotype @xmath1 . 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. @xmath20 where @xmath21 and @xmath22 . 
   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 @xmath19 and @xmath23 lead to quadratic relations for the _ mean _ fitness @xmath24 ( see equation ( [ eq : meanfitness ] ) ) of the population which becomes @xmath25 the quadratic factor originates from the fact that the linear fitness function in ( [ eq : meanfitness ] ) is multiplied by the number of individuals leading to a quadratic dependence on @xmath26 , e.g. @xmath27 . 
   this factor enters the replicator equation in evolutionary game theory , thus leading to a quadratic fitness dependence in standard game theoretic problems @xcite . 
   however , this does not imply nonlinear interactions and thus fitness effects such as presented in equation ( [ eq : fitnessfunction ] ) . 
   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 ( [ eq : fitnessfunction ] ) will  as we show in this article  lead to dynamics with more complex stability properties . usually , when analyzing the dynamics of such a model system as described above , it was assumed either that no mutations occur at all ( @xmath28 ) @xcite or that the mutation probabilities are equal ( @xmath29 ) @xcite . however , usually mutation probabilities can be highly diverse @xcite . 
   this has been considered in some examples @xcite , 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 e.g. figure [ fig : asymmetricexample ] ) . 
   in the following we assume the fitness of the individuals to be given by @xmath30 and @xmath31 with the interaction functions @xmath32 and @xmath33 . the model system is effectively one - dimensional with the variable @xmath34 completely describing the system state @xcite . at each event @xmath35 may change by at most @xmath36 depending only on the actual state of the system , so that the dynamics are markovian . 
   let the transition @xmath37 occur with the rate @xmath38 . 
   then this rate is determined by the fitnesses and mutation probabilities in the following way :    * each individual in the population receives offspring with rate @xmath39 with @xmath40 . 
   this means that genotype @xmath1 receives one offspring with rate @xmath41 . 
   * such an offspring increases the population @xmath1 with probability @xmath42 , and population @xmath2 in case it mutates with probability @xmath12 . * one individuum is chosen uniformly to die belonging to genotype @xmath2 with probability @xmath43 . 
   * equivalently the number of individuals of genotype @xmath1 can increase if genotype @xmath2 receives one offspring , which mutates to genotype @xmath1 with probability @xmath44 , and one individuum of @xmath2 is chosen to die . 
   taken these processes together the transition rate for @xmath37 is@xmath45 decreasing @xmath35 to @xmath46 is equivalent to increasing the number of individuals of genotype @xmath2 , so that the transition rate@xmath47 is obtained analogously . 
   the master equation of the markov chain is then given by@xmath48 where we define @xmath49 . 
   as also @xmath50 , the system has reflective boundaries at @xmath51 and @xmath52 . if @xmath53 ( @xmath54 ) then @xmath55 ( @xmath56 ) is an absorbing state of the dynamics 
   it is possible to derive the exact solution to such an equation similarly to the derivations in @xcite 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 @xcite describing the system in the limit of large @xmath0 . 
   as previous works have shown , this approximation already leads to very good results for moderate population sizes of @xmath57 , cf . 
   we use the transformation @xmath58,\qquad s=\frac{t}{n},\qquad{\tilde{\mu}}_{ij}=\mu_{ij}\cdot n,\\[0.1 cm ] \rho(x , s)=p_{xn}(sn)n,\qquad\tilde{g}_{j}(x)=g_{j}(xn)\cdot n \label{eq : rescaledinteraction}\end{aligned}\ ] ] yielding the fokker - planck equation . 
   the scaling factors are derived in appendix a where the derivation of the fokker - planck equation in the scaling limit @xmath59 is shown in detail ( equation ( [ eq : appfpe ] ) ) . 
   this leads to@xmath60x(1-x)-\delta{\tilde{\mu}}\right)\rho(x , s)\big]+\frac{\partial^{2}}{\partial x^{2}}\big[x(1-x)\rho(x , s)\big]\label{eq : fpe}\ ] ] with the normalization condition @xmath61 for all @xmath62 . here 
   @xmath63 and @xmath64 are the rescaled interaction functions defined in ( [ eq : rescaledinteraction ] ) , @xmath65 is the rescaled mean mutation rate and @xmath66 is the rescaled mutation rate difference . 
   these functions have to be rescaled to obtain a non - degenerate limit , the so - called weak selection regime . 
   notice , that the drift term @xmath67x(1-x)-\delta{\tilde{\mu}}\ ] ] contains three different effects . 
   the fitness difference at each point @xmath68 causes a selective drift , the mean mutation rate causes a drift directed to @xmath69 and the mutation rate difference causes a one - directional drift towards one of the boundaries @xmath70 or @xmath71 . 
   the diffusion term@xmath72 reflects the undirected genetic drift . 
   if fitness differences are on a scale of @xmath73 , then mutation , selection and genetic drift all act on the same scale @xcite . 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 @xmath74 
   genetic drift becomes negligible and the dynamics in the bulk of the system become deterministic for large @xmath0 @xcite . 
   however , as the factor @xmath75 vanishes for @xmath76 and @xmath77 the effects of mutations play an important role . clearly , without mutations @xmath70 and @xmath71 are absorbing states , but even in the strong selection regime they become non - absorbing for any positive mutation rate . 
   this is reflected in the form of the drift term @xmath78 , as only the factors reflecting mutations @xmath79 remain non - zero in the limits @xmath76 and @xmath77 .    for a one - dimensional fokker - planck equation ( [ eq : fpe ] ) the stationary distribution is @xcite@xmath80 with the potential@xmath81 which is defined up to a constant irrelevant to our calculations as it is canceled by the normalization@xmath82 thus 
   , we obtain@xmath83dx\\   & = & \ln\left(x(1-x)\right)-{\tilde{\mu}}\ln(x(1-x))+\delta{\tilde{\mu}}\ln\left[\frac{x}{1-x}\right]-\int\left[{\tilde{g}_a}(x)-{\tilde{g}_b}(x)\right]dx\end{aligned}\ ] ] where we used@xmath84.\ ] ] the stationary solution becomes@xmath85dx}\label{eq : generalsolution}\ ] ] where the constant @xmath86 defined in ( [ eq : normalization ] ) has to be calculated numerically for given parameters . 
   this stationary distribution contains contributions from the above mentioned four different effects :    1 . 
   genetic drift is reflected by @xmath87 which diverges for @xmath88 at the boundaries . 
   the mean mutation rate @xmath89 causes the balancing term @xmath90 . 
   3 .   the asymmetry in the mutation probabilities causes the term @xmath91 which diverges at one boundary and vanishes at the other . 
   all frequency - dependent selection effects are contained in the exponential factor @xmath92dx}$ ] which can take various shapes . 
   what are the possible shapes for the stationary distribution @xmath93 in the form given by equation ( [ eq : generalsolution ] ) ? 
   until now , the term @xmath91 representing the asymmetric mutation probabilities has not been considered to our knowledge and the interaction functions @xmath94 and @xmath95 in the selection term were considered to be at most linear in @xmath68 @xcite . 
   we should therefore be interested in the effects of nonlinear interaction functions and asymmetric mutation rates . 
   let us first analyze the dynamics for nonlinear interaction functions @xcite describing the effects of cooperation and limited resources as described by equation ( [ eq : fitnessfunction ] ) , so that@xmath96 already with @xmath97 such interaction functions induce dynamics stochastically switching between three metastable points . 
   an example is shown in figure [ fig : firstexample ] , where the theoretically calculated stationary distribution from equation ( [ eq : generalsolution ] ) is shown together with data from simulations with a population of @xmath98 individuals which is enough to obtain almost perfect fitting ( cf . also figure [ fig : nplot ] ) . 
   the fitness functions ( figure [ fig : firstexample]b ) both show at first an increase on increasing the number of individuals of genotype @xmath1 or @xmath2 from @xmath99 due to cooperation effects and then a strong decrease due to resource competion . 
   as the resulting fitness functions are asymmetric , also the stationary distribution is asymmetric . 
   there is a maximum at @xmath70 due to genetic drift . 
   selection drives the dynamics towards a metastable state at @xmath100 , because for @xmath101 genotype @xmath1 is fitter than @xmath2 thus increasing in frequency and for @xmath102 genotype @xmath1 is less fit than @xmath2 thus decreasing in frequency ( cf . 
   figure [ fig : firstexample]b ) . 
   the maximum of the stationary distribution is thus exactly at the point where @xmath103 . for @xmath104 
   again genotype @xmath1 is fitter than @xmath2 and thus the dynamics are driven towards @xmath71 by selection as well as genetic drift . 
   the mutational force induced by @xmath105 increases the height of the maximum at @xmath100 driving the system away from the maxima at @xmath70 and @xmath71 . 
   so , in this example genetic drift , mutation and selection all significantly influence the dynamics .    ) 
   ( red , solid ) in perfect agreement with data from simulations with @xmath98 ( blue , @xmath106 ) . 
   ( b ) shows the fitness functions of genotype @xmath1 ( blue , solid ) and @xmath2 ( red , dashed ) . 
   interaction parameters are @xmath107 , @xmath108 , @xmath109 and @xmath110 ( see equation ( [ eq : fofk ] ) and ( [ eq : interactions ] ) ) while the mutation rate is @xmath105 ( @xmath97).,width=566 ]    let us now study the influence of the mutation rates in detail . 
   interestingly , for asymmetric mutation rates ( @xmath111 ) the factor @xmath91 always diverges in the interval @xmath112 $ ] . 
   for @xmath113 it diverges for @xmath76 , otherwise for @xmath77 . 
   this can cause the emergence of a maximum of the stationary distribution at @xmath70 ( or resp . 
   @xmath71 ) . figure [ fig : asymmetricexample ] shows an example where due to asymmetric mutation rates a maximum occurs at @xmath71 , when for @xmath97 there is an absolute minimum at @xmath71 as this stable state minimizes the population s mean fitness . 
   this means , that the system dynamics are mutation dominated near @xmath71 . 
   furthermore , the maximum located at @xmath70 for @xmath97 is shifted for @xmath113 to values @xmath114 causing a minimum of the stationary distribution at @xmath70 . 
   thus , in this example both selection and the asymmetry in mutation probabilities are the driving forces of the system s dynamics . 
   ( red , solid ) and @xmath115 ( gray , dashed ) for a system with selective advantage for genotype @xmath2 as shown in ( b ) ( blue , solid : genotype @xmath1 ; red , dashed : genotype @xmath2 ) . solid and dashed lines in ( a ) show the theoretical curves from equation ( [ eq : generalsolution ] ) , crosses the data from simulations with @xmath98 ( @xmath106 : @xmath97 , @xmath116 : @xmath115 ) . 
   interaction parameters are @xmath117 , @xmath118 , @xmath119 and @xmath120 ( see equation ( [ eq : interactions ] ) ) while mean mutation rate is @xmath121.,width=566 ]    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 [ fig : mutplot]a . 
   thus , the dynamics stay close to the edges @xmath70 or @xmath71 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 ( @xmath122 ) . 
   even more , for high mutation rates the differences in the mutation rates @xmath123 can have an important influence on the dynamics . 
   for instance figure [ fig : mutplot]b illustrates a shift of an existing stable point , which completely vanishes for high @xmath124 . 
   we conclude that asymmetric mutation rates can cause the emergence of new stable states ( as in figure [ fig : asymmetricexample ] ) , the disappearance of existing ones ( as in figure [ fig : mutplot ] ) , and shift existing stable states to new positions .    ) for three different symmetric ( @xmath97 ) mutation rates : @xmath121 ( blue , solid ) , @xmath125 ( red , long dashed ) and @xmath126 ( green , short dashed ) . 
   this illustrates that for low mutation rates the dynamics stay close to the system s edges for long times . 
   here , differences in the mutation rates only slightly affect the shape of the stationary distribution . on the other hand , 
   ( b ) shows that for high mutation rates ( @xmath127 ) the dynamics are drawn stronger towards the middle . here 
   , differences in the mutation rates have a stronger impact , as the three curves with @xmath97 ( blue , solid ) , @xmath128 ( red , long dashed ) and @xmath129 ( green , short dashed ) demonstrate . 
   the increasing @xmath124 shifts the existing stable state towards the edge of the system until it vanishes . 
   system parameters were @xmath98 and a symmetric interaction function according to equation ( [ eq : interactions ] ) with one stable state at @xmath130 with @xmath131 , @xmath132.,width=566 ]    as we saw above , the fit of the theoretically calculated stationary distribution to simulation data is almost perfect for weak selection , i.e. the maximal fitness difference satisfies @xmath133 , cf . 
   e.g. @xcite . to quantify the quality of the fit we measured the stationary distribution with the same parameters as above for different population sizes @xmath0 
   this is shown in figure [ fig : nplot ] demonstrating that for increasing @xmath0 the data gets ever closer to the theoretically obtained distribution . 
   for this we define the empirical distribution @xmath134 where @xmath135 is the discrete process defined in equation ( [ eq : meq ] ) . 
   @xmath136 is a time large enough for the process to reach stationarity and depends on the system size , as does the measurement time @xmath137 which is specified in figure [ fig : nplot ] . using this definition we quantify the quality of the fit using the distance measure @xmath138 which takes the mean distance of the measured empirical distribution @xmath139 from the theoretical distribution for every point except @xmath51 and @xmath52 where the density @xmath140 diverges . here , the the theoretical distribution is calculated by integrating the theoretical density over bins of the size @xmath141 around the points @xmath142 . 
   figure [ fig : nplot]b shows that the mean distance @xmath143 decays with increasing @xmath0 . 
   the decay of this quantity is dominated by the slow convergence at the domain boundaries @xmath51 and @xmath52 . 
   therefore , we additionally defined the maximum distance measure @xmath144}\left\{\left|\pi_k-\int_{k / n-1/2n}^{k / n+1/2n}\rho^\ast(x)dx\right|\right\}\label{eq : maxdistance}\ ] ] which gives the maximum distance between the measured probability distribution from the theoretical distribution for all points except @xmath51 and @xmath52 due to the same argument as above . 
   figure [ fig : nplot]b shows that the maximum distance @xmath145 decreases with increasing @xmath0 , however slower than the mean distance @xmath143 due to the divergence of @xmath146 at the boundaries . 
   altogether we conclude , that the theoretical curve fits the data very well already for @xmath147 . 
   however , near the domain boundaries special care has to be taken , as divergences of the theoretical curve lead to larger deviations .    . 
   ( a ) shows the theoretical curve from equation ( [ eq : generalsolution ] ) ( black , solid ) for the same system as in figure [ fig : firstexample ] . 
   data points show empirical distributions as defined in equation ( [ eq : empiricaldensity ] ) obtained in simulations for @xmath148 ( blue , @xmath106 ) , @xmath149 ( red , @xmath150 ) , @xmath151 ( green , @xmath152 ) and @xmath98 ( gray , @xmath116 ) . 
   ( b ) shows the distances @xmath143 and @xmath145 of simulation data and theoretical curve for different @xmath0 for the mean distance measure @xmath143 defined in equation ( [ eq : meandistance ] ) ( blue , @xmath106 ) and the maximum distance measure @xmath145 defined in equation ( [ eq : maxdistance ] ) ( red , @xmath150 ) , demonstrating that the distance decays for increasing @xmath0 . 
   the solid line @xmath153 and the dashed line @xmath154 are added as a guide to the eye for the relation between measured distances and system size @xmath0 . 
   the measured empirical distributions for both ( a ) and ( b ) were obtained by simulating the system dynamics from an initial state drawn from @xmath140 for a mixing time @xmath155 and then recording the density for a time @xmath156.,width=566 ]    furthermore , as the fokker - planck approximation only holds for weak selection we study the quality of the solution obtained throught the fokker - planck equation in dependence of the selection strength . for this 
   we introduce a scaling factor @xmath157 to the interaction functions , so that they become @xmath158 we find that the stationary solution of the fokker - planck equation well approximates the solution of the master equation @xmath159 for selection strengths up to @xmath157 of the order of 1 ( cf . 
   figure [ fig : selectionplot ] ) . 
   we derive the solution ( [ eq : mastersolutionintext ] ) 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 [ fig : selectionplot ] illustrates , for strong selection @xmath160 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 ( [ eq : meandistance ] ) , increases with @xmath157 as a power law both in comparison with the solution from the master equation and with data from simulations , cf . 
   figure [ fig : selectionplot]d . 
   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 ( [ eq : mastersolutionintext ] ) of the master equation yields the exact distribution fitting the data for all selection strengths .    ) 
   obtained from the fokker - planck equation well fits the exact solution ( [ eq : mastersolutionintext ] ) from the master equation for weak selection , but not for strong selection . 
   ( a ) shows the stationary solution from the fokker - planck equation ( red , solid ) for the system with interaction functions defined by equation ( [ eq : scaledinteractions ] ) for very weak selection @xmath161 together with the solution from the master equation ( blue , @xmath106 ) . 
   the distributions of the same system for weak selection @xmath162 and strong selection @xmath163 are shown in ( b ) and ( c ) , respectively . 
   all curves are plotted logarithmically to better allow a comparison of the deviations over all scales . 
   ( c ) shows that the curve from the fokker - planck equation does not fit the exact master solution very well . 
   this is quantified in ( d ) showing the mean distance @xmath143 between fokker - planck and master solution ( blue , @xmath106 ) and between fokker - planck solution and an empirical density obtained from simulations ( red , @xmath150 ) . 
   the system parameters were @xmath98 , @xmath105 , @xmath131 and @xmath132 . 
   the measured empirical distributions were obtained by simulating the system dynamics from an initial state drawn from @xmath140 for a mixing time @xmath164 and then recording the density for a time @xmath165.,width=566 ]    we have shown above that the dynamics of a two - genotype 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 @xmath121 and @xmath97 , so that mutation exactly balances genetic drift . 
   then , the stationary distribution ( [ eq : generalsolution ] ) has a maximum at each point , where @xmath166 $ ] has a maximum . 
   theoretically @xmath94 and @xmath95 can be any function with an arbitrary amount of extreme points in @xmath112 $ ] , so that there is no limit to the stationary distribution s number of maxima , if selection dominates the dynamics . 
   however , for finite @xmath0 the number of possible maxima is naturally limited by @xmath167 . 
   figure [ fig : multiplestability ] shows an example , where we used periodic interaction functions@xmath168 although this is not a realistic interaction function in most applications , it demonstrates what is theoretically possible in the introduced system .    ) 
   ( red , solid ) together with simulation data with @xmath98 ( blue , @xmath106 ) for the interaction functions given in ( [ eq : perfitfunc ] ) . 
   the computed stationary distribution diverges for @xmath76 and @xmath77 while the simulation data remains finite due to the finite number of individuals . 
   ( b ) shows a sample path which exhibits switching between the different maxima of the distribution . 
   parameters are @xmath169 , @xmath170 , @xmath171 and @xmath97.,width=566 ] 
   in this article we analyzed a two - genotype system in a very general setting with ( possibly ) asymmetric mutation probabilities and nonlinear fitness functions @xcite in finite populations . 
   the underlying moran process is a well established model @xcite 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 @xcite  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 which is limited . 
   therefore , due to limited resources a population being too large can not 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 @xmath35 of genotype @xmath1 in the limit of large population sizes @xmath0 . 
   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 @xmath147 . 
   actually , if the system exhibits absorbing states then the fokker - planck method does not work to study the corresponding quasi - stationary distributions . 
   instead , wkb methods are more appropriate to describe the system dynamics as for example in @xcite , 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 ( @xmath172 ) 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 . analyzing 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 due to a higher mutational stability ( see figure [ fig : asymmetricexample ] ) . 
   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 [ fig : firstexample ] ) where three stable points exist , two caused by selection and one by genetic drift . 
   we conclude that frequency - dependent fitness together with asymmetric mutation rates induces 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 ( see figure [ fig : multiplestability ] ) . 
   all in all , we interprete 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 due to complex inter - species interactions . 
   the moran process is a standard tool to gain theoretical insights into experimental data @xcite . 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 multible 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 . 
   @xcite , 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 @xcite . 
   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 . 
   the master equation ( [ eq : meq ] ) may be transformed to a fokker - planck equation in the limit of large @xmath0 @xcite . 
   we introduce the transformation@xmath173 together with the rescaled functions@xmath174 we fix the scaling functions @xmath175 , @xmath176 and @xmath177 such that in the limit @xmath59 all terms in equation ( [ eq : meq ] ) remain finite so that mutation , selection and genetic drift all act on the same scale . 
   we further define@xmath178 and substituting all this into the master equation ( [ eq : meq ] ) , we obtain@xmath179\rho({x_{-}},s)\right.\nonumber \\   &   & \left[(1-\mu_{ba})(1+g_{b}({x_{+}}))(1-{x_{+}}){x_{+}}+\mu_{ab}(1+g_{a}({x_{+}})){x_{+}}^{2}\right]\rho({x_{+}},s)\nonumber \\   &   & -\left[(1-\mu_{ab})(1+g_{a}(x))x(1-x)+\mu_{ba}(1+g_{b}(x))(1-x)^{2}\right.\nonumber \\   &   & + \left.\left.(1-\mu_{ba})(1+g_{b}(x))x(1-x)+\mu_{ab}(1+g_{a}(x))x^{2}\right]\rho(x , s)\right\ } \label{eq : app1}\end{aligned}\ ] ] we choose @xmath180 and @xmath181 so that in the limit @xmath59 the terms stay finite . 
   further we introduce the mean mutation rate @xmath65 and the mutation rate difference @xmath66 . to not overload the notation 
   we drop the time argument @xmath182 of @xmath183 in the following calculation . 
   this leads to @xmath184\right\ } \\   &   & + n\left\ { { \tilde{g}_a}({x_{-}}){x_{-}}(1-{x_{-}})\rho({x_{-}})-{\tilde{g}_a}(x)x(1-x)\rho(x)+{\tilde{g}_b}({x_{+}}){x_{+}}(1-{x_{+}})\rho({x_{+}})-{\tilde{g}_b}(x)x(1-x)\rho(x)\right.\\   &   & + \frac{{\tilde{\mu}}}{2}\left[(1 - 2{x_{-}})\rho({x_{-}})-(1 - 2{x_{+}})\rho({x_{+}})\right]\\   &   & + \delta{\tilde{\mu}}\left[{x_{+}}\rho({x_{+}})-x\rho(x)-(1-{x_{-}})\rho({x_{-}})+(1-x)\rho(x)\right]\\   &   & + \frac{{\tilde{\mu}}}{n}\left[\left({\tilde{g}_a}({x_{+}}){x_{+}}^{2}+{\tilde{g}_b}({x_{+}})({x_{+}}^{2}-{x_{+}})\right)\rho({x_{+}})-\left({\tilde{g}_a}(x)x^{2}+{\tilde{g}_b}(x)(x^{2}-x)\right)\rho(x)\right.\\   &   & + \left.\left({\tilde{g}_a}({x_{-}})({x_{-}}^{2}-{x_{-}})+{\tilde{g}_b}({x_{-}})(1-{x_{-}})^{2}\right)\rho({x_{-}})-\left({\tilde{g}_a}(x)(x^{2}-x)+{\tilde{g}_b}(x)(1-x)^{2}\right)\rho(x)\right]\\   &   & + \frac{\delta{\tilde{\mu}}}{n}\left[\left({\tilde{g}_a}({x_{+}}){x_{+}}^{2}-{\tilde{g}_b}({x_{+}})({x_{+}}^{2}-{x_{+}})\right)\rho({x_{+}})-\left({\tilde{g}_a}(x)x^{2}-{\tilde{g}_b}(x)(x^{2}-x)\right)\rho(x)\right.\\   &   & + \left.\left.\left({\tilde{g}_a}({x_{-}})({x_{-}}^{2}-{x_{-}})-{\tilde{g}_b}({x_{-}})(1-{x_{-}})^{2}\right)\rho({x_{-}})-\left({\tilde{g}_a}(x)(x^{2}-x)-{\tilde{g}_b}(x)(1-x)^{2}\right)\rho(x)\right]\right\ } \end{aligned}\ ] ] in the limit @xmath59 the different terms with @xmath185 in front become second order derivatives with respect to @xmath68 , while the other terms become first order derivatives . 
   the terms which have a @xmath141 factor vanish in the limit @xmath59 and thus the above equation becomes @xmath186x(1-x)\rho(x)\right)-\delta{\tilde{\mu}}\rho(x)\right]+\frac{\partial^{2}}{\partial x^{2}}\left[x(1-x)\rho(x)\right]\label{eq : appfpe}\ ] ] which is the focker - planck - equation of the system . 
   we directly derive the stationary solution @xmath187 of the master equation ( [ eq : meq ] ) using the detailed balance equation @xmath188 which applies to any chain with only nearest neighbour transitions @xcite , cf . 
   also @xcite . 
   thus , using the rewritten balance equation @xmath189 iteratively , we obtain @xmath190 finally , we may use the normalization condition @xmath191 of the stationary distribution to eliminate the factor @xmath192 . 
   we then obtain the exact stationary solution of the master equation @xmath193 which can be evaluated numerically for any transition rates @xmath194 and @xmath195 . for more details on the exact solution of the master equation in the moran process 
   see for example the work by claussen and traulsen @xcite .
summary: 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 frequency - dependent interactions , derived from a game - theoretic 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 non - trivial setting of two genotypes and analyze the co - action of nonlinear frequency dependence with asymmetric mutation rates . 
   we find that their co - action 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 frequency - dependent fitness . 
   * keywords : * population dynamics ; dynamic fitness ; stochastic switching ; multistability