BACKGROUND
a phylogeny postulates shared ancestry relationships among organisms in the form of a binary tree.
phylogenies attempt to answer an important question posed in biology: what are the ancestor-descendent relationships between organisms?
at the core of every biological problem lies a phylogenetic component.
the patterns that can be observed in nature are the product of complex interactions, constrained by the template that our ancestors provide.
for example, the presence and structure of the human skull is mainly determined by its structure in our ancestors.
the relationship between the features observed in different organisms can only be understood if the phylogenetic relationships can be hypothesized.
an important method of phylogenetic inference is maximum parsimony .
under mp, the preferred hypothesis is the one that minimizes the number of evolutionary transformations required to explain the observed features  <cit> .
this optimization problem is known in computer science as the steiner tree problem, which is np-complete  <cit> .
the problem of simultaneous tree and alignment estimation under maximum parsimony is known in combinatorial optimization as the generalized tree alignment problem   <cit> .
the gtap is the steiner tree problem for the sequence edit distance.
like many biologically interesting problems, the gtap is np-hard  <cit> .
typically the steiner tree is presented under the manhattan or the hamming distances.
experimentally, the accuracy of the gtap has been subject to evaluation  <cit> .
the most recent results have shown evidence that phylogenies selected using the gtap from unaligned sequences are competitive  with the best methods and algorithms available based on coupled, but separate multiple sequence alignment and phylogeny reconstruction  <cit> .
due to its computational hardness, biologists interested in the gtap rely on heuristic procedures to find good solutions.
the simplest, and arguably the most important heuristic for the gtap is a local search.
a local search iteratively evaluates trees similar to a current solution t, where similar trees constitute the neighborhood of t. if a shorter tree s is found in the neighborhood, then t is replaced by s, and the search continues.
otherwise, t is the final solution.
local search is the work horse of most phylogenetic analysis procedures of practical use, and the core search procedures to solve the gtap in the computer programs msam  <cit> , and poy  <cit> .
it is known that the quality of a gtap analysis is heavily dependent on the fit of the local search heuristics used  <cit> , but the question of which heuristics are a better fit under what conditions remains unanswered.
in this paper, we discuss, implement, and experimentally explore existing and new local search heuristics for the gtap using simulated data.
our methods improve by more than three orders of magnitude the best local search heuristics existing to date with real data.
we begin by formally explaining the existing heuristics, and new heuristics for the gtap.
following the results of  <cit> , we use the affine-do algorithm to compute the tree length heuristically.
the algorithms
a subproblem of the gtap is the tree alignment problem  .
heuristically solving the tap with affine-do  <cit>  can be done in o, where n is the maximum sequence length and v the vertex set of the tree, and typically n≫|v|.
to simplify notation, in this section we assume that calculating the assignment of a vertex in a tree is a constant time operation  time).
existing heuristics
a local search consists of two steps: initial tree construction, and refinement .
given an initial tree t, refinement evaluates trees similar to t, in the search for a better solution.
those trees similar to t are its neighborhood.
the most commonly used neighborhood function is known as tree bisection and reconnection   <cit> .
tbr is based on two simple tree modifications: breaking an unrooted tree in two components, and joining two separate trees in one : tree breaking.
given a tree t, remove an edge  to produce two connected components, one with u, the other with v. if u  is not a leaf, then collapse it.
tree joining.
let t= and s= be two binary trees.
t and s can be joined by selecting a pair of edges ∈e and ∈e′, create subdivision vertices x in the edge  and x′ in , and add the edge .
if t  does not have edges, but only one vertex v, then take v as x .
the tbr neighborhood of t is the set of trees that can be produced by breaking t at any edge to produce two trees u and v, and then joining u and v. this neighborhood is used in the local search step of the gtap solver programs poy  <cit>  and msam  <cit> .
the most popular strategy for the initial tree construction is the wagner algorithm <cit> , a randomized, greedy strategy, of time complexity o. the wagner algorithm is used in most software packages for phylogenetic analysis under mp , including poy  <cit> .
msam takes a different approach, by using a neighbor joining tree, with time complexity o  <cit> .
deterministic algorithms are not typically used in the tree building step: for non trivial data sets, a good randomized method can be used repeatedly to initiate independent refinements resulting in different solutions.
their shared properties can give insights into the problem’s structure, and help discover better solutions.
depending on the distance function, different procedures are used to compute the score of the trees in the tbr neighborhood efficiently  <cit> .
in particular, for the hamming and manhattan distance, to calculate all of the tree scores in the tbr neighborhood has time complexity o  <cit> .
for the gtap however, it has time complexity o  <cit> , or o by increasing the hidden factor from o to o   <cit> .
exploring a neighborhood requires two additional criteria: the stopping rule, and the selection of the next candidate solution.
depending on their properties, a number of local search strategies can be described.
a classic heuristic that specifies the stopping and selection criteria is simulated annealing   <cit> .
contradictory conclusions about the applicability of sa to phylogenetic analysis can be found in the literature  <cit> .
a form of simulated annealing with better performance under the hamming and euclidean distance is known as tree-drifting  <cit> .
however, its metropolis and stopping criteria make tree-drifting inapplicable to the gtap.
the potential of simulated annealing for the gtap has remained unexplored.
sectorial search  <cit>   is a heuristic that restricts or extends the tbr neighborhood by only breaking and joining selected subtrees , or exhaustively solving such subtrees.
two variations of this scheme have been proposed: in the random-based ss, subtrees are selected uniformly at random.
in the second variation, the consensus-based ss, given a parameter 0≤n≤ <dig>  only rearrange  subtrees occurring in at least n∗m solutions found in m previous searches   <cit> .
other strategies , do not strictly belong to the set of local search heuristics.
given that local search is part of all these strategies, all of them would be more efficient if a good local search is in place.
new heuristics for the gtap
in this section, we describe four ideas to improve the local search strategies in the gtap: efficient tree length calculation during the search, better tree cost bounding, a smarter local search strategy, and initial tree building algorithms.
efficient tree updates
to apply the selection and stopping rules during tbr, it is necessary to calculate the tree length after every break, and join.
affine-do requires a directed tree as induced by its root.
if the sequence edit distance function is metric, the true tree length is independent of the root location.
given that metric distances are a common requirement under mp we assume from now on that the edit distance is metric.
it follows that, although affine-do can produce a different tree length for each possible root, there is no constraint to maintain one.
to update a tree efficiently, we do not maintain a unique rooted representation, but rather take its unrooted representation and keep all the potential roots assigned to every edge of the tree .
we call this a three directional assignment.
although we describe it for its application for the gtap, it is applicable to any algorithm that requires post-order traversal to compute  the tree length.
three directional assignment
for an unrooted binary tree, we assign to each edge  a sequence.
this sequence is the affine-do assignment to the subdivision vertex w of .
computing affine-do is dependent on the assignment to its neighbors .
in a binary tree, each interior vertex has three incident edges.
therefore, there are three possible affine-do assignments for every interior vertex .
each assignment is required to compute some subdivision vertices.
hence, we maintain the three possible assignments for each interior vertex.
these assignments can be computed with time complexity o, using first a pre-order traversal then followed by a post-order traversal, starting on any edge.
observation 1
a tree with a three directional assignment computes the length of every tree that can be produced by breaking any one edge with time complexity o.
observation 2
given two separate trees s and t with the three directional assignment, computing the length of all the trees produced by joining every pair of edges in s and t has time complexity o.
the simplest implementation of the three directions is to eagerly compute all the assignments in preparation for the first tree break, and join.
however, such an algorithm would entail overhead for greedy heuristics such as simulated annealing, where the first acceptable tree should be chosen to continue with the local search.
we solve this problem by using lazy evaluation and memoization  <cit>  as follows: eagerly assign a lazy function to each vertex and edge of the tree, but only compute its value  upon request, while memoizing the result.
in this way, we only spend time computing each vertex if used.
this technique has greater value if the tree break, and join order is carefully chosen.
in the following section, we will see how the three directional assignment can also be used to improve the estimation of each tree cost with no additional time complexity.
multiple heuristic tap solutions
the affine-do algorithm may calculate different tree length bounds depending on the root location .
nevertheless, the best of all the assignments is preferable for each tree.
computing all of the affine-do tree lengths, however, would add a o time complexity multiplicative factor to each tree break and join.
we avoid such factor and still produce better bounds for the tree cost during the search by using algorithm  <dig> on each break, and algorithm  <dig> on each join of the local search.
algorithm 1: improving the bound of a tree on each edge break
algorithm 2: improving the bound of a tree on each join
for a fixed n, the join procedure adds only a constant multiplicative factor, without increasing the time complexity.
note that if all the edges of a tree t are broken during a local search, then 2n− <dig> alignments are evaluated for the final tree, with no additional time complexity.
we call this variation of the tbr exhaustive-tbr.
smarter local searches
affine-do  <cit>  defines a compact representation of sets of sequences called a reduced alignment graph .
rags are less powerful than alignment graphs  <cit> , but are simpler and more efficient to compute and use.
it is then possible to align rag’s, find the closest sequences contained in them, and compute their rag with time complexity o, the same of a regular pairwise sequence alignment  <cit> .
ultimately, affine-do is a method to compute the distance between the closest sequences contained in a pair of rags efficiently.
rags can be used to guide a local search.
if the union of a pair of rags a and b can be efficiently computed in a new rag c, then c can be used to bound the distance between any other rag d and a or b simultaneously.
therefore, it is possible to use the union of multiple rags assigned to multiple vertices in a tree, to compute a lower bound of the closest pair of sequences contained in a pair of vertex sets .
algorithm 3: algorithm to compute merge.
the union of a single rag a is 〈ai,true〉
theorem 1
let r=merge .
all the sequences contained in x, y, and m are contained in r.
proof
at each step, either xi, yj, mk, {indel}, or any of their combinations is prepended to the result.
therefore, no element appearing in x, y, or m is missing in r. moreover, for all 0<e,f≤|x|, xe is prepended before xf if and only if e<f.
hence, the relative order of the elements in x is maintained in r. finally, for all the cases where xi is not prepended, then the indel element is included in r. it follows that that we can recover x by removing those elements in r where such indels where inserted and no element of x was.
by the definition of sequences contained in a rag  <cit> , it follows that every sequence in x is contained in r.
the analysis of y and m is symmetric.
□
theorem 2
algorithm  <dig> computes the union of x, y, and m with time complexity o where x is the longest union.
proof
the algorithm stops when i,j,k< <dig>  at each recursive step, either i or j is reduced by one, with initial values i=|x| and j=|y|.
□
the union of rags can be executed in o, on each vertex, during the affine-do computation.
affine-do is o, therefore, this method entails a small additive factor to the time complexity of affine-do.
in our implementation, we have fixed the size of the vertex sets to  <dig> vertices on all data sets experimentally.
using unions during a local search
let t be the current candidate solution during a local search, and u the set of unions of t by applying algorithm  <dig> while traversing the tree in affine-do.
if a new candidate tree s is accepted during the local search, then update u using the direction for the best subdivision vertex computed for s .
by maintaining this set of unions, we can modify the tbr local search as in algorithm  <dig>  to join only edges that are incident in unions at short distance.
we call this method union–pruning
algorithm 4: heuristic union-pruning tbr.
the threshold  <dig>  parameter was experimentally tuned
building the initial trees
the wagner algorithm is a basic procedure to compute an initial tree .
we modify this procedure in two ways.
algorithm 5: the wagner algorithm for initial tree building
union–pruning.
unions can be used to efficiently prune candidate trees during the wagner algorithm by maintaining the union set of the tree t in algorithm  <dig>  and treat each leaf to be added as a union of its own.
then use algorithm  <dig> to guide the join step in algorithm  <dig>
addition sequence
the initial sequence l in algorithm  <dig> is typically randomized, assigning equal probability to each permutation.
this algorithm is known as random addition sequence .
the randomization of l is used to obtain multiple starting points for local searches.
we have explored the following variation successfully:
<dig>  compute a minimum spanning tree  of l .
<dig>  traverse l using a bfs.
the order in which we visit the elements of l is our initial addition sequence q.
<dig>  to produce the n’th tree, produce the sequence q by flipping consecutive elements in q with probability  <dig> .
we call this procedure mst-wagner.
methods
we evaluated experimentally a number of algorithms for local searches under the gtap.
an experimental evaluation of this kind has three fundamental components: a selection of heuristics, implementation, and selection of data sets.
the overall performance is compared with the length of the trees found by each method.
algorithms compared
we compared the following heuristic local searches, in all meaningful combinations.
tap computation: using affine-do in two variations, exhaustive, and non-exhaustive.
building: wagner algorithm using ras and mst addition sequences, and the neighbor joining  algorithm.
the wagner algorithm was executed with lookahead parameters of  <dig>   <dig>   <dig>  and  <dig>  neighborhood: tbr and spr .
edge breaking order: randomized, or in length decreasing order.
join order: randomized, or in ascending order based on the distance of the union that each edge belongs to.
in the second case, the union-pruning strategy was used to filter candidates.
sector and reroot diameters:  <dig>   <dig>   <dig>  and infinity .
the rerooting order followed a breadth first search  order, around the broken edge.
the sector and reroot diameters were selected to match the simulation size .
simulated annealing: using initial temperatures of  <dig>   <dig>  and  <dig>  and coefficients of  <dig>   <dig>   <dig>  and  <dig>  the values were selected experimentally as a good sample of the performance variation observed by the authors in real gtap problems.
for the edit distance parameters we tested the following combinations of substitution, indel, and gap opening parameters : , , , .
in our experience, these parameters encompass enough variation in the gtap, while maintaining a limited number of combinations with the algorithms.
in total,  <dig> combinations of build algorithms and distance functions were tested.
for the refinement step, a total of  <dig> combinations of algorithms and edit distance functions were tested.
implementation
we implemented the algorithms under comparison in the objective caml and c programming languages.
all the algorithms are available in the author’s computer program poy version  <dig>  <cit> .
the functions are highly optimized for performance.
data sets
to generate the instance problems, we simulated sequences using dawg  <dig> .
<dig>  <cit>  with insertions and deletions following a power law distribution.
the simulations followed random binary trees of  <dig> leaves comprising all the combinations of the parameters listed in table  <dig>  this tree size was chosen to be both tractable and realistic in size without biasing trees to any particular shape.
the indel and branch parameters also were chosen to be similar to what is seen in empirical data sets.
these produced a total of  <dig> independent simulations.
each simulation was analyzed independently with  <dig> repetitions for each randomized algorithm.
nj was tested only once, as our implementation is deterministic.
an initial exploration with  <dig> repetitions showed no significant difference compared to  <dig> repetitions.
in total,  <dig>  builds, and  <dig>  refinements were performed.
due to the large number of simulations and local searches performed, we will concentrate on a reduced set of cases that represent the overall patterns observed.
all combinations of parameters were employed to generate the test data sets.
the branch length variation equals the average branch length.
RESULTS
this section begins with the difference in performance between the exhaustive  and the non-exhaustive  algorithms, which can be applied in conjunction with any other search strategy.
it continues with a comparison of the build algorithms, and the refinement algorithms.
finally, we compose the results in a simple local search heuristic which we compared with the previous best heuristic on a real dataset.
exhaustive and non-exhaustive algorithms
in the build step , the difference between e and ne is small for all equivalent algorithms with branch lengths of  <dig>  and  <dig>  .
the most striking difference, however, occurs for branch length  <dig>  , where ne shows an expected tree length 50% higher than that of e. such extreme variation shows a strong dependence on the root location when branch lengths make sequences close to random relative to each other.
for the tbr step, e significantly outperforms ne, with better minimum and expected scores .
this pattern was observed for every combination of algorithm, simulation, and edit distance parameters.
in the following two sections, we concentrate on the results obtained using the e algorithm.
the same general patterns were observed with ne, but with less competitive tree scores.
initial tree building
the initial tree building algorithms fall into two main groups: algorithms with ras, and algorithms using mst.
in all cases, mst produced significantly shorter trees .
the use of higher lookahead parameters did not produce consistent improvements in the resulting trees, while the use of the union-pruning algorithm did significantly improve the expectation, and the minimum tree cost for branch lengths  <dig>  and  <dig> .
for long branch lengths, however, no significant improvement was observed.
neighbor joining produced trees of highest score among all the algorithms for all parameters .
we do not present it in the graphs as it would make the more subtle differences between other algorithms difficult to observe.
overall, the most important improvement occurs with the mst addition sequence in first place, followed by the use of the union-pruning strategy in second.
nevertheless, we will see in the next section that the use of the mst algorithm remains limited.
refinement
to evaluate the tbr refinement experimentally, we must produce an initial tree.
although mst showed better results than ras, we found that in almost every instance tbr failed to improve the mst trees.
at the end, ras + tbr would always find better trees than mst + tbr.
for this reason, we used the second best method to construct the initial trees: ras using union-pruning.
the refinement comparison can be divided in two groups:  <dig> ) a comparison between basic tbr using union-pruning, and branch length sorting, and  <dig> ) the comparison of different algorithms using the best combination among those in  <dig>
union-pruning and branch length sorting
the behavior of tbr with union-pruning and branch length sorting is presented in figure  <dig>  the union-pruning algorithm produced significantly better trees, both in the minimum and expected scores.
this advantage disappears as sequences diverge to close to random  .
branch length sorting had a small positive impact, but not significative.
the results match our expectation: the union-pruning algorithm can positively guide the search with better taxon sampling.
we have observed this behavior in real data sets, where new terminals some times speedup the local search.
this somewhat counter-intuitive behavior is likely due to the structured nature of phylogenetic data.
the addition of new terminals increases the data support for subtrees, thereby increasing the cost penalty when these groups are violated.
union-pruning takers advantage of this in creating unions from larger sets of taxa, hence containing more information.
as the data become less structured  the effect disappears.
local search strategy
beyond the use of union-pruning, and exhaustive tap estimation, the differences among the algorithms compared are not significant .
although in general sectorial finds the shortest tree with highest frequency, the difference is typically less than two length units, compared to the second best algorithm.
in general, the algorithm with the best mean is bfs, but again, not significative.
however, due to the algorithm design, bfs is the fastest of all.
the differences observed are not significant.
all the simulations shown have branch length  <dig> , but similar patterns were observed for branch lengths  <dig>  and  <dig> .
the minima across each row is in bold.
overall performance
based on the previous experiments, we prefer a heuristic local search strategy that consists of the following steps: build initial trees using ras guided by union-pruning, followed by a refinement step consisting of tbr using the three directional heuristics, exhaustive tap, union-pruning, and cutting edges according to descending lengths.
we compared this algorithm , with that of poy version  <dig> which uses a one directional algorithm, with randomized tbr steps  <cit> .
due to limitations in poy version 3’s implementation, we only compare an edition distance with substitution parameter  <dig>  indel parameter  <dig>  and gap opening parameter  <dig>  due to the implementation limitation, msam was not included in the comparison.
for this comparison, a random subset of  <dig> published anurans  <cit>  was analyzed.
the data set includes 12s rrna, trna valine, 16s rrna, and fragments of cytochrome b, rhodopsin, tyrosinase, 28s rrna, and rag  <dig>  and a small set of  <dig> morphological, non-additive characters .
to compare the performance of poy version  <dig> and version  <dig>  we executed  <dig> independent repetitions consisting of  <dig> build, followed by refinement, and reported the resulting tree score.
this procedure can be executed in poy  <dig> with the command: poy ‐replicates  <dig> ‐seed ‐ <dig> ‐maxtrees  <dig> ‐nooneasis ‐minterminals  <dig> ‐terminalsfile rannamesph.txt ⋆.fas ⋆.ss.
the score of the trees found by each program were plotted in a density histogram .
the results show that one repetition of our new heuristic in poy version  <dig> outputs a tree which is expected to belong to the top 15% of the best trees found by this very simple search strategy.
to expect a tree within the same percentile using the old heuristic, it would be necessary to run more than  <dig> local searches.
it follows that through combination and speed and efficiency, the new heuristic is more than  <dig> times faster than the previous heuristic of poy  <dig>
discussion
we described and implemented new heuristics for the gtap.
we have shown that they find better solutions than previous approaches.
we found that a number of conditions affect the fit of the heuristic to the problem: long branch–length data sets can be better analyzed with sectorial search instead of the union-pruning, while union-pruning yields excellent results in medium, and short branch lengths.
exhaustive-tbr yields the best results overall and should always be preferred.
although the mst algorithm yields better initial results than ras, it is not preferable in the long run, and a small number of local searches should never be used to produce reliable results.
it remains to be explored the quality of the numerous meta–heuristics available in the literature.
it is now possible to explore them using a more efficient local search strategy.
CONCLUSIONS
we described new strategies that can be composed to produce a powerful local search strategy for the tree alignment problem.
the results showed that our methods improve on the best existing local search heuristics by more than three orders of magnitude.
in general, the exhaustive–tbr refinement strategy should always be used, while union-pruning should only be preferred if dense taxon sampling or short branch lengths are expected.
moreover, although the mst build strategy yields better results than the traditional wagner build, the former should not be preferred in real analyses since it tends to produce less competitive trees after the refinement step.
it is difficult to predict the performance of other high level heuristics applied to the gtap.
strategies such as sectorial search, and tree fusing should be effective.
however, divide and conquer techniques such as dcm- <dig> may have a more limited application, unless used in the spirit of sectorial search.
given that phylogenetic analysis under mp shows a simplified setting compared to other optimality criteria, it is our opinion that metaheuristics such as simulated annealing have limited applicability in the join estimation of tree and alignments for all optimality criteria, and novel strategies are needed to successfully scale to larger problem sizes.
nevertheless , all these strategies will belong to the heuristic realm, and further experimental efforts will be required.
affine-do, union–pruning, and exhaustive–tbr are some of the algorithms that we have implemented in the computer program poy version  <dig>  <cit> .
the algorithms and their implementation have had a significant impact in the biology community interested in different approaches to joint tree and phylogeny reconstruction.
by using better algorithms, algorithm engineering, and better parallel strategies, poy version  <dig> is three orders of magnitude faster than its predecessor.
the concepts, and desirable properties of this implementation should be extended to other phylogenetic inference criteria, to broaden its usability, and better serve the research purposes set for the software package.
competing interests
the authors declare that they have no competing interests.
authors’ contributions
the authors contributed equally to this work.
both authors read and approved the final manuscript.

 a phylogeny postulates shared ancestry relationships among organisms in the form of a binary tree. phylogenies attempt to answer an important question posed in biology: what are the ancestor-descendent relationships between organisms? at the core of every biological problem lies a phylogenetic component. the patterns that can be observed in nature are the product of complex interactions, constrained by the template that our ancestors provide. for example, the presence and structure of the human skull is mainly determined by its structure in our ancestors. the relationship between the features observed in different organisms can only be understood if the phylogenetic relationships can be hypothesized. an important method of phylogenetic inference is maximum parsimony. under mp, the preferred hypothesis is the one that minimizes the number of evolutionary transformations required to explain the observed features. this optimization problem is known in computer science as the steiner tree problem, which is np-complete. the problem of simultaneous tree and alignment estimation under maximum parsimony is known in combinatorial optimization as the generalized tree alignment problem. the gtap is the steiner tree problem for the sequence edit distance. like many biologically interesting problems, the gtap is np-hard. typically the steiner tree is presented under the manhattan or the hamming distances. experimentally, the accuracy of the gtap has been subject to evaluation. the most recent results have shown evidence that phylogenies selected using the gtap from unaligned sequences are competitive with the best methods and algorithms available based on coupled, but separate multiple sequence alignment and phylogeny reconstruction. due to its computational hardness, biologists interested in the gtap rely on heuristic procedures to find good solutions. the simplest, and arguably the most important heuristic for the gtap is a local search. a local search iteratively evaluates trees similar to a current solution t, where similar trees constitute the neighborhood of t. if a shorter tree s is found in the neighborhood, then t is replaced by s, and the search continues. otherwise, t is the final solution. local search is the work horse of most phylogenetic analysis procedures of practical use, and the core search procedures to solve the gtap in the computer programs msam, and poy. it is known that the quality of a gtap analysis is heavily dependent on the fit of the local search heuristics used, but the question of which heuristics are a better fit under what conditions remains unanswered. in this paper, we discuss, implement, and experimentally explore existing and new local search heuristics for the gtap using simulated data. our methods improve by more than three orders of magnitude the best local search heuristics existing to date with real data. we begin by formally explaining the existing heuristics, and new heuristics for the gtap. following the results of, we use the affine-do algorithm to compute the tree length heuristically. the algorithms a subproblem of the gtap is the tree alignment problem. heuristically solving the tap with affine-do can be done in o, where n is the maximum sequence length and v the vertex set of the tree, and typically n≫|v|. to simplify notation, in this section we assume that calculating the assignment of a vertex in a tree is a constant time operation time). existing heuristics a local search consists of two steps: initial tree construction, and refinement. given an initial tree t, refinement evaluates trees similar to t, in the search for a better solution. those trees similar to t are its neighborhood. the most commonly used neighborhood function is known as tree bisection and reconnection. tbr is based on two simple tree modifications: breaking an unrooted tree in two components, and joining two separate trees in one: tree breaking. given a tree t, remove an edge to produce two connected components, one with u, the other with v. if u is not a leaf, then collapse it. tree joining. let t= and s= be two binary trees. t and s can be joined by selecting a pair of edges ∈e and ∈e′, create subdivision vertices x in the edge and x′ in, and add the edge. if t does not have edges, but only one vertex v, then take v as x. the tbr neighborhood of t is the set of trees that can be produced by breaking t at any edge to produce two trees u and v, and then joining u and v. this neighborhood is used in the local search step of the gtap solver programs poy and msam. the most popular strategy for the initial tree construction is the wagner algorithm, a randomized, greedy strategy, of time complexity o. the wagner algorithm is used in most software packages for phylogenetic analysis under mp, including poy. msam takes a different approach, by using a neighbor joining tree, with time complexity o. deterministic algorithms are not typically used in the tree building step: for non trivial data sets, a good randomized method can be used repeatedly to initiate independent refinements resulting in different solutions. their shared properties can give insights into the problem’s structure, and help discover better solutions. depending on the distance function, different procedures are used to compute the score of the trees in the tbr neighborhood efficiently. in particular, for the hamming and manhattan distance, to calculate all of the tree scores in the tbr neighborhood has time complexity o. for the gtap however, it has time complexity o, or o by increasing the hidden factor from o to o. exploring a neighborhood requires two additional criteria: the stopping rule, and the selection of the next candidate solution. depending on their properties, a number of local search strategies can be described. a classic heuristic that specifies the stopping and selection criteria is simulated annealing. contradictory conclusions about the applicability of sa to phylogenetic analysis can be found in the literature. a form of simulated annealing with better performance under the hamming and euclidean distance is known as tree-drifting. however, its metropolis and stopping criteria make tree-drifting inapplicable to the gtap. the potential of simulated annealing for the gtap has remained unexplored. sectorial search is a heuristic that restricts or extends the tbr neighborhood by only breaking and joining selected subtrees, or exhaustively solving such subtrees. two variations of this scheme have been proposed: in the random-based ss, subtrees are selected uniformly at random. in the second variation, the consensus-based ss, given a parameter 0≤n≤ only rearrange subtrees occurring in at least n∗m solutions found in m previous searches. other strategies, do not strictly belong to the set of local search heuristics. given that local search is part of all these strategies, all of them would be more efficient if a good local search is in place. new heuristics for the gtap in this section, we describe four ideas to improve the local search strategies in the gtap: efficient tree length calculation during the search, better tree cost bounding, a smarter local search strategy, and initial tree building algorithms. efficient tree updates to apply the selection and stopping rules during tbr, it is necessary to calculate the tree length after every break, and join. affine-do requires a directed tree as induced by its root. if the sequence edit distance function is metric, the true tree length is independent of the root location. given that metric distances are a common requirement under mp we assume from now on that the edit distance is metric. it follows that, although affine-do can produce a different tree length for each possible root, there is no constraint to maintain one. to update a tree efficiently, we do not maintain a unique rooted representation, but rather take its unrooted representation and keep all the potential roots assigned to every edge of the tree. we call this a three directional assignment. although we describe it for its application for the gtap, it is applicable to any algorithm that requires post-order traversal to compute the tree length. three directional assignment for an unrooted binary tree, we assign to each edge a sequence. this sequence is the affine-do assignment to the subdivision vertex w of. computing affine-do is dependent on the assignment to its neighbors. in a binary tree, each interior vertex has three incident edges. therefore, there are three possible affine-do assignments for every interior vertex. each assignment is required to compute some subdivision vertices. hence, we maintain the three possible assignments for each interior vertex. these assignments can be computed with time complexity o, using first a pre-order traversal then followed by a post-order traversal, starting on any edge. observation 1 a tree with a three directional assignment computes the length of every tree that can be produced by breaking any one edge with time complexity o. observation 2 given two separate trees s and t with the three directional assignment, computing the length of all the trees produced by joining every pair of edges in s and t has time complexity o. the simplest implementation of the three directions is to eagerly compute all the assignments in preparation for the first tree break, and join. however, such an algorithm would entail overhead for greedy heuristics such as simulated annealing, where the first acceptable tree should be chosen to continue with the local search. we solve this problem by using lazy evaluation and memoization as follows: eagerly assign a lazy function to each vertex and edge of the tree, but only compute its value upon request, while memoizing the result. in this way, we only spend time computing each vertex if used. this technique has greater value if the tree break, and join order is carefully chosen. in the following section, we will see how the three directional assignment can also be used to improve the estimation of each tree cost with no additional time complexity. multiple heuristic tap solutions the affine-do algorithm may calculate different tree length bounds depending on the root location. nevertheless, the best of all the assignments is preferable for each tree. computing all of the affine-do tree lengths, however, would add a o time complexity multiplicative factor to each tree break and join. we avoid such factor and still produce better bounds for the tree cost during the search by using algorithm on each break, and algorithm on each join of the local search. algorithm 1: improving the bound of a tree on each edge break algorithm 2: improving the bound of a tree on each join for a fixed n, the join procedure adds only a constant multiplicative factor, without increasing the time complexity. note that if all the edges of a tree t are broken during a local search, then 2n− alignments are evaluated for the final tree, with no additional time complexity. we call this variation of the tbr exhaustive-tbr. smarter local searches affine-do defines a compact representation of sets of sequences called a reduced alignment graph. rags are less powerful than alignment graphs, but are simpler and more efficient to compute and use. it is then possible to align rag’s, find the closest sequences contained in them, and compute their rag with time complexity o, the same of a regular pairwise sequence alignment. ultimately, affine-do is a method to compute the distance between the closest sequences contained in a pair of rags efficiently. rags can be used to guide a local search. if the union of a pair of rags a and b can be efficiently computed in a new rag c, then c can be used to bound the distance between any other rag d and a or b simultaneously. therefore, it is possible to use the union of multiple rags assigned to multiple vertices in a tree, to compute a lower bound of the closest pair of sequences contained in a pair of vertex sets. algorithm 3: algorithm to compute merge. the union of a single rag a is 〈ai,true〉 theorem 1 let r=merge. all the sequences contained in x, y, and m are contained in r. proof at each step, either xi, yj, mk, {indel}, or any of their combinations is prepended to the result. therefore, no element appearing in x, y, or m is missing in r. moreover, for all 0<e,f≤|x|, xe is prepended before xf if and only if e<f. hence, the relative order of the elements in x is maintained in r. finally, for all the cases where xi is not prepended, then the indel element is included in r. it follows that that we can recover x by removing those elements in r where such indels where inserted and no element of x was. by the definition of sequences contained in a rag, it follows that every sequence in x is contained in r. the analysis of y and m is symmetric. □ theorem 2 algorithm computes the union of x, y, and m with time complexity o where x is the longest union. proof the algorithm stops when i,j,k< at each recursive step, either i or j is reduced by one, with initial values i=|x| and j=|y|. □ the union of rags can be executed in o, on each vertex, during the affine-do computation. affine-do is o, therefore, this method entails a small additive factor to the time complexity of affine-do. in our implementation, we have fixed the size of the vertex sets to vertices on all data sets experimentally. using unions during a local search let t be the current candidate solution during a local search, and u the set of unions of t by applying algorithm while traversing the tree in affine-do. if a new candidate tree s is accepted during the local search, then update u using the direction for the best subdivision vertex computed for s. by maintaining this set of unions, we can modify the tbr local search as in algorithm to join only edges that are incident in unions at short distance. we call this method union–pruning algorithm 4: heuristic union-pruning tbr. the threshold parameter was experimentally tuned building the initial trees the wagner algorithm is a basic procedure to compute an initial tree. we modify this procedure in two ways. algorithm 5: the wagner algorithm for initial tree building union–pruning. unions can be used to efficiently prune candidate trees during the wagner algorithm by maintaining the union set of the tree t in algorithm and treat each leaf to be added as a union of its own. then use algorithm to guide the join step in algorithm addition sequence the initial sequence l in algorithm is typically randomized, assigning equal probability to each permutation. this algorithm is known as random addition sequence. the randomization of l is used to obtain multiple starting points for local searches. we have explored the following variation successfully: compute a minimum spanning tree of l. traverse l using a bfs. the order in which we visit the elements of l is our initial addition sequence q. to produce the n’th tree, produce the sequence q by flipping consecutive elements in q with probability. we call this procedure mst-wagner. methods we evaluated experimentally a number of algorithms for local searches under the gtap. an experimental evaluation of this kind has three fundamental components: a selection of heuristics, implementation, and selection of data sets. the overall performance is compared with the length of the trees found by each method. algorithms compared we compared the following heuristic local searches, in all meaningful combinations. tap computation: using affine-do in two variations, exhaustive, and non-exhaustive. building: wagner algorithm using ras and mst addition sequences, and the neighbor joining algorithm. the wagner algorithm was executed with lookahead parameters of and neighborhood: tbr and spr. edge breaking order: randomized, or in length decreasing order. join order: randomized, or in ascending order based on the distance of the union that each edge belongs to. in the second case, the union-pruning strategy was used to filter candidates. sector and reroot diameters: and infinity. the rerooting order followed a breadth first search order, around the broken edge. the sector and reroot diameters were selected to match the simulation size. simulated annealing: using initial temperatures of and and coefficients of and the values were selected experimentally as a good sample of the performance variation observed by the authors in real gtap problems. for the edit distance parameters we tested the following combinations of substitution, indel, and gap opening parameters:,,,. in our experience, these parameters encompass enough variation in the gtap, while maintaining a limited number of combinations with the algorithms. in total, combinations of build algorithms and distance functions were tested. for the refinement step, a total of combinations of algorithms and edit distance functions were tested. implementation we implemented the algorithms under comparison in the objective caml and c programming languages. all the algorithms are available in the author’s computer program poy version. the functions are highly optimized for performance. data sets to generate the instance problems, we simulated sequences using dawg. with insertions and deletions following a power law distribution. the simulations followed random binary trees of leaves comprising all the combinations of the parameters listed in table this tree size was chosen to be both tractable and realistic in size without biasing trees to any particular shape. the indel and branch parameters also were chosen to be similar to what is seen in empirical data sets. these produced a total of independent simulations. each simulation was analyzed independently with repetitions for each randomized algorithm. nj was tested only once, as our implementation is deterministic. an initial exploration with repetitions showed no significant difference compared to repetitions. in total, builds, and refinements were performed. due to the large number of simulations and local searches performed, we will concentrate on a reduced set of cases that represent the overall patterns observed. all combinations of parameters were employed to generate the test data sets. the branch length variation equals the average branch length.  this section begins with the difference in performance between the exhaustive and the non-exhaustive algorithms, which can be applied in conjunction with any other search strategy. it continues with a comparison of the build algorithms, and the refinement algorithms. finally, we compose the results in a simple local search heuristic which we compared with the previous best heuristic on a real dataset. exhaustive and non-exhaustive algorithms in the build step, the difference between e and ne is small for all equivalent algorithms with branch lengths of and. the most striking difference, however, occurs for branch length, where ne shows an expected tree length 50% higher than that of e. such extreme variation shows a strong dependence on the root location when branch lengths make sequences close to random relative to each other. for the tbr step, e significantly outperforms ne, with better minimum and expected scores. this pattern was observed for every combination of algorithm, simulation, and edit distance parameters. in the following two sections, we concentrate on the results obtained using the e algorithm. the same general patterns were observed with ne, but with less competitive tree scores. initial tree building the initial tree building algorithms fall into two main groups: algorithms with ras, and algorithms using mst. in all cases, mst produced significantly shorter trees. the use of higher lookahead parameters did not produce consistent improvements in the resulting trees, while the use of the union-pruning algorithm did significantly improve the expectation, and the minimum tree cost for branch lengths and. for long branch lengths, however, no significant improvement was observed. neighbor joining produced trees of highest score among all the algorithms for all parameters. we do not present it in the graphs as it would make the more subtle differences between other algorithms difficult to observe. overall, the most important improvement occurs with the mst addition sequence in first place, followed by the use of the union-pruning strategy in second. nevertheless, we will see in the next section that the use of the mst algorithm remains limited. refinement to evaluate the tbr refinement experimentally, we must produce an initial tree. although mst showed better results than ras, we found that in almost every instance tbr failed to improve the mst trees. at the end, ras + tbr would always find better trees than mst + tbr. for this reason, we used the second best method to construct the initial trees: ras using union-pruning. the refinement comparison can be divided in two groups: ) a comparison between basic tbr using union-pruning, and branch length sorting, and ) the comparison of different algorithms using the best combination among those in union-pruning and branch length sorting the behavior of tbr with union-pruning and branch length sorting is presented in figure the union-pruning algorithm produced significantly better trees, both in the minimum and expected scores. this advantage disappears as sequences diverge to close to random. branch length sorting had a small positive impact, but not significative. the results match our expectation: the union-pruning algorithm can positively guide the search with better taxon sampling. we have observed this behavior in real data sets, where new terminals some times speedup the local search. this somewhat counter-intuitive behavior is likely due to the structured nature of phylogenetic data. the addition of new terminals increases the data support for subtrees, thereby increasing the cost penalty when these groups are violated. union-pruning takers advantage of this in creating unions from larger sets of taxa, hence containing more information. as the data become less structured the effect disappears. local search strategy beyond the use of union-pruning, and exhaustive tap estimation, the differences among the algorithms compared are not significant. although in general sectorial finds the shortest tree with highest frequency, the difference is typically less than two length units, compared to the second best algorithm. in general, the algorithm with the best mean is bfs, but again, not significative. however, due to the algorithm design, bfs is the fastest of all. the differences observed are not significant. all the simulations shown have branch length, but similar patterns were observed for branch lengths and. the minima across each row is in bold. overall performance based on the previous experiments, we prefer a heuristic local search strategy that consists of the following steps: build initial trees using ras guided by union-pruning, followed by a refinement step consisting of tbr using the three directional heuristics, exhaustive tap, union-pruning, and cutting edges according to descending lengths. we compared this algorithm, with that of poy version which uses a one directional algorithm, with randomized tbr steps. due to limitations in poy version 3’s implementation, we only compare an edition distance with substitution parameter indel parameter and gap opening parameter due to the implementation limitation, msam was not included in the comparison. for this comparison, a random subset of published anurans was analyzed. the data set includes 12s rrna, trna valine, 16s rrna, and fragments of cytochrome b, rhodopsin, tyrosinase, 28s rrna, and rag and a small set of morphological, non-additive characters. to compare the performance of poy version and version we executed independent repetitions consisting of build, followed by refinement, and reported the resulting tree score. this procedure can be executed in poy with the command: poy ‐replicates ‐seed ‐ ‐maxtrees ‐nooneasis ‐minterminals ‐terminalsfile rannamesph.txt ⋆.fas ⋆.ss. the score of the trees found by each program were plotted in a density histogram. the results show that one repetition of our new heuristic in poy version outputs a tree which is expected to belong to the top 15% of the best trees found by this very simple search strategy. to expect a tree within the same percentile using the old heuristic, it would be necessary to run more than local searches. it follows that through combination and speed and efficiency, the new heuristic is more than times faster than the previous heuristic of poy discussion we described and implemented new heuristics for the gtap. we have shown that they find better solutions than previous approaches. we found that a number of conditions affect the fit of the heuristic to the problem: long branch–length data sets can be better analyzed with sectorial search instead of the union-pruning, while union-pruning yields excellent results in medium, and short branch lengths. exhaustive-tbr yields the best results overall and should always be preferred. although the mst algorithm yields better initial results than ras, it is not preferable in the long run, and a small number of local searches should never be used to produce reliable results. it remains to be explored the quality of the numerous meta–heuristics available in the literature. it is now possible to explore them using a more efficient local search strategy.  we described new strategies that can be composed to produce a powerful local search strategy for the tree alignment problem. the results showed that our methods improve on the best existing local search heuristics by more than three orders of magnitude. in general, the exhaustive–tbr refinement strategy should always be used, while union-pruning should only be preferred if dense taxon sampling or short branch lengths are expected. moreover, although the mst build strategy yields better results than the traditional wagner build, the former should not be preferred in real analyses since it tends to produce less competitive trees after the refinement step. it is difficult to predict the performance of other high level heuristics applied to the gtap. strategies such as sectorial search, and tree fusing should be effective. however, divide and conquer techniques such as dcm- may have a more limited application, unless used in the spirit of sectorial search. given that phylogenetic analysis under mp shows a simplified setting compared to other optimality criteria, it is our opinion that metaheuristics such as simulated annealing have limited applicability in the join estimation of tree and alignments for all optimality criteria, and novel strategies are needed to successfully scale to larger problem sizes. nevertheless, all these strategies will belong to the heuristic realm, and further experimental efforts will be required. affine-do, union–pruning, and exhaustive–tbr are some of the algorithms that we have implemented in the computer program poy version. the algorithms and their implementation have had a significant impact in the biology community interested in different approaches to joint tree and phylogeny reconstruction. by using better algorithms, algorithm engineering, and better parallel strategies, poy version is three orders of magnitude faster than its predecessor. the concepts, and desirable properties of this implementation should be extended to other phylogenetic inference criteria, to broaden its usability, and better serve the research purposes set for the software package. competing interests the authors declare that they have no competing interests. authors’ contributions the authors contributed equally to this work. both authors read and approved the final manuscript. 
