The fixation probability of a beneficial mutation in a geographically structured population

One of the most fundamental concepts of evolutionary dynamics is the ‘fixation’ probability, i.e. the probability that a gene spreads through the whole population. Most natural communities are geographically structured into habitats exchanging individuals among themselves. The topology of the migration patterns is believed to influence the spread of a new mutant, but no general analytical results were known for its fixation probability. We show that, for large populations, the fixation probability of a beneficial mutation can be evaluated for any migration pattern between local communities. Specifically, we demonstrate that for large populations, in the framework of the Voter model of the Moran model, the fixation probability is always smaller than or, at best, equal to the fixation probability of a non-structured population. In the ‘invasion processes’ version of the Moran model, the fixation probability can exceed that of a non-structured population; our method allows migration patterns to be classified according to their amplification effect. The theoretical tool we have developed in order to perform these computations uses the fixed points of the probability-generating function which are obtained by a system of second-order algebraic equations.


Introduction
One of the main concepts of population genetics is the fixation probability [1]: a beneficial gene does not always spread, but has a probability π f to take over the whole population. The fixation probability can be computed for a non-structured population [2,3], where the progeny of an individual can replace any other individual in the community with equal probability, and depends monotonically on the fitness 1 + s of the new allele and the size of the community N : where c = 1/(1 + s) and one mutant has been introduced into the community. Natural communities, however, are geographically extended and subdivided into colonies that exchange individuals: the progeny of an individual in a given patch compete with the local descendants of other individuals in the same patch and also with migrants from other nearby patches. Evaluating the fixation probability of these subdivided communities is of fundamental importance (see below), but no general results analogous to that of Moran were known for them.
Evolutionary dynamics is a competition between deterministic selection pressure and stochastic events due to random sampling from one generation to the other. For small populations, stochastic events drive the evolutionary dynamics: a large beneficial mutation with fitness, say 1.1, has a probability of only 0.29 to take over a community of size N = 4, barely superior to a neutral mutation (equation (1)). For a large non-structured population, selection is the main force: the same beneficial allele in a population of size N = 10 6 has a fixation probability of 0.09, orders of magnitude higher than the fixation probability of a neutral allele. It is then of fundamental importance to identify the driving forces in a large subdivided population and determine how the fixation probability is affected by this division and the pattern of migration routes.
Evolutionary problems related to geographical subdivision have been investigated from the beginnings of the field of population genetics [4][5][6][7][8][9][10]. Maruyama [11,12] was the first author to cast the problem of fixation probability into a rigorous generalized Moran model (see figure 1) to which he applied standard techniques of branching processes and Markov chains. Under the assumption that the progeny in a patch descend from only local individuals, he obtained an The connectivity between islands is given by the matrix m ki , the probability of migration from island k to island i. In the IP version of this model, one individual is first selected for duplication and then replaces another one.
astonishing result: the migration pattern between the patches has no influence on the fixation probability of a gene and, in this respect, the population can be considered as non-structured. The simplifying assumptions of Maruyama have been analyzed and criticized by many authors (see, e.g., [1,13,14]), but the Maruyama result on fixation probability was, to a great extent, held to be true. Recently, Lieberman et al [15] were able to show that for patches of size 1 (where all the progeny descend from neighboring patches), some peculiar migration patterns can considerably amplify the effective fitness of a beneficial allele for the 'invasion' version of the Moran model. Antal et al [16,17] analyzed a similar problem and showed that the location where a new mutant is introduced can play an important role in its fixation probability. They also pointed out the fundamental difference between two versions of the generalized Moran model, called the invasion process (IP) and the Voter model (VM), in their amplifying capabilities. These new findings have spurred much interest in the field of evolutionary dynamics on graphs [14] and have applications beyond population genetics in domains such as the spread of epidemics [18] or cancer [19,20] or the appearance of collective behavior [21].
In this paper, we give the general solution of the Moran model on a graph (see figure 1) introduced by Maruyama, without his simplifying assumptions. The mathematical tool we have developed enables us to classify migration patterns according to their effect on the fixation probability and to make clear the differences between the two versions of the Moran model: the VM (die first and then be replaced) and the IP (duplicate and then replace another). We demonstrate that when the selection pressure is not too small, in the VM framework, spatial subdivision always diminishes the chances of success of a beneficial allele: the fixation probability of a beneficial mutation in a subdivided population is always smaller than or, at best, equal to that of a non-structured population of the same size. The equality happens only for a 4 'balanced' migration pattern, i.e. when the numbers of migrants sent and received are balanced for each patch.
For the IP, the fixation probability of a new mutant in a subdivided population can exceed the chances of success of the same mutant in a non-structured population. Our method allows migration patterns to be partitioned according to their effect on fixation probability. For large populations there is no absolute upper bound for the fixation probability which can approach unity. However, for a given value of imbalance between migrants sent and received by each island, there exists an upper bound to the fixation probability, which is higher than those of a non-structured population; among the patterns with the same imbalance, some topologies can amplify the fixation probability to this higher value.
The mathematical method we develop is an extension of the tool we recently developed [22] for non-structured populations and is based on the dynamics of the probabilitygenerating function (dPGF). The power of the method we present resides in its simplicity. The fixation probability is uniquely related to the fixed points of a partial differential equation and these fixed points are the roots of a system of algebraic second-order equations. In this framework, the fixation probability acquires a particularly simple geometrical interpretation. We stress that our results are not limited to migration between neighboring patches, but encompass any migration pattern.
This paper is organized as follow. We derive first the general outline of the dPGF method. The next section is dedicated to the application of the dPGF method to the VM version of the Moran process. The following section treats the IP process following the same lines of argument and highlights fundamental differences between these two models. The main points are summarized in the conclusion.

