A Continuous Method for Gene Flow | Genetics

A Continuous Method for Gene Flow

Michal Palczewski and Peter Beerli
Genetics July 1, 2013 vol. 194 no. 3 687-696; https://doi.org/10.1534/genetics.113.150904

Abstract

Most modern population genetics inference methods are based on the coalescence framework. Methods that allow estimating parameters of structured populations commonly insert migration events into the genealogies. For these methods the calculation of the coalescence probability density of a genealogy requires a product over all time periods between events. Data sets that contain populations with high rates of gene flow among them require an enormous number of calculations. A new method, transition probability-structured coalescence (TPSC), replaces the discrete migration events with probability statements. Because the speed of calculation is independent of the amount of gene flow, this method allows calculating the coalescence densities efficiently. The current implementation of TPSC uses an approximation simplifying the interaction among lineages. Simulations and coverage comparisons of TPSC vs. MIGRATE show that TPSC allows estimation of high migration rates more precisely, but because of the approximation the estimation of low migration rates is biased. The implementation of TPSC into programs that calculate quantities on phylogenetic tree structures is straightforward, so the TPSC approach will facilitate more general inferences in many computer programs.

THE estimation of population genetics parameters such as migration rates and effective population sizes is a common task for researchers in such fields as conservation biology, population biology, and biogeography. The theory of coalescence, introduced in 1982 by Kingman (1982a,b,c), is a formidable framework for describing population genetic processes.

It has changed the inference of population genetic parameters completely. We can calculate probabilities of complex interactions among individuals within and between populations, using the structured coalescent (Strobeck 1987; Notohara 1990; Wilkinson-Herbots 1998). Probabilistic inferences built on the structured coalescent (Kuhner et al. 1995; Kuhner 2006; Beerli 1998, 2006; Beerli and Felsenstein 1999; Hey 2010) are now used by many researchers. Routinely, complex population models are evaluated and, more recently, compared to each other (Beerli and Palczewski 2010). These approaches commonly integrate over many genealogies G that are augmented with migration or divergence events, using the Felsenstein equationEmbedded Image(1)(Hey 2007), where D is the data and Embedded Image is a set of model parameters, for example the effective population size Ne and immigration rates m. Beerli and Felsenstein (1999) expressed the coalescence probability density of a genealogy given the parametersEmbedded Image(2)with η number of events on the tree. The rate at which the zth event happens isλz=∑jnkjz(kjz−1)4Nj+∑in∑j,j≠inkjzmij,(3)where kjz is the number of lineages currently in population j corresponding to the time before event z, and mij is a migration rate defined as the percentage of individuals in population j that were previously in i. The variable βz is the contribution of the current event to the sum that makes λ. In other words, βz is the rate of the event considered. This rate of coalescence is 2/4Nj for a given pair of lineages and the rate of migration is mij for a given lineage.

This method allows for n2 parameters, where n is the number of populations; the parameters can be partitioned into n population sizes and n(n − 1) migration rates, thus allowing for asymmetric migration rates.

Often, we will not be able to estimate the absolute quantities of Ne and m, but only the parameter Θ, which is 4 × Ne × μ, and M, which is m/μ. For both Θ and M the mutation rate μ is the scalar.

Equation 2 is a potentially large product over all events in the genealogy, including coalescences and migration events. The state space for such augmented genealogies is potentially huge because the number of events depends on the magnitude of the parameters. For example, a low migration parameter suggests that there are few migration events in the genealogy whereas a large migration rate suggests that there are many (Figure 1).

Figure 1
Figure 1

Number of migration events in genealogies. (A) Genealogy generated with Nm = 0.400 into the population marked with open circles (○) and Nm = 0.267 into the population marked with solid circles (●). (B) Immigration rates are 10 times higher. Migration events on the genealogy are shaded according to the receiving population, looking forward in time.

The calculation of the likelihood p(D|Embedded Image) is analytically intractable and is commonly solved using Markov chain Monte Carlo (MCMC) methods (Metropolis et al. 1953; Hastings 1970). This can be very time consuming because the Markov chain needs to visit not only large numbers of probable topologies and parameter sets but also an even larger number of different configurations of migration events. Particularly, data sets that were generated by models with high migration rates among subsets of populations are difficult to analyze.

