Selection and Drift in Subdivided Populations: A Straightforward Method for Deriving Diffusion Approximations and Applications Involving Dominance, Selfing and Local Extinctions
Denis Roze, François Rousset

## 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 low-migration 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 low-migration 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 stepping-stone 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, self-fertilization, 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 self-fertilization 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 Mδp=E[δpp] (1) Vδp=E[(δp)2p], (2) where E[x] is the expectation of x. To use the diffusion method we need the following limits: Mδp=limδt0(Mδpδt) (3) Vδp=limδt0(Vδpδt). (4)

If these limits are finite, and if the higher-order moments of δp are negligible, the probability of fixation of A when present in frequency p, u(p), is given by the differential equation Mδpdudp+Vδp2d2udp2=0 (5) (Karlin and Taylor 1981) with solution u(p)=0pG(x)dx01G(x)dx, (6) where G(x)=exp[2MδxVδxdx]. (7)

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 Mδp and Vδx are the limits as N tends to infinity of N × Mδp and N × Vδp, respectively (Karlin and Taylor 1981).

In the simple case of genic selection where A has a selective advantage s over a in an ideal haploid Wright-Fisher population of size N, we have Mδp=spq+o(s) (8) Vδp=pqN+o(1N)+o(s), (9) with q = 1 – p (Ewens 1979; Karlin and Taylor 1981). For Mδp and Vδx to be finite, we need to assume that s tends to zero as N tends to infinity and that the product Ns has a finite limit, say σ. Then Equations 6 and 7 are valid, with Mδp=σpq and Vδx=pq . These expressions of Mδp and Vδx lead to the classical result (e.g., Ewens 1979) u(p)=1exp[2Nsp]1exp[2Ns], (10) where N is assumed to be large and s small. Despite all the approximations done, this formula proves very accurate for values of s as high as 0.1, even when N is as small as 10 (e.g., Carr and Nassar 1970; Gale 1990).

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 = (p1,..., pn) gives the frequencies of A in different demes. We suppose that the number of demes (n) is large and define Mδp and Vδp as the limits as n tends to infinity of the products n × Mδp and n × Vδp, where Mδp and Vδp are still the first two moments of the change of A over one generation. Using a simple argument, we show that these limits can be expressed as functions of the average frequency of A in the population (p). A rigorous demonstration that the process does converge to a diffusion (including the proof that higher-order moments vanish in the limit as n tends to infinity) would follow the general argument of Ethier and Nagylaki (1980) for Markov chains with two timescales. Wakeley (2003) explains how their equations apply to island models of structured populations. We do not repeat such arguments here, but we offer an alternative presentation of the same general idea.

## 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 Wij 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 zaa the phenotype of homozygous aa individuals and assume that the phenotype of AA individuals differs from zaa by an amount s: zAA = zaa + s. We write the phenotype of heterozygous individuals zAa = zaa + sh, h measuring the dominance of A. Defining pij as the frequency of A in individual ij (pij = 0, ½, or 1), pij1 and pij2 as the frequencies of A at each of its two homologous genes, and zij as its phenotype, we have the relation zij=zaa+s[2hpij+(12h)pij1pij2]. (11)

The expected change in frequency of A over one generation is Mδp=12nNi=1nj=1NWijpijp. (12) To the first order in s, this is Mδp=s2nNi=1nj=1NdWijdspij+o(s)=s2nNi=1nj=1N(Wijzijdzijds+Wijzidzids+Wijzdzds)pij+o(s), (13) as Wij depends on zij, zi, and z. Equations 11 and 13 give Mδp=s2Wijzij[2hpij2¯+(12h)pij1pij2¯]+s2Wijzi[2hpi2¯+(12h)pij1pij2pi¯]+s2Wijz[2hp2+(12h)pij1pij2p¯]+o(s),1 (14) where the overbar means the average among demes and individuals over the whole population.

Mδp has been defined as the limit of the product n × Mδp as n goes to infinity. For this limit to be finite, we have to assume in principle that s is of order 1/n, so that the product ns has a finite limit as n tends to infinity (for the same reason, s had to be of order 1/N in the panmictic case—see previous section). Therefore our method supposes a high number of demes (n high) and weak selection (s small), although simulations will show that these are not major constraints. As s is of order 1/n, the terms in o(s) of Mδp vanish when we take the limit of n × Mδp as n goes to infinity.

