On Computing the Coalescence Time Density in an Isolation-With-Migration Model With Few Samples

To determine the likelihood we must determine the density f(t|Θ) for coalescent of two samples in the IM model. Wang and Hey (2009) find the time to the MRCA by explicitly integrating over all possible sample paths in the system. The purpose of this letter is to demonstrate that the time to the MRCA is easy to compute from a matrix exponential. Furthermore, the matrix exponential framework has the advantage that it generalizes to more than two samples. We describe the IM model, briefly describe the solution to computing the coalescence time density from Wang and Hey (2009), and finally present an approach that computes the density through matrix exponentials. Illustration of a sample path with coalescent in population 2 at time t. In Figure 2 we illustrate the coalescent density in the two species IM model. We use the same parameters as in the simulation study in Wang and Hey (2009, Table 6): θ1 = 0.005, θ2 = 0.003, θA = 0.002, m1 = 50, m2 = 100, and T = 0.003 (the vertical line). Illustration of the density (4) and (5) for coalescent in the two species isolation with migration model. Multiple analytical approaches for computing the coalescence time density can be found in the literature. Variations of the IM model have been analyzed using Laplace transforms in, e.g.,Latter (1973), Takahata (1995), and Wilkinson-Herbots (1998). Their results can also be derived using matrix exponentiation. In Wakeley (1996), a spectral decomposition was used to obtain a continuous-time approximation to a discrete time model. Generalizations of the IM model dealing with two loci with recombination between loci are analyzed using expressions for continuous time Markov chains in Slatkin and Pollack (2006) and Simonsen and Churchill (1997). We finally emphasize that the matrix exponential framework described above generalizes to more than two samples and more than two populations. For any coalescence system that can be expressed as a finite-state homogeneous continuous time Markov chain, we can compute the density of coalescence times using matrix exponentiation. Expressing a coalescence process as such a system is straightforward although a major complication is an explosion in the number of states when the number of samples and populations increase. If the Markov chain is not homogeneous, that is, the rate matrix Q depends on the time parameter t, simple matrix exponentiation is no longer a solution to the coupled set of differential equations. The model of Innan and Watanabe (2006), for instance, consists of the same set of states and almost the same rate matrix, but has the migration rates depend linearly on the time variable t. For this system, the approach described above cannot be applied.