Here we propose a method that reduces the integration over all of these different migration events. Instead of relying on Monte Carlo methods to simulate many of these events, a one-dimensional numerical integration is proposed. This greatly simplifies the number of possible tree topologies that need to be explored. Although for any data stemming from multiple populations, there are an infinite number of possible genealogies augmented by migration events, the number of possible topologies when migration events are excluded is large but finite. Furthermore these genealogies are much simpler, since they include only coalescences. The analysis of such genealogies requires less time for situations with high migration rates where the standard methods augment the genealogies with many migration events (Figure 1).

Methods

Transition-probability structured coalescence framework

Our new framework, the transition-probability structured coalescence (TPSC), does not depend on explicit migration events, but integrates over all possible population assignments. We contrast TPSC with the event-based structured coalescence (ESC) presently incorporated into MIGRATE (Beerli 1998; Beerli and Felsenstein 1999). Although the TPSC allows for complex population structure, we describe the method using a simple two-population model with four parameters (Figure 2).

Figure 2
Figure 2

Two-population model.

Assume that there is a single stretch of nonrecombining genome L0; at the present time it is in population 1. Looking backward in time, there is an exponential distribution for the waiting time until this lineage migrates from a different population.

The probability density of the waiting time until the sample changes population one or more times during the time interval from 0 to t ism21e−m21t (4)with the immigration rate m21 from population 2 to 1; t is measured in generations and m is measured in terms of the proportion of offspring coming from a new population. A similar function can be applied to a sample from the other population.

To predict the probability of a particular lineage Li being in a particular population Zi we use a continuous-time Markov process. First construct a transition rate matrix Q of migration rates,Q=[−m21m21m12−m12] (5)and a vector of initial probabilitiesP0=[P(L∈Z1 |t0)P(L∈Z2 |t0)].(6)Now we can compute the probabilities of being in each population at time t:[P(L∈Z1 |t)P(L∈Z2 |t)]=P0eQt. (7)This framework can be extended to more than two populations: Q would still be a square matrix of migration rates, but Q would have size n, the number of populations:Q=[−∑i=1nmi1m21…mn1m12−∑i=1nmi2…mn2⋮⋮⋱⋮m1nm2n…∑i=1nmin].(8)With this framework it is possible to compute the probability density of one lineage going back in time. When looking at multiple lineages, one must also take into account coalescence events. The rate of standard coalescence for two lineages is the inverse of the population size or two times the population size for diploids. The probability density of the time t to coalescence of two independent lineages in the same population with no migration isP(t)=12Nee(−(t/2Ne)). (9)Two lineages that are not in the same population do not coalesce. Their rate of coalescence is zero. Calculating the probability that two lineages are in the same population at a specific time would require a conditional probability. This would increase the size of the Q matrix, which would include both lineages and possible coalescences. Instead we make a simplification and estimate the joint probability by assuming independence.

Thus, we can combine the probability of being in a particular population with the rate of coalescence to estimate the rate of two independent lineages coalescing in population Zp; from now on we mark Zp only by its indicator p,λ1,2,p(t)=P(L1∈p,L2∈p)2Np≈P(L1∈p)P(L2∈p)2Np, (10)where Np is the effective population size of population p. The total rate of coalescence of the two lineages is the sum over all K populations:λ1,2(t)=∑k=1Kλ1,2,k(t)=∑k=1KP(L1∈k)P(L2∈k)2Nk. (11)Expanding to multiple lineages, the total rate of coalescence isλ(t)=∑i=1n∑j≠inλi,j(t)2=∑k=1K∑i=1n∑j≠inP(Li∈k)P(Lj∈k)4Nk. (12)The 2 in the divisor offsets the double counting of the coalescence of λi,j and λj,i; n is the total number of all sampled lineages. For computational efficiency we transform toλ(t)=∑k=1K[14Nk∑i=1nP(Li∈k)(Kk−P(Li∈k))] (13)withKk=∑j=1nP(Lj∈k). (14)Disregarding the time it takes to calculate individual P(Li ∈ k), both Equations 13 and 14 can be calculated in O(nk) time.