The different averages in Equation 14 can be seen as probabilities of identity: for example, pi2¯ is the probability of obtaining two A alleles after sampling two genes with replacement from the same deme. In the diffusion limit (n tends to infinity, s tends to zero, ns finite), these probabilities converge to values that depend only on p, the overall frequency of A. Again, we take the example of pi2¯ , the probability that two genes sampled with replacement from the same deme are both A. Consider the two ancestral lineages of these genes backward in time: either these lineages stay and coalesce in the same deme or one of them (or both) migrates to another deme before coalescence occurs; in the first case, the expected coalescence time does not depend on the number of demes n, while in the latter case, coalescence time becomes infinitely long as n tends to infinity. This has been described as a separation of timescales (e.g., Wakeley 2003 and references therein), because lineages in different demes coalesce at a rate inversely related to n, while genes within demes can coalesce within a few generations, at a rate depending on deme size and migration (roughly, on m + 1/N in the island model), but essentially independent of n. Therefore in the limit as n tends to infinity, at time t independent of n in the past we can consider that the two lineages have either stayed in the same deme and coalesced (with a probability that we call r1R ) or have migrated to different demes and have not coalesced (with probability 1 − r1R ). Thus the probability that our two genes are A is equal to r1Rpt+(1r1R)pt2 , where pt is frequency of A in the whole population at time t. However, in the limit as n tends to infinity and s tends to zero, the frequency of A in the population does not change over the finite number t of generations, thus pt = p, the present frequency of A. Therefore pi2¯=r1Rp+(1r1R)p2+O(1n)=p2+r1Rpq+O(1n). (15) This expression, as well as all further expressions involving O(1/n) terms, is an average over all parental populations with the same value of p. These parental populations may differ in the frequency of demes with different copies of the A allele, so that for each possible parental population pi2¯ may depart from the above average. But the magnitude of these differences vanishes as the number of demes increases, which allows us to use diffusion equations (see Wakeley 2003, Equation 33, and also Ethier and Nagylaki 1980, Equation 1.5). In the infinite island model, one recovers the classical expression pi2¯=p2+r1Rpq for inbreeding coefficients (e.g., Crow and Kimura 1970). Here r1R is equivalent to Wright's FST (Hudson 1998; Rousset 2002).

View this table:
TABLE 1

Parameters of the model

The same reasoning is used to calculate the limit of pij2¯ , which is the probability that two genes sampled with replacement in the same individual are both A: pij2¯=p2+r0Rpq+O(1n). (16) Here r0R is the probability of coalescence for two genes sampled with replacement from the same individual, in a population with an infinite number of demes.

Denote r1D the probability of coalescence of two genes sampled from the same deme without replacement, and r0D the probability of coalescence of the homologous genes of an individual (still in the limit as n tends to infinity). We have the relations r1R=12N+12Nr0D+(11N)r1D (17) and r0R=12+12r0D. (18) r0D can be shown to be equivalent to Wright's FIT in the infinite island model. As pij1pij2¯ is the probability that the two homologous genes of an individual are both A, we have pij1pij2¯=p2+r0Dpq+O(1n). (19)

Last, we need the limit as n tends to infinity of pij1pij2pi¯ , which measures the probability that the two homologous genes of an individual and a third gene sampled from the same deme are all A. We call aR the limit as n tends to infinity of the probability that the three lineages coalesce before one of them moves to another deme, bR the limit of the probability that only two of them stay in the same deme and coalesce, and cR = 1 – aRbR the limit of the probability that two lineages move to different demes before any coalescence event has occurred. We have pij1pij2pi¯=aRp+bRp2+cRp3+O(1n)=p3+aRpq+(1cR)p2q+O(1n). (20) Denote aD, bD, and cD the same probabilities when the third gene is sampled in a different individual. We have the relations aR=1Nr0D+(11N)aD (21) and cR=(11N)cD. (22)

View this table:
TABLE 2