here, we implement and explore experimentally existing and new local search heuristics for the gtap using simulated and real data. the methods presented here improve by more than three orders of magnitude in execution time the best local search heuristics existing to date when applied to real data. at the core of every biological problem lies a phylogenetic component. phylogenies attempt to answer an important question posed in biology: what are the ancestor-descendent relationships between organisms? the problem of simultaneous tree and alignment estimation under maximum parsimony is known in combinatorial optimization as the generalized tree alignment problem. the gtap is the steiner tree problem for the sequence edit distance. typically the steiner tree is presented under the manhattan or the hamming distances. like many biologically interesting problems, the gtap is np-hard. a phylogeny postulates shared ancestry relationships among organisms in the form of a binary tree. results show that phylogenies selected using the gtap from unaligned sequences are competitive with the best methods and algorithms available. the patterns that can be observed in nature are the product of complex interactions, constrained by the template that our ancestors provide.
 a phylogeny postulates shared ancestry relationships among organisms in the form of a binary tree. phylogenies attempt to answer an important question posed in biology: what are the ancestor-descendent relationships between organisms? at the core of every biological problem lies a phylogenetic component. the patterns that can be observed in nature are the product of complex interactions, constrained by the template that our ancestors provide. the problem of simultaneous tree and alignment estimation under maximum parsimony is known in combinatorial optimization as the generalized tree alignment problem. the gtap is the steiner tree problem for the sequence edit distance. like many biologically interesting problems, the gtap is np-hard. typically the steiner tree is presented under the manhattan or the hamming distances.  experimentally, the accuracy of the gtap has been subjected to evaluation. results show that phylogenies selected using the gtap from unaligned sequences are competitive with the best methods and algorithms available. here, we implement and explore experimentally existing and new local search heuristics for the gtap using simulated and real data.  the methods presented here improve by more than three orders of magnitude in execution time the best local search heuristics existing to date when applied to real data. tree alignmenttree searchphylogenysequence alignmentdirect optimization 