Approximate Bayesian computation for Markovian binary trees in phylogenetics

Markovian binary trees (MBTs) are a flexible class of continuous-time branching processes, where each individual in a population (here, each species in a family) exists in one of $n$ phases. In this paper, we consider only the simplest case $n=2$. An individual in phase $i$ can:

• •

  transition to phase $j\neq i$ at constant rate $q_{ij}$;

• •

  give birth to a child in the same phase $i$ at constant rate $\lambda_{i}$ (remaining in the same phase); and

• •

  die at constant rate $\mu_{i}$.

The asymptotic behaviour of the population size is exponentially controlled by the growth rate $\omega$, which can be calculated from these parameters; see Supplementary Section S1.1 for more details, and for a more general formulation of the MBT process.

To infer the diversification rates of an MBT model from a set of $m$ observed phylogenetic trees, we apply the ABC-PMC method (Beaumont et al., 2009). This method proceeds over a series of iterations; at each iteration a standard ABC method is applied where parameters are proposed from a prior distribution, a dataset is simulated using the proposed parameters, and the parameters are accepted if the simulated dataset is ‘similar’ to the observed dataset when comparing various summary statistics. In the ABC-PMC method, the accepted parameters in one iteration form the basis for the prior in the next iteration, allowing the posteriors to converge towards the true values. We refer to Supplementary Section S1.2 for more details.

The effectiveness of the ABC-PMC method crucially relies on the careful selection of an appropriate suite of summary statistics to compare the observed and simulated trees. These statistics must collectively be able to capture variation in each of the parameters of the model. The statistics we use are as follows.

1. 1.

   Average branch length.

2. 2.

   Tree height (time to most recent common ancestor). Note that our simulated trees have a fixed number of extant species, but variable height. In situations where the tree height is fixed, the number of extant species may be used instead.

3. 3.

   Normalised lineage-through-time (nLTT) curve (Janzen et al., 2015). This curve shows the growth of the number of species against time, normalised by both the total time and the total number of extant species. To compute the distance between two nLTT curves, we calculate the absolute area between them:

   $$\Delta\text{nLTT}=\int_{0}^{1}|\text{nLTT}_{1}(t)-\text{nLTT}_{2}(t)|\,dt.$$

   We set the initial tolerance value for this statistic to be $0.05$. Due to the computational time required to generate a sample that can provide a close enough nLTT curve, we only use this statistic after the 20th iteration.

4. 4.

   Colless balance index (Colless, 1982). For a tree $T$, this is defined as the sum over all internal nodes $v\in V_{\text{int}}(T)$ of the absolute difference between the number of left and right descendant leaves. In other words, if $v$ has children $\{v_{l},v_{r}\}$, and $\kappa_{u}$ is the number of descendant leaves from node $u$, then the Colless index of $T$ is

   $$C(T)=\sum_{v\in V_{\text{int}}(T)}|\kappa_{v_{l}}-\kappa_{v_{r}}|.$$

5. 5.

   Balance index for each phase. To capture imbalance in the shape of the tree in each phase, we propose a novel extension of the Colless balance index to each phase; if $\kappa^{(i)}_{u}$ is the number of descendant leaves in phase $i$ from node $u$, we define the balance index for phase $i$ as

   $$C_{i}(T)=\sum_{v\in V_{\text{int}}(T)}|\kappa_{v_{l}}^{(i)}-\kappa_{v_{r}}^{(i)}|.$$

6. R1.

   Tree height for trees that contain phase-1 leaves (reducible case only). In the reducible case, species in phase 2 cannot transition to phase 1. Thus a tree with some leaves in phase 1 indicates either a large birth rate for phase 1, or a small transition rate from phase 1 to phase 2. These two situations are distinguished by the overall tree height; thus, we use the tree height for the subset of trees that contain at least one leaf in phase 1.

7. R2.

   Average branch length above phase-1 leaves (reducible case only). In the reducible case, all branches above leaves in phase 1 must themselves be in phase 1. Thus we use the average length of these branches as a summary statistic, which captures the birth and death rates for phase 1.

8. I1.

   Proportion of phase-1 leaves (irreducible case only). The proportion of leaves in phase 1 captures the variation between the growth rates of phase-1 lineages versus phase-2 lineages.

