---
abstract: |
    The goal of branch length estimation in phylogenetic inference is to estimate the divergence time between a set of sequences based on compositional differences between them. A number of software is currently available facilitating branch lengths estimation for homogeneous and stationary evolutionary models. Homogeneity of the evolutionary process imposes fixed rates of evolution throughout the tree. In complex data problems this assumption is likely to put the results of the analyses in question.

    In this work we propose an algorithm for parameter and branch lengths inference in the discrete-time Markov processes on trees. This broad class of nonhomogeneous models comprises the general Markov model and all its submodels, including both stationary and nonstationary models. Here, we adapted the well-known Expectation-Maximization algorithm and present a detailed performance study of this approach for a selection of nonhomogeneous evolutionary models. We conducted an extensive performance assessment on multiple sequence alignments simulated under a variety of settings. We demonstrated high accuracy of the tool in parameter estimation and branch lengths recovery, proving the method to be a valuable tool for phylogenetic inference in real life problems. ${\texttt{Empar}}$ is an open-source C++ implementation of the methods introduced in this paper and is the first tool designed to handle nonhomogeneous data.

    [: nucleotide substitution models; branch lengths; maximum-likelihood; expectation-maximization algorithm.]{}
author:
- 'A. M. Kedzierska$^{1,2}$'
- 'M. Casanellas$^{2,**}$'
title: '${\texttt{Empar}}$: EM-based algorithm for parameter estimation of Markov models on trees.'
---

 {#section .unnumbered}

Assuming that an evolutionary process can be represented in a phylogenetic tree, the tips of the tree are assigned operational taxonomic units (OTUs) whose composition is known. Here, the OTUs are thought of as the DNA sequences of either a single or distinct taxa. Internal vertices represent ancestral sequences and inferring the branch lengths of the tree provides information about the speciation time.

Choice of the evolutionary model and the method of inference have a direct impact on the accuracy and consistency of the results [@SulSwo97; @Fel78; @BruHal99; @Pen94; @HueHil93; @Schwartz2010]. We assume that the sites of a multiple sequence alignment (MSA) are independent and identically distributed (i.i.d. hypothesis of all sites undergoing the same process without an effect on each other), the evolution of a set of OTUs along a phylogenetic tree ${\tau}$ can be modeled by the evolution of a single character under a hidden Markov process on ${\tau}$.

Markovian evolutionary processes assign a conditional substitution (transition) matrix to every edge of ${\tau}$. Most current software packages are based on the continuous-time Markov processes where the transition matrix associated to an edge $e$ is given in the form $\exp(Q^e
t_e)$, where $Q^ee$ is an instantaneous mutation rate matrix. Although in some cases the rate matrices are allowed to vary between different lineages (cf. [@Galtier1998],[@YY99]), it is not uncommon to equate them to a *homogeneous* rate matrix $Q$, which is constant for every lineage in ${\tau}$.

Relaxing the homogeneity assumption is an important step towards increased reliability of inference (see [@mitoch]). In this work, we consider a class of processes more general than the homogeneous ones: the discrete-time Markov processes. If ${\tau}$ is rooted, these models are given by a root distribution $\pi,$ and a set of transition matrices $A^e$ (e.g. chap. 8 of @Semple2003). The transition matrices $A^e$ can freely vary for distinct edges and are not assumed to be of exponential form, thus are highly applicable in the analyses of non-homogeneous data. Among these models we find the general Markov model (${\mathtt{GMM}}$) and all its submodels, e.g. discrete-time versions of the Jukes-Cantor model (denoted as ${\mathtt{JC69}^{\ast}}$), Kimura two-parameters (${\mathtt{K80}^{\ast}}$) and Kimura 3-parameters models (${\mathtt{K81}^{\ast}}$), and the strand symmetric model ${\mathtt{SSM}}$. Though the discrete-time models provide a more realistic fit to the data [@YY99; @Ripplinger2008; @Ripplinger2010], their complexity requires a solid inferential framework for accurate parameter estimation. In continuous-time models, *maximum-likelihood estimation* (MLE) was found to outperform Bayesian methods [@Schwartz2010]. The most popular programs of phylogenetic inference (PAML [@Yang1997], PHYLIP [@Felsenstein1989], PAUP\* [@PAUP]) are restricted to the homogeneous models. Though more realistic, the use of nonhomogeneous models in phylogenetic inference is not yet an established practice. Recently, [@Jayaswal2011] proposed two new non-homogeneous models. With the objective of testing stationarity, homogeneity and inferring the proportion of invariable sites, the authors propose an iterative procedure based on the *Expectation Maximization* (${\mathtt{EM}}$) algorithm to estimate parameters of the non-homogeneous models (cf. [@barryhartigan87]). The ${\mathtt{EM}}$ algorithm was formally introduced by [@Dempster1977] (cf. @Hartley1958). It is a popular tool to handle incomplete data problems or problems that can be posed as such (e.g. missing data problems, models with latent variables, mixture or cluster learning). This iterative procedure globally optimizes all the parameters conditional on the estimates of the hidden data and computes the maximum likelihood estimate in the scenarios, where, unlike in the fully-observed model, the analytic solution to the likelihood equations are rendered intractable. An exhaustive list of references and applications can be found in [@Tanner1996], and more recently in [@Ambroise1998]. Here, we extend on the work of [@Jayaswal2011] and present ${\texttt{Empar}}$, a MLE method based on the ${\mathtt{EM}}$ algorithm which allows for estimating the parameters of the (discrete-time) Markov evolutionary models. ${\texttt{Empar}}$ is an implementation suitable for phylogenetic trees on any number of leaves and currently includes the following evolutionary models: ${\mathtt{JC69}^{\ast}}$, ${\mathtt{K80}^{\ast}}$, ${\mathtt{K81}^{\ast}}$, ${\mathtt{SSM}}$ and ${\mathtt{GMM}}.$

We test the proposed method on simulated data and analyze the accuracy of the parameter and branch length recovery. The tests are conducted in a settings analogue to that of [@Schwartz2010] and evaluate the performance of ${\texttt{Empar}}$ on the four and six-taxon trees with several sets of branch lengths, ${\mathtt{JC69}^{\ast}}$ and ${\mathtt{K81}^{\ast}}$ models under varying alignment lengths. We present an in-depth theoretical study, investigating the dependence of the performance on factors such as model complexity, size of the tree, positioning of the branches, data and total tree lengths.

Our findings suggest that the method is a reliable tool for parameter inference of small sets of taxa, best results obtained for shorter branches.

The algorithm underlying ${\texttt{Empar}}$ was implemented in C++ and is freely available to download at <http://genome.crg.es/cgi-bin/phylo_mod_sel/AlgEmpar.pl>.

METHODS {#methods .unnumbered}
=======

Models {#models .unnumbered}
------

We fix a set of $n$ taxa labeling the leaves of a rooted tree ${\tau}$. We denote by $N({\tau})$ the set of all nodes of ${\tau}$, the set of leaves as $L({\tau})$, the set of interior nodes as $Int({\tau}),$ and the set of edges as $E({\tau}).$ We are given a DNA multiple sequence alignment (MSA) associated to the taxa in ${\tau}$ and a discrete-time Markov process on ${\tau}$ associated to an evolutionary by a model ${\mathcal{M}}$, where the nodes in${\tau}$ are discrete random variables with values in the set of nucleotides $\{{\mathtt{A}},{\mathtt{C}},{\mathtt{G}},{\mathtt{T}}\}$. We assume that all sites in the alignment are i.i.d. and model evolution per site as follows: for each edge $e$ of ${\tau}$ we collect the conditional probabilities $P(y|x,e)$ (nucleotide $x$ being replaced by $y$ at the descendant node of $e$) in a transition matrix $A^e=(P(y|x,e))_{x,y}$; $\pi=(\pi_{{\mathtt{A}}},\pi_{{\mathtt{C}}},\pi_{{\mathtt{G}}},\pi_{{\mathtt{T}}})$ is the distribution of nucleotides at the root $r$ of ${\tau}$ and $\xi=\{\pi, (A^e)_{e}\}$ the set of continuous parameters of ${\mathcal{M}}$ on ${\tau}$. We denote by $X$ the set of $4^n$ possible patterns at the leaves and $Y$ the set of $4^{|Int({\tau})|}$ possible patterns at the interior nodes of ${\tau}.$ In what follows, the joint probability of observing $\textbf{x}=(x_l)_{l \in L({\tau})} \in X$ at the leaves and nucleotides $\textbf{y}=(y_v)_{v\in
Int({\tau})} \in Y$ at the interior nodes in ${\tau}$ is calculated as $$p_{\textbf{x},\textbf{y}}(\xi)=\pi_{y_r}\prod_{v \in N({\tau})\setminus\{r\}
}A^{e_{an(v),v}}_{y_{an(v)},y_{v}}$$ where $an(v)$ denotes a parent node of node $v,$ $e_{an(v),v}$ is an edge from $an(v)$ to $v$ (note: if $v$ is a leaf, then $y_v=x_v$).

In the *complete model* the states at the interior nodes are observed and the joint distribution is computed as above. On the other hand, the *observed model* assumes the variables at the interior nodes to be latent. In the latter case, the probability of observing $\textbf{x}=(x_l)_{l \in
L({\tau})}$ at the leaves of ${\tau}$ can be expressed as $$p_{\textbf{x}}(\xi)=\sum_{\textbf{y}=(y_v)_{v\in
Int({\tau})} \in Y}p_{\textbf{x},\textbf{y}}(\xi).$$ Restricting the shape of the transition matrices $A^e$ leads to different evolutionary models such as ${\mathtt{JC69}^{\ast}}$, ${\mathtt{K80}^{\ast}}$, ${\mathtt{K81}^{\ast}}$, ${\mathtt{SSM}}$ (see @CFK, @arBook, and @AllmanRhodes_chapter4 for references and background on the discrete-time models). The first three are the discrete-time versions of the widely used continuous-time JC69 [@Jukes1969], K80 [@Kimura1980] and K81 [@Kimura1981] models. The Strand Symmetric model ${\mathtt{SSM}}$ (@CS) is a discrete-time generalization to the `HKY` model [@Hasegawa1985] with equal distribution of the pairs of bases (`A,T`) and (`C,G`) at each node of the tree. It reflects the double-stranded nature of DNA and was found to be well-suited for long stretches of data [@Yap2004]. Lastly, the general Markov model ${\mathtt{GMM}}$ (@Allman2003 [@Steel1994a]) is free of restrictions on the entries of $A^e$, non-stationary, and can be thought as a non-homogeneous version of the general time reversible model ([@GTR]).

Expectation-Maximization algorithm {#expectation-maximization-algorithm .unnumbered}
----------------------------------

An algebraic approach to the Expectation-Maximization (${\mathtt{EM}}$) algorithm was first introduced in [@Pachter2005]. In this work, we adapted this approach to the context of phylogenetic trees.

Let $D$ denote a MSA recorded into a vector of $4^{|{{L({\tau})}}|}$ counts of patterns $u_D=(u_\textbf{x})_{\textbf{x}\in X},$ where each $u_\textbf{x}$ stands for the counts of a particular configuration of nucleotides **x** at the leaves, observed as columns in the alignment. We are interested in maximizing the likelihood function: $$\mathcal{L}_{obs}(\xi ; u_D)=\prod_{\textbf{x}\in X}p_{\textbf{x}}(\xi)^{u_{\textbf{x}}}$$ (up to a constant). Let $U_{cD}=(u_{\textbf{x},\textbf{y}})_{\textbf{x}\in X,\textbf{y}\in
  Y}$ be an array of counts for the complete model, where $u_{\textbf{x},\textbf{y}}$ is the number of times $\textbf{x}$ was observed at the leaves and $\textbf{y}$ at the interior nodes. The likelihood for the complete model has a multinomial form $$\label{eq:likelEM}
\mathcal{L}_c(\xi ; U_{cD})=\prod_{\textbf{x}\in X,\textbf{y}\in Y}p_{\textbf{x},\textbf{y}}(\xi)^{u_{\textbf{x},\textbf{y}}}=\prod_{\textbf{x}\in X,\textbf{y}\in Y}(\pi_{y_r}\prod_{v \in N({\tau})\setminus\{r\}
}A^{e_{an(v),v}}_{y_{an(v)},y_{v}})^{u_{\textbf{x},\textbf{y}}}$$ (up to a constant), which is guaranteed to have a global maximum given by a model-specific explicit formula (see ${\textit{Supp. mat.}}$ A).

${\mathtt{EM}}$ algorithm iteratively alternates between the expectation (${\emph{E-step}}$) and maximization step (${\emph{M-step}}$). In the ${\emph{E-step}}$ the algorithm uses the tree topology, current estimates of parameters and the observed data $u_D$ to assign a posterior probability to each of the possible $4^{|{{L({\tau})}}|}$ patterns in $X$ and the expected counts of the complete model, $u_{cD}$. This step can be efficiently performed using the peeling algorithm of [@Felsenstein2004]. In the ${\emph{M-step}}$ the updated MLE of the parameters are obtained by maximizing the likelihood of the complete model . The procedure is depicted in Fig. 1.

![Expectation-maximization algorithm.[]{data-label="fig:em_procedure"}](plots/code-crop){width="90.00000%"}

The likelihood is guaranteed to increase at each iteration of this process (e.g. [@Wu], [@Husmeier2005]). Moreover, for a compact set of parameters the algorithm converges to a critical point of the likelihood function. Although the output of the algorithm is not guaranteed to be a global maximum, multiple starting points are used for optimal solution.

Statistical tests {#statistical-tests .unnumbered}
-----------------

The substitution matrices are assumed stochastic. The number $d$ of free parameters for transition matrices in ${\mathtt{JC69}^{\ast}}$, ${\mathtt{K80}^{\ast}}$, ${\mathtt{K81}^{\ast}}$, ${\mathtt{SSM}}$ and ${\mathtt{GMM}}$ models is 1, 2, 3, 6, and 12 respectively. The root distribution under the ${\mathtt{SSM}}$ (respectively ${\mathtt{GMM}}$) model has 2 (resp. 3) free parameters, and is uniform for the remaining models considered here. For clarity of exposition, hereon the reference to the root distribution will be omitted; however, the formulas can be easily modified to include the root.

We let $\xi^e$ be the vector of free parameters defined as the off-diagonal elements of a transition matrix $A$ associated to a edge $e$ ($\xi^e_1=A_{1,2}$, $\xi^e_2$ is the next –from left to right, top down– off-diagonal entry distinct from $\xi^e_1$, etc). The procedure is repeated until $\xi^e_d$ is reached. Let $\xi=(\xi_i^e)_{i=1,\ldots,d; e\in{E({\tau})}}$ denote the vector of free parameters for an evolutionary model ${\mathcal{M}}$ as above and let $\hat{\xi}$ be its MLE. Under certain regularity conditions [@Zacks1971 Chap. 5], $\hat{\xi}$ exists, is consistent, efficient and asymptotically normal with mean $\xi$ and the covariance matrix given by the inverse of the Fisher information matrix [@Rao1973; @Efron1978]. The entries of the $d|{E({\tau})}|\times
d|{E({\tau})}|$ Fisher information matrix ${\textbf{I}}$ over free parameters are given by: $${\textbf{I}}(\xi^e_k, \xi^{e'}_{k'}) = -\mathbf{E}\left(\frac{\partial^2
      \log{\mathcal{L}}_{obs}(\xi; u_D)}{\partial
      \xi^{e}_{k}\partial \xi^{e'}_{k'}}\right)$$ (see ${\textit{Supp. mat.}}$ B for details). The Wald statistics for testing the null hypothesis $\xi_i^e=\hat{\xi}_i^e$, $e\in{E({\tau})}, i=1,\ldots,d$, is $$\label{eq:null_distr}
(\hat{\xi^e}-\xi^e)^T {\textbf{I}}^e(\hat{\xi^e}-\xi^e)\sim\chi^2_{d},$$ where ${\textbf{I}}^e$ denotes the $d\times d$ slice of ${\textbf{I}}$ corresponding to the free parameters of $e\in{E({\tau})}$. The $p-value$ can thus be easily calculated by looking at the tails of the corresponding $\chi^2$ distribution.

We tested the validity of the test statistics in our data by simulating a variety of MSAs under the complete model and compared it to the theoretical distribution . Figure I in the ${\textit{Supp. mat.}}$ C shows high fit and proves that the setting is appropriate.

Variances of the free parameters of the model and the full (observed) covariance matrix are saved in the output of ${\texttt{Empar}}$. These in turn can be used as the plug-in estimators in to calculate the $p-values$ and normal confidence intervals for the parameters.

We denote by $V^{e}_{i,i}$ the $i^{th}$ diagonal entry of the matrix $({\textbf{I}}^e)^{-1}$ corresponding to the variance of the free parameter $\xi^e_i$, $i=1,\ldots,d$. For the models with $d>1$ (i.e. all but ${\mathtt{JC69}^{\ast}}$), the variances of the free parameters can be summarized in a combined form $cV^e$ for each edge $e$: $$\label{eq:comb_var}
cV^{e}(\xi^e)=\frac{\sum^d_{j=1} \left(V^e_{j,j}+\left(\xi^e_j-\frac{\sum^d_{j=1}\xi^e_j}{d}\right)^2\right)}{d}.$$

Branch lengths {#branch-lengths .unnumbered}
--------------

The evolutionary distance between two nodes in ${\tau}$ joined by an edge $e$ with substitution matrix $A^e$ is defined as the total number of substitutions per site along $e$. This quantity is referred to as the *branch length* of edge $e$ (or of matrix $A^e$) and, following on [@barryhartigan87trans], can be approximated by: $$l({{A}^{e}})=-\frac{1}{4}\log\det({{A}^{e}}).\label{eq:brlength}$$ We denote the *total length of the tree* ${\tau}$ by ${\mathsf{L}_{{\tau}}}$, ${\mathsf{L}_{{\tau}}}=\sum_{e\in|{E({\tau})}|}l({{A}^{e}}).$ Now, let $A$ and $A'$ be two invertible $4\times 4$ matrices such that the entries of $(A')^{-1}(A-A')$ are small. From , we get $$\begin{aligned}
\label{eq:br_bound}
|l(A)-l(A')|&=&\frac{1}{4}|\log\frac{\det(A)}{\det(A')}|=\frac{1}{4}|\log\det((A')^{-1}A)|
\nonumber\\
&=&\frac{1}{4}|\log\det(\textbf{Id} +
(A')^{-1}(A-A'))|\nonumber\\&\approx& \frac{1}{4}|\log(1+Tr((A')^{-1}(A-A')))|\nonumber\\
&\approx&\frac{1}{4}|Tr((A')^{-1}(A-A'))|\leq \frac{1}{4}4||(A')^{-1}(A'-A))||_{1} \nonumber\\
&\leq& ||(A')^{-1} ||_{1} ||A-A'||_{1},\end{aligned}$$ where $||.||_1$ is the maximum absolute column sum of the matrix. Therefore if $A'$ is a good approximation to $A$, then $l(A')$ is a good approximation to $l(A)$. In what follows, we use the statistical test above to show the accurate recovery of the parameters. By the above argument, we can conclude that the estimates of the branch lengths will also be accurate (see also Results section).

Simulated data {#simulated-data .unnumbered}
==============

Performance assessment of ${\texttt{Empar}}$ was conducted on the MSAs simulated on four and six-taxon trees following [@Schwartz2010]. In the case of four taxon trees we fixed an inner node as the root and considered three types of topologies: ${{\tau}_{\mathrm{balanced}}}^4$ corresponds to the “balanced” trees with all five branches equal; the inner branch in ${{\tau}_{1:2}}$ is half the length of the exterior branches; and ${{\tau}_{2:1}}$ denotes a topology with the inner branch double the length of the external ones (see Fig. 2). In ${{\tau}_{\mathrm{balanced}}}^4$ and ${{\tau}_{2:1}}$ we let the length $l_0$ of the inner branch vary from 0.01 to 1.4, where starting from 0.05 it increases in steps of 0.05; in ${{\tau}_{1:2}}$ we let $l_0$ vary in $(0,0.7).$ For 6-taxon trees we used only balanced trees ${{\tau}_{\mathrm{balanced}}}^6$ (see Fig. 2) with $l\in (0,0.7)$.

![Unrooted trees used for simulations: ${{\tau}_{\mathrm{balanced}}}^4$, ${{\tau}_{1:2}}$, ${{\tau}_{2:1}}$ and ${{\tau}_{\mathrm{balanced}}}^6$ (*from left to right*).[]{data-label="fig_trees"}](plots/tree-eq "fig:"){width="20.00000%"} ![Unrooted trees used for simulations: ${{\tau}_{\mathrm{balanced}}}^4$, ${{\tau}_{1:2}}$, ${{\tau}_{2:1}}$ and ${{\tau}_{\mathrm{balanced}}}^6$ (*from left to right*).[]{data-label="fig_trees"}](plots/tree-short "fig:"){width="20.00000%"} ![Unrooted trees used for simulations: ${{\tau}_{\mathrm{balanced}}}^4$, ${{\tau}_{1:2}}$, ${{\tau}_{2:1}}$ and ${{\tau}_{\mathrm{balanced}}}^6$ (*from left to right*).[]{data-label="fig_trees"}](plots/tree-long "fig:"){width="20.00000%"} ![Unrooted trees used for simulations: ${{\tau}_{\mathrm{balanced}}}^4$, ${{\tau}_{1:2}}$, ${{\tau}_{2:1}}$ and ${{\tau}_{\mathrm{balanced}}}^6$ (*from left to right*).[]{data-label="fig_trees"}](plots/tree-6t "fig:"){width="20.00000%"}

We simulated multiple sequence alignments on trees with 4 and 6 leaves under ${\mathtt{JC69}^{\ast}}$ and ${\mathtt{K81}^{\ast}}$ models. We used the $\texttt{GenNon-H}$ package available from <http://genome.crg.es/cgi-bin/phylo_mod_sel/AlgGenNonH.pl>. In brief, based on an input phylogenetic tree with given branch lengths, $\texttt{GenNon-H}$ samples the substitution matrices corresponding to these lengths for all edges and uses them to generate the DNA MSAs following discrete-time Markov process on the tree. The output of this software is the alignment, the substitution matrices, root distribution (whenever non-stationary) and the variances of the continuous free parameters. We note that for the ${\mathtt{JC69}^{\ast}}$ model, there is a $1-1$ correspondence between the branch length and the free parameters of the substitution matrix. This does not hold for other models, were different substiution matrices may give the same branch length. We set the alignment length $L$ to $300nt, 500nt, 1,000nt$ and $10,000nt$ for 4-taxa and to $1,000nt$ or $10,000nt$ for 6-taxa. For the ${\mathtt{JC69}^{\ast}}$ and ${\mathtt{K81}^{\ast}}$ evolutionary models, a phylogenetic tree ${\tau}$ (with branch lengths), and a given alignment length, we run each analysis $1,000$ times. and estimated the parameters using ${\texttt{Empar}}$.

All MSAs used for the tests are accessible at the ${\texttt{Empar}}$ webpage.

Identifiability {#identifiability .unnumbered}
---------------

It is known that in certain cases the substitution parameters are not identifiable (e.g. parameters at the edges adjacent to the root of valency 2). As shown in [@Allman2003], the ${\mathtt{GMM}}$ model and its submodels, are identifiable up to a permutation of rows. [@BHident] showed that incorrect order of rows in the matrices can lead to a negative determinant of the substitution matrix from which the branch lengths cannot be calculated.

We expected this problem to arise in short data sets and large branch length, as those correspond to the substitution matrices with smaller diagonal value. For all the data sets used for tests, we calculated the percentage of cases among the 1,000 simulations for which the parameters estimated by the EM algorithm were permuted. This phenomenon was only observed in the data sets of $300nt$ and $1,000nt$. In the first case, the estimated matrices were permuted when the initial branch length was 0.55 or longer and corresponded to 0.005-0.023$\%$ of the cases; in the latter, for the branches of 0.6 or longer with at most 0.001$\%$ permuted matrices. Shorter branch lengths and longer alignments did not suffer from the above problem and recovered the underlying order in all of the cases.

As shown by [@Chang96], the entries of the Diagonal Largest in Column (DLC) substitution matrices are identifiable. Namely, there exists a unique set of substitution matrices satisfying the DLC condition and a unique root distribution that leads to a given joint distribution at the leaves. In order to ensure the reliability of the results we designed a procedure that scans the tree in the search of the permutations that maximize the number of substitution matrices with larger diagonal entries. It is not possible to maximize it for all edges, thus the goal is to find the permutations giving more weights to the lower parts of the tree, starting with the nodes corresponding to the outer branches. Given a tree ${\tau}$, we choose an interior node to be the root, directing all edges outwards. For each interior node $x$, we apply a permutation $S(x)$ of $\{{\mathtt{A}},{\mathtt{C}},{\mathtt{G}},{\mathtt{T}}\}$ that maximizes the sum of diagonal entries of the matrices assigned to the outgoing edges of $x$. Permutations $S(x)$ are applied recursively to the subtrees of ${\tau}$, moving $x$ from the outer nodes towards the root.

RESULTS AND DISCUSSION {#results-and-discussion .unnumbered}
======================

We present the results on the simulated data sets and discuss their dependence on the length of the alignments, the length of the branches and the depth of the branches in the tree–1 for the external branches and 2 for the internal branches [@Schwartz2010]. In cases with multiple branches of equal depth, we chose one of them at random.

Results of statistical tests {#results-of-statistical-tests .unnumbered}
----------------------------

Each sample gave rise to a $p-value$ based on the $\chi^2_d$ test given by . The $p-values$ are a measure of strength of evidence against the null hypothesis: for both exceptionally small or large $p-values$ one can reject the null hypothesis. We recorded the proportion of samples for which the $p-value$ lied in the interval $(0.05,0.95)$. The results are shown Table \[tab:pval\_jc\_short\] for the ${\mathtt{JC69}^{\ast}}$ model on the ${{\tau}_{1:2}}$ tree (also see Tab. I-V in the ${\textit{Supp. mat.}}$ C). We observe that even for short alignments of 300nt the null hypothesis cannot be rejected in approximately 95% of the tests.

[p[0.8cm]{}p[0.8cm]{}p[0.8cm]{}p[0.8cm]{}p[0.8cm]{}|p[0.8cm]{}p[0.8cm]{}p[0.8cm]{}p[0.8cm]{}]{} & &\
$\mathbf{l}$ $\mid$ L & & & & & & & &\
0.01 & 0.971 & 0.972 & 0.968 & 0.946 & 0.972 & 0.949 & 0.868 & 0.958\
0.05 & 0.947 & 0.951 & 0.947 & 0.948 & 0.974 & 0.943 & 0.953 & 0.952\
0.10 & 0.949 & 0.953 & 0.964 & 0.952 & 0.952 & 0.948 & 0.948 & 0.955\
0.15 & 0.952 & 0.954 & 0.958 & 0.938 & 0.946 & 0.953 & 0.940 & 0.947\
0.20 & 0.957 & 0.944 & 0.944 & 0.954 & 0.949 & 0.965 & 0.944 & 0.954\
0.25 & 0.957 & 0.955 & 0.955 & 0.956 & 0.945 & 0.939 & 0.955 & 0.936\
0.30 & 0.957 & 0.943 & 0.945 & 0.955 & 0.943 & 0.946 & 0.941 & 0.948\
0.35 & 0.952 & 0.943 & 0.958 & 0.958 & 0.948 & 0.943 & 0.950 & 0.960\
0.40 & 0.955 & 0.946 & 0.947 & 0.957 & 0.951 & 0.951 & 0.936 & 0.944\
0.45 & 0.949 & 0.944 & 0.944 & 0.947 & 0.948 & 0.955 & 0.958 & 0.958\
0.50 & 0.948 & 0.935 & 0.942 & 0.941 & 0.929 & 0.949 & 0.954 & 0.946\
0.55 & 0.954 & 0.949 & 0.946 & 0.957 & 0.936 & 0.944 & 0.944 & 0.952\
0.60 & 0.940 & 0.942 & 0.937 & 0.953 & 0.944 & 0.934 & 0.948 & 0.955\
0.65 & 0.940 & 0.934 & 0.955 & 0.952 & 0.938 & 0.938 & 0.945 & 0.948\
0.70 & 0.944 & 0.936 & 0.942 & 0.946 & 0.917 & 0.940 & 0.944 & 0.948\
0.75 & 0.922 & 0.932 & 0.947 & 0.934 & 0.922 & 0.932 & 0.943 & 0.950\
0.80 & 0.909 & 0.932 & 0.926 & 0.957 & 0.957 & 0.928 & 0.943 & 0.941\
0.85 & 0.912 & 0.912 & 0.932 & 0.948 & 0.968 & 0.930 & 0.936 & 0.947\
0.90 & 0.870 & 0.885 & 0.919 & 0.951 & 0.980 & 0.918 & 0.929 & 0.953\
0.95 & 0.852 & 0.888 & 0.939 & 0.951 & 0.981 & 0.965 & 0.908 & 0.944\
1,00 & 0.824 & 0.866 & 0.893 & 0.935 & 0.982 & 0.981 & 0.896 & 0.933\
1,05 & 0.816 & 0.853 & 0.889 & 0.930 & 0.980 & 0.981 & 0.898 & 0.937\
1.10 & 0.806 & 0.852 & 0.891 & 0.921 & 0.990 & 0.995 & 0.925 & 0.945\
1.15 & 0.784 & 0.812 & 0.867 & 0.938 & 0.980 & 0.987 & 0.982 & 0.951\
1.20 & 0.797 & 0.785 & 0.823 & 0.923 & 0.986 & 0.986 & 0.984 & 0.942\
1.25 & 0.786 & 0.803 & 0.824 & 0.938 & 0.983 & 0.981 & 0.984 & 0.941\
1.30 & 0.789 & 0.793 & 0.800 & 0.894 & 0.981 & 0.976 & 0.992 & 0.925\
1.35 & 0.755 & 0.787 & 0.786 & 0.893 & 0.973 & 0.991 & 0.989 & 0.912\
1.40 & 0.761 & 0.789 & 0.785 & 0.864 & 0.970 & 0.974 & 0.994 & 0.879\

\[tab:pval\_jc\_short\]

Error in transition matrices {#error-in-transition-matrices .unnumbered}
----------------------------

For a given branch, we quantified the divergence $D$ between the original and estimated parameters of its transition matrix $A$ using the induced $L_1$ norm: $||A-\hat{A}||_{1}$ (see ). The columns in the transition matrices of ${\mathtt{JC69}^{\ast}}$, $\mathtt{K80}^{\ast}$, and the ${\mathtt{K81}^{\ast}}$ are equal and the norm becomes: $$\label{eq:L2norm}
D=\sum_{i=1}^4\mid A_{i,1}-\hat{A}_{i,1}\mid.
$$ Figure 3 depicts the results for ${\mathtt{JC69}^{\ast}}$ and ${\mathtt{K81}^{\ast}}$ on the three 4-taxon phylogenies, different alignment lengths and depths of the branches. The shapes of the distribution of $D$ for both models are very similar. As expected, the performance is weaker for long branches and short alignments. A great improvement is observed with the increase in the alignment length, e.g. $10,000nt$ depicts very accurate estimates. The performance under the ${\mathtt{JC69}^{\ast}}$ model (Fig. 3a )is better than that of ${\mathtt{K81}^{\ast}}$ (Fig.3b) for shorter branch lengths.

Parameter dispersion {#parameter-dispersion .unnumbered}
--------------------

Figure 4 shows the variances of the estimated parameters for depth 1 and 2 branches on the ${{\tau}_{\mathrm{balanced}}}$, ${{\tau}_{1:2}}$, ${{\tau}_{2:1}}$ trees under the ${\mathtt{JC69}^{\ast}}$ model.

The variances show an exponential increase, which is most significant in the ${{\tau}_{\mathrm{balanced}}}^4$ tree for both depths, and the ${{\tau}_{1:2}}$ for depth 2. The results for the depth 1 branch in ${{\tau}_{\mathrm{balanced}}}$ and ${{\tau}_{1:2}}$ are very similar. The smallest variance was observed for the depth 2 of ${{\tau}_{2:1}}$. For alignments of length $10,000nt$ on four taxa we can say that the method is quite accurate (see also Tab.VI-VIII in the ${\textit{Supp. mat.}}$ C).

For the ${\mathtt{K81}^{\ast}}$ model we summarized the results on variances for each edge as the mean of combined variances of all samples (see formula ). The results are analogous to those of the ${\mathtt{JC69}^{\ast}}$ model, see Figure II in ${\textit{Supp. mat.}}$ C. As expected, the parameter estimates are less dispersed for shorter branches and longer alignments (see Tab. IX-XI in the ${\textit{Supp. mat.}}$ C).

Error in the branch lengths {#error-in-the-branch-lengths .unnumbered}
---------------------------

Using the formula  we calculated the actual difference $l_0-\hat{l}$ between the branch length $l_0$ computed from the original parameters $\xi$ and the branch length $\hat{l}$ computed using their MLEs $\hat{\xi^e}$. Negative values of this score imply overestimation of the branch length, while positive values indicate underestimation. The results are shown in Figures 5 and 6.

In the case of ${\mathtt{JC69}^{\ast}}$ we observe that the method presented here does not tend to underestimate or overestimate the lengths for the depth 1 branches in all the 4-taxon trees ($l_0-\hat{l}$ is centered at 0 (see Fig. 5). The depth 2 branches have a tendency towards overestimation of the length for branches longer than (approximately) $0.45$ for ${{\tau}_{1:2}}$, $0.9$ for ${{\tau}_{2:1}}$, and $0.8$ for the ${{\tau}_{\mathrm{balanced}}}^4$ trees. In the latter case, lengths longer than $1.2$ for alignments up to $1,000$nt show opposite trend of underestimating the true lengths. The values were accurate when the alignment lengths were increased in the case of ${{\tau}_{1:2}}$ and ${{\tau}_{2:1}}$. On the other hand, for ${{\tau}_{\mathrm{balanced}}}^4$ the alignments of $10,000$nt resulted in overestimation.

In the ${\mathtt{K81}^{\ast}}$ model the results are significantly more accurate (see Fig.  6). There is a trend of underestimation for branches longer than (approximately) $0.9$ for shorter alignments. That is especially noticeable for ${{\tau}_{\mathrm{balanced}}}^4$ and depth 1 branches of ${{\tau}_{1:2}}$. This trend diminishes with an increase in the alignment length. Overall, in the case of both models, the variance of the estimate is smaller for shorter lengths and both depth 1 and 2 branches of the ${{\tau}_{2:1}}$ tree.

In addition, we calculated the tree length ${\mathsf{L}_{{\tau}}}$ (i.e. the sum of its branch lengths) from the estimated parameters and compared it to the theoretical result on the original branch length $l_0$: $4.5l_0$ for ${{\tau}_{1:2}}$ (where $l_0$ is a depth 1 branch), $3l_0$ for ${{\tau}_{2:1}}$ ($l_0$ for depth 2 branch) and $5l_0$ for ${{\tau}_{\mathrm{balanced}}}^4$. The rightmost columns of Figures 5 and 6 show the results for 4-taxon trees for the ${\mathtt{JC69}^{\ast}}$ and ${\mathtt{K81}^{\ast}}$ models respectively. The length of the tree is estimated accurately for all trees, the estimates being best for ${{\tau}_{2:1}}$. The variance is small and decreasing with an increase in the data length. As the sequences get longer, the distribution is centered around the true value. This is especially visible for the ${\mathtt{K81}^{\ast}}$ model (see Fig.  6).

Results for larger trees {#results-for-larger-trees .unnumbered}
------------------------

We run the analysis on the 6-taxon balanced tree, ${{\tau}_{\mathrm{balanced}}}^6$, under the ${\mathtt{K81}^{\ast}}$ model, for alignment lengths of $1,000$nt and $10,000$nt and branch length ${\mathbf{l}}\in\{0.01, 0.1, 0.3,
0.5, 0.7, 0.9, 1.1, 1.3, 1.4\}$. The $p-values$ of the corresponding tests confirm that the performance of the algorithm is very satisfactory (see Tab. XII in the ${\textit{Supp. mat.}}$ C). We have seen in the 4-taxon study that the tree with equal branch lengths gave worse results than the unbalanced trees. Thus, we expect the results of the depth 2 branches to be similarly challenged in this case.

Figure 7 depicts the estimated tree lengths. It can be seen that the estimates are accurate and the results improve for the alignments of $10,000$nt. As expected, the variance of the estimates increases with the increase in the length of the branch. By formula , long branches correspond to small values of the determinant of the transition matrix. Thus, statistical fluctuations in the parameter estimates have greater impact on the resulting length of the tree.

Next, we calculated the difference between the original and estimated branch lengths. In Figure 8a we see that the depth 1 branches show some degree of underestimation of the length for lengths $1.1-1.4$ and alignments of $1,000$nt. In the case of $10,000$nt, the results improve and can be expected to show little bias for even longer data sets. Branches of depth 2 show higher degree of underestimation with improvement for longer data sets. The divergence of the original and estimated parameters for transition matrices given by formula is shown in Figure 8b. For branches of depth 1 and data of length $10,000$, the error is about 0.2. In the case of branches of depth 2, it is almost doubled for both alignment lengths. In both cases, branch lengths up to 0.5 give satisfactory results. The error of the estimates for longer branches seems to be approaching a plateau.

Combined variance of the estimated parameters is much decreased for the $10,000$nt data sets in comparison with the $1,000$nt, and is smaller for the depth 1 branch (see Fig. 8c ). Again, the exponential shape of the plot can be attributed to the logarithm appearing in the formula .

CONCLUSION {#conclusion .unnumbered}
==========

In order to evaluate the performance of the method proposed here under various circumstances, we conducted many tests on simulated data sets. We observed that the performance of ${\texttt{Empar}}$ is most optimal for long alignments and short branch lengths. It is worth noting that even for short alignments of 300nt or 500nt on 4 taxa, the estimated parameters approximate closely the original parameters in $\approx 95\%$ of the cases as proved by the normality test of the MLE. Moreover, the branch lengths calculated based on the parameters estimated by ${\texttt{Empar}}$ were found very accurate already for short alignments. Though the measure of divergence $D$ for the parameters of transition matrices proposed here accumulates all errors in the entries of the transition matrix, alignment length of 10,000nt showed divergence values smaller than 0.1.

In this paper, we provide the first implementation of a tool for inferring continuous parameters under the discrete-time models. The method allows for accurate estimation of branch lengths in non-homogeneous data. There are two limitations to applicability of the method. Firstly, the algorithm has an exponential computational time increase with the number of taxa. This is a restriction due to the fact that the algorithm computes large matrices of dimension that is exponential in the total number of nodes of a tree. Running time of ${\texttt{Empar}}$ on star trees with 3-8 nodes and equal branches of $0.5$ on Ubuntu 11.10, Intel Core i7 920 at 2.67 GHz with 6 Gb is given in Table \[tab:gennong\_times\]. Secondly, the memory usage of ${\texttt{Empar}}$ is approx. $8*4^{|N({\tau})|}$ and corresponds to the memory footprint of the matrix in the ${\mathtt{EM}}$ algorithm, e.g. for this matrix to fit in the memory of a 6Gb machine the bound on the number of nodes is $|N({\tau})|$ $\leq 14$.

  length $\mid$ n   3       4       5       6       7       8
  ----------------- ------- ------- ------- ------- ------- ------
  1,000             0.004   0.02    0.033   0.222   1,049   7.14
  10,000            0       0.011   0.043   0.171   1,044   6.95

  : ${\texttt{Empar}}$ performance time– estimating the parameters of ${\mathtt{K81}^{\ast}}$ on star trees with equal branch lengths of $0.5$, varying number of leaves, $L({\tau})$, for the MSAs of $1,000$ and $10,000nt$.

\[tab:gennong\_times\]

We conclude that ${\texttt{Empar}}$ is a highly reliable method for estimating branch lengths of relatively small number of taxa and trees with short branch lengths (e.g. closely related species), and achieves high accuracy even when the results are based on short sequences. In particular, ${\texttt{Empar}}$ is a reliable method to compute quartets and to be used with quartet-based methods (see [@Berry99] and [@Berry2000]) on nonhomogeneous data.

FUNDING {#funding .unnumbered}
=======

Both authors were partially supported by Generalitat de Catalunya, 2009 SGR 1284. MC is partially supported by Ministerio de Educación y Ciencia MTM2009-14163-C02-02.We thank Roderic Guigó for generously providing funding for this project under grant BIO2011-26205 from the Ministerio de Educación y Ciencia (Spain).

[45]{} natexlab\#1[\#1]{}url \#1[[\#1]{}]{}

E S Allman and J A Rhodes. , 1860 (2):0 113–144, 2003.

E S Allman and J A Rhodes. . Cambridge University Press, January 2004. ISBN 0-521-52586-1.

E S Allman and J A Rhodes. Phylogenetic invariants. In O Gascuel and MA Steel, editors, [*Reconstructing Evolution*]{}. Oxford University Press, 2007.

C Ambroise. . , 150 (1):0 154–156, January 1998. ISSN 0176-4268.

D Barry and J A Hartigan. Asynchronous distance between homologous [DNA]{} sequences. , 430 (2):0 261–276, 1987. ISSN 0006-341X.

D Barry and J A Hartigan. Statistical analysis of hominoid molecular evolution. , 20 (2):0 191–207, 1987.

V Berry and O Gascuel. Inferring evolutionary trees with strong combinatorial evidence. , 240:0 271–298, 2000.

V Berry, T Jiang, P Kearney, M Li, and T Wareham. Quartet cleaning: improved algorithms and simulations. In [*Proc. Europ. Symp. Algs. (ESA99)*]{}, volume 1643 of [*LNCS*]{}, pages 313–324, 1999.

W J Bruno and A L Halpern. Topological bias and inconsistency of maximum likelihood using wrong models. , 160 (4):0 564–566, 1999.

M Casanellas, A M Kedzierska, and J Fernández-Sánchez. The space of phylogenetic mixtures of equivariant models. arXiv:0912.1957v1, 2011.

M Casanellas and S Sullivant. The strand symmetric model. In L Pachter and B Sturmfels, editors, [*Algebraic Statistics for computational biology*]{}, chapter 16. Cambridge University Press, 2005.

J T Chang. Full reconstruction of [M]{}arkov models on evolutionary trees: identifiability and consistency. , 1370 (1):0 51–73, 1996. ISSN 0025-5564.

A Dempster, N Laird, and D Rubin. . , 39:0 1–38, 1977.

B Efron and D V Hinkley. . , 65:0 457–483, 1978.

S H Eo and J A DeWoody. Evolutionary rates of mitochondrial genomes correspond to diversification rates and to contemporary species richness in birds and reptiles. , 277:0 3587–92, Dec 2010.

J Felsenstein. Cases in which parsimony or compatibility methods will be positively misleading. , 270 (4):0 401–410, 1978.

J Felsenstein. . , pages 164–166, 1993.

J Felsenstein. . Sinauer Associates, Inc., 2004.

N Galtier and M Gouy. Inferring pattern and process: maximum likelihood implementation of a non-homogeneous model of [DNA]{} sequence evolution for phylogenetic analysis. , 1540 (4):0 871–879, 1998.

H Hartley. . , 14:0 174–194, 1958.

M Hasegawa, H Kishino, and T Yano. , 220 (2):0 160–174, 1985.

J. Huelsenbeck and D. Hillis. Success of phylogenetic methods in the four-taxon case. , 420 (2):0 247–264, 1993.

D Husmeier, R Dybowski, and S Roberts. . Springer Verlag, New York, 2005.

V Jayaswal, L S Jermiin, L Poladian, and J Robinson. , 600 (1):0 74–86, 2011.

T H Jukes and C R Cantor. . In H N Munro, editor, [*Mammalian Protein Metabolism*]{}, volume 3, pages 21–132. Academic Press, 1969.

M Kimura. . , 160 (2):0 111–120, 1980.

M Kimura. . , 780 (1):0 454–458, 1981.

P J Lockhart, M A Steel, M D Hendy, and D Penny. Recovering evolutionary trees under a more realistic model of sequence evolution. , 110 (4):0 605–612, 1994.

L Pachter and B Sturmfels, editors. . Cambridge University Press, November 2005. ISBN 0-521-85700-7.

C R Rao. . Wiley, New York, 1973.

J Ripplinger and J Sullivan. , 570 (1):0 76–85, 2008.

J Ripplinger and J Sullivan. , 270 (12):0 2790–2803, 2010.

R S Schwartz and R L Mueller. . , 100 (1):0 5, 2010.

C Semple and M Steel. , volume 24 of [*[Oxford Lecture Series in Mathematics and its Applications]{}*]{}. Oxford University Press, Oxford, 2003. ISBN 0-19-850942-1.

M A Steel, L A Székely, and M D Hendy. , 10 (2):0 153–163, 1994.

J Sullivan and D L Swofford. Are guinea pigs rodents? the importance of adequate models in molecular phylogenetics. , 40 (2):0 77–86, 1997.

D L Swofford. . , 2003. Version 4.

M A Tanner. , volume 3rd Edition. Springer Series in Statistics, New York: Springer, 1996.

S Tavar[é]{}. . In [*Some mathematical questions in biology—[DNA]{} sequence analysis ([N]{}ew [Y]{}ork, 1984)*]{}, volume 17 of [*Lectures Math. Life Sci.*]{}, pages 57–86. Amer. Math. Soc., Providence, RI, 1986.

C F J Wu. . , 110 (1):0 95–103, 1983.

Z Yang. : A program package for phylogenetic analysis by maximum likelihood. , 15:0 555–556, 1997. URL <http://abacus.gene.ucl.ac.uk/software/paml.html>.

Z Yang and A D Yoder. Estimation of the transition/transversion rate bias and species sampling. , 48:0 274–283, 1999.

V B Yap and L Pachter. , 140 (4):0 574–9, April 2004.

S Zacks. . John Wiley and Sons, Inc., New York, 1971.

L Zou, E Susko, Ch Field, and A J Roger. . , 60:0 872–875, 2011.

![Estimated tree length as a function of the initial length of a branch of ${{\tau}_{\mathrm{balanced}}}^6$ (${\mathsf{L}_{{\tau}}}=9l_0$) in $1,000$ data sets generated under the ${\mathtt{K81}^{\ast}}$ model.[]{data-label="fig:sixtaxa_len"}](plots/bran-k81-6t_all.pdf)
