## 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.

- coalescent
- migration
- gene flow
- Markov chain Monte Carlo (MCMC)
- continuous Markov model
- nonhomogeneous Poisson process

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 equation(1)(Hey 2007), where *D* is the data and is a set of model parameters, for example the effective population size *N*_{e} and immigration rates *m*. Beerli and Felsenstein (1999) expressed the coalescence probability density of a genealogy given the parameters(2)with *η* number of events on the tree. The rate at which the *z*th event happens is*j* corresponding to the time before event *z*, and *m _{ij}* is a migration rate defined as the percentage of individuals in population

*j*that were previously in

*i*. The variable

*β*is the contribution of the current event to the sum that makes

_{z}*λ*. In other words,

*β*is the rate of the event considered. This rate of coalescence is

_{z}*m*for a given lineage.

_{ij}This method allows for *n*^{2} 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 *N*_{e} and *m*, but only the parameter Θ, which is 4 × *N*_{e} × *μ*, 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).

The calculation of the likelihood p(*D*|) 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).

Assume that there is a single stretch of nonrecombining genome *L*_{0}; 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* is*m*_{21} 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 *L _{i}* being in a particular population

*Z*we use a continuous-time Markov process. First construct a transition rate matrix

_{i}*Q*of migration rates,

*t*:

*Q*would still be a square matrix of migration rates, but

*Q*would have size

*n*, the number of populations:

*t*to coalescence of two independent lineages in the same population with no migration is

*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 *Z _{p}*; from now on we mark

*Z*only by its indicator

_{p}*p*,

*N*is the effective population size of population

_{p}*p*. The total rate of coalescence of the two lineages is the sum over all

*K*populations:

*λ*

_{i}_{,}

*and*

_{j}*λ*

_{j}_{,}

*;*

_{i}*n*is the total number of all sampled lineages. For computational efficiency we transform to

*P*(

*L*∈

_{i}*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,*L _{x}* and

*L*coalescing at time

_{y}*t*is

*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 get*L _{i}*

_{,}

*and*

_{x}*L*

_{i}_{,}

*represent the*

_{y}*i*th 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) equation*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 report(19)where Θ

*is the mutation-scaled effective population size and*

_{i}*M*is the mutation-scaled immigration rate.

_{ji}### Bayesian inference using TPSC

We construct a Bayesian estimator(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:*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

This simple two-population model analyzed using TPSC leads to the transition probability matrix that takes into account only migration events:*P*_{together}). This probability is the sum of probabilities that both lineages are in population 1 and that both lineages are in population 2:*m* because *Q* depends on *m*. The rate of coalescence then becomes*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, *m*_{1→2}; the immigration rate *m*_{2→1} is zero. The population size of population 2 is inconsequential. This model is shown in Figure 3.

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:

The exact probability density of the time to coalescence can be calculated as*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).

### 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}, *M*_{21}, and *M*_{12}. The genealogy was generated using the structured coalescent with parameters Θ_{1} = 0.012, Θ_{2} = 0.01, *M*_{21} = 0, and *M*_{12} = 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.

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

### 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 4*N*_{e}*m* 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 (4*N*_{e}*m* = 100).

Table 2 summarizes standardized mean square errors (MSE) [*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

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*(*n*^{3}*k*^{2}) 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* 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:*σ*^{2} has converged to some value, then we can formulate a relation of *B*, *b*, and the acceptance rate *R*. This relation is*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.

## Footnotes

*Communicating editor: M. A. Beaumont*

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

- Copyright © 2013 by the Genetics Society of America