9. I2.

   Transition statistic for each phase (irreducible case only). We propose the transition statistic of a tree $T$ for phase $i$ as

   $$S_{i}(T)=\sum_{v\in V^{\prime}_{\text{int}}(T)}-\frac{\ln P_{v}^{i}}{d_{v}+\gamma},\qquad i=1,2,$$

   where $V^{\prime}_{\text{int}}(T)$ is the set of internal nodes with descendant leaves in both phases, $P_{v}^{i}$ is the proportion of phase-$i$ leaves in the descendant leaves of $v$, $d_{v}$ is the time of node $v$, and $\gamma$ is a small constant (here set to $0.02$).

   This statistic can be intuitively justified by considering that its dominant terms correspond to subtrees with small height (small $d_{v}$). We can approximate such a subtree by considering the paths from the root to each leaf as independent; if we assume that the root of the subtree is not in phase $i$, then the proportion of leaves in phase $i$ is approximately the probability of a transition to phase $i$ on a branch of length $d_{v}$, i.e., $e^{-q_{.i}d_{v}}$. Inverting this formula for the transition rate to phase $i$ gives the transition statistic for phase $i$. (We add a small constant $\gamma$ to the denominator to avoid extremely small $d_{v}$ values from dominating the statistic.)

We first consider the case where the MBT is reducible and the death rates are equal, i.e., $q_{21}=0,\mu:=\mu_{1}=\mu_{2}$. This results in four free parameters for the MBT model. We only consider supercritical processes (i.e., with growth rate $\omega>0$), so that the expected number of lineages grows without bound. Our observed dataset $Y_{\text{obs}}$ consists of 100 trees, each with 50 leaves. These trees are generated from the MBT model with the extinct subtrees removed.

In the ABC-PMC algorithm, we set the number of simulated trees for iteration $t$ as follows:

For earlier iterations, this enables us to quickly explore the parameter space and accept the required number of samples, which we set to $N=200$. For later iterations, we simulate (as required) datasets that are the same size as the observed dataset.

We first consider a ‘default’ set of parameters $(\lambda_{1},\lambda_{2},\mu,q_{12})=(1.5,0.51,0.15,0.5)$, which corresponds to a growth rate of $\omega=0.85$. These values are chosen to ensure that there is a noticeable difference between the birth rates of the two phases; we expect some parts of the trees will be in the fast-growing phase 1, while others will be in the slower-growing phase 2, creating a detectable imbalance.

In Supplementary Figure S11a, we show the posterior means at each iteration for one run of the algorithm. We see that the posterior means converge towards the true values as the number of iterations increases. The final posterior distributions are shown in Supplementary Figure S12; here we see that the distributions are closely concentrated around the true values. For all four parameters, the final posterior means are all within 4% of the true values. This suggests that the algorithm has achieved good convergence and can accurately recover the true parameter values.

Of particular significance is the capability to infer the true growth rate $\omega$ of the MBT process. In Supplementary Figure S11b, we illustrate (for a single run) the convergence of our estimates towards the true value as the number of iterations increases. Although the estimate is quite accurate, a minor positive bias is noticeable, consistently observed across different runs of the algorithm. This might be attributed to the fact that observed trees are conditioned on having 50 leaves (and hence are not extinct), resulting in a slightly higher effective growth rate.

We tested the accuracy of the method by varying the true parameters one at a time over a range of values, and repeated the inference (inferring all parameters) for 5 runs at each value. The results in Figure 1 show accurate inference of the true parameter values in all cases.

For most parameters, the accuracy of inference decreases as the parameter value increases (e.g., the 50% credible intervals become larger). The exception is the death rate $\mu$, where the accuracy appears relatively unaffected by the parameter value.

Next, we consider the general reducible case, where the death rates may differ between the two phases, resulting in five free parameters. The extra degree of freedom makes inference slightly more difficult. We consider the default parameters $(\lambda_{1},\lambda_{2},\mu_{1},\mu_{2},q_{12})=(1.5,0.51,0.7,0.15,0.5)$, where the parts of the tree in phase 1 both speciate and go extinct relatively quickly, while parts of the tree in phase 2 speciate and go extinct relatively slowly, and have an overall lower growth rate.

In Supplementary Figures S13 and S14, we show the trace plots and final posterior distributions for a single run. Similar to the equal death rate case, the posterior distributions show good convergence by the 30th iteration. Our inference is highly accurate for most parameters, with slight reductions in accuracy observed for $\mu_{1}$ and $q_{12}$. This decrease in accuracy could be attributed to the reducible nature of the process; as most phase-1 lineages eventually transition to phase 2, there is less information available in the data regarding the phase-1 death and transition rates. When assuming equal death rates, it is possible to infer the phase-1 death rate from phase-2 data, but this is no longer possible here. Nonetheless, our inference of the overall growth rate ($\omega=0.36$) remains accurate (Supplementary Figure S13d).

