Linear Constraints $$ \Psi $$ from Sampling Times

For any pair of nodes (i, j) (where each of i and j can either be a leaf or an internal node) with enforced divergence times (t_i, t_j), the following constraint $$ \psi (i,j)$$ must be satisfied

Let t₀ be the unknown divergence time at the root of the tree. For k calibration points $$ t_1,\dots ,t_k$$, we can setup k constraints of the form:

Optimization Criteria

Since $$ {\nu }_i=\mu {\tau }_i/{\widehat{b}}_i$$, the distribution of ν_i is influenced by both the distribution of the rates (μ_i) and the distribution of $$ {\widehat{b}}_i$$ around b_i. In traditional strict-clock models (Zuckerkandl 1962), a constant rate is assumed throughout the tree ($$ {\forall }_i{\mu }_i=\mu $$). Under this model, the distribution of ν_i is determined by deviations of $$ {\widehat{b}}_i$$ from b_i.

Langley and Fitch (1974) (LF) modeled the number of observed substitutions per sequence site on a branch i by a Poisson distribution with mean $$ \lambda =\mu {\tau }_i$$ and treated $$ s{\widehat{b}}_i$$ as if they were the total number of observed substitutions; as such, they assume $$ s{\widehat{b}}_i\sim \text{Poisson}(s\mu {\tau }_i)$$, where $$ s$$ is the sequence length. Therefore, by changing variable, we can write the log-likelihood function as:

Given $$ s$$ and $$ {\widehat{b}}_i$$, LF finds $$ \mathbf{x}$$ that maximizes the log-likelihood function and subject to the constraints $$ \Psi $$. As such,

To et al. (2016) assume $$ {\widehat{b}}_i$$ follows a Gaussian model: $$ {\widehat{b}}_i\sim \text{Gaussian}(\mu {\tau }_i,{\sigma }_i^2)$$ and approximate the variance by $$ {\widehat{b}}_i/s$$ (the method includes smoothing strategies omitted here). Then, the negative log-likelihood function can be written as:

Thus, the ML estimate can be formulated as:

Both LF and LSD have convex formulations. Langley and Fitch (1974) proved that their negative log-likelihood function is convex and thus the local minimum is also the global minimum. Our constraint-based formulation of LF also can be easily proved convex by showing its Hessian matrix is positive definite. To et al. (2016) pointed out their objective function is a weighted least squares. Using our formulation, we also see that equation (4) together with the calibration constraints form a standard convex quadratic optimization problem which has a unique analytical solution.

Motivation

LF only seeks to model the errors in $$ \widehat{\mathbf{b}}$$ and ignore true rate heterogeneity. Strict-clock assumption is now believed to be unrealistic in many settings (Schwartz and Maresca 2006; Pulquério and Nichols 2007; Ho 2014), motivating relaxed clocks, typically by assuming that μ_is are drawn i.i.d. from some distribution (Drummond et al. 2006; Akerborg et al. 2008; Volz and Frost 2017). Most methods rely on presumed parametric distributions (typically, LogNormal, Exponential, or Gamma) and estimate parameters using ML (Volz and Frost 2017), MAP (Akerborg et al. 2008), or MCMC (Drummond et al. 2006; Drummond and Rambaut 2007). The LSD method, which like LF directly models errors in $$ \widehat{\mathbf{b}}$$, is additionally justified under a Gaussian clock model. Choices of specific distributions in these methods are not motivated by the knowledge that real data follow them exactly (for example, the Normal distribution has to be misspecified as mutation rates cannot be negative).

Our goal is to avoid explicit parameter inference under a model of rate multipliers. Instead, we follow the assumption shared by existing methods like LSD and LF: we assume that given two solutions of $$ \mathbf{x}$$ both satisfying the calibration constraints, the solution with less variability in ν_i values is preferable. Thus, we prefer solutions that minimize deviations from a strict clock while allowing deviations. A natural way to minimize deviations from the clock is to minimize the variance of $$ {\tau }_i/{\widehat{b}}_i$$. This can be achieved by finding μ and all ν_i such that ν_i is centered at 1 and $$ {\sum }_{i=1}^{2n-2}{({\nu }_i-1)}^2$$ is minimized. Interestingly, the ML function used by LSD (eq. 4) is a weighted version of this approach.

