Payoff and fitness

In evolutionary game theory, the fitness of individuals is determined through games, that is, interactions with other individuals. In our case, the game is characterized by a payoff matrix. If only two types, S 1 and S 2 , interact, the payoff matrix is given by

In our model, the payoff is determined from interactions with all other individuals in the population, excluding self interactions. Thus, the payoff of type S 1 is where i is the number of S 1 individuals in the population. Equivalently, we have for the payoff of type S 2 . To avoid the complications of negative fitness, we define the fitness, f j , of type j as an exponential function of its payoff π j , Here w (0≤w<∞) controls the selective differences between players with different payoffs36. If a 11 =a 12 , and a 21 =a 22 , the payoffs are independent of interactions and are only determined by the type of the individuals. This special case corresponds to constant selection, where the fitness of one type does not depend on the frequency of the types in the population. Neutral selection corresponds to w=0 or to a 11 =a 21 =a 12 =a 22 .

In a population of n types, we use a n×n matrix to describe the payoffs. When a mutant appears, an additional column and row are added to the matrix to describe the additional interactions. In the general frequency-dependent case, 2n+1 new payoff matrix entries have to be defined. In contrast, for constant selection, only one new variable is needed to describe the fitness of a mutant. There are several ways to generate these new variables. Suppose that m is the mutant type, j is a resident type, and p is the mutant's parent type. In the simplest case, the payoffs of the mutant m against the resident types j, a mj , and the payoffs of the resident types against the new mutant, a jm , are chosen randomly and independently of the current types from some probability distributions. But it is more natural, if the mutants are similar to their parents by inheriting some aspects of their payoff entries. To this end, we randomly choose the payoff of the mutant against a certain resident from a distribution around its parent's payoff against that resident. Although we can take arbitrary distributions for the payoff entries in our model, we focus on the simplest case, where the payoff entries of the offspring are Gaussian distributed around its parent's payoff entries. In other words, the mean of the new payoff entries a mj for the mutant m is given by the payoff a pj of its parent type p against a type j individual. An equivalent rule holds for the payoff of the resident types against a mutant, the new payoff entries a jm have mean a jp . For the Gaussian distribution, a change in the variance corresponds to a change in the intensity of selection19. Thus, we always set the variance to one.

This approach implies that mutations with selective advantage are favoured, such that the average fitness increases over time. In our case, this would mean that the effective intensity of selection is also increasing, making the system nonstationary. To avoid this effect, we rescale the payoff matrix by dividing it by the largest absolute value of all payoff entries after every mutation and every extinction.

With the full information on the payoffs, we can calculate the fitness of all types in the population based on the payoff matrix. For example, type j obtains the payoff where i k is the number of type k individuals in the population and d is the number of types.

Moran dynamics

The Moran process describes the evolutionary dynamics in a finite population with overlapping generations16,36,37. We start with a homogeneous population with constant size N=1,000 and payoff a 11 =1. In every time step, one individual is chosen randomly in proportion to its fitness, and produces an identical offspring with probability 1−μ or a mutant with probability μ. The offspring then replaces a randomly chosen individual. In nature, mutation rates can range from 10−8 to 10−3 per base, per generation38. Although mutation rates are not affected by local population size, the effect of mutation rates on diversity is directly related to it. To investigate realistic mutation rates in our model, we consider it based on the population size N. As our primary interest is diversity driven by selection rather than diversity driven by mutations, we focus on low mutation rates here. When the mutation rate is high, μ>1/N, the differences between the population dynamics under frequency-dependent and constant selection are less obvious, as the diversity is mainly driven by mutation. In the case of μ=1/N, frequency-dependent selection still leads to higher diversity, compared with constant selection (Supplementary Fig. S1). In either case diversity tends to decrease for strong selection, which becomes more pronounced for higher mutation rates (Fig. 3; Supplementary Fig. S1). The reason is that there are only relatively few coexistence games and mutant types may destabilize them—and the stronger the selection, the faster this occurs.

When a mutation occurs, we generate the payoff matrix according to the method described above. We record the number of individuals of different types in every time step, which gives a straightforward picture of the population dynamics over time. As the system evolves for a long time, we record the number of types in every generation. To avoid dependence on the initial condition, we excluded the data of the first 25,000 generations (see the captions of figures) in the averages. To compare constant fitness and frequency-dependent fitness, we run simulations in both cases, which only differ in the process for generating payoff matrices.

The results under weak selection reflect the usual statistical properties of genetic data samples. The probability of m different alleles present in the population under neutral selection, P(m), can be calculated by Ewens' sampling formula, , where are the unsigned Stirling numbers of the first kind3,39. For a haploid Moran process, as in our case, the parameter θ is Nμ.

Transition probabilities between different coexistence states

Let us consider selection scenarios generated by introducing mutants. Under strong selection and low mutation rates, a population is usually in an equilibrium where different types coexist with each other. The appearance of a new mutant during a phase of coexistence can lead to establishment of a new alliance with the new mutant as an additional type, formation of a new alliance with fewer types (which may include the mutant type or not), replacement of one type from the previous alliance with the new mutant, or extinction of the mutant.