As above, we varied the true parameters one at a time over a range of values and repeated the inference 5 times for each parameter setting. The results, shown in Figure 2, indicate accurate inference, although there are some instances with noticeable errors.

As with the equal death rate case, the estimation of the phase-2 birth rate $\lambda_{2}$ and transition rate $q_{12}$ become less accurate as the parameter value increases, while the accuracy of the phase-2 death rate $\mu_{2}$ is largely unaffected by the parameter value. However, here there is a small positive bias for the phase-2 birth rate which increases with the parameter. In contrast, the estimation of the phase-1 birth rate $\lambda_{1}$ becomes more accurate as the parameter value increases, while the phase-1 death rate $\mu_{1}$ becomes less accurate. This can again be explained by the reducible nature of the process; the larger the ratio of phase-1 birth rate to death rate, the longer the tree spends in phase 1 before most lineages transition eventually to phase 2. This gives more information for phase 1, leading to better estimation of those rates.

We now consider the irreducible case, where transitions in both directions are allowed. We also allow death rates to differ between phases, resulting in six free parameters. For this case, we use the I1 and I2 summary statistics in place of the R1 and R2 statistics from Section 2.2. The remainder of the method is unchanged.

In order to verify that the summary statistics used are sensitive to the parameters of the MBT process, we visualise (in Supplementary Section S2) the variation of the summary statistics as the MBT parameters are varied. The statistics appear to have the necessary sensitivity for inference in the irreducible case.

In Supplementary Figure S15, we show the trace plots for a single run for the default parameters $(\lambda_{1},\lambda_{2},\mu_{1},\mu_{2},q_{12},q_{21})=(3,2,1,0.5,0.5,0.25)$. Again, it is clear that the process has converged by the end of the 30th iteration. In Supplementary Figure S16, we show the final posterior distribution. As in the general reducible case, the resulting posterior distributions are more dispersed for the death rates than the birth rates. Overall, our inference has high accuracy for all parameters, indicating the summary statistics used can detect the influence of different phases.

Figure 3 shows our results when each parameter is varied one at a time, with 5 replicates for each value. As with the general reducible case, the estimation of the phase-1 birth rate becomes more accurate as the rate becomes larger, while the opposite is true for the phase-2 birth rate. Again, there is a small positive bias for the phase-2 birth rate. The inference of both phase-1 and phase-2 death rates becomes less accurate as the rates increase; in addition, the phase-2 death rate is overestimated for small values and underestimated for larger values. Finally, as the transition rates increase, their accuracy decreases, and there is an underestimation of the phase 2 to phase 1 transition rate as it becomes large. Like most parameters, the accuracy of the growth rate decreases as the value increases, but there is no noticeable bias.

The general patterns observed here are mostly consistent with the reducible cases. While we can still estimate the rates satisfactorily, we note that the performance in this irreducible case is slightly less optimal compared to the simpler reducible cases, which is expected. Nonetheless, our method shows great potential in dealing with more complex scenarios.

We also studied the sensitivity of the inference for different tree sizes (number of leaves). In Supplementary Figure S17, we show the parameter inference for different tree sizes (with 25 trees in total) in the observed dataset. As expected, as the size of the trees increases (thus increasing the amount of information in the dataset), the inference becomes more accurate and less biased. This effect is less apparent, but still present, for the phase-1 rates.

In Supplementary Figure S18, we vary tree size while keeping the total number of leaves in all observed trees fixed at 5000 (thus nominally keeping the amount of information the same in the datasets). There is relatively little change compared to the previous case, although there appears to be marginally better estimation for datasets with more trees and a smaller number of leaves.

We compare the performance of our ABC method with the maximum likelihood estimates (MLEs) from the BiSSE model for the six-parameter irreducible case, as implemented in the castor package Louca and Pennell (2020).

We first consider the default parameters $(\lambda_{1},\lambda_{2},\mu_{1},\mu_{2},q_{12},q_{21})=(3,1,2,0.5,0.5,0.25)$ as in Section 3.2. As before, our observed dataset consists of 100 trees, each with 50 leaves. Because the observed dataset has multiple trees, but castor calculates the MLEs for a single tree, we average the MLEs for each parameter over the 100 trees to generate the final estimate. The results are shown in Figure 4 for 50 replicates. In all cases, it is clear that the ABC estimates are less biased than the MLEs. In Table 1, we show the relative root mean squared error (RRMSE) of the estimates. For all parameters except $\mu_{1}$, the ABC method outperforms the maximum likelihood method; the performance for $\mu_{1}$ appears to be the result of a larger variance in the ABC estimates overcoming a smaller bias.