The probability that a specific coalescent of two lineages has happened in a particular population can be calculated as the ratio of the rate that lineages coalesce in that population to the total coalescence rate,P(coalescence∈p|Li,Lj,t)=λi,j,p(t)λi,j(t). (15)With this framework it is possible to calculate the probability of an entire genealogy given the population sizes and migration rates. The probability of each coalescent event is modeled by a nonhomogeneous Poisson process. Therefore the probability of two lineages Lx and Ly coalescing at time t isP(Lx,Ly,t)=λx,y(t)e(−∫t0tλ(t)dt). (16)Here x and y are the indexes of the lineages in question.

Multiplying all coalescence probabilities results in the probability of the genealogy G given the model parameters. For our two-population model we getP(G|N1,N2,m21,m12)=∏in−1P(Li,x,Li,y,ti). (17)Here Li,x and Li,y represent the ith coalescent even on the tree where lineages x and y coalesce.

Testing the TPSC

To evaluate the merit of our approach, we evaluated the TPSC for three different situations: We calculated exact probabilities for two individuals collected in two different populations. We calculated the maximum-likelihood estimates of model parameters and compared coverage and parameter estimates of a Bayesian implementation of TPSC with MIGRATE for various simulated data sets.

Likelihood calculations

The likelihood of the genetic data D given the parameters is calculated using the Felsenstein et al. (1999) equationp(D|N1,N2,m21,m12,Mμ)=∑Gp(G|N1,N2,m21,m12) p(D|G,Mμ). (18)For the mutation model Mμ we used the F84 model (Felsenstein and Churchill 1996). Without additional information the population size parameters and the mutation rate are confounded and we express the parameters of interest as a combination of μ and a scalar, so that for diploid organisms we reportEmbedded Image(19)where Θi is the mutation-scaled effective population size and Mji is the mutation-scaled immigration rate.

Bayesian inference using TPSC

We construct a Bayesian estimatorEmbedded Image(20)The marginal posterior density for the parameters was estimated using the Metropolis–Hastings (MH) method. The implementation of such a method uses updates on the genealogy and the population genetic model parameters (Ronquist and Huelsenbeck 2003; Drummond and Rambaut 2007).

We implemented an MH algorithm, using a tree-update method similar to the one described by Nielsen (2000). The tree is updated by picking a random internal node representing a coalescence event and changing the time of the event up or down on the genealogy. In our algorithm, the probability of choosing any coalescence event is uniform, whereas in Nielsen’s algorithm the coalescence event selection is proportional to the length of a branch away from the root. The distance that each internal node is moved is a random value drawn from a normal distribution as in Nielsen’s algorithm, but unlike Nielsen’s algorithm the variance for this normal distribution is not arbitrary but is adapted to the information content of the data during the burn-in period (Appendix).

For parameter updates we use a method similar to the sliding-window proposal implemented in Mr. Bayes (Huelsenbeck et al. 2001; Ronquist and Huelsenbeck 2003). Unlike Mr. Bayes’ sliding-window proposal, which uses a uniform random number, we update the parameter by adding a normally distributed random variable. The variance of the normally distributed random variable is also adapted to information content of the data during the burn-in period. Our adaptive scheme is outlined in the Appendix.

Results

To analyze the effectiveness of our new method we have done three types of analysis. The first is an analytic treatment of two simple cases. We take a look at the probability density of time until a coalescent event. For a simple case, we can solve this analyticaly and compare the exact solution to the TPSC approximation. In the second study we simulate genealogies and use TPSC to infer the parameters used to generate these genealogies. Knowing all details of a genealogy is a rather unrealistic scenario. However, this second study tests the new model directly and without the complication of a mutation model needed to fit data to the genealogy. Finally, we did full simulation tests using DNA sequence data. We compared the ability of TPSC to the program MIGRATE, which uses a discrete coalescent method, to infer the simulated parameters.

Analysis for two lineages

Symmetric model:

First, we analyzed the structured coalescent of a two-population model with identical population sizes (N) and symmetrical migration rates (m). At the present time there are two lineages of interest, one in each population. This can be modeled by a continuous-time Markov model with the following exact transition probability matrix:Qe=[0001N−1N−2m2m02m−2m].(21)There are three states: State 3, represented by the third row, is the initial state of the lineages being in different populations. Looking backward in time, each lineage can migrate at the rate m. Either lineage migrating will result in both lineages existing in the same population. State 2, represented by the second row, is the state of both lineages existing in the same population. Either lineage can immigrate at the rate m, per lineage, or the two lineages can coalesce at the rate 1N. State 1, represented by the first row, is an absorbing state. Once the lineages are coalesced we are no longer interested in them. The probability density of time to coalescence is the derivative of the probability that the lineage is in state 1:p(tcoal=t)=ddtP(tcoal<t)=ddt(eQet)(3,1). (22)An analytic solution for this matrix exponential and derivative exists. However, the equation is very long and inconsequential. Instead of writing it out we have plotted it in Figure 4, but we have included it as a Mathematica worksheet with Supporting Information, File S1.

This simple two-population model analyzed using TPSC leads to the transition probability matrix that takes into account only migration events:Qm=[−mmm−m].(23)The first step requires the calculation of the probability that the two lineages are in the same population (Ptogether). This probability is the sum of probabilities that both lineages are in population 1 and that both lineages are in population 2:Ptogether(t,m)=e(1,1)Qmte(1,2)Qmt+e(2,1)Qmte(2,2)Qmt. (24)This is a function of m because Q depends on m. The rate of coalescence then becomesλ(t,N,m)=Ptogether(t,m)N. (25)Finally we can compute the probability density of a coalescent event. Again this is analytically tractable, but the equation is rather long, and we have included it as a Mathematica worksheet in File S1:p(tcoal)=λ(t,N,m)e−∫0xλ(x,N,m)dx. (26)We have plotted Equations 26 and 22 for various values of Nm in Figure 4. Although both N and m can vary independently, the shapes of these curves depend only on the ratio of N to m.

Asymmetric model:

In the first analytic example we created a symmetric model. In this section we explore another simplified model, one with unidirectional rather than symmetric migration.

We simplify the model from Figure 2 and consider only a two-parameter model. The parameters are the population size of population 1, Ne(1), and the immigration rate m1→2; the immigration rate m2→1 is zero. The population size of population 2 is inconsequential. This model is shown in Figure 3.

Figure 3
Figure 3

Population model with two parameters.

Just as before, two individuals were sampled, one in each population. We are interested in calculating the probability density of time until coalescence. This simple scenario can be modeled by a continuous-time Markov process. The state probabilities can be calculated exactly, using a continuous-time Markov model with a three-state Q matrix:Qf=[0001Ne−1Ne00m−m].(27)Here state 1 represents the coalesced state. This is an absorbing state. State 3 is the initial state with each sample in a different population. Since migration is a one-way state the Markov chain will go from state 3 to state 2 at the migration rate. State 2 represents both lineages being in the same population. These will coalesce at a rate that is the inverse of the population size.

The exact probability density of the time to coalescence can be calculated asp(tcoal)=ddtP(tcoal<t)=ddt(eQft)(3,1)=me−mt−me−(t/Ne)1−mNe. (28)Using TPSC, first we compute the probability that these two populations are in the same population. This is governed by a simple exponential distribution, because there is an exponential waiting time until the lineage that is able to migrates:Ptogether=1−e−mt. (29)The rate of coalescence can be computed:λ(t)=1−e−mtNe. (30)This is the probability of both lineages being in the same population. Then the probability density function becomesp(tcoal)=λ(t)e−∫0tλ(x)dx=(1−e−mt)e− (1/Ne)[(e−mt−1)/m+t]1Ne.(31)Comparisons between the exact method and TPSC, shown in Figure 4 for a symmetric migration model, reveal that the approximation works well in scenarios when the migration rate is high (Nm ≥ 1.0) and poorly when the migration rate is low (Nm ≤ 1.0). Graphs for the asymmetric case reveal the same general pattern (not shown, but included in File S2).

Figure 4
Figure 4

Graphs showing the probability density of time to coalescence of two lineages in a two-population scenario with symmetric migration. The dashed line is the exact probability density whereas the solid line is the TPSC approximation. The effective population size for each panel is Ne = 1000.

Simulated genealogies