Here we infer the probabilities of these selective consequences under the Moran process. We assume a mutation rate μ≪N−2. In this case, the average time of waiting for a new mutant is much longer than the average time a population needs to reach a new equilibrium after a mutation. We start simulations from a homogeneous population. Mutants show up at random. After a mutant appears, we wait until the population reaches a new equilibrium, and infer whether the diversity has decreased, increased or been maintained. Each state is characterized only by the number of coexisting types. For example, state one represents that the population is homogeneous, and state two represents that there are two types coexisting. We are interested in the probability that a population changes from one state to another. The transition matrix between different states T is obtained by averaging over the evolutionary trajectories. The element t ij denotes the transition probability from i to j coexisting types. For low mutation rates, t ij is very small for j>i+1. The fraction of time that the population spends in each state of diversity is then given by the stationary distribution of the Markov chain determined by the transition matrix T (Fig. 4).

Population size

For fixed selection intensity, the smaller the population size is, the larger the genetic drift. For fixed population size, the smaller the selection intensity is, the larger the genetic drift. Thus, the stochastic effect from the small population size is similar to a smaller intensity of selection. Instead of having two parameters that lead to the same effect, we focussed our discussion on the case of N=1,000 for various intensities of selection. Focussing on variable selection intensities is computationally less costly than varying the population size. For both the Moran process and the Wright–Fisher process, the required CPU time scales with the population size, but not with the intensity of selection. For comparison, we also carried out simulations for N=100, where the same patterns of difference between frequency-dependent selection and constant selection are observed (Supplementary Fig. S2).

Diploid populations

The evolutionary game dynamics for Mendelian populations has been studied in detail in the past40,41,42; the interaction of two alleles at a diploid locus can be interpreted as a special kind of two-player game, which has a symmetric payoff matrix43,44,45. Suppose there are two types of alleles, allele A and allele B. The fitness of a homozygous individual AA is w AA , the fitness of a BB individual is w BB and the fitness of a heterozygous individual AB is w AB . This can be formalized as

When w AA >w AB and w BB >w AB , it corresponds to under-dominance, where heterozygous individuals have a lower fitness than homozygous individuals. When w AA <w AB and w BB <w AB , a condition of over-dominace is described. A diploid population with more than two types of alleles at a single locus can be described by a symmetric n×n matrix, where n is the number of different alleles and matrix element w ij represents the fitness of a diploid individual with genotype ij.

We have simulated the dynamics of such a diploid population under different selection intensities based on the Moran process (Supplementary Fig. S3). In every time step, one allele is replaced, and thus the time for one generation is twice as long as the one in a haploid Moran model. Under the same mutation rate, the diversity of a diploid population (Supplementary Fig. S3a) is higher compared with the frequency-dependent case in a haploid population (Fig. 4b). This is because symmetry of the payoff matrix w ij =w ji favours coexistence of different types. In the simplest case with only two alleles, a coexistence game corresponding to over-dominance has the ordering of payoffs, w AA <w AB and w BB <w AB . Suppose allele A is a random mutant from allele B, and the payoffs of the new genotypes, w BB and w AB , are random variables with mean w BB . Thus, the probability to have a coexistence of these two alleles is 37.5%, which is larger than 25%, the probability to have a coexistence in a two-allele haploid model. In the diploid approach, the fitness of a genotype ij, w ij , is a constant number, and does not change with the composition of the frequencies of different genotypes (but the fitness of an allele is frequency dependent). Hence, this kind of frequency dependence corresponds to constant selection in a haploid population. To introduce frequency dependence on this level leads to serious mathematical intricacies40,43,45,46.

Wright-Fisher dynamics

In the Wright-Fisher Model, all individuals produce a large number of offspring proportional to their fitness. Then, all individuals from the previous generation die, and are replaced by N new individuals sampled at random from the offspring pool. This corresponds to a multinomial sampling of offspring. The expected number of offspring of a certain type, j, in the next generation is proportional to its fitness. Neglecting mutations, the expected number of type j is where, i j and f j are the number of individuals and the fitness of type j. If there is no difference in fitness between types in the population, the expected number of individuals of the different types is constant and the composition of the population will only be changed by random drift. When we consider mutations, the probability that an offspring mutates is μ. On average, there are Nμ new mutants in the population per generation.

We analyse the same quantities in the Wright–Fisher process as above in the Moran process. We see similar patterns in the differences between constant selection and frequency-dependent selection (Supplementary Fig. S4). However, for very strong selection, diversity decreases in our set-up. This can be understood as follows: consider a stable coexistence between two types. If a fluctuation leads the system away from this point, one type has a slight payoff advantage, which causes a large fitness advantage owing to our exponential payoff to fitness mapping. Such a fluctuation can lead to the immediate fixation of one type in the next generation and thus destroy the stable coexistence quickly.

Again, under weak selection and low mutation rates, we recover the diversity given by Ewens' sampling formula. Under neutral selection, random drift in a Moran process is twice as strong as in a Wright–Fisher process47. Thus, we have θ=2Nμ in Ewens' sampling formula for the Wright–Fisher process.