To compare the performance of the methods over a wider range of parameters, we also varied the true parameters so that they were simultaneously drawn from the prior distributions Unif$(0,5)$, under the constraint that the MBT be supercritical. We then performed inference as before, for datasets of 100 trees with 50 leaves each.

The results for 100 replicates are shown in Figure 5. As before, we can observe a distinct bias for the MLEs; the birth rates $\lambda_{1}$ and $\lambda_{2}$ are inferred with positive bias, as are the transition rates $q_{12}$ and $q_{21}$. Additionally, the variance of the inferred transition rates increases as the parameter values increase. The numerical approximation of MLEs can give poor estimates and lead to unreasonably high transition rates (some transition rates were estimated to be greater than 100; these outliers are not shown in the figure). In contrast, the ABC estimates are significantly more precise and show no discernible bias. For the death rates $\mu_{1}$ and $\mu_{2}$, there appears to be a positive bias at low parameter values and a negative bias at high parameter values for both methods, as was observed for the ABC method in Section 3.2. However, this effect is again much less noticeable in the ABC estimates. In Table 2, we show the RRMSE for these estimates. In all cases, it appears that the ABC method substantially outperforms the ML method, particularly for the transition rates.

To verify that similar results apply to problem sizes closer to those in the real data analysis in the following section, we repeated this analysis using observed (and simulated) datasets consisting of a single tree with 5000 leaves. The results are shown in Supplementary Figures S19–S20, and the RRMSE is shown in Supplementary Tables S1-S2. The results are consistent with our previous experiments, showing ABC estimates that have higher accuracy than ML estimates.

Markovian binary trees (MBTs) form a flexible class of continuous-time Markovian multitype branching processes [Athreya and Ney, 1972] in which the lifetime and reproduction process of each individual in the population is controlled by an underlying transient Markov chain with $n$ transient states, called phases, and one absorbing state 0. While in phase $i\in\{1,\dots,n\}$, an individual can:

• •

  transition to phase $j\neq i$ at constant rate $(D_{0})_{ij}$;

• •

  give birth to a child in phase $k$ and simultaneously transition to phase $j$ at constant rate $B_{i,kj}$;

• •

  die, that is, enter the absorbing state 0, at constant rate $d_{i}$.

Thus an MBT is parameterised by an $n\times n$ matrix $D_{0}$, an $n\times n^{2}$ matrix $B$, and an $n\times 1$ vector $\mathbf{d}$. The diagonal elements of $D_{0}$ are negative and are such that

where $\mathbf{1}$ and $\mathbf{0}$ are vectors of 1’s and 0’s, respectively, of the appropriate size.

MBTs are ‘matrix’ generalisations of birth-and-death processes: in an MBT, the lifetime of individuals is distributed according to a phase-type distribution, which is a generalisation of the exponential distribution, and the reproduction process of the individuals follows a (transient) Markovian arrival process, which is a generalisation of the Poisson process (see Latouche and Ramaswami, 1999). Note that while the MBT literature usually considers each lineage as an individual organism, in this paper they represent species, where birth events represent speciations and death events represent extinctions.

It can be shown [Hautphenne et al., 2009] that the mean number of individuals in phase $j$ in the population at time $t$, given that the population started with a single individual in phase $i$ at time 0, is given by the $(i,j)$th entry of the mean population size matrix

where

and $\otimes$ denotes the Kronecker product. The asymptotic behaviour of the MBT therefore depends on the dominant eigenvalue of $A$, which we denote by $\omega$ and call the growth rate. There are three cases:

• •

  when $\omega<0$, the mean population size goes to 0 as $t\rightarrow\infty$, and the population eventually becomes extinct with probability 1 (subcritical);

• •

  when $\omega=0$, the asymptotic mean population size is constant (critical), and the population eventually becomes extinct with probability 1;

• •

  when $\omega>0$, the mean population size grows without bounds, and the population has a positive probability of ultimate survival (supercritical).

We are primarily interested in the supercritical case, which is suitable to represent the exponential increase in the variety of species observed in real life.

In this paper, we concentrate on the special case where $n=2$, that is, there are only two possible phases. We further assume that individuals stay in the same phase when giving birth, and give birth to children in that same phase. Thus the matrices $D_{0}$, $B$, and $\mathbf{d}$ are given by

where $\lambda_{1}$ and $\lambda_{2}$ are the birth rates for each phase, $\mu_{1}$ and $\mu_{2}$ are the death rates for each phase, and $q_{12}$ and $q_{21}$ are the transition rates between phases. We denote this set of parameters by $\theta=(\lambda_{1},\lambda_{2},\mu_{1},\mu_{2},q_{12},q_{21})$. We assume $q_{12}>0$. If individuals in phase 2 cannot transition to phase 1 (i.e., $q_{21}=0$), we call the process reducible; otherwise, we call it irreducible.