To test our method we simulated genealogies from known population parameters. Using the true genealogy is equivalent to assuming that there is an infinite amount of sequence data to define the genealogy; therefore we can find the maximum-likelihood estimate of Equation 18 (cf. Felsenstein 1992). An example of such an analysis is shown in Figure 5. Each panel presents the profile-likelihood curve for each of the four parameters of a two-population model: Θ1, Θ2, M21, and M12. The genealogy was generated using the structured coalescent with parameters Θ1 = 0.012, Θ2 = 0.01, M21 = 0, and M12 = 1000. The 95% confidence intervals bracket the true parameter value for all parameters. The profile likelihood curves are strongly peaked for the mutation-scaled population sizes, but the migration parameters have wide confidence intervals.

Figure 5
Figure 5

Plots of profile-likelihood curves. Data were simulated from a two-population model with migration in one direction. Labels indicate the maximum-likelihood estimate, the “true” parameter value used to generate the genealogy, and the 95% confidence interval of the estimate.

We calculated several statistics over the maximum-likelihood estimates (MLEs) from 1000 simulated genealogies of 40 individuals, 20 per population (Table 1).

View this table:
Table 1 Accuracy of TPSC

Simulated DNA sequence data

To test the effectiveness of the TPSC we simulated DNA sequence data from two populations for a total of 40 individuals. We examined all nine combinations of three mutation-scaled population sizes Θ of 0.001, 0.01, and 0.1 and three mutation-scaled immigration rates of 10, 100, and 1000. The smallest population size is typical for nuclear data in human populations whereas the largest population size seems appropriate for species with very large effective population sizes, such as viruses or bacteria. The number of migrants 4Nem per generation ranged from 0.01 to 100, covering many potential natural scenarios. For each of the nine scenarios we simulated 100 data sets, using the simulation software MIGTREE and MIGDATA (available at http://people.sc.fsu.edu/∼pbeerli/software). DNA sequences with lengths of 500 bp were simulated using the F84 model (Hasegawa et al. 1985; Felsenstein and Churchill 1996). We chose for our simulations a DNA sequence length short enough so that even in natural populations we could expect few or no recombination events to occur.

These data sets were then run in TPSC and MIGRATE. Comparison with other programs that estimate migration rates (IMa and LAMARC) failed because of run-time constraints. Either programs did not converge within 48 hr or memory requirements were prohibitive to run 900 simulations.

TPSC and MIGRATE were run on the high-performance computing cluster at Florida State University. The run time of each separate data set was on the order of a few hours. Convergence was assessed by running TPSC multiple times from random starting genealogies on the same data to check for similar results. This procedure was then repeated using MIGRATE. Convergence of the runs of MIGRATE was assessed by repeated runs; there were potential convergence problems for data sets generated with high numbers of migrants (4Nem = 100).

Table 2 summarizes standardized mean square errors (MSE) [(1/n)∑in(x^i−xt)2/xt2] for TPSC and MIGRATE of n = 100 replicates for each set of Θ and M. Although we used a symmetric model of migration for simulation, the inference used a model with two population sizes and two migration rates that were allowed to vary independently. We report all four estimated parameters Θ^1, Θ^2, M^2→1, and M^1→2 for each combination of the true parameters, resulting in 36 comparisons of TPSC and MIGRATE. Because the true values for these parameters are symmetric, we expect that Θ^1=Θ^2=Θt and M^2→1=M^1→2=Mt.

View this table:
Table 2 Mean square errors of TPSC and MIGRATE

TPSC and MIGRATE performed similarly on the estimation of mutation-scaled effective population sizes; differences of the MSE were mostly small, although TPSC estimates usually with slightly higher MSE values. The standardized MSEs for M are larger than those for Θ for both programs. TPSC outperforms MIGRATE in the estimation of mutation-scaled migration most of the time (37 of 54). In particular TPSC’s MSE of the median of M with low true effective population size and high true migration rates is smaller than the corresponding MSEs of MIGRATE. With large population sizes (Θ = 0.1) and large migration rates (M = 1000), TPSC seems to have difficulties achieving good estimates. These are cases in which the number of migrants is so high that distinguishing large from very large values becomes difficult. The likelihood surface becomes very flat, making it difficult to get accurate estimates. Although MIGRATE seems to work better in these cases, convergence to a unique solution for a particular data set becomes difficult.

Discussion

Recently, several researchers have described similar methods to TPSC. Takahata (1988) and Hobolth et al. (2011) integrated out migration events similarly to the method described in this article, but their formulation requires much larger transition probability matrices to calculate all potential interactions among lineages. This makes it difficult to employ their methods for large numbers of individuals k. TPSC, in contrast, depends only on the number of populations n in the analysis and work increases on the order O(n3k2) instead of O((kn)k). Usually, k >> n. The program BEAST (Lemey et al. 2009) contains a phylogeographic model that may relate distantly to our method in that it presents probabilities of origin for particular pathogen strains or populations. The model of Lemey et al. (2009) similarly uses a continuous-time Markov chain to calculate location probabilities of past events. In this model, however, migration and coalescence are not intertwined. Instead a coalescent prior with a single population is used for the entire genealogy. Afterward, the locations of past states are computed on this tree. This does not take into account that individuals in small populations coalesce faster than those in large populations. In contrast, TPSC takes into account multiple population sizes, which gives information on spatial location of coalescent events.

TPSC is an approximation; it assumes independence of lineages for the calculation of the population assignment probability for the nodes in the genealogy. This leads to biased estimates for low migration rates (Figure 5, Table 1); however, TPSC outperforms event-based methods such as MIGRATE in scenarios with high immigration rates and moderate population sizes (Table 2). In such scenarios immigration events happen similarly as often as coalescence events (cf. Nordborg and Krone 2002) . With moderate immigration numbers (Nm ∼ 1) the TPSC approximation and the full solution lead to similar distributions, suggesting that TPSC can replace event-based methods for all data sets except those that include isolated populations.

We distribute our method in a stand-alone program (http://people.sc.fsu.edu/∼pbeerli/software) and will incorporate it into our program MIGRATE, allowing for switching between event-based and transition-probability structured coalescence methods.

Acknowledgments

We thank Thomas Uzzell for comments on several revisions of our text. We acknowledge the use of the high-performance computing facility at Florida State University. Our work was supported by grants DEB-0822626 and DEB-1145999 from the National Science Foundation.

Appendix

An Adaptive Scheme

The Metropolis–Hastings algorithms in this program adapt themselves to the data to ensure faster convergence. For Metropolis–Hastings the ideal acceptance rate can differ from 20% to 60% (Gelman et al. 1997; Roberts and Rosenthal 1998; Roberts and Rosenthal 2009). In a typical MCMC algorithm relatively small updates to a parameter will be accepted at a high rate. If a parameter does not change much, a likelihood and prior value will vary by only a small amount. On the other hand, a large change in a parameter when the value is already close to optimal is much more likely to be rejected.

We use the following scheme to adjust the variance of proposal distributions to adjust our acceptance ratio to a theoretical ideal. During burn-in, whenever a value is accepted for a parameter, the variance is increased by multiplying it by a value B that is slightly >1.0,σt+12=Bσt2,(A1)with t as the step number. Whenever a value is rejected, the proposal variance is decreased by a small value b that is slightly smaller than 1.0:σt+12=bσt2. (A2)If we assume that σ2 has converged to some value, then we can formulate a relation of B, b, and the acceptance rate R. This relation isB1−R=bR. (A3)In our algorithm we choose to tune our acceptance ratio as closely as possible to the ideal R = 0.44 proposed by Roberts and Rosenthal (2009). We use an arbitrary value of b = 0.99, thus ensuring that our variance is at most 1% away from the ideal variance, and solve for B. Values of b close to 1 will converge to a value closer to the ideal, although they will converge more slowly. Conversely, values of b that are farther away from 1 will converge more quickly but the final variance could be farther from the ideal. The convergence rate is exponential and thus the desired acceptance ratio can be found quickly during the burn-in.

In Figure A1 we show the convergence of a typical run to the ideal variance. We have not seen any examples where the convergence did not happen less quickly. The variance converged very early in the burn-in. It should also be noted that any errors in convergence do not result in an incorrect algorithm. Instead the result would be worse mixing and a longer run time required during the MCMC chain.

Figure A1
Figure A1

An example of the proposal variance adapting to an ideal. The acceptance rate is cumulative and has an asymptote at 0.44.

Footnotes

  • Communicating editor: M. A. Beaumont

  • Received February 26, 2013.
  • Accepted April 30, 2013.

Literature Cited

  1. ↵
    Beerli, P., 1998 Estimation of migration rates and population sizes in geographically structured populations, pp. 39–53 in Advances in Molecular Ecology, NATO Science Series A: Life Sciences, Vol. 306, edited by G. Carvalho. IOS Press, Amsterdam.
  2. ↵
    1. Beerli P.
    , 2006 Comparison of Bayesian and maximum likelihood inference of population genetic parameters. Bioinformatics 22: 341–345.
    OpenUrlAbstract/FREE Full Text
  3. ↵
    1. Beerli P.,
    2. Felsenstein J.
    , 1999 Maximum-likelihood estimation of migration rates and effective population numbers in two populations using a coalescent approach. Genetics 152: 763–773.
    OpenUrlAbstract/FREE Full Text
  4. ↵
    1. Beerli P.,
    2. Palczewski M.
    , 2010 Unified framework to evaluate panmixia and migration direction among multiple sampling locations. Genetics 185: 313–326.
    OpenUrlAbstract/FREE Full Text
  5. ↵
    1. Drummond A.,
    2. Rambaut A.
    , 2007 Beast: Bayesian evolutionary analysis by sampling trees. BMC Evol. Biol. 7: 214.
    OpenUrlCrossRefPubMed
  6. ↵
    1. Felsenstein J.
    , 1992 Estimating effective population size from sample sequences: A bootstrap Monte Carlo integration method. Genet. Res. 60: 209–220.
    OpenUrlCrossRefPubMedWeb of Science
  7. ↵
    1. Felsenstein J.,
    2. Churchill G. A.
    , 1996 A hidden Markov Model approach to variation among sites in rate of evolution. Mol. Biol. Evol. 13: 93–104.
    OpenUrlAbstract
  8. ↵
    Felsenstein, J., M. K. Kuhner, J. Yamato, and P. Beerli, 1999 IMS Lecture Notes-Monograph Series, pp. 163–185 in Statistics in Molecular Biology and Genetics: Likelihoods on coalescents: a Monte Carlo sampling approach to inferring parameters from population samples of molecular data, (Vol. 33), edited by Francoise Seillier-Moiseiwitsch. Institute of Mathematical Statistics and American Mathematical Society. Hayward, California.
  9. ↵
    1. Gelman A.,
    2. Gilks W. R.,
    3. Roberts G. O.
    , 1997 Weak convergence and optimal scaling of random walk Metropolis algorithms. Ann. Appl. Probab. 7: 110–120.
    OpenUrlCrossRef
  10. ↵
    1. Hasegawa M.,
    2. Kishino K.,
    3. Yano T.
    , 1985 Dating the human-ape splitting by a molecular clock of mitochondrial DNA. J. Mol. Evol. 22: 160–174.
    OpenUrlCrossRefPubMedWeb of Science
  11. ↵
    1. Hastings W. K.
    , 1970 Monte Carlo sampling methods using Markov chains and their applications. Biometrika 57: 97–109.
    OpenUrlAbstract/FREE Full Text
  12. ↵
    1. Hey J.
    , 2007 A model in two acts: a commentary on ‘A model of detectable alleles in a finite population’ by Timoko Ohta and Motoo Kimura. Genet. Res. 89: 365–366.
    OpenUrlCrossRefPubMed
  13. ↵
    1. Hey J.
    , 2010 Isolation with migration models for more than two populations. Mol. Biol. Evol. 27: 905–920.
    OpenUrlAbstract/FREE Full Text
  14. ↵
    1. Hobolth A.,
    2. Andersen L. N.,
    3. Mailund T.
    , 2011 On computing the coalescence time density in an isolation-with-migration model with few samples. Genetics 187: 1241–1243.
    OpenUrlFREE Full Text
  15. ↵
    1. Huelsenbeck J.,
    2. Ronquist F.,
    3. Nielsen R.,
    4. Bollback J.
    , 2001 Bayesian inference of phylogeny and it’s impact on evolutionary biology. Science 294: 2310–2314.
    OpenUrlAbstract/FREE Full Text
  16. ↵
    1. Kingman J.
    , 1982a The coalescent. Stoch. Proc. Appl. 13: 235–248.
    OpenUrlCrossRef
  17. ↵
    1. Kingman J. F. C.
    , 1982b Exchangeability and the evolution of large populations: proceedings of the international conference on exchangeability in probability and statistics, pp. 97–112 in Exchangeability in Probability and Statistics, edited by Koch G., Spizzichino F.. North-Holland Publishing, Amsterdam.
  18. ↵
    1. Kingman J. F. C.
    , 1982c On the genealogy of large populations. J. Appl. Probab. 19A: 27–43.
    OpenUrlCrossRef
  19. ↵
    1. Kuhner M.
    , 2006 Lamarc 2.0: maximum likelihood and Bayesian estimation of population parameters. Bioinformatics 22: 768–770.
    OpenUrlAbstract/FREE Full Text
  20. ↵
    1. Kuhner M. K.,
    2. Yamato J.,
    3. Felsenstein J.
    , 1995 Estimating effective population size and mutation rate from sequence data using Metropolis-Hastings sampling. Genetics 140: 1421–1430.
    OpenUrlAbstract/FREE Full Text
  21. ↵
    1. Lemey P.,
    2. Rambaut A.,
    3. Drummond A. J.,
    4. Suchard M. A.
    , 2009 Bayesian phylogeography finds its roots. PLoS Comput. Biol. 5: e1000520.
    OpenUrlCrossRefPubMed
  22. ↵
    1. Metropolis N.,
    2. Rosenbluth A. W.,
    3. Rosenbluth N.,
    4. Teller A. H.,
    5. Teller E.
    , 1953 Equation of state calculation by fast computing machines. J. Chem. Phys. 21: 1087–1092.
    OpenUrlCrossRefWeb of Science
  23. ↵
    1. Nielsen R.
    , 2000 Estimation of population parameters and recombination rates from single nucleotide polymorphisms. Genetics 154: 931–942.
    OpenUrlAbstract/FREE Full Text
  24. ↵
    Nordborg, M., and S. M. Krone, 2002 Separation of time scales and convergence to the coalescent in structured populations, pp. 194–232 in Modern Developments in Theoretical Population Genetics: The Legacy of Gustave Malécot, edited by M. Slatkin and M. Veuille. Oxford University Press, Oxford.
  25. ↵
    1. Notohara M.
    , 1990 The coalescent and the genealogical process in geographically structured population. J. Math. Biol. 29: 59–75.
    OpenUrlPubMedWeb of Science
  26. ↵
    1. Roberts G. O.,
    2. Rosenthal J. S.
    , 1998 Optimal scaling of discrete approximations to langevin diffusions. J. R. Stat. Soc. Ser. B Stat. Methodol. 60: 255–268.
    OpenUrlCrossRef
  27. ↵
    1. Roberts G. O.,
    2. Rosenthal J. S.
    , 2009 Examples of adaptive MCMC. J. Comput. Graph. Stat. 18: 349–367.
    OpenUrlCrossRef
  28. ↵
    1. Ronquist F.,
    2. Huelsenbeck J. P.
    , 2003 Mrbayes 3: Bayesian phylogenetic inference under mixed models. Bioinformatics 19: 1572–1574.
    OpenUrlAbstract/FREE Full Text
  29. ↵
    1. Strobeck C.
    , 1987 Average number of nucleotide differences in a sample from a single subpopulation: a test for population subdivision. Genetics 117: 149–153.
    OpenUrlAbstract/FREE Full Text
  30. ↵
    1. Takahata N.
    , 1988 The coalescent in two partially isolated diffusion populations. Genet. Res. 52: 213–222.
    OpenUrlCrossRefPubMedWeb of Science
  31. ↵
    1. Wilkinson-Herbots H. M.
    , 1998 Genealogy and subpopulation differentiation under various models of population structure. J. Math. Biol. 37: 535–585.
    OpenUrlCrossRefWeb of Science
Back to top

PUBLICATION INFORMATION

Volume 194 Issue 3, July 2013

Genetics: 194 (3)

ARTICLE CLASSIFICATION

INVESTIGATIONS
Population and evolutionary genetics
View this article with LENS
Email
Print
Alerts
View PDF

A Continuous Method for Gene Flow

Michal Palczewski and Peter Beerli
Genetics July 1, 2013 vol. 194 no. 3 687-696; https://doi.org/10.1534/genetics.113.150904
del.icio.us logo Digg logo Reddit logo Twitter logo CiteULike logo Facebook logo Google logo Mendeley logo
Citation