Variables of the model

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 Mδp=s2Wijzij[h(1+r0D)+(12h)[r0D+(1r0D)p]]pq+s2Wijzi[2hr1R+(12h)[aR+(1cRr0D)p]]pq+o(1n). (23)

## 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 Vδp=Var0[p]+o(1n), (24) where p′= pp and Var0 is the variance in the neutral case (s = 0). We look for an expression for Var0[p′] to the first order in 1/n. We have Var0[p]=1n2i=1nVar0[pi]+1n2i=1njiCov0[pi,pj]. (25) We show in appendix a that these variance and covariance terms can be expressed, to first order in 1/n, as functions of the frequency of A in the whole population (p) and an effective population size Ne depending on the various parameters of the model (but not on allele frequencies). By definition, Ne is the variance effective population size (Ewens 1982).

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 Vδp=2(1r1D)(1r0D)2nNpq+o(1n), (26) which is of the form pq/2Ne, with Ne=nN2(1r1D)(1r0D), (27) a result known under the form nN/(1 – FST)(1 + FIS), where FST=r1D and FIS=(r0Dr1D)(1r1D) (see Wang and Caballero 1999 for a more general result). In the case of gametic migration, monoecious individuals, and random fertilization, we have r0D=r1D , giving Ne=nN1r1D. (28)

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 Ne=n(1e)2r1R[1(1e)2(1m)2]. (29) This is also the eigenvalue effective size (Ewens 1979, 1982), derived by several independent arguments in Rousset (2003).

## 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 self-fertilization, 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) r0D=r1D11+4M, (30) (aD,bD,cD)(1(1+4M)(1+2M),6M(1+4M)(1+2M),8M2(1+4M)(1+2M)). (31) These approximations are also valid for zygotic migration and dioecious individuals (not shown). Although, as found by Whitlock (2002), these approximate expressions satisfy the relation aD=2(r1D)2(1+r1D) , this is no longer the case when the exact expressions (given in appendix b) are used.

Case ii: partial self-fertilization, no extinction: Here we assume that self-fertilization occurs at a rate α and that fertilization precedes migration. Random fertilization corresponds to α= 1/N. We obtain the approximations r0D1+2αM1+2(2α)M,r1D11+2(2α)M, (32) aD1+αM[1+2(2α)M][1+(2α)M],bD2M[3α2M+2α(M1)][1+2(2α)M][1+(2α)M],cD4(1α)(2α)M2[1+2(2α)M][1+(2α)M]. (33)

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 φ =ψ12 = 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 r0D=r1D1+eNk1+4M+2eN[1ϕ(112k)]. (34) This approximation follows from the recursion (B20) in appendix b and was also obtained by Whitlock and McCauley (1990). We could not obtain any simple approximation for aD, bD, and cD.

## 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 zaa to zero and assume that an individual of phenotype zij produces a number of gametes proportional to 1 + zij; 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 Mδp —see Equation 23.

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 Wij=2N(1e)[(1m)1+zijN(1+zi)+m1+zijN(1+zi)+e1e1+zijN(1+zi)]+O(1n). (35) Individual ij can participate to the next generation only if its deme survives (probability 1 – e). The first term in the parentheses is the frequency of gametes produced by individual ij present in deme i after migration; the second term is the frequency of gametes produced by ij in any other nonextinct deme, after migration, times the number of other nonextinct demes; and the third term is the frequency of gametes produced by ij in the pool of recolonizers, times the number of extinct demes. Simplifying Equation 35 leads to Wij=21+zij1+zi+O(1n), (36) giving ∂Wij/∂zij = 2, ∂Wij/∂zi =–2 (derivatives have to be evaluated for s = 0, which means here zij, zi, z = 0).

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 Wij=2N(1e)[(1m)(1+zij)N[(1m)(1+zi)+m(1+z)]+m1+zijN[1+z]+e1e1+zijN[1+z]]+O(1n). (37) Again, the first term in the parentheses represents the proportion of gametes produced by individual ij present in deme i after migration (provided that deme i did not go extinct); the second term is the average proportion of gametes produced by ij in the other nonextinct demes, after migration, times the number of other nonextinct demes; and the third term is the proportion of gametes produced by ij in the pool of recolonizers, times the number of extinct demes (we assume that surviving demes contribute to the pool of colonizers in proportion of the total number of gametes present in the deme). Equation 37 leads to Wij=2(1e)[(1m)(1+zij)1+(1m)zi+mz+(m+e1e)1+zij1+z]+O(1n), (38) giving ∂Wij/∂zij = 2, ∂Wij/∂zi = –2 (1 – e) (1 – m)2.