The ABC-PMC method [Beaumont et al., 2009] is an extension of the basic ABC rejection method [Beaumont et al., 2002]. We use as input to the method a set of $m$ ‘real’ observed phylogenetic trees $Y_{\text{obs}}$, with branch lengths given in units of substitutions per site (typically these would be inferred, potentially with some error, from sequence data). We assume that these trees are ultrametric and that the phases of the leaves are known (but the phases are unobserved in other parts of the tree).

The algorithm proceeds over a series of $T$ iterations. In the first iteration $t=0$, a set of parameters $\theta^{*}$ are proposed from a prior distribution $\bm{\pi}(\theta)$. In further iterations ($t\geq 1$), we adapt the posterior samples from the previous iteration to construct a new prior distribution as detailed in Algorithm 1, and then propose a set of parameters $\theta^{*}$ from this prior.

At iteration $t$, we simulate $n_{t}$ phylogenetic trees using an MBT process with parameters $\theta^{*}$, and calculate a set of summary statistics (detailed in Section 2.2 of the manuscript) for these trees. Given a set of tolerance values $\bm{\epsilon}_{t}=\{\epsilon_{tk}\}_{k=1}^{K}$, where $\epsilon_{tk}$ is the tolerance for the $k$th summary statistic, and $K$ is the total number of summary statistics, we accept the proposed parameters if the absolute difference between each summary statistic for the simulated and observed trees is less than the corresponding tolerance value; otherwise, we reject them. We continue until we have accepted $N$ samples; this results in a ‘population’ of accepted parameters $\{\theta^{(i)}_{t}\}_{i=1}^{N}$, which provides an approximation to the posterior distribution of the parameters.

To propose parameters for iteration $t\geq 1$, we draw one of the $N$ accepted parameter sets from the previous iteration, with weights $\{\omega^{(i)}_{t-1}\}_{i=1,\ldots,N}$. These weights are calculated at the end of iteration $t-1$ (details on how to calculate them at the end of the current iteration $t$ are provided below). We then perturb the sample by adding a normal perturbation $\mathcal{N}(0,\Sigma_{t-1})$, where $\Sigma_{t-1}$ is defined below, to produce the proposed parameters.

At the end of the iteration, we calculate the weights for the next iteration as follows. If $t=0$, we take $\omega_{0}^{(i)}=\frac{1}{N}$ for all $i\in\{1,2,\ldots,N\}$, so that all samples are weighted equally. If $t\geq 1$, we let

where $f(\,\cdot\,;{\mu},\Sigma)$ is the multivariate Gaussian density function, and $\Sigma_{t-1}$ is twice the weighted covariance matrix of the posterior samples from iteration $t-1$:

The weights $\{\omega^{(i)}_{t}\}_{i=1,\ldots,N}$ are then normalised so that their sum is 1.

The full formal details of the algorithm are given in Algorithm 1.

When we apply the ABC-PMC method to simulated data, we use the following values:

• •

  $T=30$ total iterations;

• •

  prior distribution Unif(0,5) for all parameters.

The initial tolerance values $\bm{\epsilon}_{0}$ are estimated from the real trees $Y_{\text{obs}}$. For each summary statistic, the corresponding tolerance value is calculated separately; here, we set them to be twice the standard deviation of the summary statistic across the trees in the observed dataset. To maintain a reasonable acceptance rate, subsequent tolerance values are a constant multiple (here we take $e^{-0.2}$) of the preceding values if the acceptance rate is above a specified threshold (we take $0.03$); otherwise, they are left unchanged.

For irreducible processes, considering the increase in the number of parameters, we use a wider acceptance region in the first iteration, which allows us to obtain the required number of accepted samples in a reasonable time. In the first iteration, if the summary statistics of the proposed sample stays within the region spanned by the minimal and maximal values of statistics from the observed trees, it will be accepted.

When proposing parameters, it is possible that the proposed parameters generate a critical or subcritical process, i.e., the growth rate $\omega\leq 0$. If this happens, we automatically reject the proposed parameters.

Likewise, it is possible that the proposed parameters generate a supercritical process, but one that is only ‘slightly’ so, i.e., has a high probability of going extinct. If we generate a tree that becomes entirely extinct before it reaches the desired number of leaves (50 in the default case), we discard it and generate another tree. If this occurs for the first 25 times in a row for a set of proposed parameters, we reject those parameters.