Calculating Evolutionary Dynamics in Structured Populations

Evolution is shaping the world around us. At the core of every evolutionary process is a population of reproducing individuals. The outcome of an evolutionary process depends on population structure. Here we provide a general formula for calculating evolutionary dynamics in a wide class of structured populations. This class includes the recently introduced “games in phenotype space” and “evolutionary set theory.” There can be local interactions for determining the relative fitness of individuals, but we require global updating, which means all individuals compete uniformly for reproduction. We study the competition of two strategies in the context of an evolutionary game and determine which strategy is favored in the limit of weak selection. We derive an intuitive formula for the structure coefficient, σ, and provide a method for efficient numerical calculation.


Introduction
Constant selection implies that the fitness of individuals does not depend on the composition of the population. In general, however, the success of individuals is affected by what others are doing. Then we are in the realm of game theory [1][2][3] or evolutionary game theory [4][5][6][7][8]. The latter is the study of frequency dependent selection; the fitness of individuals is typically assumed to be a linear function of the frequencies of strategies (or phenotypes) in the population. The population is trying to adapt on a dynamic fitness landscape; the changes in the fitness landscape are caused by the population that moves over it [9]. There is also a close relationship between evolutionary game theory and ecology [10]: the success of a species in an ecosystem depends on its own abundance and the abundance of other species.
Evolutionary graph theory is an extension of spatial games, which are normally studied on regular lattices, to general graphs [31][32][33][34]. The graph determines who meets whom and reflects physical structure or social networks. The payoff of individuals is derived from local interactions with their neighbors on the graph. Moreover, individuals compete locally with their neighbors for reproduction. These two processes can also be described by separate graphs [35].
'Games in phenotype space' [36] represent another type of spatial model for evolutionary dynamics, which is motivated by the idea of tag based cooperation [37][38][39]. In addition to behavioral strategies, individuals express other phenotypic features which serve as markers of identification. In one version of the model, individuals interact only with those who carry the same phenotypic marker. This approach can lead to a clustering in phenotype space, which can promote evolution of cooperation [36].
'Evolutionary set theory' represents another type of spatial model [40]. Each individual can belong to several sets. At a particular time, some sets have many members, while others are empty. Individuals interact with others in the same set and thereby derive a payoff. Individuals update their set memberships and strategies by global comparison with others. Successful strategies spawn imitators, and successful sets attract more members. Therefore, the population structure is described by an ever changing, dynamical graph. Evolutionary dynamics in set structured populations can favor cooperators over defectors.
In all three frameworks -evolutionary graph theory, games in phenotype space and evolutionary set theory -the fitness of individuals is a consequence of local interactions. In evolutionary graph theory there is also a local update rule: individuals learn from their neighbors on the graph or compete with nearby individuals for placing offspring. For evolutionary set theory, however, [40] assumes global updating: individuals can learn from all others in the population and adopt their strategies and set memberships. Global updating is also a feature of the model for games in phenotype space [36]. The approach that is presented in this paper requires global updating. Therefore, our result holds for evolutionary set theory and for games in phenotype space, but does not apply to evolutionary graph theory.