The minimum variance principle results in a fundamental asymmetry: multiplying or dividing the rate of a branch by the same factor are penalized differently (fig 1a). For example, the penalty for $$ {\nu }_i=4$$ is more than ten times larger than $$ {\nu }_i=1/4$$. The LF model is more symmetrical than LSD but remains asymmetrical (fig 1a). This asymmetry results from the asymmetric distribution of the Poisson distribution around its mean, especially for small mean, in log scale (fig 1b). Because of this asymmetry, methods like LSD and LF judge a very small $$ {\widehat{b}}_i/b_i$$ to be within the realm of possible outcomes, and thus penalize $$ {\nu }_i<1$$ multipliers less heavily than $$ {\nu }_i>1$$.

Fig. 1.

(a) The penalty associated with multiplying a single edge i with multiplier ν_i in LSD, LF, and LogDate approaches, as shown in equations (3–5). To allow comparison, we normalize the penalty to be zero at ν  =  1 and to be 1 at ν  =  4. (b) The CI of the ratio between estimated and true branch length using the Poisson model. For this exposition, we assume that the estimated branch length equals the number of substitutions occurring on the branch and follows a Poisson distribution (i.e., JC69 model), divided by sequence length. With these assumptions, the CI for estimate length $$ {\widehat{b}}_i$$ is between $$ 1/2{\chi }_{2s{b}_i}^2$$ and $$ 1/2{\chi }_{2s{b}_i+2}^2$$; we draw the CI for $$ \alpha /2=0.05$$ and $$ \alpha /2=0.2$$ to get 0.2–0.8 and 0.05–0.095 intervals for $$ 0.0001\le b_i\le 0.4$$. (c) and (e): Density and histograms of penalty terms (without square) used by LSD ($$ \mu {\tau }_i/{\widehat{b}}_i-1$$) and LogDate ($$ \text{log}\hspace{0.17em}\mu {\tau }_i/{\widehat{b}}_i$$) under different clock models. (c) Fixing $$ \mu {\tau }_i=0.1$$, we draw 500,000 rate multipliers (r_i) from LogNormal, Gamma, or Exponential distributions with mean 1 and variance 0.16 for LogNormal and Gamma. For strict clock, r_i = 1. We then draw estimated branch length for each replicate i from the Normal distribution with mean $$ b_i=r_i\mu {\tau }_i$$ and variance $$ b_i/s$$ for s = 200. (e) The branch lengths are estimated from the sequences using PhyML from simulated sequences of To et al. (2016), as explained in the text. Parameters of rate multiplier distributions match part (c). We omit extremely short branches (<0.001) for better visualization. (d) and (f) The penalty of LSD and LogDate versus the empirical log-likelihood of estimated length for the models described in (c) and (e), respectively. To compute the empirical likelihood, we divide estimated branch lengths into small bins and the empirical likelihood of each bin is estimated as the frequency of the data assigned to it. Ideally, increasing likelihood should monotonically decrease the penalty and two points with similar likelihood should have similar penalties. See supplementary figure S2, Supplementary Material online for extended results.

Our method is based on a principle, which we call the symmetry of ratios: the penalty for multiplying a branch by a factor of ν should be no different than dividing the branch by ν. Note that this assertion can only be applicable to true variations of the mutation rate (i.e., ignoring branch length estimation error). We further motivate this principle with more probabilistic arguments below, but here we make the following case. If one considers the distribution of rate multipliers for various branches, absent of an explicit model, it is reasonable to assume that compared with an overall rate, branches rates are as likely to increase by a factor of ν as they are to decrease by a factor of ν. When this statement is true, we shall prefer a method that penalizes ν and $$ 1/\nu $$ identically. To ensure the symmetry of ratios, we propose taking the logarithm of the multipliers ν_i before minimizing their variance. Minimizing the variance of the rates in log-scale is the essence of our method. It achieves the symmetry, and, as we show below, a better correspondence between penalty and data likelihood.

