Maximum-Likelihood Estimation of Migration Rates and Effective Population Numbers in Two Populations Using a Coalescent Approach

Corresponding author: Peter Beerli, J255 Health Sciences, Department of Genetics, University of Washington, Box 357360, Seattle, WA 98195-7360., beerli{at}genetics.washington.edu (E-mail)

Communicating editor: S. TAVARÉ

	ABSTRACT

TOP ABSTRACT MODEL IMPLEMENTATION SIMULATION RESULTS DISCUSSION APPENDIX 1 LITERATURE CITED

A new method for the estimation of migration rates and effective population sizes is described. It uses a maximum-likelihood framework based on coalescence theory. The parameters are estimated by Metropolis-Hastings importance sampling. In a two-population model this method estimates four parameters: the effective population size and the immigration rate for each population relative to the mutation rate. Summarizing over loci can be done by assuming either that the mutation rate is the same for all loci or that the mutation rates are gamma distributed among loci but the same for all sites of a locus. The estimates are as good as or better than those from an optimized F_ST-based measure. The program is available on the World Wide Web at http://evolution.genetics.washington.edu/lamarc.html.

SEVERAL methods for the estimation of migration rates between subpopulations have been proposed. We can subdivide them into two very different approaches: (1) marking individuals and tracking their individual movements and then extrapolating these individual movements to migration rates; or (2) surveying genetic markers in the populations of interest and calculating a migration rate from allele frequencies or sequence differences. Of course, one should be aware that these two approaches do not estimate the same quantity: approach 1 estimates the actual "instantaneous" migration rate of individuals, whereas approach 2 reflects an average over a time period whose length is determined by the rate of mutation per generation of the locus under study or by the time scale of genetic drift. The genetic approach tends to gives a lower estimate than the individual migration rates approach because the method is looking at changes that become established in the subpopulation gene pool.

Current estimation methods for genetic data are methods such as those related to F_ST (e.g., SLATKIN 1991 ; SLATKIN and HUDSON 1991 ), the rare allele approach of SLATKIN 1985 , a maximum-likelihood method using gene frequency distributions (RANNALA and HARTIGAN 1996 ; TUFTO et al. 1996 ), and approaches based on coalescent theory (KINGMAN 1982A , KINGMAN 1982B ), such as the cladistic approach of SLATKIN and MADDISON 1989 , the method outlined in WAKELEY 1998 , and maximum likelihood using coalescent theory (NATH and GRIFFITHS 1993 , NATH and GRIFFITHS 1996 ). We describe here a new method using a maximum-likelihood- and coalescent theory-based approach.