## 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 G(x)=exp[2nNsx], (39) while under soft selection and gametic migration, G(x)=exp[2n(N1)sx]. (40) The N – 1 factor is consistent with the fact that soft selection means no selection when N = 1. In the case of zygotic migration, N has to be replaced by N (2N)/(2N – 1) in Equation 39, while N – 1 has to be replaced by N (2N – 2)/(2N – 1) in Equation 40. These differences are significant only for small values of N.

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 G(x)=exp{n(N1)(1e)s2(1r1D)(1r0D)r1DN[1(1e)2(1m)2]}. (41)

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 G(x)exp{2nNsx1+2Nmx+4Nmh(1x)1+2Nm}, (42) which has to be integrated numerically to obtain u(p), the probability of fixation of A. For advantageous mutations (s > 0), a crude approximation can be obtained by neglecting terms in x2 in the expression of G, giving u(12nN)s1+4Nmh1+2Nm. (43) We found that this approximation often gives a good idea of the probability of fixation of a rare advantageous allele, except when both the migration rate is large and h is close to zero.

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) low-migration approximation performs better (this is discussed later).

Figure 1.

—Probability of fixation of an advantageous mutation (s = 0.01) when present as a single copy, as a function of its dominance coefficient (h). Migration is gametic, individuals are monoecious with random fertilization, and no extinction occurs (e = 0). Number of demes, n = 200; deme size, N = 100. Results from the diffusion model are: dotted line, migration rate m = 0.1; dotted/dashed line, m = 0.01; dashed line, m = 0.001. Simulation results are: solid circles, m = 0.1; open circles, m = 0.01; squares, m = 0.001. The solid line corresponds to Slatkin's (1981) low migration limit.

Effects of selfing: In case ii (monoecious individuals with rate α of self-fertilization, zygotic migration, and no extinction), approximations (32) and (33) lead to G(x)=exp{2nNsx[2h+(12h)1+(2α)Nm[α+(1α)x]1+(2α)Nm]}, (44) which has to be integrated numerically. Neglecting the terms in x2, we obtain a crude approximation for advantageous mutations: u(12nN)s1+(2α)[α+2h(1α)]Nm1+(2α)Nm. (45)

Figure 2 shows the effect of partial self-fertilization 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 low-migration 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 two-dimensional stepping-stone model (Michalakis and Olivieri 1993).

Figure 2.

—Probability of fixation of a recessive advantageous mutation (s = 0.01, h = 0), as a function of the rate of self-fertilization (α). Individuals are monoecious and migration is zygotic. Number of demes, n = 200; deme size, N = 100; no extinction occurs (e = 0). Results from the diffusion model are: dotted line, migration rate m = 0.01; solid line, migration rate m = 0.05. Simulation results are: open circles, m = 0.01; solid circles, m = 0.05.

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 (r0, r1, 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 not-too-small 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 within-deme coalescence time becomes longer, and the effect of selection during that time cannot be neglected any more.

Figure 3.

—Probability of fixation of an underdominant advantageous mutation (s = 0.02, h =–0.1) relative to a neutral mutation, as a function of the extinction rate (e). Individuals are monoecious with random fertilization, and migration and recolonization are gametic; number of demes, n = 200; deme size, N = 100; migration rate, m = 0.1. A corresponds to the “migrant pool” model and B to the “propagule pool” model. Results from the diffusion model are: solid line, number of colonizers k = 2; dashed line, k = 10; dotted line, k = 100. Simulation results are: solid circles, k = 2; open circles, k = 10; squares, k = 100.

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 r0, r1, 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 r1, 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 r0, r1, 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).

View this table:
TABLE 3

Accuracy of the diffusion results

Figure 4.