Log-transformation has long been used to reduce data skewness before applying linear regression (Stynes et al. 1986; Keene 1995; Xiao et al. 2011; Clifford et al. 2013). In molecular dating, it can be argued that log-transformation is implicitly applied in the new version of RelTime (Tamura et al. 2018) where the geometric means between sister lineages replaced the arithmetic means in its predecessor. The improvement in the accuracy of RelTime encourages a wider use of log-transformation in molecular dating. Note that log-transforming the rate multipliers before minimizing their least squares penalty is identical to applying linear least squares after log-transformation of both time and the number of substitutions. In other fields, log-transformation has been used to make the least-squares method more robust to highly skewed distributions (Aban and Meerschaert 2004; Meaney et al. 2007).

LogDate Optimization Function

We formulate the LogDate problem as follows. Given $$ \widehat{\mathbf{b}}$$ and the set of calibration constraints described earlier, we seek to find

This objective function satisfies the symmetry of ratio property (fig. 1a). Since ν_i values are multipliers of rates around μ, if we assume μ is the mean rate, the LogDate problem is equivalent to minimizing the variance of the log-transformed rate multipliers (around their mean 1). The objective function only depends on ν_i; however, note that μ is still included in the constraints and therefore is part of the optimization problem. This setting reduces the complexity of the objective function and speeds up the numerical search for the optimal solution. Since the values of ν_i close to 1 are preferred in equation (5), the optimal solution would push μ towards the mean rate.

wLogDate Optimization Function

The simple LogDate formulation, however, has a limitation: by allowing rates to vary freely in a multiplicative way, it fails to deal with the varied levels of relative branch error; that is, the ratio of the estimated branch length to the true branch length ($$ {\widehat{b}}_i/b_i$$). As $$ {\widehat{b}}_i$$ is estimated from the sequences, the error of $$ {\widehat{b}}_i$$ is directly related to the variations in the number of substitutions occurred along the branch b_i. Let us assume sequences follow the Jukes and Cantor (1969) model, and let N_i be the total number of substitutions occurred along branch i on a sequence with length s. Under Juke-Cantor model, we have $$ N_i\sim \text{Poisson}(s\mu {\tau }_i)$$ and therefore, $$ \text{var}(N_i)=s\mu {\tau }_i$$. Therefore, the variance of the expected number of substitutions around the true branch length is $$ \text{var}(N_i/s{b}_i)=s\mu {\tau }_i/s^2{b}_i^2=1/b_{i}s$$. As figure 1b shows, when b_i is small, $$ N_i/s$$ can easily vary by several orders of magnitude around b_i. Furthermore, the distribution is not symmetric: drawing values several factors smaller than the mean is more likely than drawing values above the mean by the same factor. These analyses predict that the distribution of $$ {\widehat{b}}_i/b_i$$ depends strongly on b_i—with smaller b_i gives higher variance—and is not symmetric.

The variances of the relative error $$ {\widehat{b}}_i/b_i$$ is difficult to compute analytically due to the involvement of the sequence substitution model and the method to estimate $$ {\widehat{b}}_i$$, which are both unknown. Therefore, we instead use empirical analyses of the estimated branch lengths by PhyML to demonstrate our arguments. Consistent with our prediction, supplementary figure S7a and c, Supplementary Material online illustrate that the relative error $$ {\widehat{b}}_i/b_i$$ varies more in small branches and the distribution is not symmetric. These properties of the branch length estimates are not modeled in our LogDate formulation and we seek to incorporate them in a refined version of LogDate which will be described below.

Since the true branch length b_i is unknown, a common practice is to use the estimated $$ {\widehat{b}}_i$$ in place of b_i to estimate its variance as $$ 1/{\widehat{b}}_{i}s$$. This explains why both LF and LSD objective functions (eqs. 3 and 4) have a weight of $$ {\widehat{b}}_i$$ for each term of ν_i. Following the same strategy, we propose weighting each $$ {\hspace{0.17em}\text{log}\hspace{0.17em}}^2({\nu }_i)$$ term in a way that reduces the contribution of short branches to the total penalty, and thus allows more deviations in the log space if the branch is small (and is thus subject to higher error). Since we log-transform ν_i and pursue a model-free approach, explicitly computing the weights to cancel out the variations of relative error among the branches is challenging. However, since the weights should reflect the variance of $$ {\widehat{b}}_i/b_i$$ (logarithmic scale), they should monotonically increase with $$ {\widehat{b}}_i$$ (fig. 1b) to allow more variance for the relative errors in short branches than in long branches. We use $$ \sqrt{{\widehat{b}}_i}$$ as weights, a selection driven by simplicity and empirical performance (shown in the Results section).

