Abstract
Population structure affects the relative influence of selection and drift on the change in allele frequencies. Several models have been proposed recently, using diffusion approximations to calculate fixation probabilities, fixation times, and equilibrium properties of subdivided populations. We propose here a simple method to construct diffusion approximations in structured populations; it relies on general expressions for the expectation and variance in allele frequency change over one generation, in terms of partial derivatives of a “fitness function” and probabilities of genetic identity evaluated in a neutral model. In the limit of a very large number of demes, these probabilities can be expressed as functions of average allele frequencies in the metapopulation, provided that coalescence occurs on two different timescales, which is the case in the island model. We then use the method to derive expressions for the probability of fixation of new mutations, as a function of their dominance coefficient, the rate of partial selfing, and the rate of deme extinction. We obtain more precise approximations than those derived by recent work, in particular (but not only) when deme sizes are small. Comparisons with simulations show that the method gives good results as long as migration is stronger than selection.
THE geographical distribution of a population has many consequences on its evolution. The nature of the habitat (fragmented or continuous) and the dispersal patterns of individuals across space affect the relative importance of the different processes changing allele frequencies (selection, genetic drift...) and therefore have an effect on several genetical features of populations, such as the amount of genetic diversity and the strength of inbreeding depression. Spatial structure may also affect the timescale of evolution; indeed, the rate of adaptation of a population depends critically on the probability that new advantageous mutations will go to fixation and on the time length of the fixation process.
Different models have been used to calculate the probability of fixation of an advantageous mutation in a spatially structured population, dealing with different types of population structure, and different modes of selection. For some simple models, Maruyama (1970) found that in a population subdivided into demes of constant size that never go extinct, spatial structure has no influence on fixation probabilities: a new mutation has the same probability of fixation as if the population was not fragmented. This conclusion has been supported by more recent work (Cherry and Wakeley 2003). Slatkin (1981) used a lowmigration limit to investigate the joint effects of dominance and spatial structure on fixation probabilities of advantageous mutations. In his model, the fixation of an allele within a deme occurs much faster than its transmission between demes, so that fixation is assumed to take place independently in each deme (see also Takahata 1991). He showed that the probability of fixation of recessive mutations is increased in this lowmigration limit, compared to the case of a panmictic population; indeed, population structure increases the probability that an advantageous homozygous genotype is created before the mutation is lost. However, the probability of fixation of dominant mutations is decreased. Barton (1993) worked on another model to study the effect of catastrophic extinctions of demes and showed that extinctions can substantially reduce the fixation probability of beneficial alleles. More recently, several authors used diffusion approximations to derive fixation probabilities and fixation times in metapopulations. Cherry (2003a,b), Cherry and Wakeley (2003), and Wakeley (2003) have shown that diffusion methods can be used in island models of population structure and studied the joint effects of dominance, population structure, and local extinctions on fixation probabilities and times of mutant alleles. Whitlock (2003) has presented another approach, which he also applies to the steppingstone model of population structure.
In this article we provide a different method to construct diffusion approximations, which is in several ways more general than previous ones. The method consists of expressing the mean of allele frequency change over one generation in terms of partial derivatives of a “fitness function” and probabilities of genetic identity computed in a neutral model. We make no assumption on the particular form of the fitness function; indeed, the method can be applied to many different selection models (including multilevel selection). We then take the example, considered by previous authors, of selection through effects on fecundity (different genotypes producing different numbers of gametes) and consider the effects of dominance, selffertilization, and local extinctions on the probability of fixation of new mutations. We confirm some previous approximations; however, our method gives more precise results when migration is strong, when deme size is small, and also when extinctions occur. As discussed later, the diffusion method gives satisfying results as long as migration rates are not too small; with dominance, Slatkin's method performs better when migration is very low.
LIFE CYCLE
Throughout this article, we use an island model of population structure, with possible local extinctions. We assume that the population is subdivided into n demes. At the beginning of the life cycle, each deme contains N adult individuals and has a probability e of becoming extinct. In the surviving demes, individuals produce a very large number of gametes and die; in practice we treat the number of gametes produced as if it were infinite. Gamete fusion may precede or follow migration (zygotic or gametic migration); we first do not make any assumption about the breeding system, but later we consider the cases of random union of gametes and of partial selffertilization at a rate α. Migrants produced by a deme can reach any other deme with the same probability. The migration rate m is defined as the proportion of migrant zygotes (or gametes) present in any surviving deme once migration has occurred (“backward” migration rate). We assume that migrants arriving in extinct demes do not survive, each extinct deme being recolonized later by a finite number k of individuals (Slatkin 1977). These recolonizers come from “propagules” that formed in each surviving deme after the first round of dispersal. Again, recolonization can be zygotic (if it occurs after fertilization) or gametic (if it occurs before; in this case each extinct deme is recolonized by 2k gametes). We also define a parameter φ measuring the probability that two colonizers come from the same deme (this is specified later). We assume that colonizers reproduce immediately, producing a very large number of zygotes. Finally, in all demes N individuals are sampled randomly among all the zygotes to form the new adult generation; exactly N adults are always in each deme, recolonized or not. Thus, recolonizers experience two rounds of migration and reproduction within one “generation.” This way of representing the dynamics of recolonization is obviously simplistic, but may be seen as a first step in understanding the consequences of extinction and recolonization (Slatkin 1977; Wade and McCauley 1988; Whitlock and Barton 1997; Pannell and Charlesworth 1999); more complex models may be tractable within the framework presented here (Rousset 2004, Chaps. 10 and 11).
THE DIFFUSION APPROXIMATION
Diffusion theory has provided good approximations for probabilities of fixation of mutant alleles (Kimura 1964). In this section we briefly review the methodology in the case of a panmictic population containing two alleles a and A. We call p the frequency of A and δp the change in frequency of A over a short time interval δt. We first suppose that time and frequency are continuous and define
If these limits are finite, and if the higherorder moments of δp are negligible, the probability of fixation of A when present in frequency p, u(p), is given by the differential equation
To apply this continuous model to the case of a population with discrete generations, one states that time is measured in units of N generations, where N is the population size, and that the short time interval δt corresponds to one generation. Then δt = 1/N, and
In the simple case of genic selection where A has a selective advantage s over a in an ideal haploid WrightFisher population of size N, we have
In the following, we consider an island model of population structure and look for an expression for the probability of fixation of a mutant allele A in a population initially fixed for allele a; p is now the frequency of A in the whole population, while the vector p = (p_{1},..., p_{n}) gives the frequencies of A in different demes. We suppose that the number of demes (n) is large and define
EXPECTED CHANGE IN ALLELE FREQUENCY
We use the direct fitness method (Taylor and Frank 1996; Rousset and Billiard 2000) to express M_{δ}_{p} as a function of probabilities of identity and partial derivatives of a fitness function. The parameters and variables of the model are summarized in Tables 1 and 2. In all the following, the subscript ij refers to the jth adult individual of the ith deme. The fitness W_{ij} of this individual is defined as the expected number of its gametes that will participate to the next adult generation. Fitness may depend on the phenotype of the individual, on the mean phenotype of the individuals present in its deme, and on the mean phenotype in the whole population. We call z_{aa} the phenotype of homozygous aa individuals and assume that the phenotype of AA individuals differs from z_{aa} by an amount s: z_{AA} = z_{aa} + s. We write the phenotype of heterozygous individuals z_{Aa} = z_{aa} + sh, h measuring the dominance of A. Defining p_{ij} as the frequency of A in individual ij (p_{ij} = 0, ½, or 1), p_{ij}_{1} and p_{ij}_{2} as the frequencies of A at each of its two homologous genes, and z_{ij} as its phenotype, we have the relation
The expected change in frequency of A over one generation is
The different averages in Equation 14 can be seen as probabilities of identity: for example,
The same reasoning is used to calculate the limit of
Denote
Last, we need the limit as n tends to infinity of
Because adult population size remains constant after regulation, the partial derivatives of the fitness function sum to zero (Rousset and Billiard 2000); therefore we can express M_{δ}_{p} as a function of the two first partial derivatives only. After replacing the average products of allele frequencies by the expressions calculated above, we obtain
VARIANCE IN ALLELE FREQUENCY CHANGE
To use the diffusion method, we also need the limit as n tends to infinity of n × V_{δ}_{p}, where V_{δ}_{p} = E [(δp)^{2}p] and δp is the change in frequency of A over one generation. We have
When e = 0 (no extinction), covariances in the change in allele frequencies in different demes equal zero (because sampling occurs independently in each deme). Equations A1, A2, A3, A4, A5 of appendix a lead to
Extinctions add an additional complication because they introduce correlations in the change in allele frequency in different demes (because the frequency of A in the migrant pool depends on which demes did survive). We show in appendix a that when e > 0, the variance effective population size is
PROBABILITIES OF COALESCENCE
The expressions for M_{δ}_{p} and V_{δ}_{p} derived in the previous sections involve variables representing various probabilities of coalescence in a neutral model, evaluated in the limit as n tends to infinity. These variables can be depend on the life cycle considered) and calculating values at equilibrium. These recurrence equations are given in appendix b for three different cases: (i) monoecious individuals with random fertilization, gametic migration, and no extinction; (ii) monoecious individuals with rate α of selffertilization, zygotic migration, and no extinction; and (iii) monoecious individuals with random fertilization, gametic migration, extinction rate e, and gametic recolonization. Here we give only approximate solutions assuming that migration rate is small and deme size large, so that N = O(1/m), and neglecting terms in O(m).
Case i: random fertilization, no extinction: Approximating Equations B5, B6, B7, B8, B9 for small m and large N, one obtains (with M = Nm)
Case ii: partial selffertilization, no extinction: Here we assume that selffertilization occurs at a rate α and that fertilization precedes migration. Random fertilization corresponds to α= 1/N. We obtain the approximations
Case iii: random fertilization, extinctions: We assume here that both migration and recolonization are gametic. In the life cycle section, we introduced the parameters k (number of colonizers) and φ (probability of common origin). More precisely, φ measures the probability that two genes sampled randomly among the 2k gametes colonizing an extinct deme were in the same deme before migration. We also need the parameters ψ_{1} and ψ_{2}, defined as the probabilities that three genes sampled randomly among the 2k colonizers were in the same deme before migration (ψ_{1}) or in two different demes (ψ_{2}). Following Slatkin (1977), we distinguish two modes of recolonization: the “migrant pool,” where the 2k gametes recolonizing a deme are taken randomly in the whole population after migration (in this case φ =ψ_{1} =ψ_{2} = 0), and the “propagule pool,” where the 2k gametes were all in the same deme after migration [in this case φ = (1 – m)^{2}, ψ_{1} = (1 – m)^{3}, and ψ_{2} = 3m(1 – m)^{2}]. When the extinction rate e is O(m) we obtain
EXAMPLES OF FITNESS FUNCTIONS
We consider two cases of individual selection: “soft” and “hard” selection. The first case corresponds to a regulation of the number of gametes (or zygotes) produced by the different surviving demes just before migration, while in the second case there is no such regulation. However, in both cases there is a phase of population regulation, as N adults are sampled in each deme at the end of every generation. This definition of hard selection is the same as in Barton (1993), but different from in Whitlock (2002), which explains some discrepancies between our results and results in Whitlock (2003). In Whitlock's hard selection model, each deme contributes to the next adult generation in proportion to the number of juveniles produced in the deme before migration; although this definition sounds simple, we do not know any realistic life cycle maintaining a constant population size in this case.
In both cases, we fix z_{aa} to zero and assume that an individual of phenotype z_{ij} produces a number of gametes proportional to 1 + z_{ij}; the numbers of gametes produced by Aa and AA individuals, relative to aa individuals, are then 1 + sh and 1 + s, respectively.
We note at this point that we need only the limit as n tends to infinity of the fitness function, as terms of order 1/n and higher will vanish when we will calculate
Soft selection: In this case, we assume that the number of gametes produced in each deme is regulated before migration to a constant value that is the same for all demes; we suppose that this value is still very large. Recall that the fitness of individual ij is the expected number of its gametes that will participate to the next adult generation. We obtain
Hard selection: Here we assume that the number of gametes or juveniles produced in a deme before migration is not regulated: a deme fixed for A produces 1 + s times the number of gametes (or juveniles) produced in a deme fixed for a. Here the fitness function is given by
RESULTS
Additivity: In the additive case (h = ½) and when no extinction occurs (e = 0), simple results are obtained. Under hard selection and gametic migration one arrives
Therefore, in the additive case with no extinction, probabilities of fixation depend neither on population structure nor on inbreeding, as found by Maruyama (1970), except for effects of order 1/N that have been ignored in previous analyses. Previous simulation studies have also missed them, often because large values of N were assumed. When extinctions occur (e > 0), the expression of G is more complex and depends on population structure. For example, under soft selection and haploid migration, we find
Effects of dominance: When h ≠ ½, the function G is a complicated function of the parameters. Here we present some results relative to case i (monoecious individuals with random fertilization, gametic migration, and no extinction). If we assume that m is small and N large, so that we can use approximations 30 and 31, we obtain
Figure 1 shows the probability of fixation of a rare advantageous allele as a function of its dominance coefficient, for different values of the migration rate. Here we used the exact expressions of the probabilities of identity for case i (given in appendix b); however, Equation 43 gives very close results (not shown). We used the soft selection model, but under hard selection the results are qualitatively and quantitatively very similar for large N. Figure 1 illustrates the fact that population structure increases the probability of fixation of recessive advantageous mutants, while it decreases the probability of fixation of dominant advantageous mutants. It also shows that the diffusion approximation gives poor results for small values of m; in the last case Slatkin's (1981) lowmigration approximation performs better (this is discussed later).
Effects of selfing: In case ii (monoecious individuals with rate α of selffertilization, zygotic migration, and no extinction), approximations (32) and (33) lead to
Figure 2 shows the effect of partial selffertilization on the probability of fixation of a recessive advantageous mutation (h = 0). Again, we used the exact expressions of the probabilities of identity and integrated numerically the G function, but Equation 45 gives very close results. We still observe a discrepancy between analytical and simulation results at small values of m.
Effects of the extinction rate: Figure 3 shows some results concerning case iii (monoecious individuals with random fertilization, gametic migration, extinction rate e, and gametic recolonization), taking the example of an advantageous underdominant mutant. Underdominance corresponds to the case where heterozygotes are disadvantaged (h < 0); it can result from chromosomal rearrangements (Lande 1985) and possibly from mutations affecting transcriptional regulation during the development (Gibson 1996). In a large panmictic population, such mutations have a very low probability of fixation (e.g., Kimura 1962). Using a lowmigration limit, Lande (1985) showed that local extinctions can increase the probability of fixation of underdominant mutations, in an island model of population structure. This effect has also been observed by simulation in a twodimensional steppingstone model (Michalakis and Olivieri 1993).
Figure 3 represents the probability of fixation of an underdominant mutation, relative to a neutral mutation, as a function of the extinction rate e, and for different values of k, the number of colonizers. Figure 3A corresponds to the migrant pool model of colonization, while Figure 3B corresponds to the propagule pool model. We used the soft selection model, but again results are very similar with hard selection (not shown). The effect of the rate of extinction on the fixation probability of the mutant is nonmonotonic. This may result from two opposite forces: extinctions increase the probability that an advantageous AA genotype will appear before allele A is lost, in particular when the number of colonizers is small, but also decrease the chance of fixation of A once AA genotypes have been created (by increasing random drift, extinctions reduce the probability of fixation of advantageous alleles).
DISCUSSION
Figures 1 and 2 show that the analytical model gives good predictions when the migration rate is not too small, while for small values of m predictions are not as good (except when h = ½ and e = 0, since in that case we obtain an expression independent on m for the probability of fixation). This comes from the fact that we calculated the probabilities of coalescence before migration (r_{0}, r_{1}, a, b, c) in the neutral case (s = 0); as s increases, these probabilities become more and more different from the neutral values, and this difference depends on the migration rate. As previously noted, what matters in the argument based on a separation of timescales is that changes in allele frequencies within demes, due to drift and migration, are fast relative to selection (see also Wakeley 2003), hence that s is small relative to m + 1/N. In practice, the validity of the diffusion approximation depends on weak selection (small s) and nottoosmall migration (roughly, m > s). This can also be expressed in terms of our coalescent argument: for relatively high values of m, two lineages from the same deme either separate to different demes (in the past) or coalesce rapidly within the deme, and during this short time selection does not have much effect. For small values of the migration rate, however, the mean withindeme coalescence time becomes longer, and the effect of selection during that time cannot be neglected any more.
Although the validity of the diffusion approximation depends on a large number of demes (large n) and weak selection (small s), these two conditions are in fact not very restrictive: Table 3 shows that our solution gives accurate results for values of n as low as 10 and values of s as high as 0.1. When the exact expressions of the probabilities of identity are used, no other assumption needs to be made; for example, deme size (N) and rate of extinction (e) can be arbitrarily large or small. This is a difference between the approach used in Cherry and Wakeley (2003), Cherry (2003a), and Whitlock (2003) and our approach. Indeed, these authors calculate various quantities that are similar to the probabilities of identity in Equation 14 by using the fact that the distribution of allele frequencies among demes can be approximated by a βdistribution when N is large and m is small. From our coalescent argument, we obtain the first moments of the distribution of allele frequencies among demes without having to assume that N is large or m small. In fact, Cherry (2003a) and Whitlock's (2003) approximations for the G function are the same as our Equation 42, which assumes small m and large N. We found that Equation 42 gives satisfying results as long as deme size is not too small; for small values of N, however, using the exact expressions of r_{0}, r_{1}, a, and c gives more precise results. For example, with n = 50, N = 15, and m = 0.3, Equation 42 predicts that the probability of fixation of a deleterious allele with s = –0.005 and h = 0 (relative to a neutral allele) is 0.0072, while using the exact expressions of r_{1}, a, and c gives a value of 0.01128, and simulations give 0.01121 (for soft selection). Although assuming either hard or soft selection does not make much difference when N is large and m small (indeed, Equation 42 holds in both cases), this is not the case when N is small, as can be shown by using expressions of the G function that do not make any assumption on N and m—obtained from equations B5, B6, B7, B8, B9 and 36 and 38. This is illustrated by Figure 4, which shows the probability of fixation of a recessive advantageous allele obtained from these expressions and from Equation 42, for low values of N. These examples show that using exact values of the probabilities of coalescence leads to more precise approximations; moreover, exact values of r_{0}, r_{1}, a, and c can be calculated without too much difficulty for many different life cycles, while it is not always easy to obtain approximations for the distribution of allele frequencies among demes (for example, when extinctions occur).
Barton (1993) obtained approximations for the probability of fixation of a single mutant in a model with extinction and additive selection, assuming that deme size is large. We observed that his approximation for the case φ = 0 (Equation 11b in Barton 1993) gives similar results as our solution when migration and extinction rates are small (not shown); however, in our model recolonization does not occur exactly as in Barton (1993); therefore, one must be cautious when comparing our results.
An important consequence of population subdivision is an increase of the fixation probability of recessive advantageous mutations. Figure 5 shows this effect for a population of size 10,000; the different curves correspond to different modes of subdivision (10 demes of 1000 individuals, 100 demes of 100 individuals, and 1000 demes of 10 individuals). One can see that the increase in fixation probability of the recessive mutation is substantial only with “strong” population structure (small migration rate and/or small deme size); a moderate structure such as N = 100, m = 0.1 has virtually no effect on the fixation probability of recessive alleles.
We have demonstrated a simple method to construct diffusion approximations for populations subdivided according to an island model. It is based on the separation of timescales and on the additional fact that for the different life cycles, expressions for M_{δ}_{p} can be deduced mechanically from the fitness function with minimal risk of error by the “direct fitness” method, where the fitness function itself is a direct expression of the life cycle considered. This method has been useful for analyzing selection on various traits under different population structures (e.g., Taylor and Frank 1996; Rousset and Gandon 2002). In this context, fitness has to be defined as the expected number of “successful” gametes of an individual, i.e., those that effectively participate to the gene pool of the next adult generation. Many different uses of the word “fitness” can be found in the literature; it often means “relative fecundity” or “relative survival.” Defining a fitness function as we do proves to be useful since it leads to general expressions like Equation 23.
The direct fitness method used here has also proven convenient to study kin selection (Taylor and Frank 1996; Rousset and Billiard 2000). A substantial part of kin selection theory has been devoted to the evolution of altruistic traits; in some cases, it can be shown that such traits always increase in frequency if a condition of the form –c + rb > 0 holds, where c and b are the cost and benefit of the altruistic act, and r the “relatedness” between interactants, assumed to be independent of the frequency of the allele coding for altruism (Hamilton 1964). Whether such a condition exists is not always obvious, and the present application of the direct fitness formalism may clarify such problems, since from Equation 23 one easily obtains conditions under which an allele increases in frequency in a structured population. For example, when h = ½, Equation 23 shows that M_{δ}_{p} is always positive provided that
Many other applications are outside the scope of this article, but are worth mentioning. For example, one can compute the distribution of allele frequency of recessive deleterious alleles in a metapopulation and consider the implication of different demographic regimes on the purge of such alleles. Average times to fixation or loss of new alleles are easily obtained once the expressions of M_{δ}_{p} and V_{δ}_{p} are known (Crow and Kimura 1970); although in this article we focused on fixation probabilities, fixation times are equally important in determining the rate of evolution of populations. Cherry and Wakeley (2003), Cherry (2003b), and Whitlock (2003) showed that the diffusion method provides good approximations for fixation times in structured populations. In general, a lower rate of migration slows down fixation, while extinction/recolonization can speed it up.
One may wonder about the generality of the present approach with respect to the type of population structure. Although in this article we stayed in the framework of the island model, the methodology presented here can be used with other spatial models, as long as coalescence occurs at two different timescales: a fast one during which selection can be neglected and a slow one that becomes infinitely long as population size tends to infinity. This is the case, for example, for cytoplasmic genes: coalescence is fast if lineages are sampled within the same organism and slow if they are sampled in different organisms. The two scales can also be family vs. population (e.g., Caballero and Hill 1992). In these examples many questions arise about conflicting levels of selection; multilevel selection is easily incorporated into this framework by redefining the fitness function.
Acknowledgments
We thank Joshua Cherry, John Wakeley, and Michael Whitlock for providing articles in press and Joshua Cherry, Montgomery Slatkin, John Wakeley, Michael Whitlock, and an anonymous reviewer for helpful comments on the manuscript. This is paper ISEM 200379.
APPENDIX A
No extinction: When e = 0, V_{δ}_{p} is given by
Extinctionrecolonization: Extinctions introduce correlations in the change in allele frequency in different demes (because the frequency of A in the migrant pool depends on which demes did survive). V_{δ}_{p} is expressed as
We then need to calculate the average covariance of the change in frequency of A in different demes, to the first order in 1/n. We start by noting that at the beginning of the life cycle, the probability of sampling two A alleles with replacement is p^{2} if we sample them from the whole population,
The first term of the covariance is the probability of sampling two A alleles from two different demes, at the next generation. With a probability of order 1/n, these two lineages come from the same deme in the previous generation and are both A with probability
APPENDIX B
Case i: monoecious individuals with random fertilization, gametic migration, and no extinction: The recurrence equations are
Case ii: monoecious individuals with rate α selffertilization, zygotic migration, and no extinction: Because we assume zygotic migration, we need to calculate the probabilities that three genes sampled in three different individuals from the same deme all have a common ancestor (d^{D}), that two of them have a common ancestor (e^{D}), and that they have no common ancestor (f
^{D}), in the infinite island model. The recurrence equations are
The exact expressions of a^{D}, b^{D}, and c^{D} are complicated and not given here for space reasons, but are available by request to the authors.
Case iii: monoecious individuals with random fertilization, gametic migration, extinction rate e, and gametic recolonization: When extinctions occur, the recurrence equations become
Footnotes

Communicating editor: L. Excoffier
 Received February 11, 2003.
 Accepted May 19, 2003.
 Copyright © 2003 by the Genetics Society of America