—Probability of fixation of an advantageous recessive mutation (h = 0, s = 0.02), relative to a neutral mutation, as a function of deme size (N). Migration is gametic, individuals are monoecious with random fertilization, and no extinction occurs (e = 0). Number of demes, n = 100; migration rate, m = 0.3. Dashed line, soft selection; dotted line, hard selection; solid line, approximation (42) for N large and m small.

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.

Figure 5.

—Effect of population structure on the probability of fixation of a recessive advantageous mutation (s = 0.01, h = 0). x-axis, migration rate; y-axis, probability of fixation of the mutation when present as a single copy. Migration is gametic and fertilization random; extinctions do not occur (e = 0). Solid line, n = 1000, N = 10; dashed line, n = 100, N = 100; dotted line, n = 10, N = 1000.

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 Wijzij+2r1R1+r0DWijzi>0. (46) This is the correct form of Hamilton's rule that one obtains by modern inclusive fitness techniques; it shows that –c, b, and r have to be defined as ∂Wij/∂zij, ∂Wij/∂zi, and 2r R1/(1 + r D0), respectively. When h ≠ ½ we can also define a relatedness coefficient from Equation 23, but now this coefficient depends on p, the frequency of A. Conversely, our model can be used to include the effects of genetic drift in kin selection models and derive fixation probabilities and fixation times of alleles coding for altruistic behaviors.

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 2003-79.

## APPENDIX A

No extinction: When e = 0, Vδp is given by Vδp=1n2i=1nVar0[pi]+o(1n) (A1) since the covariance in the change in allele frequency in two different demes equals zero (sampling occurs independently in each deme). We have Var0[pi]=E0[pi]E0[qi]E0[piqi] (A2) and E0[pi]=(1m)pi+mp (A3) and since the averages over i of E0[piqi] and piqi both tend to (1r1R)pq as n tends to infinity, we obtain Vδp=1n[1(1m)2]r1Rpq+o(1n)=(r1Rr1D)pqn+o(1n) (A4) as r1D=(1m)2r1R . Using Equation 17, one arrives at Vδp=2(1r1D)(1r0D)2nNpq+o(1n) (A5) and thus we can define a variance effective population size as Ne=nN2(1r1D)(1r0D) (A6) so that Vδp=pq2Ne+o(1n). (A7)

Extinction-recolonization: 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 Vδp=1n2i=1nVar0[pi]+1n2i=1njiCov0[pi,pj]+o(1n). (A8) The first term is obtained as above. Var0[pi] is still given by Equation A2, with E0[pi]=(1e)[(1m)pi+mp]+ep, (A9) giving 1n2i=1nVar0[pi]=1n[1(1e)2(1m)2]r1Rpq+o(1n). (A10)

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 p2 if we sample them from the whole population, pi2¯ if we sample them from the same deme, and pipj¯ if we sample them from different demes (ij). This gives p2=1npi2¯+(11n)pipj¯ (A11) and since pi2¯ tends to p2+r1Rpq as n tends to infinity, pipj¯=p21nr1Rpq+o(1n). (A12) We have Cov[pi,pj]=E[pipj]E[pi]E[pj]. (A13) From Equations A9 and A12 we obtain the second term of the covariance: 1n2ijE[pi]E[pj]=p21n(1e)2(1m)2r1Rpq+o(1n). (A14)

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 p2+r1Rpq+O(1n) , while with a probability of order 1 they come from different demes, in which case the probability that they are A is given by Equation A12. We calculate the probability that they come from the same deme as follows: if the two demes where the genes are sampled have not become extinct in the previous generation [probability (1 – e)2], the probability that the two lineages come from a single deme in the previous generation is 1(1m)2n(1e)+o(1n), (A15) while if one or both demes became extinct [probability 1 – (1 – e)2], this probability is 1/[n (1 – e)]. Taking all cases into account, one arrives at 1n2ijE[pipj]=p2+[1(1e)2(1m)2n(1e)1n]r1Rpq+o(1n), (A16) which, combined with Equations A10 and A14, finally gives Vδp=1(1e)2(1m)2n(1e)r1Rpq+o(1n). (A17) Therefore, when extinctions occur, the variance effective population size is Ne=n(1e)2r1R[1(1e)2(1m)2]. (A18)