The shortest branches require even more care. When the branch is very short, for a limited-size alignment, the evolution produces zero mutations with high probability. For these no-event branches, tree estimation tools report arbitrary small lengths (see supplementary fig. S7, Supplementary Material online), rendering $$ {\widehat{b}}_i$$ values meaningless for very small branches. To deal with this challenge, the r8s’s implementation of LF Sanderson (2003) collapses all branches with length $$ {\widehat{b}}_i<1/s$$. To et al. (2016) proposed adding a smoothing constant $$ c/s$$ to each $$ {\widehat{b}}_i$$ to estimate the variance of $$ {\widehat{b}}_i$$, where c is a parameter that the user can tune. Following a similar strategy, we propose adding a small constant $$ \stackrel{\sim }{b}$$ to each $$ {\widehat{b}}_i$$. We choose $$ \stackrel{\sim }{b}$$ to be the maximum branch length that produces no substitutions with probability at least $$ 1-\alpha $$ for $$ \alpha \in [0,1]$$. Recall that N is the total number of actual substitutions on a branch. Under the Jukes and Cantor (1969) model, it is easy to show that $$ \arg {\max }_{\stackrel{\sim }{b}}\mathrm{Pr}(N=0|b=\stackrel{\sim }{b})\ge 1-\alpha =-\frac{1}{s}\text{log}(1-\alpha )$$. We choose this value as $$ \stackrel{\sim }{b}$$ and set $$ \alpha =0.01$$ by default. Thus, we define the wLogDate as follows:

Solving the Optimization Problem

Both LogDate and wLogDate problems (eqs. 5 and 6) are nonconvex, and hence solving them is nontrivial. The problem is convex if $$ 0\le {\nu }_i\le e$$. For small clock deviation and small estimation error in $$ {\widehat{b}}_i$$, the ν_i values should be small so that the problem becomes convex with one local minimum. However, as $$ {\nu }_i\le e$$ convexity is not guaranteed, we have to rely on gradient-based numerical methods to search for multiple local minima and select the best solution we can find. To search for local minima, we use the SciPy solver with trust-constr Lalee et al. (1998) method. To help the solver work efficiently, we incorporate three techniques that we next describe.

Computing Jacobian and Hessian matrices analytically helps speedup the search. By taking the partial derivative of each ν_i, we can compute the Jacobian, J, of equation (6). Also, since equation (6) is separable, its Hessian H is a $$ (2n-2)\times (2n-2)$$ diagonal matrix. Simple derivations give us:

Sparse matrix representation further saves space and computational time. The Hessian matrix is diagonal, allowing us to store only the diagonal elements. In addition, the constraint matrix defined by $$ \Psi $$ is highly sparse. If all sampling times are given at the leaves, the number of nonzero elements in our $$ (n-1)\times (2n-1)$$ matrix is $$ O(n\hspace{0.17em}\text{log}\hspace{0.17em}n)$$ (supplementary claim 3, Supplementary Material online). If the tree is either caterpillar or balanced, the number of nonzeroes reduced to $$ \Theta (n)$$. Thus, we use sparse matrix representation implemented in the SciPy package. This significantly reduces the running time of LogDate.

Starting from multiple feasible initial points is necessary given that our optimization problem is nonconvex. Providing initial points that are feasible (i.e., satisfied the calibration constraints) helps the SciPy solver work efficiently. We designed a heuristic strategy to find multiple initial points given sampling times $$ t_1,\dots ,t_n$$ of all the leaves (as is common in phylodynamics).

We first describe the process to get a single initial point. We compute the root age t₀ and μ using root-to-tip regression (RTT) (Shankarappa et al. 1999). Next, we scale all branches of T to conform with $$ \Psi $$ as follow: let $$ m=\text{arg}{\text{min}}_i{t}_i$$ (breaking ties arbitrarily). Let d(r, i) denote the distance from the root r to node i and P(r, m) denote the path from r to m. For each node i in P(r, m), we set $$ {\tau }_i={\widehat{b}}_i(t_m-t_0)/d(r,m)$$. Then going upward from m to r following P(m, r), for each edge (i, j) we compute $$ t_j=t_i-{\tau }_i$$ and recursively apply the process on the clade i. At the root, we set t_m to the second oldest (second minimum) sampling time and repeat the process on a new path until all leaves are processed. Since RTT gives us μ, to find ν we simply set $$ {\nu }_i=\mu {\tau }_i/{\widehat{b}}_i$$.