Results
Consider a game between two strategies, A and B. If two A players interact, both get payoff a; if A interacts with B, then A gets b and B gets c; if two B players interact, both get d. These interactions are represented by the payoff matrix We consider a population of finite size N. Each individual uses either strategy A or B. In the framework that we investigate here, the population structure specifies how people interact to derive their payoff. It could be that some individuals interact while others do not, or that some individuals interact stronger or more frequently than others. For example, in evolutionary set theory individuals interact with others who are in the same set and two individuals interact as many times as they have sets in common; in games in phenotype space, individuals interact with others who share the same phenotype. Based on these interactions, individuals derive a cumulative payoff, p. The fitness of an individual is given by 1zwp where the parameter w characterizes the intensity of selection. In this paper we consider the limit of weak selection, w?0.
Reproduction is proportional to fitness but subject to mutation. With probability 1{u the offspring adopts the strategy of the parent. With probability u a random strategy is chosen (which is either A or B).
A state of the population contains all information that can affect the payoffs of players. It assigns to each player a strategy (A or B) and a 'location' (in space, phenotype space etc). Thus, one can think of a state as a binary vector which specifies the strategy of each individual, together with a real N|N matrix whose ij-th entry specifies the weight of the interaction of individual i with j. For example, in evolutionary set theory, the ij-th entry of this matrix gives the number of sets i and j have in common [40]. Note that this matrix is not necessarily symmetric: the weight of i's interaction with j might be different from the weight of j's interaction with i. In this paper, whenever we refer to the number of interactions between individuals, we always count them with their weights or multiplicities.
For our proof we assume a finite state space and we study the Markov process defined by gameplay together with the update rule on this state space. The Markov process has a unique stationary distribution defined over all states.
It is shown in [41] that for weak selection, the condition that A is more abundant than B in the stationary distribution of the mutation-selection process described above can be written as Therefore, the crucial condition specifying which strategy is more abundant is a linear inequality in the payoff values, a,b,c,d. The structure coefficient, s, can depend on the population structure, the update rule, the population size and the mutation rate, but not on the payoff values, a,b,c and d. This 'structural dominance' condition (2) holds for a wide variety of population structures and update rules, including games in well mixed populations [12,13], games on graphs [32][33][34], games in phenotype space [36] and games in set structured populations [40]. For a large well-mixed population we obtain s~1. Therefore, the standard risk-dominance type condition, azbwczd, specifies if A is more abundant than B. Spatial structure leads to s values that are greater than 1. The larger s the greater is the deviation from the well mixed population. For very large s strategy A is more abundant than B if awd. Therefore, spatial structure promotes Pareto efficiency over risk dominance [41]. If a spatial model generates sw1 then it is a mechanism for the evolution of cooperation [42].
Here we derive a formula for s that holds for all processes satisfying two conditions: (i) global updating, which means individuals compete uniformly with all others for reproduction and (ii) constant birth or death rate which means the payoff from the game can affect either the birth rate or the death rate but not both.
These assumptions are fulfilled, for example, by games in phenotype space [36] and by games on sets [40]. They do not hold, however, for games on graphs [32]. The first assumption is necessary because our calculation requires that the update rule depends only on fitness, and not on locality. Local update rules are less well-behaved and can even lead to negative values of s. The second assumption insures that the change in the frequency of players is due only to a change in selection. Without this second assumption the conditions would be more complicated.
For each state of the system, let N A be the number of individuals using strategy A; the number of individuals using strategy B is N B~N {N A . Furthermore, let I AA denote the total number of encounters that A individuals have with other A individuals. Note that every AA pair is counted twice because each A individual in the pair has an encounter with another A individual. As specified before, whenever we say 'number of interactions' we count the interactions together with their weights (if such weights occur in the model). Let I AB denote the total number of interactions that an A individual has with B individuals. Our main result is that the structure coefficient, s, can be written as The notation S:T 0 means that the quantity is averaged over all states of the stochastic process under neutral drift, w~0; each term of the average is weighted by the frequency of the corresponding state in the stationary distribution. Intuitively, s captures how much more likely it is, on average, for an individual to play with

Author Summary
At the center of any evolutionary process is a population of reproducing individuals. The structure of this population can greatly affect the outcome of evolution. If the fitness of an individual is determined by its interactions with others, then we are in the world of evolutionary game theory. The population structure specifies who interacts with whom. We derive a simple formula that holds for a wide class of such evolutionary processes. This formula provides an efficient computational method for studying evolutionary dynamics in structured populations.
his own kind rather than with the other kind. An illustration of this formula is shown in Figure 1. This formula suggests a simple numerical algorithm for calculating the s-factor for any spatial process with global updating. We let the process run for a very long time assuming that all individuals have the same fitness. Thus, we simulate mutation and neutral drift on a spatial structure. For each state we evaluate N B , I AA , and I AB . We add up all I AA N B terms to get the numerator in eq (3). We add up all I AB N B terms to get the denominator. The resulting s can be used for any game given by the payoff matrix (1) to determine if strategy A is more frequent than strategy B in the limit of weak selection.
The rigorous proof of eq (3) is given in Appendix A; here we provide an intuition for it. For symmetry reasons, at neutrality, we have the following identities SI AA N B T 0~S I BB N A T 0 and SI AB N B T 0~S I BA N A T 0 . Using these symmetries together with our formula (3), we rewrite condition (2) as Denoting by H XY~IXY =N X the average number of interactions of X individuals with Y individuals, we can further rewrite eq. (4) as Here x A is the frequency of A individuals, p A is the average payoff of an A-individual and p B is the average payoff of a B-individual. These are p A~a H AA zbH AB and p B~c H BA zdH BB . A standard replicator equation for deterministic evolutionary game dynamics of two strategies in a well-mixed population can be written as _ where _ x x A is the time derivative of the change due to selection and p A~a x A zb(1{x A ) and p B~c x A zd(1{x A ) denote the average payoffs for A and B if the frequency of A is x A . This equation describes how selection alone changes the frequency of strategy A over time. Hence, the condition that strategy A is favored by selection is S _ x x A Tw0 where the average is now taken over all states of the mutation-selection process, in the presence of game (w=0). In the limit of weak selection, one can write the first-order Taylor expansion of this inequality to obtain S _ Since at neutrality the average change in the frequency of A is zero, our condition for strategy A to be favored over strategy B becomes S L Lw _ x x A T 0 w0 which is precisely inequality (5). Therefore inequality (5) has a very intuitive interpretation.

Evolution of cooperation
As a particular game we can study the evolution of cooperation. Consider the simplified Prisoner's Dilemma payoff matrix: This means cooperators, C, pay a cost c for others to receive a benefit, b. Defectors, D, pay no cost and distribute no benefits. The game is a Prisoner's Dilemma if bwcw0. As shown in [41], if we use equation (2) we can always write the critical benefit-to-cost ratio as b c provided sw1. If the benefit-to-cost ratio exceeds this critical value, then cooperators are more abundant than defectors in the mutation-selection equilibrium of the stochastic process for weak selection. A higher s corresponds to a lower benefit-to-cost ratio and is thus better for the evolution of cooperation. From eqs (3) and (7) we can write b c This formula is very useful for finding the critical benefit-to-cost ratio numerically. Moreover, we can rewrite the critical benefit-tocost ratio in terms of average number of interactions rather than total number of interactions as b c These equations provide intuitive formulations of the critical benefit-to-cost ratio for processes with global updating.

Computational example: Evolutionary dynamics on sets
Our new formula for s (eq. 3) gives a simple numerical algorithm for calculating this quantity in any spatial process with global updating and constant birth or death rate. We simulate this Figure 1. Calculation of s for a very simple example with population size N~3. Suppose there is a 'spatial' process which has two mixed states. These two states must have the same frequency in the stationary distribution at neutrality, because the process cannot introduce asymmetries between A and B at neutrality. Each mixed state can be described by a weighted, directed graph: in a state with i A players, let p i be the probability that an A plays with another A and let q i be the probability that an A plays with a B. These probabilities are enough since for the calculation of s we only need the AA edges and the AB edges. Note also that the pure states, all-A and all-B, do not contribute to the calculation. We obtain s~p 2 =(2q 1 zq 2 ). doi:10.1371/journal.pcbi.1000615.g001 process under neutral drift for many generations. For each state we evaluate N B , I AA , and I AB . We add up all N B I AA products to get the numerator in eq (3), and then we add up all N B I AB products to get the denominator. The resulting s can be used for any game given by the payoff matrix (1) to determine if strategy A is more frequent than strategy B in the limit of weak selection.
In this section we use the simple numerical algorithm suggested by our formula (3) to find s for evolutionary dynamics on sets [40]. In that paper, the authors compute an exact analytic formula for s that depends on the parameters of their model. We compare our simulated estimates for s with their theoretical values and find perfect agreement (Figure 2). Furthermore, we use our computational method to calculate s in an extension of the original model. An analytic solution for this extended model has not yet been found. Thus our simulated estimates constitute the first ''solution'' of this extended model (Figure 3).
The original set-structured model describes a population of N individuals distributed over M sets. Individuals interact with others who belong to the same set. Two individuals interact as many times as they have sets in common, and these interactions lead to payoffs from a game as described in general in Section 2. Reproductive updating follows a Wright-Fisher process, where N individuals are selected with replacement to seed the next generation. The more fit an individual, the more likely it is to be chosen as a parent. An offspring adopts the parent's strategy with probability 1{u, as described in Section 2. The offspring adopts the parent's set memberships, but this inheritance is also subject to mutation; with probability v, an offspring adopts a random list of set memberships. This updating process can be thought of as imitation-based dynamics where both strategies and set memberships are subject to selection [40].
To obtain exact analytical calculations, it is assumed that each individual belongs to exactly KƒM sets. In Figure 2, we pick values for N,M,K, and u and plot s as a function of the set mutation rate, v. The continuous curves are based on the analytic formula for s derived in [40]. The new numerical algorithm generates the data points. There is perfect agreement between these two methods.
In Figure 3, we consider a variant of this model. Instead of belonging to exactly K sets, individuals now belong to at most K sets. With probability v, an offspring adopts a random list of at most K memberships, the length of which is uniformly random. So far there exists no analytical solution for this model but we can use eq. Figure 2. Agreement of simulations with analytic results. We test our simulation procedure against the analytic results of the set model of [40]. Parameters used are N~100 and M~10. K~1,2 or 3 is the number of sets an individual is in, u is the strategy mutation, and v is the set mutation. We run simulations for 10 7 generations. We use a low strategy mutation (u~0:2) in (A) and a high strategy mutation (u~0:002) in (B). doi:10.1371/journal.pcbi.1000615.g002 (3) to compute s numerically. We interpolate the numerical results with smooth curves. We observe that for low mutation, Fig. 3(A), the case Kƒ3 gives a s which is smaller than the K~3 case. Hence, for low mutation, allowing people to be in at most K sets turns out to be worse for cooperation than restricting them to be in exactly K sets. However, for high strategy mutation, Fig. 3(B), the s for Kƒ3 is greater than the one for K~3. Hence, for high strategy mutation, allowing individuals to be in at most K sets seems to be better for cooperation than restricting them to be in exactly K sets. This suggests that there exists an intermediate strategy mutation rate where the two cases are similar.

Discussion
It has been shown that evolutionary dynamics in a structured population can be described by a single parameter, s, if we are merely interested in the question, which of the two competing strategies, A or B, is more abundant in the limit of weak selection [41]. Payoff matrix (1) describes the interaction between the two strategies A and B and the inequality sazbwczsd specifies that A is more abundant than B in the mutation-selection equilibrium. In general the parameter s can depend on the population structure (which specifies who interacts with whom for accumulating payoff and for evolutionary updating), the population size and the mutation rates; but it does not depend on the entries of the payoff matrix. The s parameter has been explicitly calculated for a number of models including games on graphs, games in phenotype space, games in set structured populations and a simple model of multi-level selection [42].
Here we provide a general formula for the s factor, which holds for the case of global updating. Global updating means that all members of the population compete globally (as opposed to locally) for reproduction. For example, global updating arises in the following way: one individual reproduces and another random individual dies (in order to maintain constant population size); the offspring of the first individual might inherit (up to mutation) the strategy and the 'location' of the parent. Global updating is a feature of models for games in phenotype space [36] and for games on sets [40].
Our main result, eq (3), provides both an intuitive description of what the s factor is and an efficient way for numerical computation.

Materials and Methods
Here we give the proof of equation (3). It is based on the following three claims which we prove in the next subsection: Claim 1. First, we show that for structures and update rules with either constant death rate or constant birth rate the condition for strategy A to be favored over strategy B is equivalent to where birth A and death A are the total birth and death rates of A players and Sbirth A {death A T~P S (birth A {death A ) S p S is the change due to selection averaged over all states of the system, weighted by the probability p S that the system is in each state. The change due to selection in the frequency of A in each state is the difference between the number of A's that are born and the number of A's that die.
Claim 2. We show that for global updating, condition (11) is equivalent to Here S:T 0 denotes the average over the stationary distribution in the neutral process, w~0. Claim 3. Finally we claim that, in the limit of weak selection, for structures satisfying global updating and constant death or birth, the difference between the birth rate and death rate of an individual i in state S can be written in terms of the payoff of individual i as: where p tot is the total payoff of players in the given state S.
Combining the three claims, we conclude that condition (10) is equivalent to Using the weighted number of interactions between players, we can rewrite the total payoffs in any given state as p tot A~a I AA zbI AB p tot~a I AA zbI AB zcI BA zdI BB Thus, condition (14) is equivalent to However, since 1{x A~xB , by symmetry at neutrality we have that SI AA x B T 0~S I BB x A T 0 and SI AB x B T 0~S I BA x A T 0 . Hence (16) is equivalent to sazbwczsd ð17Þ where s~S This concludes the proof of the main result. Below we give the proofs for the three claims made above.

Proofs of Claims
Proof of Claim 1. By assumption, either birth or death has a fixed rate; assume without loss of generality that death is constant with rate d. In a given state, the expected change in the frequency of A individuals is We simplify this equation using the following three relations: birth A zbirth B~d eath A zdeath B since the population size is fixed; death A~d x A and death B~d x B since the death rate is constant and, finally x B~1 {x A . Moreover, we know that on average selection and mutation balance each other, so the average total change in the frequency of A individuals is zero, i.e. SDx A T~0. Using all these into (19) we conclude that This proves the claim. Note that this claim holds for any intensity of selection. Proof of Claim 2. As in [41], we are assuming that the transition probabilities are differentiable functions of w at w~0. Then, in the limit of weak selection, we can write the first-order Taylor expansion of Sbirth A {death A T at w~0 For global updating, the average change due to selection in the neutral process is zero, i.e. Sbirth A {death A T 0~0 . Moreover, using the product rule, we write: Here we used the fact that for neutrality, under global updating in a fixed population size, individuals have equal birth and death rates; hence, (birth A {death A ) S j w~0~0 for all states S. This gives the desired result.
Proof of Claim 3. Again, we assume without loss of generality that the death rate is constant, equal to d. In neutrality, all individuals have effective payoff 1. As noted in the proof of Claim 2, an individual has equal birth and death rates at neutrality, w~0. Thus, in the limit of weak selection, we can write the first-order Taylor expansion at w~0 and obtain birth i~d zw Lbirth i Lw w~0 : ð23Þ When w=0, the birth rate of each individual depends on the effective payoff of any other individual, which itself is a function of w: f j~1 zwp j . Hence (23) can be rewritten using the chain rule as Because the population size is fixed, we have P i birth i~P i death i~d . Hence, summing (24) we obtain When w~0 all individuals have the same fitness. Therefore, by the symmetry imposed by global updating, we have: Lbirth i =Lf i j w~0~L birth j =Lf j j w~0 for all i and j and Lbirth i =Lf j j w~0~L birth k =Lf l j w~0 for all i=j and k=l. It thus follows from (25) that for each j=i Thus, we can rewrite (24) as birth i~d zw N N{1 which gives the desired result.