Our method integrates over all possible genealogies and migration events. This should be superior to pairwise estimators such as those using F_ST or related statistics (cf. FELSENSTEIN 1992B ) and also more powerful than the cladistic approach of SLATKIN and MADDISON 1989 , which needs to know the true genealogy. Our approach estimates similar parameters as the program of Bahlo and Griffiths (http://www.maths.monash.edu.au/~mbahlo/mpg/gtree.html), which is based on the work of GRIFFITHS and TAVARE 1994 and NATH and GRIFFITHS 1996 . It differs in the way we search the genealogy space, and our methods support mutation models for different types of data: the infinite allele model for electrophoretic markers, a one-step model for microsatellites, and a finite-sites model for nucleotide sequences. The sequence model is more useful than the infinite sites model, which forces the researcher to discard data when there are multiple mutations at the same sites.

	MODEL

TOP ABSTRACT MODEL IMPLEMENTATION SIMULATION RESULTS DISCUSSION APPENDIX 1 LITERATURE CITED

We propose a method to make a maximum-likelihood estimate of population parameters for geographically subdivided populations. The general outline of such estimates involves extending coalescence theory (KINGMAN 1982A , KINGMAN 1982B ) to include migration events. Migration models with coalescents were first developed by TAKAHATA 1988 and TAKAHATA and SLATKIN 1990 for two gene copies and discussed more generally for n gene copies in two populations by HUDSON 1990 , with generalization to multiple populations and different models of migration by NOTOHARA 1990 . NATH and GRIFFITHS 1993 , NATH and GRIFFITHS 1996 used this migration-coalescent process for maximum-likelihood estimation of one parameter, the effective number of migrants = 4N_em, where N_e is the effective population size and m is the migration rate per generation.

Our migration model consists of two populations. These have effective population sizes N⁽¹⁾_e and N⁽²⁾_e, which need not be equal. The populations exchange migrants at rates m₁ and m₂ per generation (Figure 1).

View larger version (8K): In this window In a new window Download PPT slide   Figure 1. Two-population model with four parameters: m1 is the migration rate per generation from population 2 to 1, m2 from population 1 to 2, and N(i)e is the effective population size.

Because we observe only genetic differences and do not know the mutation rate µ, we must absorb µ into the parameters so that we use as our parameters

and

and occasionally write

We make the following assumptions: that diploid individuals in each population reproduce according to any standard population genetic model that converges to the same diffusion equation as a Wright-Fisher model, that populations have a constant population size through time so that they do not grow or shrink, and that the mutation rate µ is constant. These are also the assumptions used by KINGMAN 1982A for coalescence theory. We use the convention with sequence data that µ is the rate of mutation per generation per site, whereas with microsatellite or enzyme electrophoretic data we take µ to be the rate of mutation per generation per locus. To emphasize this distinction between sites and loci, we use a capital in the first case and a small for the latter ones. The mutation rate µ can vary among loci according to a gamma distribution. In addition we also need to assume that the two populations exchange migrants with constant rates per generation.

The likelihood formula for this family of methods (FELSENSTEIN 1988 , FELSENSTEIN 1992A ; KUHNER et al. 1995 , KUHNER et al. 1998 ) is based on the product between the genealogy likelihood Prob(D|G) (e.g., FELSENSTEIN 1973 ; SWOFFORD et al. 1996 ) and the prior probability of the coalescent genealogy Prob(G|) (KINGMAN 1982A , KINGMAN 1982B )

(1)

where the parameters are

TAKAHATA and SLATKIN 1990 found that Kingman's results could not be extended to cases with geographical structure. They were unable to obtain the distribution of times to coalescence in the presence of migration. We have avoided this difficulty by having the genealogy G specify not only the coalescences but also the times and places of migration events. With this information in G, its probability density becomes for two populations a product of terms for the intervals between events in the genealogy

(2)

where

(3)

Prob(G|) is obtained by computing for each interval the probability that nothing happens during the interval times the probability of the particular event at the bottom (beginning) of the interval. The variable v_i is the total number of coalescences in subpopulation i and w_.i is the sum of the numbers of migration events to subpopulation i over all time intervals T. Furthermore, p_j is the probability that no event occurs during time interval j, u_j is the length of time interval j, and k_ji is the number of lineages in subpopulation i during time interval j.

If we extend the parameter estimation from a single locus to many unlinked loci, we face the additional problem that, although µ may be constant per locus, it can vary between loci. Neglecting variation of µ, combining loci using the likelihood framework is simple:

(4)

On the other hand, if we assume that µ follows a gamma distribution and therefore can vary between loci, we need to include an additional parameter, the shape parameter of this distribution. The other parameter of this gamma distribution is a scale parameter that can be conveniently chosen to be . By parameterizing in this way, 1/ is the squared coefficient of variation of . We cannot estimate the mean mutation rate of this gamma distribution directly. As N_e is the same for all loci, we can use an arbitrary (we have chosen the effective population size of the first population) and scale the other parameters relative to it. We replace µ in the parameters _i and _i with a new variable, the mutation rate x. is the mutation rate x scaled by 4N⁽¹⁾_e. We integrate over all possible values of instead of x. This does not change the shape of our likelihood surface, because N⁽¹⁾_e in is constant, so that we have

(5)

where

and

This integral over all mutation rate values is implemented using a discrete gamma approximation with 100 intervals (cf. YANG 1994 ).

Note that the gamma distribution described here is the distribution of µ across loci, not the more commonly used gamma distribution of rates across sites within a locus. If needed the latter can be approximated by Hidden Markov model methods in the case of DNA sequences.

	IMPLEMENTATION

TOP ABSTRACT MODEL IMPLEMENTATION SIMULATION RESULTS DISCUSSION APPENDIX 1 LITERATURE CITED

To calculate the likelihood L() for the parameters using (1) we ought to sum over all possible genealogies, including all topologies with all possible branch lengths and with all possible migration scenarios. This is not possible, as even with two tips there are an infinite number of possible sequences of migration events with different branch length. We have to resort to an approximation of the likelihood function using a Metropolis-Hastings Markov chain Monte Carlo approach (METROPOLIS et al. 1953 ; HAMMERSLEY and HANDSCOMB 1964 ; HASTINGS 1970 ; KUHNER et al. 1995 ). The derivation of the importance sampling function is shown in the Appendix 1.

The parameter estimation is simple in principle: we have to find the maximum of the likelihood function, which has five or six dimensions, depending on whether µ is allowed to vary between loci. This is done by a modified damped Newton-Raphson method (DAHLQUIST and BJORK 1974 ) using explicit equations for the first and second derivatives.

The genealogy sampler is of more interest. A large sample of genealogies G₁, G₂, ... ,G_m is drawn from a distribution whose density is proportional to Prob(D|G) Prob(G|₀). This sampling scheme is biased toward genealogies that contribute more to the likelihood with a given parameter set ₀. We find the relative likelihood at other values by dividing by the density from which we sampled the genealogies (see also Appendix 1). If m is the number of sampled genealogies, we get

(6)

In the program Coalesce (KUHNER et al. 1995 ) a region of internal nodes in the genealogy is erased and rebuilt according to the coalescent model. This scheme is not applicable to migration models. We adopt a new scheme that is much more general and that can be used in any extension to the coalescent model. The first genealogy for each locus is constructed with a UPGMA method (as implemented by M. K. Kuhner and J. Yamato in FELSENSTEIN 1993 ), and the minimal necessary migration events are inserted using a Fitch parsimony algorithm (FITCH 1971 ). The times for coalescent events or migration events on this first genealogy are calculated with (3) using a uniform random number for p_j and solving for u_j, and parameter values, which are either guessed by the researcher or calculated using an F_ST-based method (see Appendix 1).

The genealogies are sampled as follows (cf. Figure 2): A coalescent node or tip z on this genealogy is chosen at random, the lineage below it is dissolved, and the node is used as the starting point to simulate the ancestry using a migration-coalescent process through a series of time intervals (Figure 2). The other lineages do not change and are taken as given and form the partial genealogy G_p. For the further description of the rearrangement process we use "up" and "top" for being or moving closer to the tips of the genealogy and "down" and "bottom" for being or moving toward the root of the genealogy. Starting at the chosen node z (see Figure 2), we draw a new time interval by solving for u_j in

(7)

after substituting a uniform deviate for p^'_j, where k^'_ji is the number of lineages of population i in the partial genealogy G_p during the time interval j. If this new time is further back in time than the bottom of the time interval j on the partial genealogy G_p, the active lineage advances down to that time and the simulation process starts again. When the newly drawn time is in the active time interval, an event, either a migration or coalescence, is drawn according to their relative probabilities of occurrence, e.g.,

(8)

View larger version (21K): In this window In a new window Download PPT slide   Figure 2. Transition from an old genealogy Go to a new genealogy Gn. Dotted lines are the times of coalescences or migrations. (A) Genealogy Go with migration. The black bar marks a migration from the white population to the black or vice versa; z is the node to be picked. (B) Partial genealogy Gp after drawing the coalescent node z at random and dissolving the branch to the next coalescent below. (C) Simulation of the coalescent with migration. One possible outcome with three consecutive steps is shown: (1) using Equation 7 a new time interval is drawn and a migration event from white to black is also drawn (Equation 8) and the lineage is extended down to that event; (2) a new time interval is drawn: it extends too far back, so the lineage advances down to the time j; (3) a new time interval is drawn with a coalescent event at its bottom end. The process stops at time k. (D) the final configuration Gn.

If a migration event occurs the process moves down to the time when the migration event happened. Below a migration event the active lineage is in the other population and the migration-coalescent process continues (see also Figure 2C). After a coalescent event occurs the process stops. The transformation from one genealogy into another allows change from any genealogy over several steps into any other and therefore allows the process to search the whole genealogy space. On the other hand, the resimulation of a single line guarantees that newly found genealogy is similar to the old genealogy.

In cases where all the lineages have not coalesced by the time of the root of the partial genealogy, simulation on the active lineage and the bottom lineage of G_p continues until the lineages coalesce (Figure 3). When simulating on two lineages we draw new time intervals using (3) instead of (7). There is one exception to this rule at the bottom of the genealogy: if by chance the dissolved lineage is the bottommost lineage of one of the populations on G_o, we keep simulating with a single active lineage below the root of G_p until the former root of G_o is reached and then start to simulate on both of the two remaining lines. The lineage between the root of G_p and the old root on G_o on the partial genealogy has a known history, so there is no need to simulate this history again.

View larger version (11K): In this window In a new window Download PPT slide   Figure 3. Simulation on one or two lineages at the bottom of a genealogy. Striped lineages are active. On these lineages the migration-coalescent process is used to find new time intervals and events. (A) old genealogy Go, with 1 indicating the branch to be cut to get genealogy B, and 2 indicating the branch to be cut to get C. (B) The root R is below the cut-point D, and so below the root we need to calculate the migrations and coalescences for two lineages, because there is no previous information present. (C) In this genealogy a single lineage simulation is needed, using the information present between R and D. Below D we need to simulate on both remaining lineages.

A change to the newly found genealogy G_n is accepted with probability

(9)

where r_M is the Metropolis term

(10)

and r_H is the Hastings term

(11)

This acceptance probability r is based on the work of HASTINGS 1970 . In our specific implementation it can be simplified. The scheme for changing the genealogy creates a new genealogy from an old one by first randomly removing a branch with all its possible migration events, which creates a partial genealogy G_p. The probability of a particular change is only governed by the number of coalescent nodes in the genealogy, and this number is constant so that the probabilities Prob(G_p|G_n) and Prob(G_p|G_o) are equal. For a specific triplet of original, partial, and new genealogies G_o, G_p, and G_n

(12)

holds. The resimulation process (7)–(8) is driven only by the parameters and so Prob(G|G_p) is proportional to Prob(G|) so that

(13)

This results in cancellation of the Prob(G|) terms in r, so that we have

(14)

If the genealogy is accepted, it is the starting point for a new cycle; otherwise the previous one continues to be used. After having sampled a large number of genealogies we estimate the parameters by maximizing (6). In theory a single long Markov chain would be enough to find the genealogies that contribute most to the likelihood, but if we start from a very bad ₀ we can need an excessively long time to find this region. If the likelihood of a new genealogy is only slightly better than that of the old one, the search is moving nearly randomly and can spend a long time in regions that do not contribute much to the final likelihood. The sampler will move more quickly if the likelihood surface is steep. If we start with bad parameters there is some chance that the likelihood surface is locally very flat. To overcome this, we restart a new Markov chain with the improved and the last genealogy of the most recent chain. Our implementation uses a number of short chains in which the parameter estimates are based only on a few genealogies to nd regions with good parameter values. It then uses an arbitrary number of long chains to get good estimates of the parameters. All but the last chain are discarded and not used for the final result. KUHNER et al. 1998 have not found a decrease in performance compared to the scheme of taking several chains into account for the final result. If we run this genealogy sampling process long enough and with enough chains, the sampler will find its way to the region with the biggest contribution to the likelihood and therefore approximate the maximum-likelihood estimate of our parameters . The main problem with importance sampling schemes is that little is known about how we can be sure that we have sampled enough genealogies from the region that contributes most to our likelihood function. It is common practice in our group to run 10 "short" chains and then 2 "long" chains. Using random values for the first parameter set, we have found that for a given data set the results converge to a value that is not dependent on these initial parameter values when the total numbers of sampled genealogies exceed 30,000 genealogies (data not shown). These are 10 short chains with 1000 genealogies and 2 long chains with 10,000 genealogies sampled. It seems to be sufficient to deliver acceptable estimates for simulated data and averages of those, although we would recommend choosing a better start parameter than a random guess and running the program at least twice as long. The number of genealogies needed for an accurate estimate is certainly dependent on many factors, such as starting parameters, first genealogy, number of sampled individuals, variability in the data set, and number of migration events.

	SIMULATION RESULTS

TOP ABSTRACT MODEL IMPLEMENTATION SIMULATION RESULTS DISCUSSION APPENDIX 1 LITERATURE CITED

We have simulated genealogies with known true parameters _T using a coalescent with migration and evolved data according to our data models along the branches starting from the root of the genealogy (HUDSON 1983 ). The data that resulted at the tips were used as the input data for our method. The mutation rate was held constant for all simulations except one, in which we used a gamma-distributed mutation rate with shape parameter = 1.0. For all other simulations we used 10 short chains sampling 1000 genealogies, keeping 100 of these genealogies for the parameter estimation and 3 long chains sampling 10,000 genealogies, using 1000 of them. The last chain delivered the parameters shown in the tables, and all other chains were discarded.

For the comparison with an F_ST-based estimator (see Appendix 1) we used a symmetrical parameter setup for the simulations, ₁ = ₂ and ₁ = ₂. However, the parameters that were estimated were not constrained to be symmetric.

Our simulations generally recover the true _T very well (Table 1). With 1000 bp the means are better than with 500 bp. The means for are not close to the true _T. With low _T values is grossly overestimated. The upward bias in is huge with 500-bp data and shrinks when longer sequences are used (Table 1 and Table 2).

Sequences in two populations can be similar to each other due to migration or due to low variation caused by a small . The simulated datasets with low sometimes show by chance no variation at all in short sequences. In these cases one would expect the estimate of to be 0.0 and the estimate of to be indeterminate. The likelihood surface for such data is very flat and some of the runs (Table 1) failed to find a maximum. These cases were omitted from the tables. If the migration parameter is high, it is consistently underestimated with high . This is caused by our allowing only a limited number of migration events per genealogy in the Markov chain Monte Carlo runs, which were 1000 migration events per genealogy, so that solutions with higher numbers of migration events were not proposed. Our method delivers similar estimates for and similar or better estimates for than an F_ST-based estimator (Table 3).

View this table: In this window In a new window   Table 3. FST-based estimator with three parameters ( 1 2, 1, 2, see also Appendix 1)

The number of cases in which the F_ST-based estimator (Table 3) is usable is often dramatically lower than that with the ML estimator (Table 1; see also Appendix 1). For example, with _T = 0.001 and _T = 10 the F_ST-based estimator could be used in only 45 out of 100 runs. This behavior seems not to improve with longer sequences. Increasing the number of sites and the number of loci helps to reduce the variance of the ML estimates (Table 2 and Figure 4). Estimates of from datasets with very long sequences are biased and are smaller than _T. This is the result of low acceptance rates during a run, when the starting genealogy is so well defined by the data that in our simulation runs the chains were too short and too few different migration scenarios were tried. The estimation of nonsymmetrical true parameters works quite well (Table 4). We encounter similar biases as with symmetrical parameters: with low the estimates are biased upward and with high they are biased downward.

View larger version (20K): In this window In a new window Download PPT slide   Figure 4. Example of estimation of population parameters for two populations with 1, 3, and 10 loci. The graphs are cross-sections through the parameter space passing through the peak of the likelihood surface. The axes are on a logarithmic scale. The dashed lines give the true values of the parameter values: 1 = 0.05, 2 = 0.005, 1 = 10.0, 2 = 1.0. The gray area is an approximate 95% confidence set based on a likelihood-ratio test criterion defined as the range of parameter values with log-likelihood values equal to or higher than the maximum -0.5 x 9.4877 (2d.f. = 4; = 0.05).

The assumption that mutation rate varies according to a gamma distribution adds more noise to the estimation. All biases are similar to the ones shown in the tables. The estimation of the shape parameter 1/ is certainly not very precise with only a few loci (Figure 5). One of the runs shown in Figure 5 obviously contains almost no information about 1/ and its log-likelihood curve is nearly flat, so that almost any value of 1/ can be accepted for this dataset. The increase of the log-likelihood values with very small 1/ and far from the maximum is due to imprecisions in the calculation using the discrete gamma approximation. In evaluations using a more exact but much slower quadrature scheme (SIKORSKI et al. 1984 ) the peaks of the log-likelihood curves are at similar values, but log-likelihood for very low 1/ values are monotonically decreasing and not, as with the discrete gamma approximation, slightly increasing.

View larger version (14K): In this window In a new window Download PPT slide   Figure 5. Log-likelihood curves of shape parameter 1/. The true 1/ is 1.0. The curves were evaluated with and parameters from the maximum-likelihood estimate. Left, 10 runs with 10 loci each. Right, sum of the log-likelihood curves over 10 runs with 10 loci each.

	DISCUSSION

TOP ABSTRACT MODEL IMPLEMENTATION SIMULATION RESULTS DISCUSSION APPENDIX 1 LITERATURE CITED

The estimation of the migration parameter seems to be rather imprecise, and this is true for both migrate and F_ST (e.g., Table 1 and Table 3). The standard deviations for are large for single-locus estimates. More sites per locus give better estimates but the improvement of the variance is small, and given the obstacles of sequencing long stretches of DNA for several individuals one might better invest in the investigation of an additional unlinked locus. Each new locus has its own genealogy that is not correlated with that of the other loci and so increases the power of estimating parameters more, as we can see in Table 2.

There is considerable bias in the estimates of when datasets either have no or little genetic variation or have very high variation. Datasets with almost no variation inflate and produce an upward bias. The likelihood surfaces become very flat because of the lack of variation, and almost any migration parameter value is possible. These problems can be overcome by adding more variation, which for low true population size _T means either adding more base pairs or adding more individuals. With high _T and high true migration parameter _T, the estimates are biased downward because of our need to truncate the number of migration events in the reconstructed genealogies at some high arbitrary number. If we did not bias against high numbers of migration events, the program could add more and more migration events and would eventually crash the computer by running out of memory or would run too long.

In single runs we can also see a "fatal attraction" to 0.0 if the true parameters are very close to 0.0. The genealogy sampler is using the current parameter values to propose coalescences and migration events. If one or more of these parameters are very close to 0.0, then it can be seen by inspecting (2) that as a result our procedure will either most often propose an immediate coalescent event if one of the is close to 0.0 or rarely propose an immigration event with very low . If a chain never proposes any instances of an event, its rate of occurrence will be estimated to 0.0, and this situation may persist indefinitely, a fatal attraction. Even if the parameters are prevented from becoming exactly zero this means that it may take large amounts of simulation to escape these values.

Because long sequences define the genealogy very well, our Markov chain Monte Carlo sampler will often reject newly proposed genealogies and only a few different genealogies are sampled for the parameter estimation. Therefore, the program does not readily explore the whole migration-genealogy space and tends to stick to the good starting genealogy. This generates a bias in because we have started with a genealogy that has a minimal number of migration events on it (Table 2). We may be able to overcome this mixing problem by adding a "heating" scheme to our importance sampling. Currently the only way to overcome this is to run the chains much longer than the defaults (for practical guidance see BEERLI 1997 ).

Summing over loci assuming that the mutation rate is constant between loci delivers much better estimates of the parameters than single-locus estimation. This model is rather unnatural for real data, especially microsatellite data or electrophoretic data. Summing over loci using a gamma-distributed mutation rate and estimating the shape parameter 1/ increases the variation of the parameters but allows more realistic estimation of the parameters with several loci.

Conclusion:
Our method based on coalescents with migration delivers similar or better estimates than the F_ST method based on the expectations derived by MAYNARD SMITH 1970 , which for certain data also produces fine results. Qualitative comparison with other methods than the F_ST measure, for example RANNALA and HARTIGAN 1996 , or the cladistic approach of SLATKIN and MADDISON 1989 (see also HUDSON et al. 1992 , their Table 1) shows that all these methods have a tendency to overestimate when is not very high and all methods have quite large variances. The maximum-likelihood estimators have smaller variances than the other approaches. But the biggest drawback of the F_ST methods compared to ML methods seems to be their inability to estimate population size and migration rate for data where the homozygosity between populations happens to be equal to or less than the homozygosity within a population (see also HUDSON et al. 1992 ).

Direct comparison of the different methods is currently difficult because each one estimates a different number of parameters or uses a different migration model. We expand our model to accept more populations and different migration schemes. This should then facilitate direct comparison. We believe that taking the history of mutation into account and summing over all possible genealogies gains the most information from the data. We believe that our approach and the approaches of NATH and GRIFFITHS 1996 and Bahlo and Griffiths (http://www.maths.monash.edu.au/~mbahlo/mpg/gtree.html) will produce better estimates than other approaches.

The methods described here are available in the program Migrate over the World Wide Web at http://evolution.genetics.washington.edu/lamarc.html. We distribute it as C source code and also as binaries for several platforms including PowerMacintoshes and Windows 95/98 for Intel-compatible processors.

	ACKNOWLEDGMENTS

We thank Mary K. Kuhner and Jon Yamato for many discussions, much help in debugging, and comments on the manuscript. This work was funded by National Science Foundation grant no. BIR 9527687 and National Institutes of Health grant no. GM51929, both to J.F., and in part by a Swiss National Funds fellowship to P.B.

Manuscript received January 5, 1998; Accepted for publication February 23, 1999.

	APPENDIX 1

TOP ABSTRACT MODEL IMPLEMENTATION SIMULATION RESULTS DISCUSSION APPENDIX 1 LITERATURE CITED

Derivation of the importance sampling function:
We want to calculate Equation 1 but would need to sum over all possible genealogies. This function can be transformed into an importance sampling function by assuming that L() is an expectation and we sample from a distribution whose density is g instead of the correct density, f,

(A1)

(A2)

(A3)

Suppose that

(A4)

(A5)

and

(A6)

because

(A7)

then we have

(A8)

so that

(A10)

where

The expectation can be estimated by its average over the simulation

(A11)

In Markov chain Monte Carlo approaches, the goal is to sample from the posterior and concentrate the sampling on regions with higher probabilities (HAMMERSLEY and HANDSCOMB 1964 ; CHIB and GREENBERG 1995 ; KASS et al. 1998 ). In our scheme we are approximating the target function Prob(G|) Prob(D|G) with Prob (G|₀) Prob(D|G); this may be a rather crude approximation when ₀ is very different from , but after several updating chains the target and sampling distribution are very similar. When and ₀ are very similar, our approach is nearly optimal if it is run long enough to avoid problems of autocorrelation. Sampling from the prior alone, Prob(G|) for example, is very inefficient as was shown by J.F. in 1988 (J. FELSENSTEIN, unpublished data).

Data models:
Our migration estimation divides naturally into two parts: calculation of Prob(G|) and Prob(D|G). This makes it easy to implement models for different kinds of data: any data model influences only the latter calculation, which is the genealogy likelihood calculation.

For sequences we are using the model of change originated by one of us (J.F.) as implemented in PHYLIP 3.2 in 1984, described in KISHINO and HASEGAWA 1989 and described as F84 by SWOFFORD et al. 1996 . It is a variant of the Kimura two-parameter model, which allows for different transition and transversion ratios and variable base frequencies. In migrate the model is the same as in PHYLIP 3.5 and includes also rate variation among sites (FELSENSTEIN and CHURCHILL 1996 ), but the inclusion of rate variation between sites was not tested.

For microsatellite data, we have implemented a one-step mutation model in which the probability of making a net change of i steps in time interval u is

(A12)

In addition to the stepwise mutation model a Brownian motion approximation is available; this fast approximation to the exact model is described elsewhere.

Our "infinite" allele model is approximated with a (k + 1)-allele model. The observed alleles are = (a₁, a₂, ... , a_k). All unobserved alleles are pooled into a_k+1.

(A13)

F_ST calculations:
We use a calculation based on F⁽ⁱ⁾_W, the homozygosity within the population i and F_B, the homozygosity between populations. For sequences, F_W and F_B can be calculated using the heterozygosity calculations described in SLATKIN and HUDSON 1991 . The original recurrence equations F_W and F_B were developed by MAYNARD SMITH 1970 and MARUYAMA 1970 . We use an equation based on NEI and FELDMAN 1972 . The exact equations are simplified by removing quadratic terms like µ², m², and µm and divisions by number of individuals in a population (e.g., m/N). For two populations in equilibrium we get the equation system

(A14)

where µ is the mutation rate, m_i is the immigration rate into population i, and N_i is the subpopulation size. To fit these calculations into our framework we replace 4Nµ with and m/µ with . We cannot solve for the four parameters we would need for an accurate comparison with our likelihood method. We can solve the equation system (A14) by assuming either that the population sizes are the same in both subpopulations (SN) or that the migration rates are symmetric (SM). With the SN approach solving for _{1 2}, ₁, and ₂ we get

(A15)

The SM approach with ₁, ₂, and _{1 2} results in

(A16)

Comparing (A15) and (A16) one can recognize that the estimation of _i in (A15) breaks down when F_B F⁽ⁱ⁾_W. In (A16) the can be reasonably estimated only if 2F_B F⁽¹⁾_W + F⁽²⁾_W. For the comparison with our maximum-likelihood estimator (Table 3) we used SN because of our emphasis on estimation of migration rates. There are differences between the two F_ST methods even when the true parameters are symmetrical (data not shown), but we do not have the impression that one F_ST method works better than the other over the range we are using for the comparison with our migrate method (Table 5).

	LITERATURE CITED

TOP ABSTRACT MODEL IMPLEMENTATION SIMULATION RESULTS DISCUSSION APPENDIX 1 LITERATURE CITED

BEERLI, P., 1997 Migrate documentation version 0.3. Distributed over the Internet: http://evolution.genetics.washington.edu/lamarc.html.

CHIB, S. and E. GREENBERG, 1995  Understanding the Metropolis-Hastings algorithm. Am. Stat. 49:327-335.

DAHLQUIST, G., and A. BJÖRK, 1974 Numerical Methods. Prentice-Hall, New York.

FELSENSTEIN, J., 1973  Maximum likelihood estimation of evolutionary trees from continuous characters. Am. J. Hum. Genet. 25:471-492[Medline].

FELSENSTEIN, J., 1988  Phylogenies from molecular sequences: inference and reliability. Annu. Rev. Genet. 22:521-565[Medline].

FELSENSTEIN, J., 1992a  Estimating effective population size from sample sequences: a bootstrap Monte Carlo integration method. Genet. Res. 60:209-220[Medline].

FELSENSTEIN, J., 1992b  Estimating effective population size from samples of sequences: inefficiency of pairwise and segregating sites as compared to phylogenetic estimates. Genet. Res. 59:139-147[Medline].

FELSENSTEIN, J., 1993 PHYLIP: phylogenetic inference package version 3.5c. Distributed over the Internet: http://evolution.genetics.washington.edu/phylip.html.

FELSENSTEIN, J. and G. A. CHURCHILL, 1996  A hidden Markov model approach to variation among sites in rate of evolution. Mol. Biol. Evol. 13:93-104[Abstract].

FITCH, W. M., 1971  Toward defining the course of evolution: minimum change for a specified tree topology. Syst. Zool. 20:406-416.

GRIFFITHS, R. and S. TAVARÉ, 1994  Sampling theory for neutral alleles in a varying environment. Philos. Trans. R. Soc. Lond. Ser. B Biol. Sci. 344:403-410[Medline].

HAMMERSLEY, J., and D. HANDSCOMB, 1964 Monte Carlo Methods. Methuen and Co., London.

HASTINGS, W., 1970  Monte Carlo sampling methods using Markov chains and their applications. Biometrika 57:97-109[Abstract/Free Full Text].

HUDSON, R., 1983  Properties of a natural allele model with intragenic recombination. Theor. Popul. Biol. 23:183-201[Medline].

HUDSON, R., M. SLATKIN, and W. MADDISON, 1992  Estimation of levels of gene flow from DNA sequence data. Genetics 132:583-589[Abstract].

HUDSON, R. R., 1990 Gene genealogies and the coalescent process, pp. 1–44 in Oxford Surveys in Evolutionary Biology, Vol. 7, edited by R. DAWKINS and M. RIDLEY. Oxford University Press, Oxford, U.K.

KASS, R. E., B. P. CARLIN, A. GELMAN, and R. M. NEAL, 1998  Markov chain Monte Carlo in practice: a roundtable discussion. Am. Stat. 52:93-100.

KINGMAN, J., 1982a  The coalescent. Stoch. Proc. Appl. 13:235-248.

KINGMAN, J., 1982b On the genealogy of large populations, pp. 27–43 in Essays in Statistical Science, edited by J. GANI and E. HANNAN. Applied Probability Trust, London.

KISHINO, H. and M. HASEGAWA, 1989  Evaluation of the maximum likelihood estimate of the evolutionary tree topologies from DNA sequence data, and the branching order in hominoidea. J. Mol. Evol. 29:170-179[Medline].

KUHNER, M., J. YAMATO, and J. FELSENSTEIN, 1995  Estimating effective population size and mutation rate from sequence data using Metropolis-Hastings sampling. Genetics 140:1421-1430[Abstract].

KUHNER, M., J. YAMATO, and J. FELSENSTEIN, 1998  Maximum likelihood estimation of population growth rates based on the coalescent. Genetics 149:429-434[Abstract/Free Full Text].

MARUYAMA, T., 1970  Effective number of alleles in a subdivided population. Theor. Popul. Biol. 1:273-306[Medline].

MAYNARD SMITH, J., 1970  Population size, polymorphism, and the rate of non-darwinian evolution. Am. Nat. 104:231-237.

METROPOLIS, N., A. ROSENBLUTH, M. ROSENBLUTH, A. TELLER, and E. TELLER, 1953  Equation of state calculations by fast computing machines. J. Chem. Phys. 21:1087-1092.

NATH, H. and R. GRIFFITHS, 1993  The coalescent in two colonies with symmetric migration. J. Math. Biol. 31:841-851[Medline].

NATH, H. and R. GRIFFITHS, 1996  Estimation in an island model using simulation. Theor. Popul. Biol. 50:227-253[Medline].

NEI, M. and M. W. FELDMAN, 1972  Identity of genes by descent within and between populations under mutation and migration pressures. Theor. Popul. Biol. 3:460-465[Medline].

NOTOHARA, M., 1990  The coalescent and the genealogical process in geographically structured population. J. Math. Biol. 29:59-75[Medline].

RANNALA, B. and J. HARTIGAN, 1996  Estimating gene flow in island populations. Genet. Res. 67:147-158[Medline].

SIKORSKI, K., F. STENGER, and J. SCHWING, 1984  Algorithm 614: A FORTRAN subroutine for numerical integration in H_p.. ACM Transact. Math. Soft. 10:152-160.

SLATKIN, M., 1985  Rare alleles as indicators of gene flow. Evolution 39:53-65.

SLATKIN, M., 1991  Inbreeding coefficients and coalescence times. Genet. Res. 58:167-175[Medline].

SLATKIN, M. and R. HUDSON, 1991  Pairwise comparisons of mitochondrial DNA sequences in stable and exponentially growing populations. Genetics 129:555-562[Abstract].

SLATKIN, M. and W. MADDISON, 1989  A cladistic measure of gene flow inferred from the phylogenies of alleles. Genetics 123:603-613[Abstract/Free Full Text].

SWOFFORD, D., G. OLSEN, P. WADDELL and D. HILLIS, 1996 Phylogenetic inference, pp. 407–514 in Molecular Systematics, edited by D. HILLIS, C. MORITZ and B. MABLE. Sinauer Associates, Sunderland, MA.

TAKAHATA, N., 1988  The coalescent in two partially isolated diffusion populations. Genet. Res. 52:213-222[Medline].

TAKAHATA, N. and M. SLATKIN, 1990  Genealogy of neutral genes in two partially isolated populations. Theor. Popul. Biol. 38:331-350[Medline].

TUFTO, J., S. ENGEN, and K. HINDAR, 1996  Inferring patterns of migration from gene frequencies under equilibrium conditions. Genetics 144:1911-1921[Abstract].

WAKELEY, J., 1998  Segregating sites in Wright's island model. Theor. Popul. Biol. 53:166-174[Medline].

YANG, Z., 1994  Maximum likelihood phylogenetic estimation from DNA sequences with variable rates over sites: approximative methods. J. Mol. Evol. 39:306-314[Medline].

This article has been cited by other articles:

				C. R. Bonvicino, P. R. Goncalves, J. A. de Oliveira, L. F. B. de Oliveira, and M. S. Mattevi Divergence in Zygodontomys (Rodentia: Sigmodontinae) and Distribution of Amazonian Savannas J. Hered., December 9, 2008; (2008) esn105v1. [Abstract] [Full Text] [PDF]

				M. Bonhomme, S. Cuartero, A. Blancher, and B. Crouau-roy Assessing Natural Introgression in 2 Biomedical Model Species, the Rhesus Macaque (Macaca mulatta) and the Long-Tailed Macaque (Macaca fascicularis) J. Hered., October 30, 2008; (2008) esn093v1. [Abstract] [Full Text] [PDF]

				J. I. Schmidt, K. J. Hundertmark, R. T. Bowyer, and K. G. McCracken Population Structure and Genetic Diversity of Moose in Alaska J. Hered., October 3, 2008; (2008) esn076v1. [Abstract] [Full Text] [PDF]

				C. P. Burridge, D. Craw, D. Fletcher, and J. M. Waters Geological Dates and Molecular Rates: Fish DNA Sheds Light on Time Dependency Mol. Biol. Evol., April 1, 2008; 25(4): 624 - 633. [Abstract] [Full Text] [PDF]

				D. Garrigan, S. B. Kingan, M. M. Pilkington, J. A. Wilder, M. P. Cox, H. Soodyall, B. Strassmann, G. Destro-Bisol, P. de Knijff, A. Novelletto, et al. Inferring Human Population Sizes, Divergence Times and Rates of Gene Flow From Mitochondrial, X and Y Chromosome Resequencing Data Genetics, December 1, 2007; 177(4): 2195 - 2207. [Abstract] [Full Text] [PDF]