Email updates

Keep up to date with the latest news and content from Algorithms for Molecular Biology and BioMed Central.

Open Access Research

Stochastic errors vs. modeling errors in distance based phylogenetic reconstructions

Daniel Doerr1, Ilan Gronau2, Shlomo Moran3* and Irad Yavneh3

Author Affiliations

1 Center for Biotechnology, Bielefeld University, Bielefeld, Germany

2 Department of Biological Statistics and Computational Biology, Cornell University, Ithaca, USA

3 Computer Science Department, Technion - Israel Institute of Technology, Haifa, Israel

For all author emails, please log on.

Algorithms for Molecular Biology 2012, 7:22  doi:10.1186/1748-7188-7-22

The electronic version of this article is the complete one and can be found online at:

Received:23 December 2011
Accepted:28 June 2012
Published:31 August 2012

© 2012 Doerr et al.; licensee BioMed Central Ltd.

This is an Open Access article distributed under the terms of the Creative Commons Attribution License (, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.



Distance-based phylogenetic reconstruction methods use evolutionary distances between species in order to reconstruct the phylogenetic tree spanning them. There are many different methods for estimating distances from sequence data. These methods assume different substitution models and have different statistical properties. Since the true substitution model is typically unknown, it is important to consider the effect of model misspecification on the performance of a distance estimation method.


This paper continues the line of research which attempts to adjust to each given set of input sequences a distance function which maximizes the expected topological accuracy of the reconstructed tree. We focus here on the effect of systematic error caused by assuming an inadequate model, but consider also the stochastic error caused by using short sequences. We introduce a theoretical framework for analyzing both sources of error based on the notion of deviation from additivity, which quantifies the contribution of model misspecification to the estimation error. We demonstrate this framework by studying the behavior of the Jukes-Cantor distance function when applied to data generated according to Kimura’s two-parameter model with a transition-transversion bias. We provide both a theoretical derivation for this case, and a detailed simulation study on quartet trees.


We demonstrate both analytically and experimentally that by deliberately assuming an oversimplified evolutionary model, it is possible to increase the topological accuracy of reconstruction. Our theoretical framework provides new insights into the mechanisms that enables statistically inconsistent reconstruction methods to outperform consistent methods.

Phylogenetic reconstructions; Substitution models; Additive substitution rate functions


Phylogenetic reconstruction is the task of determining the topology of an evolutionary tree underlying a given set of samples (species) using sequence data extracted from them. This is typically done by assuming some simplified model for DNA sequence evolution, in most cases modeling it as a homogeneous continuous-time Markov process [1-3]. Distance-based reconstruction algorithms tackle this task by first computing a set of n 2 pairwise distances between the n input samples and then finding a tree which fits these distances. The distance measures used for this purpose typically reflect the rates of certain substitution events along the evolutionary paths in question. We thus refer to these distance measures as substitution rate (SR) functions. The distance-based approach is based on the fact that if the SR function used is additive for the underlying substitution model, and the input sequences are sufficiently long, then the topology of the true tree can be efficiently recovered with high probability. However, since the underlying evolutionary model is usually unknown, this assumption is rarely satisfied in practice.

Substitution models used for phylogenetic reconstruction range from the simplest Jukes-Cantor (JC) model [4], through slightly more complex and flexible models, such as Kimura’s two-parameter (K2P) model [5] and the Hasegawa-Kishino-Yano model (HKY) [6], to the General Time-Reversible (GTR) model [7,8]. In previous works [9,10] we observed that substitution models which are not too restrictive or too general have many inherently different additive SR functions. We used this basic observation to demonstrate that it is possible to adjust for each given set of DNA sequences a “good” additive SR function, which leads to significantly increased phylogenetic reconstruction accuracy, compared to other additive SR functions. This exploits our ability to predict the stochastic noise associated with each SR function. When the SR function used for distance estimation is additive for the underlying substitution model, this stochastic noise is the only cause for inaccurate reconstruction. However, in the scenario, which is very common in practice, where the SR function in use is not additive for the model, an additional systematic bias is introduced in the distance estimates. This systematic bias in distance estimation results in a phylogenetic reconstruction method that might be statistically inconsistent in some cases. In this papera, we extend our previous line of research to this scenario, by removing the constraint of additivity. We do this by considering both the stochastic noise and systematic error.

Several previous studies have demonstrated the utility of phylogenetic reconstruction methods that are not generally statistically consistent. The maximum parsimony method has been long known to be inconsistent in some cases [11,12]. However, in other cases it was shown to be more likely to produce accurate reconstructions, compared with the maximum likelihood method [13-15]. More recently, it has been demonstrated that reconstruction accuracy can be improved by deliberately assuming an oversimplified substitution model, when reconstructing a tree using maximum likelihood [16,17]. In the context of distance-based reconstruction, non-additive distance measures have been shown in several cases to lead to improved accuracy when compared with additive measures [18,19]. Overall, these studies provide convincing evidence for the need to consider inconsistent phylogenetic reconstruction methods. However, none of them provide a rigorous framework for characterizing the cases in which inconsistent methods outperform consistent ones.

In this paper we develop a theoretical framework which provides a practical and systematic way to quantify the effect of distance-estimation-bias on the accuracy of distance-based reconstruction. This framework is based on a novel method for measuring the deviation from additivity of SR functions. Coupled with the results in [9], this method enables evaluation of both the systematic bias and stochastic noise of SR functions. Such evaluation is important, because there is often a tradeoff between these two sources of error, stemming from the fact that simpler models with fewer parameters (such as JC) have smaller stochastic noise at the expense of greater estimation bias. Our framework allows us to consider this tradeoff when deciding which SR function to use for a given data set. This allows us to characterize a wide range of cases in which an SR function associated with an oversimplified evolutionary model results in increased reconstruction accuracy.

This finding falls in line with previous studies demonstrating the usefulness of phylogenetic reconstruction methods that are not generally consistent. Previous studies have attributed the increased accuracy of inconsistent methods mainly to the fact that these methods have a bias toward reconstructing certain topologies, leading to increased accuracy in cases where the phylogeny being reconstructed has the “favored topology”. We notice a similar behavior using our theoretical characterization of non-additive SR functions. However, somewhat surprisingly, we find that non-additive SR functions often have an advantage even when the phylogeny being reconstructed has an “unfavorable topology”. This is due to the reduced stochastic noise of the non-additive SR function (compared with it additive alternatives), which compensates for its topological bias.

Our paper is organized as follows. Section “Background” outlines some of the required background and introduces several new concepts that are central in our analysis. Section “Deviation from additivity in homogeneous substitution models” provides the main analytic results in the paper, and introduces deviation from additivity as a measure of distance estimation bias. In that section we prove a general upper bound for this deviation and establish a connection with reconstruction accuracy. We then study deviation from additivity and stochastic error of the JC distance formula when applied to data generated under the K2P model. In Section “Performance of Non affine-additive SR functions in quartet resolution” we study the effect of deviation from additivity and stochastic error on the accuracy of quartet reconstruction. In the case of quartets we can draw a tight connection between the different sources of error in distance estimation and inaccuracy of reconstruction. We present a useful heuristic, based on the so-called Fisher criterion ( [20,21]), for comparing the expected accuracy of two SR functions in this context. In Section “Simulations on Hasegawa’s Tree” we extend our study to larger trees using experiments on simulated data based on the tree obtained by Hasegawa in [6]. Finally, In Section “Inferring trees from genomic sequences” we demonstrate our approach through a series of experiments reconstructing trees from bacterial gene sequences.


In this section we provide a brief exposition of DNA substitution models and substitution rate functions used for distance estimation. We concentrate on details essential to this study and refer the reader to a previous paper [9] and standard textbooks [1,2] for a more complete survey.

Substitution Models

In this work, a DNA substitution model is simply a set of stochastic 4×4 transition matrices closed under matrix product (i.e., PQ PQ ). These matrices serve to describe the substitution process along evolutionary paths in a phylogenetic tree. All substitution models addressed in this paper are time-reversible [7]. A model tree in a time reversible substitution model , or an tree, is an undirected tree t =(VE ) in which each edge eE is associated with a transition matrix P e . An -tree t implies an inter-leaf transition matrix P ij for each pair of leaves { i , j } L ( T ) , namely P ij = e pat h T ( i , j ) P e . Most common models are defined using rate matrices, which are 4×4 matrices whose off-diagonal elements are non-negative substitution rates, and whose rows sum to 0. A stochastic transition matrix P is obtained from a rate matrix R through matrix exponentiation: P=eR.

A common assumption made on the substitution process is that it is homogeneous throughout time. This means that all rate matrices in the model are proportional to each other. Such a substitution model is thus termed homogeneous, and it is defined by a unit rate matrixR as follows: R = { e t R : t > 0 } . Note that the definition of the unit rate matrix associated with a given homogeneous model is somewhat arbitrary b, but once the unit R is defined, it implies a bijection (or equivalence) between rate matrices in R and the parameter t, which corresponds to evolutionary time. We will make use of this equivalence extensively throughout this paper.

We use the Kimura’s two-parameter (K2P) model [5] as a concrete example for demonstrating our approach. A rate matrix in this model is defined by two rate parameters: α, which is the rate of transition -type (ti) substitutions ( A G , C T ), and β, which is the rate of transversion -type (tv) substitutions ( { A,G } { C,T } ). Each K2P rate matrix can be represented as a product of a unit rate matrix, in which α + 2β =1, and a scalar t corresponding to evolutionary time .

K2P = e t R α , β | t > 0 , α β > 0 , α + 2 β = 1 ; R α , β = α β β α β β β β α β β α (1)

Each unit rate matrix of the K2P model defines a homogeneous sub-model, which is identified by its unique transition-transversion (ti-tv) ratio R = α 2 β 1 2 . The Jukes-Cantor (JC) model [4] is a special homogeneous sub-model of K2P, in which R = 1 2 (i.e., α =β ). Although the K2P model is defined in (1) as a union of its homogeneous sub-models, it is important to note that this union is closed under matrix product, implying that K2P adheres to our definition of a proper substitution model. Conversely, some commonly used substitution models, such as GTR and HKY, are defined as a union of homogeneous models, but are not themselves closed under matrix product [22].

Transition matrices in the K2P model have the same symmetric structure as the underlying rate matrices, with two distinct transition parameters: pα – the probability of a transition-type substitution; pβ– the probability of a transversion-type substitution. The transformations between (α,β,t ) and (pα,pβ) are given by the following equations:

αt = 1 2 ln ( 1 2 p β 2 p α ) + 1 4 ln ( 1 4 p β ) βt = 1 4 ln ( 1 4 p β ) . (2)

p α = 1 4 1 + e 4 βt 2 e 2 αt 2 βt p β = 1 4 1 e 4 βt . (3)

Substitution rate functions

A substitution rate (SR) function for a model is a non-negative continuous function Δ : R + that maps each transition matrix onto a numerical value of “substitution rate”. An SR function Δ induces the following dissimilarity mapping over the leaves of an -tree t : D Δ T ( i , j ) = Δ ( P ij ) , for all { i , j } L ( T ) . Of particular interest in phylogenetic reconstruction are additive SR functions.

Definition 2.1 (Additive SR function)

An SR function Δ is said to be additive for a substitution model if for all P , Q , Δ (PQ)=Δ (P) + Δ (Q).

It is often explicitly required that an SR function be additive for the assumed model (see [9]). The evolutionary time, t, typically serves as the standard additive measure in most common substitution models. Throughout this study we follow the special case of K2P, focusing on the two SR functions defined below.

Δ K2P ( p α , p β ) = 1 2 ln ( 1 2 p β 2 p α ) 1 4 ln ( 1 4 p β ) = αt + 2 βt = t . (4)

Δ JC ( p α , p β ) = 3 4 ln 1 4 3 ( p α + 2 p β ) = 3 4 ln 1 3 ( e 4 βt + 2 e 2 αt 2 βt ) . (5)

The first SR function, ΔK2P, is the common SR function suggested for the K2P model in [5], and it is clearly additive, as it maps the transition probabilities onto evolutionary time t . The second SR function, ΔJC, maps the transition probabilities onto evolutionary time only in the special case of the JC model where α =β . Under other homogeneous sub-models of K2P, it is non-additive. This non-additivity is analyzed in details in section Deviation from additivity in homogeneous substitution models.

Additive metrics, Affine-additive mappings, and Near-additivity

The core idea behind distance-based phylogenetic reconstruction is that a phylogenetic tree t can be accurately and efficiently reconstructed from pairwise distances which are additive with respect to t[23,24].

Definition 2.2 (Additive metric)

A metric d defined over the leaf-set L (t ) of a tree t is t -additive (or additive w.r.t t ), if there exists a positive edge-weighting function w : E ( T ) R + , such that for each i,jL (t ), D ( i , j ) = e pat h T ( i , j ) w ( e ) . d is additive for a set S if it is t -additive for some tree t where L (t )=S .

It is well known that additive SR functions imply additive metrics: if Δ is an additive SR function for a model , then for any -tree t, D Δ T (the dissimilarity mapping induced by Δ on t ) is a t -additive metric. The inherent difficulty in reconstructing phylogenies using additive SR functions is that computing the implied t -additive metric requires the exact values of the inter-taxon transition matrices {Pij}, and getting these exact values from alignments of finite length is practically impossible. Therefore, a distance-based reconstruction algorithm is useful in a realistic setting only if it has some robustness to error in distance estimation. In [25], Atteson observed that the topology of a phylogenetic tree t can be accurately (and efficiently) reconstructed from any dissimilarity mapping d which is sufficiently close to a t -additive metric, using certain “robust” distance-based algorithmsc. Formally, “sufficiently close” is defined by the following relation:

Definition 2.3 (Near-additive mapping)

A dissimilarity mapping d on L (t ) is said to be near-additive w.r.t. t iff there exists a t -additive mapping d s.t.

| | D , D | | = max { i , j } L ( T ) { | D ( i , j ) D ( i , j ) | } < 1 2 w min ( D ) , (6)

where w min ( D ) is the minimal weight assigned to an internal edged by the edge weighting function corresponding to the additive metric d.

For our results we will be using a generalization of this criterion, in which the mapping d can be any affine-additive mapping, defined below.

Definition 2.4 (Affine-additive mapping)

A dissimilarity mapping D is said to be affine-additive w.r.t. a phylogenetic tree t, if there is a t -additive metric d, and scalars a >0,B s.t. D = aD + b (i.e., D ( i , j ) = aD ( i , j ) + b for all { i , j } L ( T ) ).

As with additive metrics, affine-additive mapping are also associated with edge weights. Let d be a t -additive mapping corresponding to the edge-weighting function w (·). Then the edge weighting function w ( · ) corresponding to the affine additive mapping D = aD + b is given by: w ( e ) = aw ( e ) for all internal edges, and w ( e ) = aw ( e ) + 1 2 b for all external edges. When B is positive, D is actually an additive metric, but when B is negative, the weights of external edges implied by w ( · ) might be negative, and D might even yield negative dissimilarities. The generalization of Atteson’s theorem to cases where d is affine-additive follows from the observation that the robust distance-based reconstruct algorithms considered by Atteson are invariant to affine transformations of their input distances. From this point on, when we say a dissimilarity mapping d is near additive, we mean it satisfies (6) with respect to some affine-additive mapping d.

Local consistency

Atteson’s result plays a central role in arguing the statistical consistency of distance-based phylogenetic reconstruction. Typically, this is done by assuming that the inter-leaf distances are computed using an SR function Δ which is additive for the underlying substitution model , as follows:

1. If Δ is additive for , then for each -tree t the mapping D Δ T defined by D Δ T ( i , j ) = Δ ( P ij ) for all i,jL (t ), is a t -additive metric.

2. As the length of the input sequences grows, the estimated transition matrices { P ̂ ij } converge (w.h.p.) to the true matrices {Pij}.

3. When { P ̂ ij } are sufficiently close to {Pij}, the estimated dissimilarity map D ̂ defined by D ̂ ( i , j ) = Δ ( P ̂ ij ) is sufficiently close to D Δ T , and is thus near-additive.

4. The near-additivity of the estimated dissimilarity map D ̂ implies accurate topological reconstruction, assuming a robust distance-based algorithm is used.

This line of argument has been used in numerous works studying statistical consistency of distance-based algorithms (e.g., [25-27]), and in all these cases an additive SR function is assumed. Notice, however, that this line of argument remains valid when D Δ T is near additive w.r.t. t . For instance, consistent reconstruction of any -tree is guaranteed by using an affine-additive SR function Δ , which is an affine transformation of some additive SR function Δ : Δ = + b (with a >0). An SR function that is not affine-additive in a given substitution model does not guarantee consistency across all -trees, but it still can be consistent for specific -trees.

Definition 2.5 (Consistent SR function)

An SR function Δ of a substitution model is said to be consistent w.r.t. an -tree t if D Δ T is near-additive w.r.t t .

The main idea endorsed in this paper is that if an SR function only deviates slightly from some SR function which is affine-additive for , then it might be consistent with respect to many -trees of interest, and as such should be considered for use in distance based reconstructions.

Deviation from additivity in homogeneous substitution models

In order to assess whether a given SR function Δ is consistent w.r.t. a given model tree t, one has to find an affine-additive mapping d which minimizes the ratio | | D Δ T , D | | w min ( D ) (see Definition 2.3). This task seems hard in a general setting, but in the special case of homogeneous substitution models it is tractable. Consider a homogeneous substitution model R . The unit rate matrix Rimplies a 1-1 mapping between evolutionary time t and rate matrices in R . It is thus useful to view an SR function for R as a function Δ : R + R + which maps the evolutionary timet to a dissimilarity measure Δ (t ).

It can be shown that such Δ is affine-additive in the model if and only if Δ (t )=at + B for some a R + , b R . We define the deviation of an SR function Δ from a given affine-additive function at + B in an interval [t0,t1] as 1 a max { | Δ ( t ) at b | : t [ t 0 , t 1 ] } (the factor 1 a normalizes the deviation to units of evolutionary time). The deviation from additivity of Δ within [t0,t1] is defined as the minimum deviation of Δ from any affine-additive function in that interval.

Definition 2.6 (Deviation from additivity)

Let Δ : R + R + be an SR function in a homogeneous substitution model. The deviation from additivity of Δ in an interval [t0,t1] is defined by:

dev ( Δ , [ t 0 , t 1 ] ) = inf a R + , b R max t [ t 0 , t 1 ] | Δ ( t ) at b | a . (7)

Lemma 2.7 below presents the basic relation between deviation from additivity and consistency. In Section Performance of Non affine-additive SR functions in quartet resolution we demonstrate the tightness of this relation.

Lemma 2.7

Let be a homogeneous model, and let t be an -tree with edge lengths (measured in time units) denoted by {te}. Let tmin=min{te:et }, and assume that all inter-leaf distances in t fall within the interval [t0,t1]. Then any SR function Δ in for which dev ( Δ , [ t 0 , t 1 ] ) < 1 2 t min is consistent w.r.t. t .


We need to show that D Δ T is near-additive w.r.t. t . Since dev ( Δ , [ t 0 , t 1 ] ) < 1 2 t min , there are a R + , b R which satisfy

max t [ t 0 , t 1 ] | Δ ( t ) at b | a < 1 2 t min .

For all i,jL (t ), denote t ij = e pat h T ( i , j ) t e , and let d be the dissimilarity map associated with evolutionary time: d (i,j )=tij. Clearly, d is an additive metric, and the dissimilarity mapping D = aD + b is an affine-additive mapping. The internal-edge-weights associated with D are given by w ( e ) = at ( e ) (see discussion following Definition 2.4), implying that w min ( D ) = a t min . We thus have:

| | D , D Δ T | | max t [ t 0 , t 1 ] | Δ ( t ) at b | < 1 2 a t min = 1 2 w min ( D ) .

An upper bound on the deviation of an SR function Δ from additivity in a given interval [t0,t1] is implied from the error associated with its linear interpolation At + B within that interval ( A = Δ ( t 1 ) Δ ( t 0 ) t 1 t 0 and B = t 1 Δ ( t 0 ) t 0 Δ ( t 1 ) t 1 t 0 ). Figure 1a demonstrates this for ΔJC under a homogeneous sub-model of K2P, and Lemma 2.8 below presents a general upper bound on the deviation from additivity. For this purpose, we assume that the SR function Δ is a monotone increasing continuous function of t with continuous first and second derivatives.

thumbnailFigure 1. Deviation from additivity and stochastic error. (a) ΔJC is portrayed (green) in the homogeneous sub-model of K2P with R =10 in the interval t ∈[0.8,2]. Its linear interpolation in that interval, Δint=At + B, is plotted in blue, and the maximum difference between the two functions is designated by X. The deviation of ΔJC from additivity in this setting is X 2 A . (b) The affine-additive SR function minimizing its deviation from ΔJC is Δ = Δ int + 1 2 X . The two SR functions ΔJC and Δ are shown with their stochastic error margins, assuming sequence length of 500 bp.

Lemma 2.8

Let Δ : R + R + be an SR function in a homogeneous substitution model, and let [t0,t1] be an interval. Let Δint(t )=At + B be the linear interpolation of Δ in [t0,t1] defined above, and let F = max t [ t 0 , t 1 ] { | Δ ( t ) | } . Then

dev ( Δ , [ t 0 , t 1 ] ) ( t 1 t 0 ) 2 F 16 A . (8)


Let us start by introducing a couple of auxiliary notations:

ψ ( a , b , t ) = Δ ( t ) at b ψ ( a , b ) = max t [ t 0 , t 1 ] | ψ ( a , b , t ) | .

We are looking for a R + and b R which minimize 1 a ψ ( a , b ) . Let ψmin=mint ∈[t0,t1]{ψ (A,B,t )}, ψmax=maxt ∈[t0,t1]{ψ (A,B,t )}, and let b = B + 1 2 ( ψ max + ψ min ) . Then ψ ( A , b ) = 1 2 ( ψ max ψ min ) . A bound for dev (Δ,[t0,t1]) will thus follow by showing that ψ max ψ min ( t 1 t 0 ) 2 F 8 .

Since Δint(t )=At + B is a linear interpolation of Δ in t0t1, we have ψ (ABt0)=ψ (ABt1)=0. Let tmin be an arbitrary point in the interval t0t1s.t. ψ (ABtmin)=ψmin≤0 and let (t2t3) be the maximal open interval in t0t1 containing tmin in which ψ (ABt )<0 (this interval can be empty if ψmin=0). We define a similar interval (t4t5) in which ψ (ABt )>0 around some arbitrary tmaxs.t. ψ (ABtmax)=ψmax. Note that the intervals (t2t3) and (t4t5) are disjoint, and that Δint is the linear interpolation of Δ in both these intervals (since ψ (ABt2)=ψ (ABt3)=ψ (ABt4)=ψ (ABt5)=0). Therefore, the bound on the error of polynomial interpolation (see, e.g., [28], p. 187) implies that

ψ min ( t 3 t 2 ) 2 F 8 and ψ max ( t 5 t 4 ) 2 F 8 ,

Combining these, we get

dev ( Δ , [ t 0 , t 1 ] ) 1 A ψ ( A , b ) = 1 2 A ( ψ max ψ min ) ( t 5 t 4 ) 2 + ( t 3 t 2 ) 2 F 16 A ( t 1 t 0 ) 2 F 16 A . (9)


In Appendix 3 we prove that if Δ does not intersect its linear interpolation Δint=At + B within the interval (t0,t1), then the function At + Bmentioned in the proof above is, in fact, the affine-additive function which minimizes the deviation from additivity of Δ in [t0,t1]. This means that, in such cases, the first inequality in (9) holds in equality. The last inequality in (9) also holds in equality in such cases, because we are guaranteed to have either [t2,t3]=[t0,t1] (when Δ is bounded from above by its linear interpolation) or [t4,t5]=[t0,t1] (when Δ is bounded from below by its linear interpolation). Thus, in such a case, the bound of Lemma 2.8 is reduced to the bound on interpolation error (middle inequality in (9)). Cases where Δ does not intersect its linear interpolation are frequent among many SR functions of interest, as this condition holds when Δ is either convex or concave.

Deviation of ΔJCfrom Additivity in K2P

We now turn to study the deviation of ΔJC from additivity in homogeneous sub-models of K2P with ti-tv ratio R > 1 2 . First, we express ΔJCas a function of the ti-tv ratio R and the time t, using (5) and the relations α 2 β = R and α + 2β =1.

Δ JC ( R , t ) = 3 4 ln 1 3 ( e 4 βt + 2 e 2 αt 2 βt ) = 3 4 ln 1 3 ( e 2 t R + 1 + 2 e t 2 R + 1 R + 1 ) = 3 4 ln 1 3 e 2 t R + 1 1 + 2 e t 2 R 1 R + 1 = 3 2 ( R + 1 ) t 3 4 ln 1 3 1 + 2 e t 2 R 1 R + 1 . (10)

Note that the homogeneous K2P sub-model with R = 1 2 is the JC model; in this case the second term of (10) vanishes, leaving Δ JC ( 1 2 , t ) = t . For other homogeneous sub-models of K2P, where R > 1 2 , ΔJC is not affine-additive (i.e., not of the form at + B for a >0), and we can use the result in Lemma 2.8 to bound the deviation of ΔJCfrom additivity. Denoting ρ = 2 R 1 R + 1 , we get

Δ JC ( R , t ) ∂t = 3 2 ( R + 1 ) + 3 2 ρ e ρt 1 + 2 e ρt > 0 . (11)

2 Δ JC ( R , t ) t 2 = 3 2 ρ 2 e ρt ( 1 + 2 e ρt ) 2 < 0 . (12)

3 Δ JC ( R , t ) t 3 = 3 2 ρ 3 ( 1 2 e ρt ) e ρt ( 1 + 2 e ρt ) 3 . (13)

We get that for any given ti-tv ratio R > 1 2 , ΔJC(R,t ) is a concave monotone increasing function, and its second derivative attains a global minimum of 3 16 ρ 2 at t = ln ( 2 ) ρ . By the note following Lemma 2.8, the deviation of ΔJC from additivity in an interval [t0,t1] can be evaluated by computing the linear interpolation Δint=At + B of ΔJC in [t0,t1], and finding t ∈[t0,t1] which maximizes ΔJC(t )−Δint(t ) (see Figure 1a). A bound on this deviation from additivity can be obtained through Lemma 2.8 by plugging in the slope of the linear interpolation, A, and the maximum value, F, attained by the second derivative of ΔJC in [t0,t1]. Using Lemma 2.7 and an expression for dev (ΔJC(R,t ),[t0,t1]), it is possible to map out coherent collections of homogeneous K2P-trees for which ΔJC is guaranteed to be consistent. Each collection is defined by a range of ti-tv ratios [0.5,Rmax], a range of inter-leaf distances [t0,t1], and a lower bound on the weights of internal edges in the tree, given by tmin=2dev (ΔJC(Rmax,t ),[t0,t1]).

After determining a collection of trees for which a given non-affine-additive SR function, Δ, is consistent, one can compare the performance of Δ with additive alternatives. In our case, we compare Δ =ΔJC, which is not affine additive when R > 1 2 , to the standard additive SR function ΔK2P. The potential advantage of ΔJCover ΔK2P lies in its reduced stochastic noise . Informally, this occurs because JC relies on the accuracy of estimating a single parameter - the sum p =pα + 2pβ, while ΔK2P relies on the accuracy of estimating each of the two parameters pαand pβ separately. The stochastic noise of an SR function is measured by the standard deviation of the statistical estimator associated with it, denoted σ (ΔJC) and σ (ΔK2P), respectively. We use the result in [9] to get a first order approximation (based on the delta method [29]) of σ (ΔK2P) for sequences of length k and model parameters Rt :

σ ( Δ K2P ) ( e 4 t R + 1 1 ) + 4 ( e 2 t R + 1 1 ) + 2 ( e 4 Rt R + 1 ( e 4 t R + 1 + 1 ) 2 ) 16 k . (14)

By a similar application of the delta method to ΔJC, we obtain:

σ ( Δ JC ) p ( t , R ) ( 1 p ( t , R ) ) k ( 1 4 3 p ( t , R ) ) 2 , (15)

where k is the sequence length and p ( t , R ) = p α + 2 p β = 3 4 1 4 e 2 t R + 1 1 2 e ( 2 R + 1 ) t R + 1 (see (3)).

Figure 1 provides an illustrative comparison of ΔJCand ΔK2P under the homogeneous sub-model of K2P with ti-tv ratio R =10, and within the inter-leaf time interval of [0.8,2]. Figure 1a shows the deviation of ΔJCfrom additivity in that setting, using its linear interpolation Δint=At + B . Note that Lemma 2.8 and the subsequent note imply that dev ( Δ JC , [ 0 . 8 , 2 ] ) = X 2 A , where X =maxt∈[0.8,2]{ΔJC(t )−Δint(t )}. Figure 1b depicts ΔJC in the same setting with its stochastic error margins (ΔJC±σ (ΔJC)), alongside its closest affine-additive function Δ = Δ int + 1 2 X and its stochastic error margins (Δ±σ (Δ)). These stochastic error margins are determined by assuming a sequence length of 500 bp in the first-order approximations given in (14) and (15), where σ (Δ) is given by scaling σ (ΔK2P) by the slope A of the linear interpolation. Note how the margins of ΔJCare actually more tightly concentrated around its affine-additive approximation Δ than the margins of Δ. This implies that, despite its deviation from additivity in this setting, distances obtained using ΔJCare actually more likely to be near-additive than distances obtained using ΔK2P.

Performance of Non affine-additive SR functions in quartet resolution

The quartet tree is the smallest phylogenetic tree with non-trivial topology. Focusing on quartets enables a close study of the effects of deviation from additivity and stochastic noise on reconstruction accuracy. The topology of a quartet spanning four taxa {1,2,3,4} can be represented by the split notation (ij |kl ) (where {ijkl }={1,2,3,4}), indicating that the internal edge of the quartet separates ij from kl . All distance based quartet resolution algorithms essentially reduce to the four-point method (FPM) [26,30], which resolves this split using the six observed pairwise distances { d ̂ ij : { i , j } { 1 , 2 , 3 , 4 } } : it first partitions the six observed distances into three sums d ̂ 12 + d ̂ 34 , d ̂ 13 + d ̂ 24 , and d ̂ 14 + d ̂ 23 , and then determines the quartet split according to the minimal sum (the sum d ̂ ij + d ̂ kl corresponds to the split (ij |kl )). We will focus on the task of reconstructing homogeneous K2P quartets using FPM with distances { d ̂ ij } estimated using either ΔJCor ΔK2P. We note that most of our findings easily generalize to more sophisticated homogeneous substitution models, replacing ΔJC by any concave distance function and ΔK2P by some SR function corresponding to the evolutionary time t .

For concreteness, we assume henceforth that the quartet split is (12|34), meaning that the sum of the exact evolutionary times t12 + t34 is minimal. We start by analyzing the impact of the deviation from additivity of ΔJC on the consistency of quartet resolutions. First, observe that any monotone distance function is consistent for quartets in which t12 and t34 are the smallest interleaf distances - as is the case with symmetric quartets, in which all external edges are of the same length. Therefore, we study two prototypes of asymmetric quartets. The length of the internal edge in both types is ti, and each type has two long external edges of length tl, and two short external edges of length ts. In type A quartets (Figure 2a), the short edges are on one side of the split and the long edges are on the other side. In this case d12and d34 re the smallest and largest interleaf distances (resp.). Hence, the concavity of ΔJC increases the separation between the sum d12 + d34 and the other two competing sums, leading to an expected improvement in reconstruction accuracy. The other quartet configuration (type B; Figure 2b) has a short edge and a long edge on both sides of the split. In this case, the interval of interpolation is [d13,d24], and the distance d12=d34 is near the center of this interval. Thus the concavity of ΔJC decreases the separation between the sums d13 + d24 and d12 + d34 by approximately twice the deviation from additivity of ΔJC in that range.

thumbnailFigure 2. Performance of the Four Point Method using ΔJC on K2P quartets with ti-tv ratio R = 2. The concave non affine-additive SR function ΔJC is shown (dashed green line) in the interval [t0,t1], where t0 and t1 are the smallest and largest of the six pairwise distances (resp.). The dashed blue line shows the linear interpolation Δint=At + B of ΔJC in the interval [t0,t1]. Horizontal dotted lines correspond to half of the two competing sums computed by FPM under the two SR functions (see legend). (a) In quartets of type A, t0=t12 and t1=t34, and so Δint(1,2) + Δint(3,4)=ΔJC(1,2) + ΔJC(3,4). However, for i ∈{1,2} and j ∈{3,4}, Δint(i,j )<ΔJC(i,j ). Therefore, the deviation from additivity of ΔJC increases its FPM separation, denoted SEPJC, compared to the FPM separation SEPint of Δint. (b) In quartets of type B, t0=t13 and t1=t24, and so Δint(1,3) + Δint(2,4)=ΔJC(1,3) + ΔJC(2,4). However, Δint(1,2)=Δint(3,4)<ΔJC(1,2)=ΔJC(3,4), and so Δint(1,2) + Δint(3,4)<ΔJC(1,2) + ΔJC(3,4). Therefore, the deviation from additivity of ΔJC decreases its FPM separation, denoted SEPJC, compared to the FPM separation SEPint of Δint. Note that SEPint remains invariant in both types of quartets under fixed ti whereas SEPJC changes, depending on the type of quartet and the ts/tl ratio.

When the deviation from additivity exceeds half the length of the internal edge, the sum d13 + d24 becomes the minimal sum, and ΔJC becomes inconsistent. Note that this demonstrates the tightness of the condition stated in Lemma 2.7, and in this sense, type B quartets provide a worst case scenario for quartet resolution by a concave SR functione.

Next we turn to compare the accuracy of ΔJC with that of ΔK2P when used to reconstruct its “worst case scenario” quartets of type B. Interestingly, ΔJC ends up outperforming ΔK2P on many of these quartets, due to its reduced stochastic noise (as predicted in our discussion revolving around Figure 1b). For example, consider a series of homogeneous K2P quartets of type B with ti-tv ratio R =5, whose edge lengths were set as follows: ti=0.2, tl=1.0, and ts∈[0.2,1.0]. We assessed reconstruction accuracy for both SR functions (ΔJCand ΔK2P) across this series of quartets, by generating 100,000 simulations of the substitution process using 1,000 bp long sequences for each quartet (Figure 3a). Despite its deviation from additivity, ΔJC outperforms the additive SR function ΔK2P on many of these quartets (as long as tl/ts<3.6) . Note that as ts shrinks, the deviation of ΔJC from additivity increases, since the interval [t0,t1] expands. This experiment appears to indicate that the deviation of ΔJC from additivity has to be quite large for ΔK2P to outperform it.

thumbnailFigure 3. Performance of ΔJC and ΔK2P on a series of quartets of type B. A series of homogeneous K2P quartets is considered (left illustration), with ti-tv ratio of R =5, and edge lengths ti=0.2, tl=1, and ts∈[0.2,1]. (a) Reconstruction accuracy using FPM and either ΔJC (dashed green) or ΔK2P (solid red) plotted against tl/ts. Accuracy ratio is estimated using 100,000 independent replicates for each value of ts in the interval [0.2,1] (in steps of 0.01), with sequence length 1,000 bp. (b) Fisher’s Criterion (FC) for the sums corresponding to splits (12|34) and (13|24) under either ΔJC (dashed green) or ΔK2P (solid red) plotted against tl/ts.

Fisher’s criterion for separability

We now present a simple and general framework based on the so-called Fisher Criterion (FC) for predicting the relative accuracy of two SR functions in resolving quartets. FC measures the effective separation between normal random variables Xn (μ1σ1) and Yn (μ2σ2) using the following measuref ( [20,21]):

FC ( X , Y ) = | μ 1 μ 2 | σ 1 2 + σ 2 2 . (16)

We use FC to measure the separability of the distance sum corresponding to the true split (which should be the minimal sum for consistent SR functions) from the two remaining sums. For the expectation μ of each sum we use the true distances as computed by the SR function on the actual model parameters. For the variance σ 2, we use the sum of the approximate variances of the two distances involved in the sum. We expect that an SR function which provides a larger separation of the smallest sum from the two other sums will imply a better reconstruction probability.

We note that FC is not an exact indicator of the separability in our case, because the necessary criteria for this are not satisfied in our model. Namely, the two distance sums are not normally distributed, and they are correlated through the substitution process along the external edges of the quartet. Nevertheless, as Figure 3b suggests, FC turns out to provide a quite reliable comparison of the expected performance of ΔJC and ΔK2P for the quartet series considered in the aforementioned experiment. Figure 3b exhibits for each quartet the FC of ΔJC alongside that of ΔK2P, both associated with the comparison of the true split (12|34) and the “ΔJC favored split” (13|24). As shown, the trends observed in both FC plots closely resemble the trends observed in the reconstruction accuracy plot (Figure 3a), and the the equilibrium point of the FC values of ΔJC and ΔK2P is very close to the equilibrium point of the accuracy of reconstructions of these two functions (near tl/ts=3.6).

A useful feature of this framework is the natural way in which it teases apart the stochastic noise from the deviation from additivity. If we denote the numerator of FC by SEP (for “separation”) and its denominator by NOISE, then a comparison of FC estimates between two SR function Δ1,Δ2 can be represented as a ratio of ratios:

FC ( Δ 1 ) FC ( Δ 2 ) = SEP ( Δ 1 ) SEP ( Δ 2 ) / NOISE ( Δ 1 ) NOISE ( Δ 2 ) . (17)

Figure 4 illustrates how a comparison between the expected performance of ΔJC and that of ΔK2P can be carried out by tracing the SEP and NOISE ratios along four series of homogeneous K2P quartet: the bottom-left plot corresponds to the quartet series considered in Figure 3; the plot above it corresponds to the same series with ti-tv ratio R =2; the two plots on the right describe two quartet series in which the weight of the short edges is constant ts=0.2, and the weight of the long edges ranges in [0.2,1]. These four series demonstrate several typical trends in the behavior of the SEP and NOISE ratios. First, we observe that the NOISE ratio decreases (favoring ΔJC) as the diameter of the quartet (t24) increases (it is almost constant in the two series on the left, and monotone decreasing in the series on the right). This is because the diameter provides the major contribution to the stochastic noise (for both ΔJC and ΔK2P), and as it increases, the ratio between the stochastic noise of ΔK2P and ΔJC increases as well. We also observe a natural decrease in the NOISE ratio with an increase in the ti-tv ratio (the NOISE ratio for R =5 is consistently smaller than for R =2). Concerning the SEP ratio, we see it becomes smaller (favoring ΔK2P) as the quartet becomes more unbalanced (the SEP ratio decreases along the X axis in each of the four plots). This is because the deviation of ΔJC from additivity increases as the inter-leaf distance interval [t0,t1]=[t13,t24] expands. Deviation of ΔJC from additivity also increases with the ti-tv ratio, as the substitution model further departs from the assumptions of JC (the SEP ratio for R =5 is consistently smaller than for R =2).

thumbnailFigure 4. SEP and NOISE ratios.SEP (ΔJC)/SEP (ΔK2P) (dashed) and NOISE (ΔJC)/NOISE (ΔK2P) (dotted) plotted against tl/ts for four series of homogeneous K2P quartets of type B. Top two series have ti-tv ratio of R =2, and bottom two series have ti-tv ratio of R =5. Left two series have external edge lengths tl=1 and ts∈[0.2,1], and right two series have external edge lengths tl∈[0.2,1] and ts=0.2. The length of the internal edge is constant ti=0.2 in all four series.

The two series on the right side of Figure 4 demonstrate well the tradeoff between the effects of stochastic noise and deviation from additivity. In both series, the SEP and NOISE ratios decrease as the quartets become more unbalanced (due to the trends listed above). However, the rates of decrease of these two ratios are different due to the different ti-tv ratios, and this determines the expected relative performance of the two SR functions across the series. When R = 2, the SEP ratio decreases at a slower rate than the NOISE ratio, and ΔJC is expected to outperform ΔK2P across the entire series. When R = 5, the SEP ratio decreases at a faster rate than the NOISE ratio, and when the quartets are sufficiently unbalanced (tl/ts>4) ΔK2P is expected to outperform ΔJC.

Simulations on Hasegawa’s Tree

In this section we describe experiments done on simulated data sets generated along the seven-taxon tree assembled by Hasegawa, Kishino, and Yano in 1985 [1,6]. This tree, spanning seven eutherian mammals (Figure 5a), was reconstructed originally using mitochondrial DNA sequences. It has a caterpillar topology (meaning that every internal node is incident to an external edge), and it has long external edges and short internal edges, making it a suitable representative of small phylogenetic trees spanning moderately distant species. These features also make it particularly challenging for distance-based reconstruction.

thumbnailFigure 5. Simulations on Hasegawa’s Tree. (a) Reconstruction accuracy of four different SR functions on different scaled versions of Hasegawa’s tree [6]. The tree with scaled edge weights is depicted (left) next to the graph (right) plotting reconstruction accuracy of four SR functions. Different scales of the tree are considered, indicated by the diameter of the tree (X axis). Reconstruction accuracy (Y axis) is measured for each scaled tree by the average normalized RF distance between the reconstructed tree and the true tree across 10,000 simulated data sets. Simulations were carried out assuming a ti=tv ratio of R = 2 and sequence length of 500 bp. (b) A similar plot is shown for a semi-symmetric caterpillar tree.

In our study we used the tree structure and edge lengths to generate simulated data sets. We considered the tree in various scales, by setting the tree diameter (largest inter-taxon path length) to values in the interval [0.1,2.0]. For each scale considered, 10,000 simulations were carried out, where in each simulation 500 bp sequences were evolved along the tree according to a homogeneous K2P substitution model with ti-tv ratio of R =2. For each simulated data set, estimated values of the K2P statistics pα and pβ, denoted by p ̂ α and p ̂ β , were extracted for all 7 2 pairs of taxa. Subsequently, several distance matrices were computed for each data set by applying different SR functions to these estimated statistics. Reconstruction accuracy was evaluated by applying the Neighbor Joining (NJ) algorithm [31,32] to these distance matrices and recording the Robinson-Foulds topological distance (RF) [33] between the reconstructed tree and the Hasegawa tree. Sequence simulation was performed using SeqGen [34] (by choosing the HKY model with uniform base frequencies), and tree reconstruction was performed using the version of NJ implemented in the PHYLIP package [35].

We studied the reconstruction accuracy associated with four different SR functions: ΔJC, ΔK2P, Δtv, and ΔR=2. The first two are as described in Equations (5) and (4), respectively. The third SR function, Δtv, considers only tv-type substitutions: Δ tv ( p α , p β ) = 1 4 log 1 4 p β ( t ) = βt , and the fourth SR function, ΔR=2, is based on a maximum likelihood (ML) estimatorg of the time t from the estimated transition probabilities p ̂ α , p ̂ β , given that R =2. Informally, this function, which uses knowledge of the true value of R (which is typically unknown to the user), is optimal in our setting, because it has similar stochastic noise as ΔJC, and it is additive since it coincides with ΔK2P when applied to transition probabilities p ̂ α , p ̂ β that are consistent with a ti-tv ratio of R =2.

The performance of these four SR functions is traced across the different tree scales in Figure 5a. For each SR function Δ and scale s, we recorded the average normalized RF distance from the true tree to each of the 10,000 trees reconstructed using Δ . The RF distance was normalized by its maximum value which is twice the number of internal edges in the tree (in our case 2×4=8). As observed previously in [9], ΔK2P performed well in shorter scales, and Δtv performed well in longer scales. However, both additive SR functions were significantly outperformed in nearly all cases by ΔJC. Surprisingly, ΔJC even slightly outperformed ΔR=2. We speculate that this happened due to a bias similar to the one observed in type A quartets in Section Performance of Non affine-additive SR functions in quartet resolution, improving the performance of concave SR functions such as ΔJC on certain K2P-trees.

To test this hypothesis, we went through a similar experiment with a more symmetric seven-taxon caterpillar tree, with internal edges of uniform length tint, and external edges of uniform length text=5tint (Figure 5b). The symmetry of this tree was expected to reduce the effect of the reconstruction bias observed in Hasegawa’s tree, and indeed, ΔJC performed much more poorly on this tree. Despite this fact, ΔJC still outperformed ΔK2P in all scales and Δtv in the smaller scales (s <1.1).

Inferring trees from genomic sequences

In this section we describe our study comparing various SR functions on genomic DNA sequences. Next to ΔJC and ΔK2P we also considered the well known LogDet SR function [36,37], denoted here as ΔLogDet. Extending our study to this setting is challenging in two respects. First of all, unlike the simulated case, the true tree is not known with complete confidence, and accuracy of reconstruction can only be determined by using a well-accepted reference tree that may contain some errors. Secondly, the true substitution model is also unknown and is likely to violate the assumptions of both JC and K2P models and even the relaxed assumptions of the general time-reversible model (in which ΔLogDet is additive). Hence, we have to assume in this case that ΔJC, ΔK2P, and ΔLogDet are all non affine-additive, where ΔJC and ΔK2P are still likely to exhibit higher deviation from additivity than ΔLogDet, since they make stronger assumptions on the substitution model.

The genomic data set

In building the genomic data set, we made use of a set of 31 clusters of orthologous groups (COGs) which was compiled by Ciccarelli et al. and used for inferring phylogenetic relationships amongst a large number of species in [38,39]. These 31 gene families were selected to capture the evolutionary history of the species containing them. This was done in [38] by making sure that the genes in these families have the following properties: (1) they are highly conserved across species, (2) they have a small number of paralogs, and (3) they are weakly affected by horizontal gene transfer. We scanned the NCBI genome database and found 199 bacterial genomes that contained all annotated COGs. For each of the 31 COGs, we extracted the appropriate protein sequence in each of the 199 bacterial species, choosing an arbitrary paralog in cases of multiple hits. We followed a procedure similar to the one described in [38,39] to obtain reliable multiple-sequence alignments for each COG: we computed a 199-way multiple alignment of the protein sequences of each COG using HMMalign [40] and then mapped each protein sequence back to its coding DNA sequence. The conserved parts of each of the 31 DNA alignments were extracted using GBLOCKS [41] to filter out alignment columns with 50% or more gap symbols. The alignments were manually scanned, and 36 species which contributed a large number of gaps to the alignments were removed from the subsequent analysis. The 31 different alignments were concatenated to form one long 163-way multiple sequence DNA alignment.

For the reference tree we used the phylogenetic tree of microbial species provided by the ARB-SILVA Living Tree Project [42]. This tree, spanning 8,029 species at the time of writing, is based on a widely accepted analysis of the small subunit (SSU) 16S RNA. A subtree spanning our 163 bacterial species was extracted from this tree and treated as the true phylogenetic tree in our analysis.

Reconstruction accuracy for ten-species subsets

We used the base set of 163 species to generate 40,000 random 10-species sub-alignments. The random selection process was guided to generate species subsets corresponding to a wide range of diameter scales (a blind random selection process is biased toward subsets with large diameters). For each of the 40,000 subsets, a 10-way subalignment was extracted from the original 163-way alignment, and in this alignment we extracted only columns corresponding to four-fold degenerate sites that do not have any gap symbol. This is done to make sure the sites used for distance estimation have undergone a substitution process that is as uniform as possible along the different lineages and across the different sites. Each sub-alignment was used to compute three distance matrices – one under ΔJC, one under ΔK2P, and one under ΔLogDet. The latter was calculated by the version that is implemented in the PHYLIP package. The NJ algorithm was then applied to the three matrices and the resulting trees were compared to the true tree (as depicted by the appropriate LTP subtree) according to the RF distance.

As an additional comparison, we used a fourth reconstruction technique. This method (termed BIONJ-GTR) used the BIONJ reconstruction algorithm [43] on distances obtained under the general time-reversible model with invariant sites and Gamma distribution of rates across variant sites (GTR+Γ +I) [8,44].

The PhyML package [45] was used to infer this tree for each of the 40,000 subsets. We selected the GTR+Γ +I model since it was found by the MEGA5 software [46] to provide the best fit to the sequence data. The 40,000 sampled instances were partitioned into eight bins according to the RF distance observed between the BIONJ-GTR tree and the true (LTP) tree, and average RF distances were recorded for each of the three SR functions in each bin. This allowed us to observe trends throughout these 40,000 samples (Figure 6). Of the 40,000 trees inferred under ΔJC, 83.1% showed an equal or lower RF distance than those reconstructed by the BIONJ-GTR method. Moreover, ΔJC outperformed ΔK2P and ΔLogDet on average in all partitions, and ΔLogDet showed by far the worst performance with 48.7% of all reconstructed trees achieving higher RF distances to the reference tree than those inferred by BIONJ-GTR. As with our results on simulated data sets, we see that the SR functions with lower stochastic error but inferior model fit performed best. Unsurprisingly, the GTR+G+I model itself, which was predicted to have the best fit to the sequence data, was often outperformed by the simpler JC and K2P models. Note that the difference in performance between ΔJC and the two other SR functions is greater for subsets that are more accurately reconstructed by the BIONJ-GTR approach (the lower bins). This appears to indicate that over-simplified distance methods are particularly beneficial when the sequence data conveys a stronger phylogenetic signal.

thumbnailFigure 6. Evaluation against BIONJ-GTR tree. The 40,000 subsets of size 10 were partitioned according to the the RF-distance between the reference LTP tree and the tree reconstructed using BIONJ-GTR (X axis). The (left) Y axis describes the mean difference between the RF-distance associated with a tree reconstructed using a particular SR function (ΔK2P, ΔJC, or ΔLogDet) and the RF-distance associated with the BIONJ-GTR tree. The bar plot in the background depicts the number of subsets in each bin.


In this paper we explored the basic properties of methods for estimating evolutionary distances, and studied how these properties affect the accuracy of distance-based phylogenetic reconstruction. We considered both the systematic bias and the stochastic noise (variance) of the distance estimators, and examined the tradeoff between these two factors. We focused on the common task of phylogenetic reconstruction under homogeneous substitution models. Assuming homogeneous models simplifies the analytical framework, since in such models each SR function is reduced to a univariate function of the evolutionary time t . However, obtaining accurate estimates of t is still a hard task in this setting, since the unit rate matrix is unknown. An SR function Δ is guaranteed to yield consistent reconstruction across all trees in a homogeneous model only if it is additive, meaning that it is a linear function of t . When Δ is not additive, it introduces a systematic bias in distance estimates, which we denoted here as deviation from additivity . Some SR functions are only additive in one homogeneous model, whereas others are additive across a wider collection of homogeneous models. This less constrained additivity is typically achieved at a price of increased estimation noise. We studied the tradeoff between “deviation from additivity” and “estimation noise” via a case study where the model tree is a homogeneous K2P-tree with an unknown ti-tv ratio R . In this case, Kimura’s distance formula ΔK2P is always additive, while the less noisy Jukes Cantor’s formula, ΔJC, is additive only when R = 1 2 .

A study of this type requires a way to measure the deviation from additivity of a non-additive SR function Δ in a given range of distances [t0,t1]. To this end, we introduced the concept of affine-additive distance functions, and defined the deviation from additivity of Δ in [t0,t1] as the distance of Δ from its closest affine-additive function in [t0,t1]. We established a tight connection between this measure and statistical consistency of reconstruction (Lemma 2.7) and derived an upper bound for deviation from additivity in homogeneous models (Lemma 2.8). We applied these results in analyzing the deviation from additivity of ΔJC, and its effect on the accuracy of reconstructing homogeneous K2P-trees. We then showed, both analytically (in Section Deviation from additivity in homogeneous substitution models) and through experiments on simulated data sets (in Sections Performance of Non affine-additive SR functions in quartet resolution and Simulations on Hasegawa’s Tree), that, compared to ΔK2P, it is often better to use the non-additive but less noisy estimates of ΔJC, even when R is quite high. Somewhat surprisingly, we found this to be the case even when the tree being reconstructed has an “unfavorable” topology. Our experiments on bacterial gene sequences (Section Inferring trees from genomic sequences) also indicate that the simple and less noisy SR functions perform better on average than ones that are expected to better fit the true substitution process.

The framework presented in this paper implies a practical way for selecting SR functions which are likely to increase the accuracy of distance estimation. The practicality of the method is drawn from the fact that the criteria by which we select an SR function depend only a relatively crude information about the tree being reconstructed. For instance, in the case of a homogeneous K2P-tree, one can easily obtain from the input sequences rough estimates of both the ti-tv ratio R and the range of inter-leaf times [t0,t1]. These estimates can then be used to compare the expected accuracies of ΔJC and ΔK2P on the given input, and determine which of them is more likely to yield an accurate phylogeny. For quartets, a tight comparison can be made using the FC-based approach suggested in Section Fisher’s Criterion for Separability, and for larger trees, a cruder comparison can be made using a plot like the one presented in Figure 1b. A promising avenue of further research is to extend the FC-based approach to allow tighter prediction of reconstruction accuracy of trees spanning more than four taxa.


aThis is a WABI 2011 special issue invited paper. Extended abstract of this paper appeared in [47]. bTypically, the unit rate matrix is assumed to be the one corresponding to one substitution per site. cMany common distance-based algorithms, such as the Neighbor Joining (NJ) algorithm [31,32], are known to be robust in this sense. dIn a tree, edges which touch leaves are external, and all other edges are internal. eTypes A and B quartets represent the Farris zone and Felsenstein zone, resp. (see, e.g., [1], Chapter 9). fWe use here the square root of the criterion commonly used in the literature, because we prefer to think in terms of distances rather than squares of distances. This has no practical influence, since we use FC only for comparing between different choices, not for assessing the quality of a give choice.g This ML estimate is obtained by a simple numerical method for maximizing the likelihood function (see, e.g., [1]).


Tightness of Lemma 2.8

Let f (t ) be a (continuous) function on some interval [t0,t1]. We prove below that if f does not intersect its linear interpolation At + B in that interval, then dev ( f , [ t 0 , t 1 ] ) = 1 A max t [ t 0 , t 1 ] | f ( t ) At b | . We use the following notations, conforming to the notations in the proof of Lemma 2.8:

ψ ( a , b , t ) = f ( t ) at b ψ ( a , b ) = max t [ t 0 , t 1 ] | ψ ( a , b , t ) | ψ ( a ) = min b R ψ ( a , b ) .

Lemma 2.9

Let f (t ) be a monotone increasing function in the interval [t0,t1] and let At + B be its linear interpolation in [t0,t1]. If either f (t )≥At + B for all t ∈[t0,t1] or f (t )≤At + B for all t ∈[t0,t1], then for all a >0, we have 1 a ψ ( a ) 1 A ψ ( A ) .


We prove the minimality of 1 A ψ ( A ) in the case where f (t )≥At + B for all t ∈[t0,t1]. The other case (where f (t )≤At + B for all t ∈[t0,t1]) can be proven in an identical fashion.

For a >0, let Babe the maximum value of Bs.t. ψ ( a , b , t ) 0 for all t ∈[t0,t1]. Evidently, ψ ( a ) = 1 2 ψ ( a , b a ) . If the linear interpolation of f (t ) in [t0,t1] is given by At + B, then BA=B . We need to show that for every a >0, it holds that (a,Ba)> (A,BA). Let tAbe a point in [t0,t1] s.t. ψ (A,BA,tA)=ψ (A,BA). Note that if a <A, then the two linear functions At + BAand at + Baintersect at (t0,f (t0)), and if a >A, then they intersect at (t1,f (t1)) (see Figure 7).

thumbnailFigure 7. Proof of Lemma 2.9. A function f(t) is portrayed (bold) with its linear interpolation At + B =At + BA(green) in the interval [t0,t1], s.t.f (t)≥At + B for allt ∈[t0,t1]. Equation (18) is illustrated for a <A on the right, and equation (19) is illustrated for a >A on the left.

For a <A, we get the following equality (Figure 7; right):

ψ ( A , b A , t A ) + A ( t A t 0 ) = f ( t A ) f ( t 0 ) = ψ ( a , b a , t A ) + a ( t A t 0 ) . (18)

Hence, since ψ (a,Ba)≥ψ (a,Ba,t ) for every t ∈[t0,t1], and since a <A, we get

( A , b A , t A ) + aA ( t A t 0 ) < ( a , b a , t A ) + Aa ( t A t 0 ) ( A , b A ) < ( a , b a ) .

Similarly, if a >A, we get the following equality (Figure 7; left)

A ( t 1 t A ) ψ ( A , b A , t A ) = f ( t 1 ) f ( t A ) = a ( t 1 t A ) ψ ( a , b a , t A ) , (19)

and a >A implies that

aA ( t 1 t A ) ( A , b A ) > Aa ( t 1 t A ) ( a , b a ) ( A , b A ) < ( a , b a ) .

Competing interests

The authors declare that they have no competing interests.

Authors’ contributions

All authors participated in discussing, formulating, and modulating the research. DD performed the simulations and experiments of Sections Simulations on Hasegawa’s Tree and Section Inferring trees from genomic sequences. IG and SM initiated and directed the research and drafted the manuscript. IY performed the analysis in Sections Deviation from additivity in homogeneous substitution models and Section Performance of Non affine-additive SR functions in quartet resolution and contributed to the ideas of the project. All authors contributed to the writing and editing of the manuscript, and all authors read and approved the final manuscript.


This research was supported by the Israel Science Foundation (ISF) grant No. 509/11. We also acknowledge the support for the publication fee by the Deutsche Forschungsgemeinschaft and the Open Access Publication Funds of Bielefeld University. The third author would also like to thanks the Max Planck Institute for Informatics for supporting a visit under which part of this research was carried out.


  1. Felsenstein J: Inferring Phylogenies. Sunderland: MA Sinauer Associated Inc; 2004. OpenURL

  2. Semple C, Steel M: Phylogenetics. Oxford University Press; 2003. OpenURL

  3. Papoulis A, Pillali SU: Probability, Random Variables and Stochastic Processes. New York: McGraw Hill Higher Education; 2002. OpenURL

  4. Jukes T, Cantor C: Evolution of Protein Molecules. In Mammalian Protein Metab. Edited by Munro H. New York: Academic Press; 1969:21-132. OpenURL

  5. Kimura M: A simple method for estimating evolutionary rates of base substitutions through comparative studies of nucleotide sequences.

    J Mol Evol 1980, 16(2):111-120. PubMed Abstract | Publisher Full Text OpenURL

  6. Hasegawa M, Kishino H, Yano T: Dating of the human-ape splitting by a molecular clock of mitochondrial DNA.

    J Mol Evol 1985, 22(2):160-174. PubMed Abstract | Publisher Full Text OpenURL

  7. Tavaré S: Some Probabilistic and Statistical Problems in the Analysis of DNA Sequences.

    Lectures on Mathematics in the Life Sci 1986, 17:57-86. OpenURL

  8. Lanave C, Preparata G, Saccone C, Serio G: A new method for calculating evolutionary substitution rates.

    J Mol Evol 1984, 20:86-93. PubMed Abstract | Publisher Full Text OpenURL

  9. Gronau I, Moran S, Yavneh I: Towards Optimal Distance Functions for Stochastic Substitution Models.

    J Theor Biol 2009, 260(2):294-307. PubMed Abstract | Publisher Full Text OpenURL

  10. Gronau I, Moran S, Yavneh I: Adaptive Distance Measures for Resolving K2P Quartets: Metric Separation versus Stochastic Noise.

    J Comp Biol 2010, 17(11):1391-1400. OpenURL

  11. Felsenstein J: Cases in which parsimony or compatability methods will be positively misleading.

    Syst Zool 1978, 27:401-410. Publisher Full Text OpenURL

  12. Cavender J: Taxonomy with confidence.

    Math Biosci 1978, 40:271-280. Publisher Full Text OpenURL

  13. Steel M, Penny D: Parsimony, likelihood, and the role of models in molecular phylogenetics.

    Mol Biol Evol 2000, 17:839-850. PubMed Abstract | Publisher Full Text OpenURL

  14. Sober E: A likelihood justification of parsimony.

    Cladistics 1985, 1:209-233. Publisher Full Text OpenURL

  15. Felstenstein J, Sober E: Parsimony and likelihood: an exchange.

    Syst Zool 1986, 35:617-626. Publisher Full Text OpenURL

  16. Yang Z: How often do wrong models produce better phylogenies?

    Mol Biol Evol 1997, 14:105-108. PubMed Abstract | Publisher Full Text OpenURL

  17. Bruno WJ, Halpern AL: Topological bias and inconsistency of maximum likelihood using wrong models. [] webcite

    Mol Biol Evol 1999, 16(4):564-566. PubMed Abstract | Publisher Full Text OpenURL

  18. Zharkikh A: Estimation of evolutionary distances between nucleotide sequences.

    J Mol Evol 1994, 39(3):315-329. PubMed Abstract | Publisher Full Text OpenURL

  19. Gascuel O, Guindon S: Efficient Biased Estimation of Evolutionary Distances When Substitution Rates Vary Across Sites.

    Mol Biol Evol 2002, 19(4):534-543. PubMed Abstract | Publisher Full Text OpenURL

  20. Fisher R: The use of multiple measurements in taxonomic problems.

    Ann of Eugenics 1936, 7:177-188. OpenURL

  21. Duda R, Hart P: Pattern Classification and Scene Analysis. Hoboken: John Wiley and Sons; 1973. OpenURL

  22. Sumner J, Fernandez-Sanchez J, Jarvis P: Lie Markov Models.

    J Theor Biol 2012, 298:16-31. PubMed Abstract | Publisher Full Text OpenURL

  23. Buneman P: The recovery of trees from measures of dissimilarity. In Mathematics in the Archeological and Historical Sciences. Edited by Hodson F, Kendall D, Tautu P. Edinburgh University Press; 1971:387-395. OpenURL

  24. Sattath S, Tversky A: Additive similarity trees.

    Psychometrica 1977, 42(3):319-345. Publisher Full Text OpenURL

  25. Atteson K: The Performance of Neighbor-Joining Methods of Phylogenetic Reconstruction.

    Algorithmica 1999, 25:251-278. Publisher Full Text OpenURL

  26. Erdos P, Steel M, Szekely L, Warnow T: A few logs suffice to build (almost) all trees (I).

    Random Struct Algorithms 1999, 14:153-184. Publisher Full Text OpenURL

  27. Erdos P, Steel M, Szekely L, Warnow T: A few logs suffice to build (almost) all trees (II).

    Theoret Comput Sci 1999, 221:77-118. Publisher Full Text OpenURL

  28. Johnson L, Riess R: Numerical Analysis. Boston: Addison Wesley; 1977. OpenURL

  29. Oehlert G: A note on the delta method.

    Am Statistician 1992, 46:27-29. OpenURL

  30. Zaretskii K: Constructing a tree on the basis of a set of distances between the hanging vertices.

    Uspekhi Mat Nauk 1965, 20(6):90-92.

    [In Russian]


  31. Saitou N, Nei M: The neighbor-joining method: a new method for reconstructing phylogenetic trees.

    Mol Biol Evol 1987, 4:406-425. PubMed Abstract | Publisher Full Text OpenURL

  32. Studier J, Keppler K: A note on the neighbor-joining algorithm of Saitou and Nei.

    Mol Biol Evol 1988, 5(6):729-731. PubMed Abstract | Publisher Full Text OpenURL

  33. Robinson F, Foulds R: Comparison of phylogenetic trees.

    Math Biosci 1981, 53:131-147. Publisher Full Text OpenURL

  34. Rambaut A, Grass NC: Seq-Gen: an application for the Monte Carlo simulation of DNA sequence evolution along phylogenetic trees.

    Comput Appl Biosci 1997, 13(3):235-238. PubMed Abstract OpenURL

  35. Felsenstein J: PHYLIP - Phylogeny Inference Package (Version 3.2).

    Cladistics 1989, 5:164-166. OpenURL

  36. Steel M: Recovering a tree from the leaf colourations it generates under a Markov model.

    Appl Math Lett 1994, 7(2):19-24. Publisher Full Text OpenURL

  37. Lockhart P, Steel M, Hendy M, Penny D: Recovering evolutionary trees under a more realistic model of sequence evolution.

    Mol Biol Evol 1994, 11(4):605-612. PubMed Abstract | Publisher Full Text OpenURL

  38. Ciccarelli FD, Doerks T, von Mering C, Creevey CJ, Snel B, Bork P: Toward Automatic Reconstruction of a Highly Resolved Tree of Life.

    Science 2006, 311(5765):1283-1287. PubMed Abstract | Publisher Full Text OpenURL

  39. von Mering C, Hugenholtz P, Raes J, Tringe SG, Doerks T, Jensen LJ, Ward N, Bork P: Quantitative Phylogenetic Assessment of Microbial Communities in Diverse Environments.

    Science 2007, 315(5815):1126-1130. PubMed Abstract | Publisher Full Text OpenURL

  40. Durbin R, Eddy SR, Krogh A, Mitchison G: Biological Sequence Analysis: Probabilistic Models of Proteins and Nucleic Acids. Cambridge University Press; 1999. OpenURL

  41. Talavera G, Castresana J: Improvement of phylogenies after removing divergent and ambiguously aligned blocks from protein sequence alignments.

    Syst Bio 2007, 56:564-577. Publisher Full Text OpenURL

  42. Yarza P, Ludwig W, Euzeby J, Amann R, Schleifer KH, Glockner FO, Rossello-Mora R: Update of the All-Species Living Tree Project based on 16S and 23S rRNA sequence analyses.

    Syst Appl Microbiol 2010, 33:291-299. PubMed Abstract | Publisher Full Text OpenURL

  43. Gascuel O: BIONJ: an improved version of the NJ algorithm based on a simple model of sequence data.

    Mol Biol Evol 1997, 14(7):685-695. PubMed Abstract | Publisher Full Text OpenURL

  44. Rodriguez F, Oliver JL, Marin A, Medina JR: The general stochastic model of nucleotide substitution.

    J Theor Biol 1990, 142:485-501. PubMed Abstract | Publisher Full Text OpenURL

  45. Guindon S, Gascuel O: A simple, fast and accurate algorithm to estimate large phylogenies by maximum likelihood.

    Syst Biol 2003, 52:696-704. PubMed Abstract | Publisher Full Text OpenURL

  46. Tamura K, Peterson D, Peterson N, Stecher G, Nei M, Kumar S: MEGA5: Molecular Evolutionary Genetics Analysis using Maximum Likelihood, Evolutionary Distance, and Maximum Parsimony Methods.

    Mol Biol Evol 2011, 28:2731-2739. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  47. Doerr D, Gronau I, Moran S, Yavneh I: Stochastic Errors vs. Modeling Errors in Distance Based Phylogenetic Reconstructions. In Algorithms in Bioinformatics, Volume 6833 of Lecture Notes in Computer Science. Edited by Przytycka T, Sagot MF. Berlin / Heidelberg: Springer; 2011:49-60. OpenURL