BACKGROUND
phylogenetic trees are used in biology to represent the evolutionary history of a set x of species   <cit> .
they are trees whose leaves are bijectively labeled by x and whose internal vertices represent the ancestors of the species set; they can be rooted or unrooted.
since in a rooted tree edges have a direction, the concepts of indegree and outdegree of a vertex are well defined.
binary rooted  trees are rooted  trees whose internal vertices have outdegree  <dig>  nonbinary rooted  trees have no restriction on the outdegree of inner vertices.
biological events in which a species derives its genes from different ancestors, such as hybridization, recombination and horizontal gene transfer events, cannot be modelled by a tree.
to be able to represent such events, a generalization of trees is considered which allows vertices with indegree two or higher, known as reticulations.
this model, which is called a rooted phylogenetic network, is of growing importance to biologists  <cit> .
for detailed background information we refer the reader to  <cit> .
although phylogenetic networks are more general than phylogenetic trees, trees are still often the basic building blocks from which phylogenetic networks are constructed.
specifically, there are many techniques available for constructing gene trees.
however, when more genes are analyzed, topological conflicts between individual gene phylogenies can arise for methodological or biological reasons .
this has led computational biologists to try and quantify the amount of reticulation that is needed to simultaneously explain two trees.
to state this problem more formally, we have that a phylogenetic tree t on  is a refinement of a phylogenetic tree t′ on the same set  if t can be obtained from t′ by deleting edges and identifying their incident vertices.
then, we say that a phylogenetic network n on displays a phylogenetic tree t on  if t can be obtained from a subgraph of n by contracting edges.
informally, this means that  t can be obtained from n by, for each reticulation vertex of n, “switching off” all but one of its incoming edges and then suppressing all indegree- <dig> outdegree- <dig> vertices .
given two rooted phylogenetic trees t <dig> and t <dig> on , the problem then becomes to determine the minimum number of reticulation events contained in a phylogenetic network n on  displaying both trees .
the value we are minimizing is often called the hybridization number and instead of the term phylogenetic network, the term hybridization network is often used.
it is known that the problem of computing hybridization numbers is both np-hard and apx-hard  <cit> , but it is not known whether it is in apx .
until recently, most research on the hybridization number of two phylogenetic trees had focused on the question of how to exactly compute this value using fixed parameter tractable  algorithms, where the parameter in question is the hybridization number r of the two trees.
for an introduction to fpt we refer to  <cit> .
for binary trees, algorithmic progress has been considerable in this area, with various authors reporting increasingly sophisticated fpt algorithms  <cit> .
the fastest algorithms currently implemented are the algorithm available inside the package dendroscope <cit> , based on  <cit> , and the sequence of progressively faster algorithms in the hybridnet family  <cit> .
the fastest theoretical fpt algorithm has running time o  <cit> , where n is the number of taxa in the trees.
even though in practice it rarely happens that trees are binary, the nonbinary variant of the problem has been less studied.
the nonbinary version is also fpt  <cit>  and a  algorithm has recently been implemented in dendroscope <cit> .
such  algorithms do, however, have their limits.
the running time still grows exponentially in r, albeit usually at a slower rate than algorithms that have a running time of the form n
f, where f is some function of r. in practice this means that existing algorithms can only handle instances of binary trees when r is at most 40- <dig> and instances of nonbinary trees when r is at most 5- <dig>
these limitations are problematic.
due to ongoing advances in dna sequencing, more and more species and strains are being sequenced.
consequently, biologists use trees with more and more taxa and software that can handle large trees is required.
for such large and/or difficult trees one can try to generate heuristic or approximate solutions, but how far are such solutions from optimality?
in  <cit>  we showed that the news is worrying.
indeed, we showed that polynomial-time constant-ratio approximation algorithms exist if and only if such algorithms exist for the problem directed feedback vertex set .
however, dfvs is a well-studied problem in combinatorial optimization and to this day it is unknown if it permits such an algorithm.
pending a major breakthrough in computer science, it therefore seems difficult to build polynomial-time algorithms which approximate hybridization number well.
on the positive side, we showed that in polynomial time an algorithm with approximation ratio o is possible.
however, this algorithm is purely of theoretical interest and is not useful in practice.
new algorithms: cyclekiller and nonbinarycyclekiller
in this article we extend the theoretical work of  <cit>  slightly and give it a practical twist to yield a fast approximation algorithm which we have made publicly available as the program cyclekiller.
furthermore, we give an implementation of the algorithm presented in  <cit> , available as nonbinarycyclekiller.
the worst-case running time of these approximation algorithms is exponential.
however, as we demonstrate with experiments, the running time of our algorithms is in practice extremely fast.
for large and/or massively discordant binary trees, cyclekiller is typically orders of magnitude faster than the hybridnet algorithms and the algorithm in dendroscope.
the performance gap between nonbinarycyclekiller and its exact counterparts is less pronounced, but still significant, especially in its fastest mode of operation.
of course, exact algorithms attempt to compute optimum solutions, whereas our algorithms only give approximate solutions.
nevertheless, our experiments show that when cyclekiller and nonbinarycyclekiller are run in their most accurate mode of operation, an approximation ratio very close to  <dig> is not unusual, suggesting that the algorithms often produce solutions close to optimality and well within the worst-case approximation guarantee.
the idea behind the binary and nonbinary algorithm is similar.
specifically, we describe an algorithm with approximation ratio d for the hybridization number problem on two binary trees and an algorithm with approximation ratio d for the hybridization number problem on two nonbinary trees by combining a c-approximation for the problem maf  with a d-approximation for the problem dfvs.
both these problems are np-hard so polynomial-time algorithms attaining c= <dig> or d= <dig> are not realistic.
nevertheless, there exist extremely fast fpt algorithms for solving maf on binary trees exactly , the fastest is rspr by whidden, beiko and zeh  <cit>  although the maf algorithm inside  <cit>  is also competitive.
moreover, we observe that the type of dfvs instances that arise in practice can easily be solved using integer linear programming  , so d= <dig> is also often possible.
combining these two exact approaches gives us, in the binary case, an exponential-time approximation algorithm with worst-case approximation ratio  <dig> that for large instances still runs extremely quickly; this is the 2-approx option of cyclekiller.
in practice, we have observed that the upper bound of  <dig> is often pessimistic, with much better approximation ratios observed in experiments .
we find that this algorithm already allows us to cope with much bigger trees than the hybridnet algorithms or the algorithm in dendroscope.
nevertheless, for truly massive trees it is often not feasible to have c= <dig>  fortunately there exist linear-time algorithms which achieve c= <dig>  <cit> .
this, coupled with the fact that  it remains feasible to use an exact  solver for dfvs, means that in practice we achieve a 4-approximation for gigantic binary trees; this is the 4-approx option of cyclekiller.
again, the ratio of  <dig> is a worst-case bound and we suspect that in practice we are doing much better than  <dig>  however, this cannot be experimentally verified due to the lack of good lower bounds for such massive instances.
in any case, the main advantage of this option is that it can, without too much effort, cope with trees with hundreds or thousands of taxa and hybridization number of a similar order of magnitude.
an implementation of cyclekiller and accompanying documentation can be downloaded from http://skelk.sdf-eu.org/cyclekiller.
networks created by the algorithm can be viewed in dendroscope.
for the nonbinary case, there also exist exact and approximation algorithms for maf  <cit> .
in case when one of the input trees is binary we can still use the exact  and approximate  algorithms given in  <cit>   to obtain respectively a 4-approximation and a 6-approximation of the hybridization number problem for nonbinary trees.
when both input trees are nonbinary, then we must use the somewhat less optimized exact  and approximate  algorithms described in  <cit> .
we then obtain 4- and 7-approximations .
to measure the approximation ratios attained by nonbinarycyclekiller in practice we have also implemented and made publicly available the exact nonbinary algorithm terminusest, based on the theoretical results in  <cit> .
terminusest will be of independent interest because it is currently the fastest exact nonbinary solver available.
cyclekiller, nonbinarycyclekiller and terminusest can be downloaded respectively from http://skelk.sdf-eu.org/cyclekiller <cit> , from http://homepages.cwi.nl/~iersel/cyclekiller <cit> , and from http://skelk.sdf-eu.org/terminusest <cit> .
theoretical and practical significance
we have described, implemented and made publicly available two algorithms with two desirable qualities: they terminate quickly even for massive instances of hybridization number and give a non-trivial guarantee of proximity to optimality.
these are the first algorithms with such properties.
both algorithms are based on a non-trivial marriage of maf and dfvs solvers , meaning that further advances in solving maf and dfvs will directly lead to improvements in cyclekiller and nonbinarycyclekiller.
this article also improves the theoretical work given in  <cit> , which also proposed using dfvs but beginning from a trivial agreement forest  known as a chain forest.
here we use a smarter starting point: an  maf, and it is this insight which makes a 2-approximation  possible when using an exact dfvs solver.
other articles have also had the idea of cycle-breaking in afs: the advanced fpt algorithm of whidden et al.
<cit>  – which has not been implemented – and the algorithms in the aforementioned hybridnet family.
however, both algorithms start the cycle-breaking from many starting points.
in contrast, our algorithm requires only a single starting point, i.e.
a single  solution to maf.
here, we only present the theory behind the binary algorithm.
the nonbinary case is more involved and we refer the reader to  <cit>  in which we introduce it.
note that our results for the binary case do not follow from the results for the nonbinary case in  <cit>  because here we obtain a better constant in the approximation ratio.
after a presentation of the binary algorithm in section “?
?”, we will show the results of some experiments with binary trees in section “practical experiments with binary trees” and nonbinary trees in section “practical experiments with nonbinary trees”.
finally, in section “practical experiments on biologically relevant trees” we demonstrate that both terminusest and nonbinarycyclekiller are easily capable of generating optimal  solutions on a real biological dataset originally obtained from the greenphyldb database.
technical note
at the time the experiments on binary trees were conducted  hybridnet was the fastest algorithm available in its family.
it has recently been superceded by the faster ultranet <cit> .
we believe, however, that it is neither necessary nor desirable to re-run the binary experiments, for the following reasons.
in the same period the solver rspr has also increased dramatically in speed , leading to a corresponding speed-up in cyclekiller.
in fact, both rspr and the algorithms in the hybridnet family are constantly in flux and are always being improved, so any experimental setup is prone to age extremely quickly.
however, the conclusions that we can derive from these experiments are unlikely to change much over time.
given that the algorithms in the hybridnet family  implicitly have to explore exponentially many optimal and sub-optimal solutions to the maf problem, the running time of maf solvers  is likely for the foreseeable future to remain much better than the running time of solvers for hybridization number.
the central message is stable: approximating hybridization number by splitting it into maf and dfvs instances yields extremely competitive approximation ratios for instances that exact hybridization number solvers will probably never be able to cope with.
methods
preliminaries
let  be a finite set .
a rooted phylogenetic-tree is a rooted tree with no vertices with indegree  <dig> and outdegree  <dig>  a root with indegree  <dig> and outdegree at least  <dig>  and leaves bijectively labelled by the elements of .
we identify each leaf with its label and use l to refer to the leaf set  of t. a rooted phylogenetic -tree is called binary if each nonleaf vertex has outdegree two.
we henceforth call a rooted, binary phylogenetic -tree a tree for short.
for a tree t and a set x′⊂x, we use the notation t to denote the minimal subtree of t that contains all elements of x′ and t|x′ denotes the result of suppressing all indegree- <dig> outdegree- <dig> vertices in t.
the following definitions apply only to binary trees.
definitions for nonbinary trees are analogous but slightly more technical  <cit> .
we define a forest as a set of trees.
each element of a forest is called a component.
let t be a tree and  a forest.
we say that  is a forest fort if, for all f∈f, t|l is isomorphic to f and the trees {t),f∈f} are vertex-disjoint subtrees of t whose leaf-set union equals l. if t <dig> and t <dig> are two trees, then a forest  is an agreement forest of t <dig> and t <dig> if it is a forest for t <dig> and t <dig>  the number of components of  is denoted |f|.
we define cleaning up a directed graph as repeatedly suppressing indegree- <dig> outdegree- <dig> vertices, removing indegree- <dig> outdegree- <dig> vertices and removing unlabelled outdegree- <dig> vertices until no such operation is possible.
observe that, if  is a forest for t,  can be obtained from t by removing |f|− <dig> edges and cleaning up.
from now on we consider t <dig> t <dig> as trees on the same taxon set.
problem: maximum agreement forest instance: two rooted, binary phylogenetic trees t <dig> and t <dig> solution: an agreement forest  of t <dig> and t <dig> objective: minimize |f|− <dig>
the directed graph ig, called the inheritance graph, is the directed graph whose vertices are the components of  and which has an edge  precisely if either
● there is a directed path in t <dig> from the root of t1) to the root of t1) or;
● there is a directed path in t <dig> from the root of t2) to the root of t2).
an agreement forest  of t <dig> and t <dig> is called an acyclic agreement forest if the graph ig is acyclic.
a maximum acyclic agreement forest  of t <dig> and t <dig> is an acyclic agreement forest of t <dig> and t <dig> with a minimum number of components.
problem: maximum acyclic agreement forest instance: two rooted, binary phylogenetic trees t <dig> and t <dig>  solution: an acyclic agreement forest  of t <dig> and t <dig> objective: minimize |f|− <dig>
we use maf  and maaf  to denote the optimal solution value of the problem maf and maaf respectively, for an instance t <dig> t <dig>
a rooted phylogenetic network on  is a directed acyclic graph with no vertices with indegree  <dig> and outdegree  <dig> and leaves bijectively labelled by the elements of .
rooted phylogenetic networks, which are sometimes also called hybridization networks, will henceforth be called networks for short in this paper.
a tree t on  is displayed by a network n if t can be obtained from a subtree of n by contracting edges.
a reticulation is a vertex v with δ−≥ <dig>  denoting the indegree of v).
the reticulation number  of a network n with root ρ is given by
r=∑v≠ρ−1).
it was shown that the optimum to maaf is equal to the optimum of the following problem  <cit> .
problem: minimumhybridizationinstance: two rooted binary phylogenetic trees t <dig> and t <dig>
solution: a rooted phylogenetic network n that displays t <dig> and t <dig> objective: minimize r.
moreover, it was shown that, for two trees t <dig> t <dig>  any acyclic agreement forest for t <dig> and t <dig> with k+ <dig> components can be turned into a phylogenetic network that displays t <dig> and t <dig> and has reticulation number k, and vice versa.
thus, any approximation for maaf gives an approximation for minimumhybridization.
finally, a feedback vertex set of a directed graph is a subset of the vertices that contains at least one vertex of each directed cycle.
equivalently, a subset of the vertices of a directed graph is a feedback vertex set if removing these vertices from the graph makes it acyclic.
problem: directed feedback vertex set instance: a directed graph d. goal: find a feedback vertex set of d of minimum size.
we note that the definition of minimumhybridization easily generalises to nonbinary trees, since the definition of display allows the image of each input tree in the network to be more “resolved” than the original tree.
however, the definitions of  agreement forests are different in the nonbinary case  <cit> .
the algorithm for binary trees
we show how maaf can be approximated by combining algorithms for maf and dfvs.
in particular, we will prove the following theorem.
theorem 1
if there exists a c-approximation for maf and a d-approx-imation for dfvs, then there exists a d-approximation for maaf .
note that this theorem does not follow from theorem  <dig>  of  <cit> , since there the approximation ratio for maaf is a d-approximation.
to prove the theorem, suppose there exists a c-approximation for maf.
let t <dig> and t <dig> be two trees and let m be an agreement forest returned by the algorithm.
then,
|m|−1≤c·maf≤c·maaf.
an m-splitting is an acyclic agreement forest that can be obtained from m by removing edges and cleaning up.
lemma 2
let t <dig> and t <dig> be two trees and m an agreement forest of t <dig> and t <dig>  then, there exists an m-splitting of size at most maaf+|m|.
proof
consider a maximum acyclic agreement forest f of t <dig> and t <dig>  for i∈{ <dig> }, f can be obtained from t
i
by removing a set of edges, say efi, and cleaning up.
moreover, also m can be obtained from t
i
by removing a set of edges, say emi, and cleaning up.
now consider the forest s obtained from t <dig> by removing em1∪ef <dig> and cleaning up.
then,
● s is an agreement forest of t <dig> and t <dig> because it can be obtained from t <dig> by removing edges em2∪ef <dig> and cleaning up;
● s is acyclic because it can be obtained by removing edges from f, which is acyclic, and cleaning up;
● s can be obtained from m by removing edges and cleaning up.
hence, s is an m-splitting.
furthermore, |s|≤|ef1|+|em1|+ <dig>  the lemma follows since |ef1|=maaf and |m|=|em1|+ <dig>
let optsplittingt <dig> t <dig> denote the size of a minimum-size m-splitting.
combining lemma  <dig> and eq.
<dig>  we obtain
optsplittingt <dig> t2−1≤maaf
we will now show how to find an approximation for the problem of finding an optimal m-splitting.
we do so by reducing the problem to dfvs.
we construct an input graph d for dfvs  as follows.
for every vertex of m that has outdegree  <dig> , we create a vertex in d. there is an edge in d from a vertex u to a vertex v precisely if in either t <dig> or t <dig>  there is a directed path from u to v. an example is in figure  <dig>  we claim the following.
lemma  <dig>
a subset v′ of the vertices of d is a feedback vertex set of d if and only if removing v′ from m makes it an acyclic agreement forest.
proof.
we show that d∖v′ has a directed cycle if and only if the inheritance graph of m∖v′ has a directed cycle.
to prove this, first suppose that there is a cycle v <dig> v <dig> …,v
k
=v <dig> in the inheritance graph of m∖v′.
the vertices in the inheritance graph of m∖v′ correspond to the roots of the components of m∖v′.
since these roots have outdegree  <dig> in m∖v′, they had outdegree  <dig> in m, and are thus vertices of d. so the vertices v <dig> v <dig> …,v
k
that form the cycle are vertices of d. since these vertices are in the inheritance graph of m∖v′, they can not be in v′ and so they are vertices of d∖v′.
the reachability relation between these vertices in d∖v′ is the same as in the inheritance graph of m∖v′.
so, the vertices v <dig> v <dig> …,v
k
form a cycle in d∖v′.
now suppose that there is a cycle w <dig> w <dig> …,w
k
=w <dig> in d∖v′.
each of the vertices w <dig> w <dig> …,w
k
is a vertex with outdegree- <dig> in m. some of them might be roots of components, while others are not.
however, observe that if there is a directed path from a vertex u to a vertex v in t <dig>  then there is also a directed path from the root of the component of m∖v′ that contains u to the root of the component of m∖v′ that contains v. hence, there is a directed cycle in the inheritance graph of m∖v′, formed by the roots of the components of m∖v′ that contain w <dig> w <dig> …,w
k
.
proof of theorem 1
suppose that there exists a d-approximation for dfvs.
let fvs be a feedback vertex set returned by this algorithm and let mfvs be a minimum feedback vertex set.
then, removing the vertices of mfvs from m gives an optimal m-splitting.
furthermore, optsplittingt <dig> t2=|m|+|mfvs|.
this is because for every vertex in a cycle c, its parent in m must participate in some cycle that contains elements of c. so if we start by removing the root of the component we are splitting and subsequently remove only those vertices whose parents have already been removed we see that we add at most one component per vertex.
in fact, because vertices of d all have out-degree  <dig> in m, we add exactly one component per vertex.
by removing the vertices of fvs from m, we obtain an acyclic agreement forest  such that
() |f|−1=|m|+|fvs|−1≤|m|+d·|mfvs|−1≤d=d−1)≤dmaaf,
where the last inequality follows from eq.
<dig>  thus,  is a d-approximation to maaf, which concludes the proof of theorem  <dig>
theorem  <dig> implies that a solution to the maaf problem for a given instance can be constructed by  finding a solution  to the maf problem for the same instance  constructing the extended inheritance graph d for   finding a solution v for the dfvs problem on the graph d and  modifying  accordantly to v.
RESULTS
practical experiments with binary trees
to assess the performance of cyclekiller, a simulation study was undertaken.
we generated  <dig> synthetic datasets, an easy, a medium and a hard one, containing respectively  <dig>   <dig> and  <dig> pairs of rooted binary phylogenetic trees.
the easy data set was created by varying two parameters, namely the number of taxa n and the number of rspr-moves k used to obtain the second tree from the first .
the  <dig> pairs of rooted binary phylogenetic trees were created by varying n in { <dig> , <dig> } and k in { <dig> ,...,25}, and then creating  <dig> different instances per each combination of parameters.
each pair  of rooted binary phylogenetic trees for a given set of parameters n and k is created as follows: the first tree t <dig> on x={x <dig> …,xn} is generated by first creating a set of n leaf vertices bijectively labeled by the set .
then, two vertices u and v, both with indegree  <dig>  are randomly picked and a new vertex w, along with two new edges  and , is created.
this is done until only one vertex with no ancestor, the root, is present.
the second tree t <dig> is obtained from t <dig> by applying k rspr-moves.
the medium and the hard data sets were generated in the same way as the easy one, but for different choices of the parameters: n in { <dig> , <dig> } and k in { <dig> , <dig> } for the medium one and n in { <dig> , <dig> } and k in { <dig> , <dig> } for the hard one.
the exact hybridization number has been computed by hybridnet <cit> , available from http://www.cs.cityu.edu.hk/~lwang/software/hn/treecomp.html or with dendroscope <cit> , available from http://dendroscope.org/.
we will refer to these algorithms as the exact algorithms.
each instance has been run on a single core of an intel xeon e <dig> processor.
each run that took more than one hour was aborted.
for each instance, we ran our program with the option 2-approx, and, in case the latter did not finish within one hour, we ran it again, this time using the option 4-approx, always with a one-hour limit .
we used the program rspr v <dig>   <cit>  to solve or approximate maf and glpk v <dig>   to solve the following simple polynomial-size ilp formulation of dfvs:
min∑v∈vxvs.t.0≤ℓv≤|v|−1for allv∈vℓv≥ℓu+1−|v|xu−|v|xvfor alle=∈eℓv∈ℤfor allv∈vxv∈{ <dig> }for allv∈v
given a directed graph d=, the binary variables x
v
model whether a vertex is in the feedback vertex set, and the integer variables ℓ
v
model the positions of the surviving vertices in the induced topological order.
the edge constraints enforce the topological order.
note that an edge constraint is essentially eliminated if one or both endpoints of the edge are in the feedback vertex set.
for all instances of the easy data set, cyclekiller finished with the 2-approx option within the one hour limit, while for  <dig> instances the exact algorithms were unable to compute the hybridization number.
note that, even for “easy” instances, computing the exact hybridization number can take a very long time.
to give the reader an idea, for  <dig> runs of the easy data, dendroscope and hybridnet did not complete within  <dig> days.
table  <dig> shows a summary of the results.
it can be seen that cyclekiller was much faster than the exact algorithms.
moreover, for  <dig> % of the instances for which an exact algorithm could find a solution, cyclekiller also found an optimal solution.
while the theoretical worst-case approximation ratio of the 2-approx option of cyclekiller is  <dig>  in our experiments it performed very close to a 1-approximation.
cyclekiller
the third column indicates for how many instances at least one exact algorithm finished within one hour.
the fifth column indicates for how many instances the 2-approx option of cyclekiller finished within one hour.
for the remaining instances, the 4-approx option finished within one hour, as can be seen from the seventh column.
the average running time for the 2-approx and the 4-approx in seconds are reported respectively in the sixth and eighth column.
the average approximation ratio  is taken over all instances for which at least one exact method finished.
the last column indicates the percentage of those instances for which cyclekiller found an optimal solution.
for the medium data set, cyclekiller finished with the 2-approx option for  <dig> instances, and for the remaining ones with the 4-approx option.
the exact algorithms could compute the hybridization number for only  <dig> instances .
for  <dig> % of these instances, cyclekiller also found an optimal solution, but with a much better running time.
regarding the hard data set,  <dig> runs were completed with the 2- approx option and for the remaining ones we were able to use the 4-approx option within the given time constraint.
unfortunately, the exact algorithms were unable to compute the hybridization number for any tree-pair of this data set and hence we could not compute the average approximation ratios.
over all our experiments, the maximum hybridization number that the exact algorithms could handle was 25a.
in contrast, the 2-approx option of cyclekiller could be used for instances for which the size of a maf was up to  <dig>  and thus for instances for which the hybridization number was at least  <dig>
to find the limits of the 4-approx option of cyclekiller, we also tested it on randomly generated trees.
on a normal laptop, it could construct networks with up to  <dig>  leaves and up to  <dig>  reticulations within  <dig> minutes.
since the number of reticulations found is at most four times the optimal hybridization number, this implies that the 4-approx option of cyclekiller can handle hybridization numbers up to at least  <dig> .
these randomly generated trees are, however, biologically meaningless and, therefore, we conducted the extensive experiment described above on trees generated by rspr moves.
finally we note that over all experiments the worst approximation ratio we encountered was  <dig> .
practical experiments with nonbinary trees
to run the simulations with nonbinarycyclekiller, we used a subset of the trees from the easy set of binary experiments.
we then applied random edge contractions in order to obtain nonbinary trees.
hence, we have the same two parameters as before, namely the number of taxa n∈{ <dig> ,100} and the number of rspr-moves k∈{ <dig> , <dig> }, and an additional parameter ρ∈{ <dig> ,75} which measures the percentage of the edges of an original binary tree that were contracted in order to obtain a nonbinary tree.
we could only use smaller values of n and k from the easy set of experiments because exact solvers for nonbinary maf  and exact solvers for nonbinary minimumhybridization  are slower than their binary counterparts.
we performed two runs of experimentsb.
one run with instances consisting of one binary and one nonbinary tree, and one run with instances consisting of two nonbinary trees.
for the experiments with one binary and one nonbinary tree, we were still able to use the rspr algorithm  <cit> , which has a better running time and approximation ratio compared to the available algorithm for two nonbinary trees.
when rspr is used in exact mode, nonbinarycyclekiller yields a theoretical worst-case approximation ratio of  <dig>  when rspr is used in its 3-approximation mode, nonbinarycyclekiller yields a theoretical worst-case approximation ratio of  <dig> .
the results of this run are summarized in table  <dig>
terminusest
nonbinarycyclekiller
opt
r
r
we list the average hybridization number found ), the average running time in seconds  and where applicable the average approximation ratio  for the three algorithms.
for the experiments with two nonbinary trees, the rspr software can no longer be used, and instead we used the exact and 4-approximate maf algorithm described in  <cit> .
this makes nonbinarycyclekiller behave as a 4-approximation and 7-approximation respectively .
note that the exact algorithm  <cit>  is considerably slower than rspr, meaning that in practice nonbinarycyclekiller struggles with two nonbinary trees more than when at most one of the trees is nonbinary.
the results for this run are summarized in table  <dig>
terminusest
nonbinarycyclekiller
opt
r
r
the layout of the table is the same as that of table  <dig>
the exact hybridization number in both runs was computed by terminusest <cit> .
each instance that took longer than  <dig> minutes to compute was aborted and the running time was set to  <dig> seconds.
the averages of the running-times are taken over all instances, with running-time taken to be  <dig> if the program timed out for that instance.
note that we did not compare the performance of nonbinarycyclekiller to dendroscope because terminusest has better running times than the exact nonbinary minimumhybridization solver inside dendroscope .
to enable a clearer analysis we divided the trees into representative “simple” and “tricky” ones based on two parameters, n and k. parameter values for the simple set were n∈{ <dig> }, k∈{ <dig> ,15} and for the tricky set n∈{ <dig> }, k= <dig>  in addition we varied the percentage of contracted edges .
in table  <dig> we show running times and solution quality of our algorithm when one of the input trees is binary.
for the simple set of instances  we see that the more accurate version of our algorithm, the 4-approximation, had a better running time than the exact algorithm, and at the same time had an average approximation ratio very close to  <dig>  far more interesting is to see what happens with tricky instances.
as predicted, the running time of the exact algorithm is much higher for tricky instances due to the higher hybridization numbers.
on the other hand, the running time of the 4-approximation does not rise significantly at all, whilst still attaining an approximation ratio again very close to  <dig>  another thing to note is that the percentage of contraction only seems to affect the running time of the exact algorithm.
the practical worst-case approximation ratio observed in these experiments was  <dig>  for the 4-approximation and  <dig> for the 6-approximation.
it is worth noting that, for the 4-approximation, the running time for the 75%-contraction trees is considerably lower than the one for the 50%-contraction trees.
this is due to the fact that a high contraction in both trees causes the hybridization number of the instance to drop, and a lower hybridization number leads to a better running time.
also note that the exact solver terminusest seems more able to cope with the tricky 25%-contraction instances than the tricky 50%-contraction instances.
this is probably because, although low contraction rates yield a higher hybridization number, the trees remain “relatively binary” and this can induce more efficient branching in the underlying fpt algorithm  <cit> .
it is plausible that with 50%-contraction the instances suffer from the disadvantage of relatively high hybridization number without the branching advantages associated with  binary trees.
to find the limits of the 7-approx option of nbck, we also tested it on huge, biologically meaningless, randomly generated trees.
below some results:
●  <dig> leaves, 25%-contraction, on average  <dig> reticulations in  <dig> sec.
●  <dig> leaves, 50%-contraction, on average  <dig> reticulations in  <dig> sec.
●  <dig> leaves, 75%-contraction, on average  <dig> reticulations, in  <dig> sec.
computation times of this last run of experiments do not include the network construction.
practical experiments on biologically relevant trees
finally, we tested our methods on phylogenetic trees obtained from greenphyldb  <cit>  – version  <dig>  a database containing twenty-two full genomes of members of the plantae kingdom, ranging from algae to angiosperms.
we were able to retrieve from the database the  <dig> rooted phylogenetic trees associated to the gene families contained in the database , along with the rooted phylogenetic tree describing the history of the twenty-two species contained in greenphyldb .
note that the species tree for these species is not completely resolved, i.e.
it is nonbinary.
among the gene trees,  <dig> contain less than  <dig> species and they were discarded.
of the remaining  <dig> trees, only  <dig> were directly usable for testing our methods.
indeed, because of gene duplication events arising in genomes, some species host several copies of the same gene, hence individual gene trees usually have several leaves labeled with identical species names.
unfortunately, our methods do not handle such multi-labeled gene trees .
we thus transformed the mul trees into trees containing single copies of labels, applying the tools described in  <cit>  to the forest f of  <dig> trees.
as in section  <dig>  of  <cit> , we obtained four data sets: f <dig>  f <dig>  f3p and f3s, respectively containing  <dig>   <dig>   <dig> and  <dig> trees.
note that only f3s contains nonbinary trees.
finally, for each single labeled tree g∈, we restricted the species tree s  to the leaves of g and we applied our methods to all so obtained pairs .
the results are presented in tables  <dig>   <dig>   <dig> and  <dig>  for all four datasets both terminusest and nonbinarycyclekiller ran extremely quickly, rarely taking more than a couple of seconds for each species-gene tree pair.
moreover, the clear conclusion with this dataset is that, although the species-gene pairs are often incompatible, there are rarely many cycles to kill and optimum solutions to the hybridization number problem are generally extremely close to optimal solutions to maf.
f
1
originally obtained from greenphyldb database
common taxa is the number of taxa after restricting the gene tree and the species tree to common taxa.
opt is the exact hybridization number, as computed by terminusest.
ratio 4-approx  is the ratio of the solution obtained by nonbinarycyclekiller  to the solution obtained by terminusest.
gap  is the absolute gap between the optimum maf solution  and the exact hybridization number, as computed by terminusest.
gap  is the absolute gap between the optimum maf solution and the reticulation number of the solution generated by nonbinarycyclekiller running in its 4-approx mode.
time t-est is the running time  of terminusest, and time 4-approx is the running time  of nonbinarycyclekiller running in its 4-approx mode.
in  <dig> instances terminusest returned the same size solution as rspr, in  <dig> cases terminusest returned the same size solution as nonbinarycyclekiller , and in  <dig> cases nonbinarycyclekiller  returned the same size solution as rspr.
f
2
originally obtained from greenphyldb database
in  <dig> instances terminusest returned the same size solution as rspr, in  <dig> cases terminusest returned the same size solution as nonbinarycyclekiller , and in  <dig> cases nonbinarycyclekiller  returned the same size solution as rspr.
in  <dig> instances terminusest returned the same size solution as rspr, in  <dig> cases terminusest returned the same size solution as nonbinarycyclekiller , and in  <dig> cases nonbinarycyclekiller  returned the same size solution as rspr.
in  <dig> instances terminusest returned the same size solution as rspr, in  <dig> cases terminusest returned the same size solution as nonbinarycyclekiller , and in  <dig> cases nonbinarycyclekiller  returned the same size solution as maf.
in this dataset the gene trees were also nonbinary, meaning that nonbinarycyclekiller had to use the maf algorithm described in  <cit>  instead of rspr.
CONCLUSIONS
our experiments with binary trees show that cyclekiller is much faster than available exact methods once the input trees become sufficiently large and/or discordant.
in over 96% of the cases cyclekiller finds the optimal solution and in the remaining cases it finds a solution very close to the optimum.
we have shown that the most accurate mode of the program produces solutions that are at most a factor  <dig> from the optimum.
in practice, the average-case approximation ratio that we observed was  <dig> .
the fastest mode of the algorithm can be used on trees with thousands of leaves and provably constructs networks that are at most a factor of  <dig> from the optimum.
our experiments with nonbinary trees highlight once again that the cycle-breaking technique described in this article is intrinsically linked to the current state-of-the-art in maf algorithms.
terminusest is faster than the most accurate mode of nonbinarycyclekiller when both trees are nonbinary due to the fact that maf solvers for two nonbinary trees have not yet been optimized to the same extent as their binary counterparts.
in fact, terminusest is the best avaible exact method for nonbinary trees and can handle instances for which the optimum is up to 15- <dig>  for other instances, nonbinarycyclekiller in its fastest mode is much faster than terminusest and produces solutions that are at most a factor  <dig> from the optimum .
finally, for instances with one binary and one nonbinary tree, the most accurate mode of nonbinarycyclekiller is again much faster than terminusest and produces solutions that are at most a factor  <dig>  from the optimum .
endnotes
a in  <cit> , it has been shown that this number can go up to  <dig> when running dendroscope on a similar processor but allocating all cores for one instance, i.e.
exploiting the possibilities of parallel computation of this implementation.
b we note that nonbinarycyclekiller uses a row-generation ilp formulation - based on  <cit>  - to solve dfvs, rather than the polynomial-size formulation used by cyclekiller.
ilp is in neither case a bottleneck for the running time.
competing interests
the authors declare that they have no competing interests.
authors’ contributions
all the authors conceived the ideas, designed and conducted the experiments, and wrote and approved the paper.
authors’ information
leo van iersel was supported by a veni grant of the netherlands organisation for scientific research .
nela lekić was supported by a vrije competitie grant of of the netherlands organisation for scientific research .

 phylogenetic trees are used in biology to represent the evolutionary history of a set x of species. they are trees whose leaves are bijectively labeled by x and whose internal vertices represent the ancestors of the species set; they can be rooted or unrooted. since in a rooted tree edges have a direction, the concepts of indegree and outdegree of a vertex are well defined. binary rooted trees are rooted trees whose internal vertices have outdegree nonbinary rooted trees have no restriction on the outdegree of inner vertices. biological events in which a species derives its genes from different ancestors, such as hybridization, recombination and horizontal gene transfer events, cannot be modelled by a tree. to be able to represent such events, a generalization of trees is considered which allows vertices with indegree two or higher, known as reticulations. this model, which is called a rooted phylogenetic network, is of growing importance to biologists. for detailed background information we refer the reader to. although phylogenetic networks are more general than phylogenetic trees, trees are still often the basic building blocks from which phylogenetic networks are constructed. specifically, there are many techniques available for constructing gene trees. however, when more genes are analyzed, topological conflicts between individual gene phylogenies can arise for methodological or biological reasons. this has led computational biologists to try and quantify the amount of reticulation that is needed to simultaneously explain two trees. to state this problem more formally, we have that a phylogenetic tree t on is a refinement of a phylogenetic tree t′ on the same set if t can be obtained from t′ by deleting edges and identifying their incident vertices. then, we say that a phylogenetic network n on displays a phylogenetic tree t on if t can be obtained from a subgraph of n by contracting edges. informally, this means that t can be obtained from n by, for each reticulation vertex of n, “switching off” all but one of its incoming edges and then suppressing all indegree- outdegree- vertices. given two rooted phylogenetic trees t and t on, the problem then becomes to determine the minimum number of reticulation events contained in a phylogenetic network n on displaying both trees. the value we are minimizing is often called the hybridization number and instead of the term phylogenetic network, the term hybridization network is often used. it is known that the problem of computing hybridization numbers is both np-hard and apx-hard, but it is not known whether it is in apx. until recently, most research on the hybridization number of two phylogenetic trees had focused on the question of how to exactly compute this value using fixed parameter tractable algorithms, where the parameter in question is the hybridization number r of the two trees. for an introduction to fpt we refer to. for binary trees, algorithmic progress has been considerable in this area, with various authors reporting increasingly sophisticated fpt algorithms. the fastest algorithms currently implemented are the algorithm available inside the package dendroscope, based on, and the sequence of progressively faster algorithms in the hybridnet family. the fastest theoretical fpt algorithm has running time o, where n is the number of taxa in the trees. even though in practice it rarely happens that trees are binary, the nonbinary variant of the problem has been less studied. the nonbinary version is also fpt and a algorithm has recently been implemented in dendroscope. such algorithms do, however, have their limits. the running time still grows exponentially in r, albeit usually at a slower rate than algorithms that have a running time of the form n f, where f is some function of r. in practice this means that existing algorithms can only handle instances of binary trees when r is at most 40- and instances of nonbinary trees when r is at most 5- these limitations are problematic. due to ongoing advances in dna sequencing, more and more species and strains are being sequenced. consequently, biologists use trees with more and more taxa and software that can handle large trees is required. for such large and/or difficult trees one can try to generate heuristic or approximate solutions, but how far are such solutions from optimality? in we showed that the news is worrying. indeed, we showed that polynomial-time constant-ratio approximation algorithms exist if and only if such algorithms exist for the problem directed feedback vertex set. however, dfvs is a well-studied problem in combinatorial optimization and to this day it is unknown if it permits such an algorithm. pending a major breakthrough in computer science, it therefore seems difficult to build polynomial-time algorithms which approximate hybridization number well. on the positive side, we showed that in polynomial time an algorithm with approximation ratio o is possible. however, this algorithm is purely of theoretical interest and is not useful in practice. new algorithms: cyclekiller and nonbinarycyclekiller in this article we extend the theoretical work of slightly and give it a practical twist to yield a fast approximation algorithm which we have made publicly available as the program cyclekiller. furthermore, we give an implementation of the algorithm presented in, available as nonbinarycyclekiller. the worst-case running time of these approximation algorithms is exponential. however, as we demonstrate with experiments, the running time of our algorithms is in practice extremely fast. for large and/or massively discordant binary trees, cyclekiller is typically orders of magnitude faster than the hybridnet algorithms and the algorithm in dendroscope. the performance gap between nonbinarycyclekiller and its exact counterparts is less pronounced, but still significant, especially in its fastest mode of operation. of course, exact algorithms attempt to compute optimum solutions, whereas our algorithms only give approximate solutions. nevertheless, our experiments show that when cyclekiller and nonbinarycyclekiller are run in their most accurate mode of operation, an approximation ratio very close to is not unusual, suggesting that the algorithms often produce solutions close to optimality and well within the worst-case approximation guarantee. the idea behind the binary and nonbinary algorithm is similar. specifically, we describe an algorithm with approximation ratio d for the hybridization number problem on two binary trees and an algorithm with approximation ratio d for the hybridization number problem on two nonbinary trees by combining a c-approximation for the problem maf with a d-approximation for the problem dfvs. both these problems are np-hard so polynomial-time algorithms attaining c= or d= are not realistic. nevertheless, there exist extremely fast fpt algorithms for solving maf on binary trees exactly, the fastest is rspr by whidden, beiko and zeh although the maf algorithm inside is also competitive. moreover, we observe that the type of dfvs instances that arise in practice can easily be solved using integer linear programming, so d= is also often possible. combining these two exact approaches gives us, in the binary case, an exponential-time approximation algorithm with worst-case approximation ratio that for large instances still runs extremely quickly; this is the 2-approx option of cyclekiller. in practice, we have observed that the upper bound of is often pessimistic, with much better approximation ratios observed in experiments. we find that this algorithm already allows us to cope with much bigger trees than the hybridnet algorithms or the algorithm in dendroscope. nevertheless, for truly massive trees it is often not feasible to have c= fortunately there exist linear-time algorithms which achieve c=. this, coupled with the fact that it remains feasible to use an exact solver for dfvs, means that in practice we achieve a 4-approximation for gigantic binary trees; this is the 4-approx option of cyclekiller. again, the ratio of is a worst-case bound and we suspect that in practice we are doing much better than however, this cannot be experimentally verified due to the lack of good lower bounds for such massive instances. in any case, the main advantage of this option is that it can, without too much effort, cope with trees with hundreds or thousands of taxa and hybridization number of a similar order of magnitude. an implementation of cyclekiller and accompanying documentation can be downloaded from http://skelk.sdf-eu.org/cyclekiller. networks created by the algorithm can be viewed in dendroscope. for the nonbinary case, there also exist exact and approximation algorithms for maf. in case when one of the input trees is binary we can still use the exact and approximate algorithms given in to obtain respectively a 4-approximation and a 6-approximation of the hybridization number problem for nonbinary trees. when both input trees are nonbinary, then we must use the somewhat less optimized exact and approximate algorithms described in. we then obtain 4- and 7-approximations. to measure the approximation ratios attained by nonbinarycyclekiller in practice we have also implemented and made publicly available the exact nonbinary algorithm terminusest, based on the theoretical results in. terminusest will be of independent interest because it is currently the fastest exact nonbinary solver available. cyclekiller, nonbinarycyclekiller and terminusest can be downloaded respectively from http://skelk.sdf-eu.org/cyclekiller, from http://homepages.cwi.nl/~iersel/cyclekiller, and from http://skelk.sdf-eu.org/terminusest. theoretical and practical significance we have described, implemented and made publicly available two algorithms with two desirable qualities: they terminate quickly even for massive instances of hybridization number and give a non-trivial guarantee of proximity to optimality. these are the first algorithms with such properties. both algorithms are based on a non-trivial marriage of maf and dfvs solvers, meaning that further advances in solving maf and dfvs will directly lead to improvements in cyclekiller and nonbinarycyclekiller. this article also improves the theoretical work given in, which also proposed using dfvs but beginning from a trivial agreement forest known as a chain forest. here we use a smarter starting point: an maf, and it is this insight which makes a 2-approximation possible when using an exact dfvs solver. other articles have also had the idea of cycle-breaking in afs: the advanced fpt algorithm of whidden et al. – which has not been implemented – and the algorithms in the aforementioned hybridnet family. however, both algorithms start the cycle-breaking from many starting points. in contrast, our algorithm requires only a single starting point, i.e. a single solution to maf. here, we only present the theory behind the binary algorithm. the nonbinary case is more involved and we refer the reader to in which we introduce it. note that our results for the binary case do not follow from the results for the nonbinary case in because here we obtain a better constant in the approximation ratio. after a presentation of the binary algorithm in section “??”, we will show the results of some experiments with binary trees in section “practical experiments with binary trees” and nonbinary trees in section “practical experiments with nonbinary trees”. finally, in section “practical experiments on biologically relevant trees” we demonstrate that both terminusest and nonbinarycyclekiller are easily capable of generating optimal solutions on a real biological dataset originally obtained from the greenphyldb database. technical note at the time the experiments on binary trees were conducted hybridnet was the fastest algorithm available in its family. it has recently been superceded by the faster ultranet. we believe, however, that it is neither necessary nor desirable to re-run the binary experiments, for the following reasons. in the same period the solver rspr has also increased dramatically in speed, leading to a corresponding speed-up in cyclekiller. in fact, both rspr and the algorithms in the hybridnet family are constantly in flux and are always being improved, so any experimental setup is prone to age extremely quickly. however, the conclusions that we can derive from these experiments are unlikely to change much over time. given that the algorithms in the hybridnet family implicitly have to explore exponentially many optimal and sub-optimal solutions to the maf problem, the running time of maf solvers is likely for the foreseeable future to remain much better than the running time of solvers for hybridization number. the central message is stable: approximating hybridization number by splitting it into maf and dfvs instances yields extremely competitive approximation ratios for instances that exact hybridization number solvers will probably never be able to cope with. methods preliminaries let be a finite set. a rooted phylogenetic-tree is a rooted tree with no vertices with indegree and outdegree a root with indegree and outdegree at least and leaves bijectively labelled by the elements of. we identify each leaf with its label and use l to refer to the leaf set of t. a rooted phylogenetic -tree is called binary if each nonleaf vertex has outdegree two. we henceforth call a rooted, binary phylogenetic -tree a tree for short. for a tree t and a set x′⊂x, we use the notation t to denote the minimal subtree of t that contains all elements of x′ and t|x′ denotes the result of suppressing all indegree- outdegree- vertices in t. the following definitions apply only to binary trees. definitions for nonbinary trees are analogous but slightly more technical. we define a forest as a set of trees. each element of a forest is called a component. let t be a tree and a forest. we say that is a forest fort if, for all f∈f, t|l is isomorphic to f and the trees {t),f∈f} are vertex-disjoint subtrees of t whose leaf-set union equals l. if t and t are two trees, then a forest is an agreement forest of t and t if it is a forest for t and t the number of components of is denoted |f|. we define cleaning up a directed graph as repeatedly suppressing indegree- outdegree- vertices, removing indegree- outdegree- vertices and removing unlabelled outdegree- vertices until no such operation is possible. observe that, if is a forest for t, can be obtained from t by removing |f|− edges and cleaning up. from now on we consider t t as trees on the same taxon set. problem: maximum agreement forest instance: two rooted, binary phylogenetic trees t and t solution: an agreement forest of t and t objective: minimize |f|− the directed graph ig, called the inheritance graph, is the directed graph whose vertices are the components of and which has an edge precisely if either ● there is a directed path in t from the root of t1) to the root of t1) or; ● there is a directed path in t from the root of t2) to the root of t2). an agreement forest of t and t is called an acyclic agreement forest if the graph ig is acyclic. a maximum acyclic agreement forest of t and t is an acyclic agreement forest of t and t with a minimum number of components. problem: maximum acyclic agreement forest instance: two rooted, binary phylogenetic trees t and t solution: an acyclic agreement forest of t and t objective: minimize |f|− we use maf and maaf to denote the optimal solution value of the problem maf and maaf respectively, for an instance t t a rooted phylogenetic network on is a directed acyclic graph with no vertices with indegree and outdegree and leaves bijectively labelled by the elements of. rooted phylogenetic networks, which are sometimes also called hybridization networks, will henceforth be called networks for short in this paper. a tree t on is displayed by a network n if t can be obtained from a subtree of n by contracting edges. a reticulation is a vertex v with δ−≥ denoting the indegree of v). the reticulation number of a network n with root ρ is given by r=∑v≠ρ−1). it was shown that the optimum to maaf is equal to the optimum of the following problem. problem: minimumhybridizationinstance: two rooted binary phylogenetic trees t and t solution: a rooted phylogenetic network n that displays t and t objective: minimize r. moreover, it was shown that, for two trees t t any acyclic agreement forest for t and t with k+ components can be turned into a phylogenetic network that displays t and t and has reticulation number k, and vice versa. thus, any approximation for maaf gives an approximation for minimumhybridization. finally, a feedback vertex set of a directed graph is a subset of the vertices that contains at least one vertex of each directed cycle. equivalently, a subset of the vertices of a directed graph is a feedback vertex set if removing these vertices from the graph makes it acyclic. problem: directed feedback vertex set instance: a directed graph d. goal: find a feedback vertex set of d of minimum size. we note that the definition of minimumhybridization easily generalises to nonbinary trees, since the definition of display allows the image of each input tree in the network to be more “resolved” than the original tree. however, the definitions of agreement forests are different in the nonbinary case. the algorithm for binary trees we show how maaf can be approximated by combining algorithms for maf and dfvs. in particular, we will prove the following theorem. theorem 1 if there exists a c-approximation for maf and a d-approx-imation for dfvs, then there exists a d-approximation for maaf. note that this theorem does not follow from theorem of, since there the approximation ratio for maaf is a d-approximation. to prove the theorem, suppose there exists a c-approximation for maf. let t and t be two trees and let m be an agreement forest returned by the algorithm. then, |m|−1≤c·maf≤c·maaf. an m-splitting is an acyclic agreement forest that can be obtained from m by removing edges and cleaning up. lemma 2 let t and t be two trees and m an agreement forest of t and t then, there exists an m-splitting of size at most maaf+|m|. proof consider a maximum acyclic agreement forest f of t and t for i∈{ }, f can be obtained from t i by removing a set of edges, say efi, and cleaning up. moreover, also m can be obtained from t i by removing a set of edges, say emi, and cleaning up. now consider the forest s obtained from t by removing em1∪ef and cleaning up. then, ● s is an agreement forest of t and t because it can be obtained from t by removing edges em2∪ef and cleaning up; ● s is acyclic because it can be obtained by removing edges from f, which is acyclic, and cleaning up; ● s can be obtained from m by removing edges and cleaning up. hence, s is an m-splitting. furthermore, |s|≤|ef1|+|em1|+ the lemma follows since |ef1|=maaf and |m|=|em1|+ let optsplittingt t denote the size of a minimum-size m-splitting. combining lemma and eq. we obtain optsplittingt t2−1≤maaf we will now show how to find an approximation for the problem of finding an optimal m-splitting. we do so by reducing the problem to dfvs. we construct an input graph d for dfvs as follows. for every vertex of m that has outdegree, we create a vertex in d. there is an edge in d from a vertex u to a vertex v precisely if in either t or t there is a directed path from u to v. an example is in figure we claim the following. lemma a subset v′ of the vertices of d is a feedback vertex set of d if and only if removing v′ from m makes it an acyclic agreement forest. proof. we show that d∖v′ has a directed cycle if and only if the inheritance graph of m∖v′ has a directed cycle. to prove this, first suppose that there is a cycle v v …,v k =v in the inheritance graph of m∖v′. the vertices in the inheritance graph of m∖v′ correspond to the roots of the components of m∖v′. since these roots have outdegree in m∖v′, they had outdegree in m, and are thus vertices of d. so the vertices v v …,v k that form the cycle are vertices of d. since these vertices are in the inheritance graph of m∖v′, they can not be in v′ and so they are vertices of d∖v′. the reachability relation between these vertices in d∖v′ is the same as in the inheritance graph of m∖v′. so, the vertices v v …,v k form a cycle in d∖v′. now suppose that there is a cycle w w …,w k =w in d∖v′. each of the vertices w w …,w k is a vertex with outdegree- in m. some of them might be roots of components, while others are not. however, observe that if there is a directed path from a vertex u to a vertex v in t then there is also a directed path from the root of the component of m∖v′ that contains u to the root of the component of m∖v′ that contains v. hence, there is a directed cycle in the inheritance graph of m∖v′, formed by the roots of the components of m∖v′ that contain w w …,w k. proof of theorem 1 suppose that there exists a d-approximation for dfvs. let fvs be a feedback vertex set returned by this algorithm and let mfvs be a minimum feedback vertex set. then, removing the vertices of mfvs from m gives an optimal m-splitting. furthermore, optsplittingt t2=|m|+|mfvs|. this is because for every vertex in a cycle c, its parent in m must participate in some cycle that contains elements of c. so if we start by removing the root of the component we are splitting and subsequently remove only those vertices whose parents have already been removed we see that we add at most one component per vertex. in fact, because vertices of d all have out-degree in m, we add exactly one component per vertex. by removing the vertices of fvs from m, we obtain an acyclic agreement forest such that () |f|−1=|m|+|fvs|−1≤|m|+d·|mfvs|−1≤d=d−1)≤dmaaf, where the last inequality follows from eq. thus, is a d-approximation to maaf, which concludes the proof of theorem theorem implies that a solution to the maaf problem for a given instance can be constructed by finding a solution to the maf problem for the same instance constructing the extended inheritance graph d for finding a solution v for the dfvs problem on the graph d and modifying accordantly to v.  practical experiments with binary trees to assess the performance of cyclekiller, a simulation study was undertaken. we generated synthetic datasets, an easy, a medium and a hard one, containing respectively and pairs of rooted binary phylogenetic trees. the easy data set was created by varying two parameters, namely the number of taxa n and the number of rspr-moves k used to obtain the second tree from the first. the pairs of rooted binary phylogenetic trees were created by varying n in {, } and k in {,...,25}, and then creating different instances per each combination of parameters. each pair of rooted binary phylogenetic trees for a given set of parameters n and k is created as follows: the first tree t on x={x …,xn} is generated by first creating a set of n leaf vertices bijectively labeled by the set. then, two vertices u and v, both with indegree are randomly picked and a new vertex w, along with two new edges and, is created. this is done until only one vertex with no ancestor, the root, is present. the second tree t is obtained from t by applying k rspr-moves. the medium and the hard data sets were generated in the same way as the easy one, but for different choices of the parameters: n in {, } and k in {, } for the medium one and n in {, } and k in {, } for the hard one. the exact hybridization number has been computed by hybridnet, available from http://www.cs.cityu.edu.hk/~lwang/software/hn/treecomp.html or with dendroscope, available from http://dendroscope.org/. we will refer to these algorithms as the exact algorithms. each instance has been run on a single core of an intel xeon e processor. each run that took more than one hour was aborted. for each instance, we ran our program with the option 2-approx, and, in case the latter did not finish within one hour, we ran it again, this time using the option 4-approx, always with a one-hour limit. we used the program rspr v to solve or approximate maf and glpk v to solve the following simple polynomial-size ilp formulation of dfvs: min∑v∈vxvs.t.0≤ℓv≤|v|−1for allv∈vℓv≥ℓu+1−|v|xu−|v|xvfor alle=∈eℓv∈ℤfor allv∈vxv∈{ }for allv∈v given a directed graph d=, the binary variables x v model whether a vertex is in the feedback vertex set, and the integer variables ℓ v model the positions of the surviving vertices in the induced topological order. the edge constraints enforce the topological order. note that an edge constraint is essentially eliminated if one or both endpoints of the edge are in the feedback vertex set. for all instances of the easy data set, cyclekiller finished with the 2-approx option within the one hour limit, while for instances the exact algorithms were unable to compute the hybridization number. note that, even for “easy” instances, computing the exact hybridization number can take a very long time. to give the reader an idea, for runs of the easy data, dendroscope and hybridnet did not complete within days. table shows a summary of the results. it can be seen that cyclekiller was much faster than the exact algorithms. moreover, for % of the instances for which an exact algorithm could find a solution, cyclekiller also found an optimal solution. while the theoretical worst-case approximation ratio of the 2-approx option of cyclekiller is in our experiments it performed very close to a 1-approximation.        cyclekiller the third column indicates for how many instances at least one exact algorithm finished within one hour. the fifth column indicates for how many instances the 2-approx option of cyclekiller finished within one hour. for the remaining instances, the 4-approx option finished within one hour, as can be seen from the seventh column. the average running time for the 2-approx and the 4-approx in seconds are reported respectively in the sixth and eighth column. the average approximation ratio is taken over all instances for which at least one exact method finished. the last column indicates the percentage of those instances for which cyclekiller found an optimal solution. for the medium data set, cyclekiller finished with the 2-approx option for instances, and for the remaining ones with the 4-approx option. the exact algorithms could compute the hybridization number for only instances. for % of these instances, cyclekiller also found an optimal solution, but with a much better running time. regarding the hard data set, runs were completed with the 2- approx option and for the remaining ones we were able to use the 4-approx option within the given time constraint. unfortunately, the exact algorithms were unable to compute the hybridization number for any tree-pair of this data set and hence we could not compute the average approximation ratios. over all our experiments, the maximum hybridization number that the exact algorithms could handle was 25a. in contrast, the 2-approx option of cyclekiller could be used for instances for which the size of a maf was up to and thus for instances for which the hybridization number was at least to find the limits of the 4-approx option of cyclekiller, we also tested it on randomly generated trees. on a normal laptop, it could construct networks with up to leaves and up to reticulations within minutes. since the number of reticulations found is at most four times the optimal hybridization number, this implies that the 4-approx option of cyclekiller can handle hybridization numbers up to at least. these randomly generated trees are, however, biologically meaningless and, therefore, we conducted the extensive experiment described above on trees generated by rspr moves. finally we note that over all experiments the worst approximation ratio we encountered was. practical experiments with nonbinary trees to run the simulations with nonbinarycyclekiller, we used a subset of the trees from the easy set of binary experiments. we then applied random edge contractions in order to obtain nonbinary trees. hence, we have the same two parameters as before, namely the number of taxa n∈{,100} and the number of rspr-moves k∈{, }, and an additional parameter ρ∈{,75} which measures the percentage of the edges of an original binary tree that were contracted in order to obtain a nonbinary tree. we could only use smaller values of n and k from the easy set of experiments because exact solvers for nonbinary maf and exact solvers for nonbinary minimumhybridization are slower than their binary counterparts. we performed two runs of experimentsb. one run with instances consisting of one binary and one nonbinary tree, and one run with instances consisting of two nonbinary trees. for the experiments with one binary and one nonbinary tree, we were still able to use the rspr algorithm, which has a better running time and approximation ratio compared to the available algorithm for two nonbinary trees. when rspr is used in exact mode, nonbinarycyclekiller yields a theoretical worst-case approximation ratio of when rspr is used in its 3-approximation mode, nonbinarycyclekiller yields a theoretical worst-case approximation ratio of. the results of this run are summarized in table   terminusest nonbinarycyclekiller opt r r we list the average hybridization number found ), the average running time in seconds and where applicable the average approximation ratio for the three algorithms. for the experiments with two nonbinary trees, the rspr software can no longer be used, and instead we used the exact and 4-approximate maf algorithm described in. this makes nonbinarycyclekiller behave as a 4-approximation and 7-approximation respectively. note that the exact algorithm is considerably slower than rspr, meaning that in practice nonbinarycyclekiller struggles with two nonbinary trees more than when at most one of the trees is nonbinary. the results for this run are summarized in table   terminusest nonbinarycyclekiller opt r r the layout of the table is the same as that of table the exact hybridization number in both runs was computed by terminusest. each instance that took longer than minutes to compute was aborted and the running time was set to seconds. the averages of the running-times are taken over all instances, with running-time taken to be if the program timed out for that instance. note that we did not compare the performance of nonbinarycyclekiller to dendroscope because terminusest has better running times than the exact nonbinary minimumhybridization solver inside dendroscope. to enable a clearer analysis we divided the trees into representative “simple” and “tricky” ones based on two parameters, n and k. parameter values for the simple set were n∈{ }, k∈{,15} and for the tricky set n∈{ }, k= in addition we varied the percentage of contracted edges. in table we show running times and solution quality of our algorithm when one of the input trees is binary. for the simple set of instances we see that the more accurate version of our algorithm, the 4-approximation, had a better running time than the exact algorithm, and at the same time had an average approximation ratio very close to far more interesting is to see what happens with tricky instances. as predicted, the running time of the exact algorithm is much higher for tricky instances due to the higher hybridization numbers. on the other hand, the running time of the 4-approximation does not rise significantly at all, whilst still attaining an approximation ratio again very close to another thing to note is that the percentage of contraction only seems to affect the running time of the exact algorithm. the practical worst-case approximation ratio observed in these experiments was for the 4-approximation and for the 6-approximation. it is worth noting that, for the 4-approximation, the running time for the 75%-contraction trees is considerably lower than the one for the 50%-contraction trees. this is due to the fact that a high contraction in both trees causes the hybridization number of the instance to drop, and a lower hybridization number leads to a better running time. also note that the exact solver terminusest seems more able to cope with the tricky 25%-contraction instances than the tricky 50%-contraction instances. this is probably because, although low contraction rates yield a higher hybridization number, the trees remain “relatively binary” and this can induce more efficient branching in the underlying fpt algorithm. it is plausible that with 50%-contraction the instances suffer from the disadvantage of relatively high hybridization number without the branching advantages associated with binary trees. to find the limits of the 7-approx option of nbck, we also tested it on huge, biologically meaningless, randomly generated trees. below some results: ●  leaves, 25%-contraction, on average reticulations in sec. ●  leaves, 50%-contraction, on average reticulations in sec. ●  leaves, 75%-contraction, on average reticulations, in sec. computation times of this last run of experiments do not include the network construction. practical experiments on biologically relevant trees finally, we tested our methods on phylogenetic trees obtained from greenphyldb – version a database containing twenty-two full genomes of members of the plantae kingdom, ranging from algae to angiosperms. we were able to retrieve from the database the rooted phylogenetic trees associated to the gene families contained in the database, along with the rooted phylogenetic tree describing the history of the twenty-two species contained in greenphyldb. note that the species tree for these species is not completely resolved, i.e. it is nonbinary. among the gene trees, contain less than species and they were discarded. of the remaining trees, only were directly usable for testing our methods. indeed, because of gene duplication events arising in genomes, some species host several copies of the same gene, hence individual gene trees usually have several leaves labeled with identical species names. unfortunately, our methods do not handle such multi-labeled gene trees. we thus transformed the mul trees into trees containing single copies of labels, applying the tools described in to the forest f of trees. as in section of, we obtained four data sets: f f f3p and f3s, respectively containing and trees. note that only f3s contains nonbinary trees. finally, for each single labeled tree g∈, we restricted the species tree s to the leaves of g and we applied our methods to all so obtained pairs. the results are presented in tables and for all four datasets both terminusest and nonbinarycyclekiller ran extremely quickly, rarely taking more than a couple of seconds for each species-gene tree pair. moreover, the clear conclusion with this dataset is that, although the species-gene pairs are often incompatible, there are rarely many cycles to kill and optimum solutions to the hybridization number problem are generally extremely close to optimal solutions to maf. f 1 originally obtained from greenphyldb database common taxa is the number of taxa after restricting the gene tree and the species tree to common taxa. opt is the exact hybridization number, as computed by terminusest. ratio 4-approx is the ratio of the solution obtained by nonbinarycyclekiller to the solution obtained by terminusest. gap is the absolute gap between the optimum maf solution and the exact hybridization number, as computed by terminusest. gap is the absolute gap between the optimum maf solution and the reticulation number of the solution generated by nonbinarycyclekiller running in its 4-approx mode. time t-est is the running time of terminusest, and time 4-approx is the running time of nonbinarycyclekiller running in its 4-approx mode. in instances terminusest returned the same size solution as rspr, in cases terminusest returned the same size solution as nonbinarycyclekiller, and in cases nonbinarycyclekiller returned the same size solution as rspr. f 2 originally obtained from greenphyldb database in instances terminusest returned the same size solution as rspr, in cases terminusest returned the same size solution as nonbinarycyclekiller, and in cases nonbinarycyclekiller returned the same size solution as rspr. in instances terminusest returned the same size solution as rspr, in cases terminusest returned the same size solution as nonbinarycyclekiller, and in cases nonbinarycyclekiller returned the same size solution as rspr. in instances terminusest returned the same size solution as rspr, in cases terminusest returned the same size solution as nonbinarycyclekiller, and in cases nonbinarycyclekiller returned the same size solution as maf. in this dataset the gene trees were also nonbinary, meaning that nonbinarycyclekiller had to use the maf algorithm described in instead of rspr.  our experiments with binary trees show that cyclekiller is much faster than available exact methods once the input trees become sufficiently large and/or discordant. in over 96% of the cases cyclekiller finds the optimal solution and in the remaining cases it finds a solution very close to the optimum. we have shown that the most accurate mode of the program produces solutions that are at most a factor from the optimum. in practice, the average-case approximation ratio that we observed was. the fastest mode of the algorithm can be used on trees with thousands of leaves and provably constructs networks that are at most a factor of from the optimum. our experiments with nonbinary trees highlight once again that the cycle-breaking technique described in this article is intrinsically linked to the current state-of-the-art in maf algorithms. terminusest is faster than the most accurate mode of nonbinarycyclekiller when both trees are nonbinary due to the fact that maf solvers for two nonbinary trees have not yet been optimized to the same extent as their binary counterparts. in fact, terminusest is the best avaible exact method for nonbinary trees and can handle instances for which the optimum is up to 15- for other instances, nonbinarycyclekiller in its fastest mode is much faster than terminusest and produces solutions that are at most a factor from the optimum. finally, for instances with one binary and one nonbinary tree, the most accurate mode of nonbinarycyclekiller is again much faster than terminusest and produces solutions that are at most a factor from the optimum. endnotes a in, it has been shown that this number can go up to when running dendroscope on a similar processor but allocating all cores for one instance, i.e. exploiting the possibilities of parallel computation of this implementation. b we note that nonbinarycyclekiller uses a row-generation ilp formulation - based on - to solve dfvs, rather than the polynomial-size formulation used by cyclekiller. ilp is in neither case a bottleneck for the running time. competing interests the authors declare that they have no competing interests. authors’ contributions all the authors conceived the ideas, designed and conducted the experiments, and wrote and approved the paper. authors’ information leo van iersel was supported by a veni grant of the netherlands organisation for scientific research. nela lekić was supported by a vrije competitie grant of of the netherlands organisation for scientific research. 