The dynamics of the probability-generating function equation and its fixed points
We have shown recently [22] that the differential equation governing the probability-generating function (PGF) can be used to compute efficiently the behavior of the Moran process for a nonstructured population. We recall this method, which we will generalize to the Moran process on a graph. Consider a community of fixed size N , where n haploid individuals carry the allele A with the fitness 1 + s (compared with 1 for the others). Death events occur continuously, and when an individual dies, it is immediately replaced by the duplicate of another individual, chosen at random among the others, the randomness of this choice being weighted by the parent's fitness. The transition probabilities for the A allele to decrease or increase its population number by one unit during a short time dt read where µ is the death rate. Let us denote by P(n, t) the probability density of observing n A individuals at time t, and set φ(z, t) = n P(n, t)z n as its PGF. Then, in units of generation time (µ/N = 1), φ obeys the following equation, where c = 1/(1 + s). Note that for a beneficial mutation, which we consider in this paper, s > 0 and c < 1. The solution of this equation is given elsewhere [22]. The stationary solution φ s (z) = φ(z, t → ∞), however, can easily be used to compute the fixation probability. We observe from expression (4) that φ s (z) obeys the first-order ODE N φ s − dφ s /dz = K , the solution of which is Note that as φ(1, t) = 1, π 0 + π f = 1; moreover, as π 0 = φ s (0) and π f = (1/N !)φ (N ) s (0), they represent, respectively, the loss of allele and the fixation probability.
The heart of our method lies in the fact that z = c is a fixed point of equation (4); if at time t = 0, there were exactly n 0 A alleles present, for all subsequent times φ(c, t) = φ(c, 0) = c n 0 . Inserting this relation into equation (5) yields which is the celebrated Moran result [2], usually obtained by more complicated means. Consider now the Moran process on a graph: the community is subdivided into M colonies (or islands in Maruyama's terms), the size of the ith island being N i . When an individual dies on the island i, it has a probability m ki of being replaced by the progeny of an individual on the island k. The matrix of m ki , the migration probability from island k to island i, uniquely specifies the migration pattern between islands. Let n i be the number of A on island i. The transition probabilities for A to decrease or increase in number by one on this island is a small modification of expressions (2) and (3): Let us make two remarks before going further. Firstly, there are two possibilities for the Moran model on a graph. In the VM, first an individual dies and then is replaced by the progeny of another individual. This is the case in for example a tropical forest, where space is the limiting resource and a seed can grow only when an adult tree dies [23]. In the IP, used in for example cancer propagation [19], first an individual duplicates and then the progeny replace another individual. The N k in the denominator of equations (7) and (8) must be replaced by N i ; moreover, the normalization constraint is different for these two cases: The difference between the two versions becomes important for unbalanced graphs (see below). Secondly, in order to simplify the notation, we set all islands to the same population size 6 N i = N . All the results obtained here can be trivially generalized to variable island size (see appendix A.4). We denote by P(n 1 , n 2 , . . . , n M ; t) the probability of observing n i individuals on island i at time t, and let φ(z 1 , . . . , z M ; t) be its PGF. It can be shown (see appendix in A.1) that φ obeys the following equation for the VM: where ∂ i = ∂/∂z i . The above equation can be treated along the same lines as before. We first note that the stationary solution for a strongly connected graph reads where again π 0 and π f stand for the loss and fixation probabilities. If we could find a fixed point ζ =(ζ 1 , . . . , ζ M ) of equation (11), then we would immediately deduce the fixation probability We see here how the general problem of finding the fixation probability is reduced to finding the fixed point of a partial differential equation (PDE). A priori, the existence of such a fixed point ζ cannot be taken for granted. We shall see in the next section that a true fixed point actually exists for certain classes of connectivity and that a 'quasi-fixed point' exists in the general case, when the global population is large enough. In this paper, we are interested in the fixation probability of one mutant appearing on an island chosen at random. Therefore the initial condition is and φ(z 1 , . . . , z M ; t = 0) = (1/M) i z i . The fixation probability is then given by We call attention to the versatility of this method. The appearance of a mutant on a specific island, for example, can be treated with the same ease by using (13). Also note that for the initial condition (14), the solution (15) remains valid for weakly connected graphs (see the dominated islands example in appendix A.2).

Exact results
We first concentrate on the VM; the IP will be discussed in the next section and treated as a natural extension. We first investigate for what kind of connectivity matrix the fixation probability is the same as that for a non-structured population. Obviously ζ i = c yields the nonstructured fixation probability. Inserting this specific value into equation (11) and arranging the summation orders, we obtain where is called the temperature of node i by Lieberman et al [15]. T i > 1 represents a node that sends out more migrants than it receives. If all the nodes are balanced, i.e. T i = 1 ∀i, then ζ i = c is indeed a fixed point, and according to equation (15) This generalizes the Maruyama finding to balanced islands. This result has already been obtained for the IP process with a different method by Lieberman et al, who call graphs obeying this condition 'isothermal'. Other particular solutions for simple symmetries such as star connectivity (see figure 2 and appendix A.2) can be readily obtained.

Approximate solutions
Our aim, however, is to evaluate the fixation probability for any given connectivity, i.e. to determine the function π f (m i j ). We have already shown that [22] for the classical Moran model, and for not too small a selection pressure N s 1, second-order terms not containing the factor N in the dPGF (4) can be neglected. Following the same approximation here, we find that there exist quasi-fixed points ζ = (ζ 1, . . . , ζ M ) which obey the following system of second-order algebraic equations: The above equations can be directly solved by standard numerical means. They are open, however, to a geometrical interpretation. Summing the above expression over the index k yields 8 This shows that the quasi-fixed points corresponding to the various connectivities lie on a particular hypersurface. Note that the exact fixed points found previously naturally belong to this hypersurface. An important consequence of the above equation is that we can now sort migration patterns according to their effect on the fixation probability. This also allows us to demonstrate that π f has an upper bound corresponding to π n.s. f . To illustrate these concepts and the geometrical meaning of these equations, we first consider the case of two islands, which readily generalizes to arbitrary M. There are only two independent migration parameters in this case: m 21 = 1 − m 11 and m 12 = 1 − m 22 . Rearranging the terms, the above set of equations can be written as These are the equations of two conics crossing at ζ 1 = ζ 2 = 1 and therefore have one more solution. We can also regard the above equation as a linear system in migration probabilities m i j , which has a solution only if ζ = (ζ 1 , ζ 2 ) belongs to the curve Choosing a point along this curve uniquely determines the fixation probability π f (equation (15)). On the other hand, once the point ζ is chosen, equation (21) determines the level curves of the function π f (m 11 , m 22 ), which are straight lines. This construction is detailed in figure 3.
The important point we wish to stress is that π f is bounded between two extremum values, the maximum corresponding to π n.s. f and the minimum to the neutral case 1/M N . This can be observed directly from the geometrical construction of figure 3(a). Alternatively, this maximum principle can be demonstrated by searching for the constrained extremum of the function π f (ζ 1 , ζ 2 ) (equation (15)) under the constraint of equation (22), using the Lagrange multipliers (see appendix A.3). Figure 4 presents a comparison of the above theoretical computations and direct numerical resolution of the Moran stochastic process for two islands by a Gillespie algorithm [24] (see appendix A.5). As can be observed, there is excellent agreement between the theoretical results and 'true' fixation probabilities obtained by numerical simulations.
Except for the ease of visualization when M = 2, nothing is specific to this case, and all of the above arguments can be extended to an arbitrary number of islands. Every quasi-fixed point ζ = (ζ 1 , . . . , ζ M ) belonging to the hypersurface (20) uniquely determines a fixation probability. In the M(M − 1) space of connectivity coefficients, the level curves of π f (m i j ) are hyper-planes given by equation (19).
The reverse procedure to explicitly compute the function π f (m i j ) is simple: once the m i j are specified, the fixed point ζ is found by solving the system of M second-order algebraic equations. Finding ζ thus allows for the computation of π f through the algebraic expression (15). Figure 5 presents a comparison of theoretical and numerical fixation probability, for five islands of population size N = 50 and for various selection pressures; for each selection pressure, 1200 random graphs have been generated. The relative deviation between theory and simulation is ∼10 −3 , which is due to the precision of numerical simulation and the number of stochastic trajectories used for evaluating the fixation probability.
The deviation of π f from its maximum π n.s. f is related to the deviation from 'isothermy'. The 'nonisothermy' of a graph can be measured by = Var(T i − 1). Figure 5(c) shows that even though cannot be used as a single parameter to evaluate the fixation probability (a) (b) Figure 3. Fixation probability as a function of connectivities for two islands.
(or equivalently the deviation of π f from the maximum value), these two quantities are highly correlated.
The fixed point approximation developed above was derived for large selection pressure, N s 1. The domain of validity of our method is, however, much broader. In fact, it is not the population size of an island N , but the population size of the whole community M N that sets the limit of the validity: M N s 1, as can be seen in figure 6. It is not the fixation probability, but the limit of validity of our approximation that seems to follow the Maruyama conjecture that 'certain quantities are independent of the geographical structure of population' [12] . The limit M N s 1 is in fact what makes our method relevant: natural populations are always large and spatially subdivided. Therefore even small additive fitnesses are amenable to our analytical treatment.

Fixed points of the invasion process
For a non-structured population, there is no difference between VM (first die and then be replaced) and IP (first duplicate and then replace) versions of the Moran model. For a subdivided population however, this difference has profound consequences when the islands are not balanced [16]. We can understand this intuitively by considering an island that sends more The deviation from isothermal = Var(T i − 1) versus the fixation probability π num f of graphs. There is a high positive correlation (C ≈ 0.9) between and π f − π n.s. f migrants than it receives. In the IP version, when one A (beneficial) allele duplicates on this island, it has a high probability of spreading to other islands. In the VM version, as death occurs first, this advantage disappears. It is, however, less obvious why the IP advantage for one island could persist upon averaging over all islands.  Mathematically, the treatment of the IP is similar to VM, except that migration coefficients should obey row normalization i m ki = 1 instead of column normalization (equations (9) and (10)). The dPGF equation is therefore slightly different from expression (11): but as before, we can find the fixation probability of a graph by the fixed point of the dPGF through the same equations (13) and (15). For example, the isothermal theorem holds for IP and is demonstrated along the same lines as for the VM: for the point ζ = (c, c, . . . , c), the second-order terms of the dPGF (11) where the temperature is now T i = k m ki . (T i − 1) has the meaning of a balance between migrants received and sent by island i. All other exact results such as the 'star' graph treated in the preceding section have an analogue for the IP (see appendix A.2). In the limit M N s 1, second-order terms of the dPGF can be neglected and quasi-fixed points ζ = (ζ 1, . . . , ζ M ) now obey the following system of second-order algebraic equations: The process of estimating the fixation probability is similar to the previous discussion of the VM version. Figure 7 presents, for 500 random graphs and various selection pressures, a comparison between the fixation probability found from expression (24) and that found by direct numerical simulations. The same high degree of precision is observed as in the previous case. It can be seen in figure 7 that many migration patterns show fixation probabilities higher than the non-structured value. The difference in the evaluation of the fixed point of the invasion process is the presence of the term T k on the right-hand side of equations (24) compared to equations (19). This modification excludes a constraint independent of the graph topology similar to equation (20), and the fixed points no longer belong to a universal hypersurface. An important consequence is that for the IP, the fixation probability does not obey a global 'maximum principle'. We can, however, give a restricted version of this principle. Each migration pattern {m ki } can be considered as a point in an M(M − 1) space. We can partition this space into subspaces E T of dimension (M − 1), each subspace labeled by the vector T = (T 1 , . . . , T M ): the subspace E T contains all migration patterns which have temperature distribution T. Using the same arguments as in the previous section, we can show that in each subspace E T , the fixed points belong to the (M − 1) hypersurface On each hypersurface, the fixation probability has a maximum value π max f (T) > π n.s. f . For a given temperature distribution T, the set of all graphs with fixation probability ∈ [π max f (T), π n.s. f ] are fitness amplifiers. For some particular classes of graph, π f can greatly exceed the nonstructured value and approach unity, as noted was by other authors [15,18] (see appendix A.2).

Conclusion
The evolutionary dynamics of a non-structured population is the interplay between deterministic selection pressure and stochastic sampling between generations. In a geographically subdivided population, migration between colonies is the third ingredient to be taken into account. A migration pattern is a more complex quantity than a number such as the fitness, and it has to be specified by an M × M matrix, where the coefficient m ki weights migration importance from patch k to patch i. In the preceding sections, in the framework of the Moran model, we have developed a mathematical method that allows for the computation of the fixation probability of a beneficial allele in a subdivided population for arbitrary migration patterns.
This mathematical method allows us to partition the migration patterns according to their effect on the fixation probability and gives rise to important fundamental results. For example, the celebrated result of Maruyama [11,12] was that the fixation probability does not depend on population subdivision and is equal to that of a non-structured population π n.s. f . In order to compute the fixation probability, Maruyama used a severe approximation. Our investigation shows (in the VM version of the Moran model) that the Maruyama solution corresponds to the upper bound of the fixation probabilities and the effective fitness of a beneficial allele is degraded by the population subdivision.
The IP version of the Moran model is more often used for epidemic or cancer propagation. In this framework, Lieberman et al [15] were able to build particular migration patterns capable of amplifying the fixation probability of a mutant in a subdivided population compared to a nonstructured one. The method we develop in this paper gives the general conditions under which a pattern can amplify the fitness of a mutant. We show that the important parameter in this case is the (vector) of migration imbalance T, that sets the upper bound for the fixation probability π max f (T). One can then find all topologies corresponding to this imbalance which have a fixation probability in the range of [π n.s.
f , π max f (T)] and hence are fitness amplifiers. The mathematical tool we have developed is based on finding the fixed points of a partial differential equation; for large selection pressure, M N s 1, the fixed points are solutions of a system of second degree algebraic equations. The problem of finding the fixation probability is therefore mapped into a simple problem of geometry and of conics crossing in an Mdimensional space. To date, no general analytical tool has been available for estimating the fixation probability of a subdivided population. Some solutions have been found for simple symmetries, but it was necessary to resort to numerical simulations in order to investigate a particular migration pattern [14]. Although faster algorithms have been developed [25], the numerical simulations remain extremely costly for large populations. The method we have developed in this paper requires negligible computational time and allows the precise investigation of wide classes of topologies and problems. For example, one could envisage problems in which the migration patterns themselves, i.e. the coefficients m ki , are subject to selection in order to prevent epidemic propagation or help the emergence of cooperative behavior.
All species are distributed through space and are therefore subdivided into colonies. Understanding the evolutionary process for these populations is therefore of fundamental importance in our comprehension of evolution in general. We believe that the dPGF method we propose here is a valuable tool in advancing our comprehension of this process. 16 where the symbol · · · stands for summation over states n. The computation rules are simple enough: which result in equations (11) and (23) depending on the choice of row or column normalization of m ki . Note that the stationary solution (12) of the dPGF is valid for connected graphs.

A.2. Some particular solutions of the Moran process on a graph
This paper contains the general solution of fixation probability for the Moran process. Some graphs have simple symmetries and their π f can be given in closed form. Here we explore a few such cases. These graphs, which have already been solved by other authors (at least for IP), have been selected to illustrate the simplicity of the fixed point method.

A.2.1. Dominated islands.
A trivial extension of the isothermal case is the 'dominated islands', where a group of islands can only receive migrants from another group. This is the case, for example, along a river, where downstream patches do not send migrants upstream. The point ζ = (1, . . . , 1) is a trivial fixed point of the dPGF (11) due to the normalization constraint on probabilities n P(n; t) = 1. We saw that for the isothermal graph, ζ = (c, . . . , c) is also a fixed point. We could then explore under which connectivity the point ζ = (c, . . . , c, 1, . . . , 1), where J components are c and M − J components 1, constitutes a solution. Following the same methodology as that for the isothermal case, we obtain the condition As m ki 0, condition (A.4) is equivalent to m ki = 0, for i J and k > J.
This condition partitions the islands into two groups of J and M − J islands when there is no migration from the second group toward the first one; moreover, the islands in the first group are isothermal among themselves. According to (15), the fixation probability reads where π n.s. f is the non-structured fixation probability (18) A.2.2. Star connectivity. Consider figure 2(b), where one central island, labeled 0, is connected to P peripheral islands and where self-migration is prohibited, m ii = 0. Peripheral islands do not communicate between each other, m i j = 0 if i, j > 0. For the VM, this configuration imposes m i0 = 1/P and m 0i = 1. The symmetric disposition of the islands imposes a fixed point of the form ζ = (ζ 0 , ζ 1 , ζ 1 , . . . , ζ 1 ). Plugging this solution into the VM dPGF (11), we obtain ζ 0 = c (Pc + 1) / (P + c) , For the IP, migrations are m i0 = 1 and m 0i = 1/P (i > 0). Using the dPGF (23), we obtain ζ 0 = c (P + c) / (Pc + 1) , We emphasize that these are exact solutions of equations (11) and (23), respectively. In both cases, the fixation probability is When P 1 for IP, this reduces to A.2.3. Superstar. This configuration (figure 2(c)) is an extension of the star configuration, where a central island A is connected to P secondary islands B; each secondary island is connected to Q tertiary islands C. For the VM, the migration coefficients are m AC = 1, m C B = 1/Q and m B A = 1/P, and all other coefficients are zero. There is no longer an exact fixed point for this case; the quasi-fixed points are For the IP, m AC = 1/P Q, m C B = 1, m B A = 1, and the ζ i are the same as the above expressions up to a circular permutation of indices A, B and C. In both cases, the fixation probability is For IP, in the limit P 1, Q 1, this reduces to Note that expressions (A.5) and (A.6) for the star and superstar are identical to those found by Lieberman et al [15]. f . For non-isothermal distributions, α T > α 1 . Points of C T between the two curves π f = α 1 and π f = α T correspond to fitness amplifier graphs.
This 'maximum principle' holds as long as our approximate expression of the fixation probability is valid, that is, for N Ms > 1.
For the IP case, the maximum principle has a restricted version. Consider all graphs corresponding to a given (IP) temperature distribution T = {T 1 , . . . , T M }. Then the quasi-fixed point belongs to the hypersurface Given T, all the above arguments for the maximum can be repeated, except that the extrema of π f are now given by the crossing of C T and the curve T k 1 − c/ζ 2 k = constant, k = 1, . . . , M. As a consequence, the maximum of the fixation probability π max f (T) > π n.s. f . Figure

A.4. Generalization to islands of variable size
The generalization to islands of variable size is straightforward. We describe this generalization for VM, and a similar generalization can be provided for the invasion process. The dPGF reads The temperature of an island is defined as

20
The isothermal condition always corresponds to T k = 1 for all the islands and it corresponds to the fixed point ζ k = c ∀k.
In the case of the appearance of a mutant at a random site, the fixation probability is related to the coordinates ζ i of the fixed point by One can readily show that the 'maximum principle' for VM also holds in this more general case. We setâ Using the Lagrange multipliers method it is easy to show thatĝ is never larger thanâ and that the range of variation ofâ under the constraint (A.11) is c â 1. Following the same reasoning as in section A.4, one then deduces The fixation probability given by (A.10) is maximum for a non-structured population.

A.5. Numerical resolution
A.5.1. Numerical simulation. The stochastic equations given by the rates (7) and (8) can be seen as 2M chemical reactions (two for each island) for the species A: which we solve by using the classical Gillespie algorithm [24] written in C++. Here, we are interested only in the fixation probability and not in the fixation time; the program can therefore be accelerated by computing only the nature of the event that occurs at each turn (and not its time of occurrence). In general, to solve for the fixation probability of a given graph, R =