## APPENDIX B

Case i: monoecious individuals with random fertilization, gametic migration, and no extinction: The recurrence equations are r0D=r1D=(1m)2(12N+(112N)r1D) (B1) aD=(1m)3[1(2N)2+32N(112N)r1D+(112N)(122N)aD] (B2) bD=3m(1m)2[12N+(112N)r1D]+(1m)3[32N(112N)(1r1D)+(112N)(122N)bD] (B3) cD=m3+3m2(1m)+3m(1m)2(112N)(1r1D)+(1m)3(112N)(122N)cD, (B4) which gives r0D=r1D=(1m)22N(1m)2(2N1) (B5) aD=(1m)3Y[N+(2N1)(1m)2] (B6) bD=3m(1m)2NY[2N+(1m)(2m)(2N1)] (B7) cD=2m2N2Y[2N(32m)+(3m)(1m)2(2N1)] (B8) with Y=[2N(2N1)(1m)2][2N2(2N1)(N1)(1m)3]. (B9)

Case ii: monoecious individuals with rate α self-fertilization, 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 (dD), that two of them have a common ancestor (eD), and that they have no common ancestor (f D), in the infinite island model. The recurrence equations are r0D=(1α)r1D+α(12+12r0D) (B10) r1D=(1m)2[12N+12Nr0D+(11N)r1D] (B11) aD=(1m)2(1α)(1Nr1D+1NaD+(12N)dD)+(1m)2α(14N+34Nr0D+12(11N)(r1D+aD)) (B12) bD=[1(1m)2][(1α)r1D+α(12+12r0D)]+(1m)2(1α)[1N(1r1D)+1NbD+(12N)eD]+(1m)2α[34N(1r0D)+12(11N)(1r1D+bD)] (B13) cD=[1(1m)2][(1α)(1r1D)+α12(1r0D)]+(1m)2[(1α)(1NcD+(12N)fD)+α12(11N)cD] (B14) dD=(1m)3[1(2N)2+3(2N)2r0D+32N(11N)(r1D+aD)+(11N)(12N)dD] (B15) eD=3m(1m)2[12N+12Nr0D+(11N)r1D]+(1m)3[3(2N)2(1r0D)+32N(11N)(1r1D+bD)+(11N)(12N)eD] (B16) fD=m3+3m2(1m)+3m(1m)2[12N(1r0D)+(11N)(1r1D)]+(1m)3[32N(11N)cD+(11N)(12N)fD], (B17) which gives for r0D and r1D r0D=1[1(1m)2][1αN]1[1(1m)2][1(2α)N], (B18) r1D=(1m)21[1(1m)2][1(2α)N]. (B19)

The exact expressions of aD, bD, and cD 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 r0D=r1D=[(1e)(1m)2+e(112k)ϕ]×[12N+(112N)r1D]+e2k (B20) (see Rousset 2003). aD=[(1e)(1m)3+e(112k)(11k)ψ1]×[1(2N)2+(32N)(112N)r1D+(112N)(11N)aD]+e[1(2k)2+32k(112k)ϕ(12N+(112N)r1D)] (B21) bD=[(1e)(1m)3+e(112k)(11k)ψ1]×[32N(112N)(1r1D)+(112N)(11N)bD]+[3(1e)m(1m)3+e(112k)(11k)ψ2][12N+(112N)r1D]+e32k(112k)[(1ϕ)+ϕ(112N)(1r1D)] (B22) cD=(1e)[m3+3m2(1m)+3m(1m)2(112N)(1r1D)+(1m)3(112N)(11N)cD]+e(112k)(11k)[1ψ1ψ2+ψ2(112N)(1r1D)+ψ1(112N)(11N)cD]. (B23) Here again, we give only the solution for r0D=r1D , while the solutions for aD, bD, and cD are available on request: r0D=r1D=2k(1m)2(1e)+2Ne+(2k1)eϕ2k[2N(2N1)(1e)(1m)2](2k1)(2N1)eϕ. (B24)

## Footnotes

• Communicating editor: L. Excoffier