hybridization numberphylogenetic networksapproximation algorithmsdirected feedback vertex set the problem of computing the minimum number of such events to explain discordance between two phylogenetic trees is a hard computational problem. as a spin-off from our simulations we also present terminusest, which is the fastest exact method currently available that can handle nonbinary trees: this is used to measure the accuracy of the nonbinarycyclekiller algorithm. even for binary trees, exact solvers struggle to solve instances with reticulation number larger than 40- 
 reticulate events play an important role in determining evolutionary relationships. the problem of computing the minimum number of such events to explain discordance between two phylogenetic trees is a hard computational problem. even for binary trees, exact solvers struggle to solve instances with reticulation number larger than 40-  here we present cyclekiller and nonbinarycyclekiller, the first methods to produce solutions verifiably close to optimality for instances with hundreds or even thousands of reticulations.  using simulations, we demonstrate that these algorithms run quickly for large and difficult instances, producing solutions that are very close to optimality. as a spin-off from our simulations we also present terminusest, which is the fastest exact method currently available that can handle nonbinary trees: this is used to measure the accuracy of the nonbinarycyclekiller algorithm. all three methods are based on extensions of previous theoretical work:1635- tcbb 10:18- sidma 28:49-66) and are publicly available. we also apply our methods to real data. hybridization numberphylogenetic networksapproximation algorithmsdirected feedback vertex set 