To find multiple initial points, we repeatedly apply RTT to a set of randomly selected clades of T and scale each clade using the aforementioned strategy. Specifically, we randomly select a set S of some internal nodes in the tree and add the root to S. Then, by a postorder traversal, we visit each node $$ u\in S$$ and date the clade under u using the scaling strategy described above. We then remove the entire clade u from the tree but keep the node u as a leaf (note that the age of u is already computed) and repeat the process for the next node in S. The root will be the last node to be visited. After visiting the root, we have all the τ_i for all i. After having all the branches in time unit, we find $$ \mathbf{x}$$ to minimize either equation (5) or (6), depending on whether LogDate or wLogDate is chosen. In a tree of n leaves, we have $$ 2^{(n-1)}-1$$ ways to select the initial nonempty set S, giving us enough room for randomization.

Computing Confidence Interval

With the ability of wLogDate to work on any combination of sampling times/calibration points on both leaves and internal nodes (as long as at least two time points are provided), we design a method to estimate the CIs for the estimates of wLogDate. We subsample the sampling times/calibration points given to us repeatedly to create N replicate data sets (where N is 100 by default, but can be adjusted). Note that our subsampling is not a bootstrapping procedure as node sampling times cannot be resampled with replacement. We then compute the time tree for each replicate to obtain N different estimates for the divergence time of each node, from which we can compute their CIs (95% as default). This sampling would work best when we have a fairly large number of calibration points, which is the case in phylodynamic settings where all (or nearly all) sampling times for the leaves are given, or in large phylogenies where abundant calibration points can be obtained from fossils. Although we refer to the resulting intervals as CIs, it is important to recognize that the resampling procedure is not strictly justified via bootstrap theory because subsampling is necessarily without replacement and sampled nodes are not independent of each other.

Phylodynamics Setting

To et al. (2016) simulated a data set of HIV env gene. Their time trees were generated based on a birth–death model with periodic sampling times. There are four tree models, namely D995_11_10 (M1), D995_3_25 (M2), D750_11_10 (M3), and D750_3_25 (M4), each of which has 100 replicates for a total of 400 different tree topologies. M1 and M2 simulate intra-host HIV evolution and are ladder-like whereas M3 and M4 simulate inter-host evolution and are balanced. Also, M4 has much higher root-to-tip distance (mean: 57) compared with M1–M3 (22, 33, and 29). Starting from conditions simulated by To et al. (2016), we use the provided time tree to simulate the clock deviations. Using an uncorrelated model of the rates, we draw each rate from one of three different distributions, each of which is centered at the value $$ \mu =0.006$$ as in To et al. (2016). Thus, we set each μ_i to $$ x_i\mu $$ where x_i is drawn from one of three distributions: LogNormal (mean : 1.0, std: 0.4), Gamma ($$ \alpha =\beta =6.05$$), and Exponential (λ = 1). Sequences of length 1,000 were simulated for each of the model conditions using SeqGen Rambaut and Grass (1997) under the same settings as To et al. (2016).

H1N1 2009 Pandemic

We reanalyze the H1N1 biological data provided by To et al. (2016) which includes 892 H1N1pdm09 sequences collected worldwide between March 13, 2009 and June 9, 2011. We reuse the estimated PhyML (Guindon et al. 2010) trees, 100 bootstrap replicates, and all the results of the dating methods other than wLogDate that are provided by To et al. (2016).

San Diego HIV

We study a data set of 926 HIV-1 subtype B pol sequences obtained in San Diego between 1996 and 2018 as part of the PIRC study. We use IQTree (Nguyen et al. 2015) to infer a tree under the GTR + Γ model, root the tree on 22 outgroups, then remove the outgroups. Because of the size, we could not run BEAST.

West African Ebola Epidemic

