Distribution of coalescent histories under the coalescent model with gene flow

We propose a coalescent model for three species that allows gene flow between both pairs of sister populations. The model is designed for multilocus genomic sequence alignments, with one sequence sampled from each of the three species, and is formulated using a Markov chain representation that allows use of matrix exponentiation to compute analytical expressions for the probability density of coalescent histories. The coalescent history distribution as well as the gene tree topology distribution under this coalescent model with gene flow are then calculated via numerical integration. We analyze the model to compare the distributions of gene tree topologies and coalescent histories for species trees with differing effective population sizes and gene flow rates. Our results suggest conditions under which the species tree and associated parameters are not identifiable from the gene tree topology distribution when gene flow is present, but indicate that the coalescent history distribution may identify the species tree and associated parameters. Thus, the coalescent history distribution can be used to infer parameters such as the ancestral effective population sizes and the rates of gene flow in a maximum likelihood (ML) framework. We conduct computer simulations to evaluate the performance of our method in estimating these parameters, and we apply our method to an Afrotropical mosquito data set (Fontaine et al., 2015).