We study the data set of Zaire Ebola virus from Africa, which includes 1,610 near-full length genomes sampled between March 17, 2014 and October 24, 2015. The data were collected and analyzed by Dudas et al. (2017) using BEAST and reanalyzed by Volz and Frost (2017) using IQTree to estimate the ML tree and treedater to infer node ages. We run LSD, LF, and wLogDate on the IQTree from Volz and Frost (2017) and use the BEAST trees from Dudas et al. (2017), which include 1,000 sampled trees (BEAST-1000) and the Maximum clade credibility tree (BEAST-MCC). To root the IQTree, we search for the rooting position that minimizes the triplet distance (Sand et al. 2013) between the IQTree and the BEAST-MCC tree.

Methods Compared

For the phylodynamics data, we compared wLogDate with three other methods: LSD (To et al. 2016), LF (Langley and Fitch 1974), and BEAST (Drummond and Rambaut 2007). For all methods, we fixed the true rooted tree topology and only inferred branch lengths. For LSD, LF, and wLogDate, we used phyML (Guindon et al. 2010) to estimate the branch lengths in substitution units from sequence alignments and used each of them to infer the time tree. LSD was run in the same settings as the QPD* mode described in the original paper (To et al. 2016). LF was run using the implementation in r8s (Sanderson 2003). wLogDate was run with ten feasible starting points. For the Bayesian method BEAST, we also fixed the true rooted tree topology and only inferred node ages. Following To et al. (2016), we ran BEAST using HKY+$$ \Gamma 8$$ and coalescent with constant population size tree prior. We used two clock models on the rate parameter: the strict-clock (i.e., fixed rate) model and the LogNormal model. For the strict-clock prior, we set clock rate prior to a uniform distribution between 0 and 1. For the LogNormal prior, we set the ucld.mean prior to a uniform distribution between 0 and 1, and ucld.stdev prior to an exponential distribution with parameter $$ 1/3$$ (default). We always set the length of the MCMC chain to 10⁷ generations, burn-in to 10%, and sampling to every 10⁴ generations (identical to To et al. [2016]).

For the autocorrelated rate model, we compared wLogDate with LF and RelTime (Tamura et al. 2018), which is one of the state-of-the-art model-free dating methods. We randomly chose subsets of the internal nodes (10% on average) as calibration points and created 20 tests for each of the 10 replicates (for a total of 200 tests). We also compared wLogDate with DAMBE using this data set. Because DAMBE can only be run in interactive mode where each calibration point has to be manually placed onto the tree, we could not run DAMBE on the 200 tests with hundreds of calibration points in total. Therefore, we instead ran DAMBE only once on each of the ten trees and infer a unit time tree for each of them (i.e., calibrate the root to be at 1 unit time backward) and compared the results to that of wLogDate. DAMBE does not accept identical sequences so we removed identical sequences from the simulated alignments and trees before running DAMBE and ran wLogDate using these reduced trees to have a fair comparison.

Evaluation Criteria

On the simulated phylodynamics data set where the ground truth is known, we compare the accuracy of the methods using several metrics. We compute the RMSE of the true and estimated vector of the divergence times (τs) and normalize it by tree height. We also rank methods by RMSE rounded to two decimal digits (to avoid different ranks when errors are similar). In addition, we examine the inferred divergence tMRCA and mutation rate. The comparison of methods mostly focuses on point-estimates of these parameters and the accuracy of the estimates (as opposed to their variance). In one analysis, we also compare the CIs produced by wLogDate and BEAST on one model condition (M3 with LogNormal rate distribution). Finally, we examine the correlation between variance of the error in wLogDate and divergence times and branch lengths.

On the simulated data with autocorrelated rate, we show the distributions of the divergence times estimated by wLogDate, LF, and RelTime and report the RMSE normalized by tree height for each replicate. To compare with DAMBE in inferring unit time trees, we report the average relative error of the inferred to the true divergence times. After removing identical sequences, there are 438 internal nodes in total across the 10 tree replicates. For each internal nodes, we compute the relative error of its divergence time inferred by either DAMBE or wLogDate to its true divergence time in the normalized true time tree, which is $$ |\widehat{t_i}-t_i|/t_i$$ where $$ \widehat{t_i}$$ and t_i are the inferred and true divergence times of node i, respectively. We report the average relative error per tree replicate and the average of all 438 nodes for DAMBE and wLogDate.

On real data, we show LTT plots (Nee et al. 1994), which trace the number of lineages at any point in time and compare tMRCA times to the values reported in the literature. We also compare the runtime of wLogDate to all other methods in all analyses.