A practical approximation algorithm for solving massive instances of hybridization number for binary and nonbinary trees

Background

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-50.

Results

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.

Conclusions

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 (SIDMA 26(4):1635-1656, TCBB 10(1):18-25, SIDMA 28(1):49-66) and are publicly available. We also apply our methods to real data.

Phylogenetic trees are used in biology to represent the evolutionary history of a set $\mathcal{X}$ of species (or taxa) [1, 2]. They are trees whose leaves are bijectively labeled by $\mathcal{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 (phylogenetic) trees are rooted (phylogenetic) trees whose internal vertices have outdegree 2. Nonbinary rooted (phylogenetic) 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 [3]. For detailed background information we refer the reader to [4–6].

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 (e.g. aforementioned reticulate phenomena such as hybridization). 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 (a refinement of) 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-1 outdegree-1 vertices (i.e. replacing paths of these vertices by one edge). 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 (where an indegree-d reticulation counts as d−1 reticulation events). 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 [7], but it is not known whether it is in APX (i.e. whether it admits a polynomial-time approximation algorithm that achieves a constant approximation ratio).

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 (FPT) algorithms, where the parameter in question is the hybridization number r of the two trees. For an introduction to FPT we refer to [8, 9].

For binary trees, algorithmic progress has been considerable in this area, with various authors reporting increasingly sophisticated FPT algorithms [10–13]. The fastest algorithms currently implemented are the algorithm available inside the package DENDROSCOPE[14], based on [15], and the sequence of progressively faster algorithms in the HYBRIDNET family [11, 16, 17]. The fastest theoretical FPT algorithm has running time O(3.18^rn) [13], 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 [18, 19] and a (non-FPT) algorithm has recently been implemented in DENDROSCOPE[14].

Such (FPT) 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(r), 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-50 and instances of nonbinary trees when r is at most 5-10.

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 [20] 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 (DFVS). 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(logr log logr) 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 [20] 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 [21], 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 1 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(c+1) for the hybridization number problem on two binary trees and an algorithm with approximation ratio d(c+3) for the hybridization number problem on two nonbinary trees by combining a c-approximation for the problem MAF (Maximum Agreement Forest) with a d-approximation for the problem DFVS. Both these problems are NP-hard so polynomial-time algorithms attaining c=1 or d=1 are not realistic. Nevertheless, there exist extremely fast FPT algorithms for solving MAF on binary trees exactly (i.e. c=1), the fastest is RSPR by Whidden, Beiko and Zeh [22, 23] although the MAF algorithm inside [17] is also competitive. Moreover, we observe that the type of DFVS instances that arise in practice can easily be solved using Integer Linear Programming (ILP) (and freely-available ILP solver technology such as GLPK), so d=1 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 2 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 2 is often pessimistic, with much better approximation ratios observed in experiments (1.003 on average for the simulations presented in this article). 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=1. Fortunately there exist linear-time algorithms which achieve c=3 [13]. This, coupled with the fact that (even for such trees) it remains feasible to use an exact (d=1) 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 4 is a worst-case bound and we suspect that in practice we are doing much better than 4. 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 [13, 21, 24]. In case when one of the input trees is binary we can still use the exact (thus c=1) and approximate (c=3) algorithms given in [13] (referred to as RSPR) 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 (c=1) and approximate (c=4) algorithms described in [21]. We then obtain 4- and 7-approximations (because in the nonbinary case d=1 is still easily attainable using ILP).

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 [18]. 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[25], from http://homepages.cwi.nl/~iersel/cyclekiller[26], and from http://skelk.sdf-eu.org/terminusest[27].

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 (both exact and approximate), 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 [20], which also proposed using DFVS but beginning from a trivial Agreement Forest (AF) known as a chain forest. Here we use a smarter starting point: an (approximate) MAF, and it is this insight which makes a 2-approximation (rather than the 6-approximation implied by [20]) 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. [13] – 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 (approximate) 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 [21] in which we introduce it. Note that our results for the binary case do not follow from the results for the nonbinary case in [21] 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 (respectively, nearly 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 (i.e. for the preliminary version of this article [28]) HYBRIDNET was the fastest algorithm available in its family. It has recently been superceded by the faster ULTRANET[17]. 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 (it is now at v1.2), 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 (and the theoretical algorithm in [13]) implicitly have to explore exponentially many optimal and sub-optimal solutions to the MAF problem, the running time of MAF solvers (and thus also CYCLEKILLER) 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.

Preliminaries

Let be a finite set (e.g. of species). A rooted phylogenetic-tree is a rooted tree with no vertices with indegree 1 and outdegree 1, a root with indegree 0 and outdegree at least 2, and leaves bijectively labelled by the elements of . We identify each leaf with its label and use L(T) to refer to the leaf set (or label 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 ${\mathcal{X}}^{\prime }\subset \mathcal{X}$, we use the notation $T\left({\mathcal{X}}^{\prime }\right)$ to denote the minimal subtree of T that contains all elements of ${\mathcal{X}}^{\prime }$ and $T|{\mathcal{X}}^{\prime }$ denotes the result of suppressing all indegree-1 outdegree-1 vertices in $T\left({\mathcal{X}}^{\prime }\right)$.

The following definitions apply only to binary trees. Definitions for nonbinary trees are analogous but slightly more technical [21].

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 for T if, for all $F\in \mathcal{F}$, T|L(F) is isomorphic to F and the trees $\left\{T\right(L\left(F\right)),F\in \mathcal{F}\}$ are vertex-disjoint subtrees of T whose leaf-set union equals L(T). 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 $\left|\mathcal{F}\right|$.

We define cleaning up a directed graph as repeatedly suppressing indegree-1 outdegree-1 vertices, removing indegree-0 outdegree-1 vertices and removing unlabelled outdegree-0 vertices until no such operation is possible. Observe that, if is a forest for T, can be obtained from T by removing $\left|\mathcal{F}\right|-1$ edges and cleaning up. From now on we consider T₁,T₂ as trees on the same taxon set.

Problem: Maximum Agreement Forest (MAF)Instance: Two rooted, binary phylogenetic trees T₁ and T₂.Solution: An agreement forest of T₁ and T₂.Objective: Minimize $\left|\mathcal{F}\right|-1$.

The directed graph $\mathit{\text{IG}}(T_1,T_2,\mathcal{F})$, called the inheritance graph, is the directed graph whose vertices are the components of and which has an edge (F,F^′) precisely if either

• there is a directed path in T₁ from the root of T₁(L(F)) to the root of T₁(L(F^′)) or;

• there is a directed path in T₂ from the root of T₂(L(F)) to the root of T₂(L(F^′)).

An agreement forest of T₁ and T₂ is called an acyclic agreement forest if the graph $\mathit{\text{IG}}(T_1,T_2,\mathcal{F})$ is acyclic. A maximum acyclic agreement forest (MAAF) of T₁ and T₂ is an acyclic agreement forest of T₁ and T₂ with a minimum number of components.

Problem: Maximum Acyclic Agreement Forest (MAAF)Instance: Two rooted, binary phylogenetic trees T₁ and T₂. Solution: An acyclic agreement forest of T₁ and T₂.Objective: Minimize $\left|\mathcal{F}\right|-1$.

We use MAF (T₁,T₂) and MAAF (T₁,T₂) 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 1 and outdegree 1 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 δ⁻(v)≥2 (with δ⁻(v) denoting the indegree of v). The reticulation number (sometimes also called hybridization number) of a network N with root ρ is given by

It was shown that the optimum to MAAF is equal to the optimum of the following problem [29].

Problem: MINIMUMHYBRIDIZATIONInstance: Two rooted binary phylogenetic trees T₁ and T₂.

Solution: A rooted phylogenetic network N that displays T₁ and T₂.Objective: Minimize r(N).

Moreover, it was shown that, for two trees T₁,T₂, any acyclic agreement forest for T₁ and T₂ with k+1 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 (DFVS)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 (acyclic) agreement forests are different in the nonbinary case [21].

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(c+1)-approximation for MAAF (and thus for MINIMUMHYBRIDIZATION).

Note that this theorem does not follow from Theorem 2.1 of [21], since there the approximation ratio for MAAF is a d(c+3)-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,

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(T₁,T₂)+|M|.

Proof

Consider a maximum acyclic agreement forest F of T₁ and T₂. For i∈{1,2}, F can be obtained from T _i by removing a set of edges, say $E_F^i$, and cleaning up. Moreover, also M can be obtained from T _i by removing a set of edges, say $E_M^i$, and cleaning up.

Now consider the forest S obtained from T₁ by removing $E_M^1\cup E_F^1$and cleaning up. Then,

• S is an agreement forest of T₁ and T₂ because it can be obtained from T₂ by removing edges $E_M^2\cup E_F^2$ 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, $\left|S\right|\le \left|E_F^1\right|+\left|E_M^1\right|+1$. The lemma follows since $\left|E_F^1\right|=\text{MAAF}(T_1,T_2)$ and $\left|M\right|=\left|E_M^1\right|+1$.

Let OptSplitting${}_{T_1,T_2}\left(M\right)$ denote the size of a minimum-size M-splitting. Combining Lemma 2 and Eq. 1, we obtain

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 (called the extended inheritance graph) as follows. For every vertex of M that has outdegree 2 (in M), 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₂ (or in both) there is a directed path from u to v. An example is in Figure 1. We claim the following.

Two binary trees T ₁ and T ₂ and the auxiliary graph D . A maximum agreement forest M of T₁ and T₂ is obtained by deleting the dashed edges. Graph D can be made acyclic by deleting either both filled or both unfilled vertices. Hence, removing either v₁ and v₂ or v₃ and v₄ from M makes it an acyclic agreement forest for T₁ and T₂, see Lemma 3. The acyclic agreement forest M∖{v₁,v₂} obtained by removing v₁ and v₂ from M is depicted on the right.

Full size image

Lemma 3.

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 2 in M∖V^′, they had outdegree 2 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-2 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₁ (or 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, OptSplitting${}_{T_1,T_2}\left(M\right)=\left|M\right|+\left|\text{MFVS}\right|$. 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 2 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

where the last inequality follows from Eq. 2. Thus, is a d(c+1)-approximation to MAAF, which concludes the proof of Theorem 1.

Theorem 1 implies that a solution to the MAAF problem for a given instance can be constructed by (i) finding a solution to the MAF problem for the same instance (ii) constructing the extended inheritance graph D for (iii) finding a solution V for the DFVS problem on the graph D and (iv) modifying accordantly to V.

Practical experiments with binary trees

To assess the performance of CYCLEKILLER, a simulation study was undertaken. We generated 3 synthetic datasets, an easy, a medium and a hard one, containing respectively 800, 640 and 640 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 (note that this number is an upper bound on the actual rSPR distance). The 800 pairs of rooted binary phylogenetic trees were created by varying n in {20,50,100,200} and k in {5,10,...,25}, and then creating 40 different instances per each combination of parameters. Each pair (T₁,T₂) of rooted binary phylogenetic trees for a given set of parameters n and k is created as follows: The first tree T₁ on $\mathcal{X}=\{x_1,\dots ,x_n\}$ 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 0, are randomly picked and a new vertex w, along with two new edges (w,u) and (w,v), 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 {50,100,200,300} and k in {15,25,40,55} for the medium one and n in {100,200,400,500} and k in {40,60,80,100} for the hard one.

The exact hybridization number has been computed by HYBRIDNET[11], available from http://www.cs.cityu.edu.hk/~lwang/software/Hn/treeComp.html or with DENDROSCOPE[14], 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 E5506 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 (see Section “New algorithms: CYCLEKILLER and NONBINARYCYCLEKILLER”). We used the program RSPR v1.03 [22, 23] to solve or approximate MAF and GLPK v4.47 (http://www.gnu.org/software/glpk/) to solve the following simple polynomial-size ILP formulation of DFVS:

Given a directed graph D=(V,E), 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 33 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 9 runs of the easy data, DENDROSCOPE and HYBRIDNET did not complete within 10 days. Table 1 shows a summary of the results. It can be seen that CYCLEKILLER was much faster than the exact algorithms. Moreover, for 96.6% 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 2, in our experiments it performed very close to a 1-approximation.

For the medium data set, CYCLEKILLER finished with the 2-approx option for 613 instances, and for the remaining ones with the 4-approx option. The exact algorithms could compute the hybridization number for only 199 instances (out of 640). For 97.5% of these instances, CYCLEKILLER also found an optimal solution, but with a much better running time. Regarding the hard data set, 444 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 25^a. In contrast, the 2-approx option of CYCLEKILLER could be used for instances for which the size of a MAF was up to 97, and thus for instances for which the hybridization number was at least 97.

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 10,000 leaves and up to 10,000 reticulations within 10 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 2,500. 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 1.2.

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∈{20,50,100} and the number of rSPR-moves k∈{5,10,15,20}, and an additional parameter ρ∈{25,50,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 (upon which NONBINARYCYCLEKILLER is built) and exact solvers for nonbinary MINIMUMHYBRIDIZATION (which is important to measure the accuracy of NONBINARYCYCLEKILLER in practice) are slower than their binary counterparts.

We performed two runs of experiments^b. 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 [13, 22], 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 4. When RSPR is used in its 3-approximation mode, NONBINARYCYCLEKILLER yields a theoretical worst-case approximation ratio of 6 (see Section “New algorithms: CYCLEKILLER and NONBINARYCYCLEKILLER”). The results of this run are summarized in Table 2.

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 [21]. This makes NONBINARYCYCLEKILLER behave as a 4-approximation and 7-approximation respectively (see Section “New algorithms: CYCLEKILLER and NONBINARYCYCLEKILLER”). Note that the exact algorithm [21] 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 3.

The exact hybridization number in both runs was computed by TERMINUSEST[27].

Each instance that took longer than 10 minutes to compute was aborted and the running time was set to 600 seconds. The averages of the running-times are taken over all instances, with running-time taken to be 600 if the program timed out for that instance. (We used a shorter time-out than in the binary experiments because of the observation that, in the nonbinary case, exact algorithms running longer than 10 minutes almost always took longer than 60 minutes too.)

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 (data not shown).

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∈{20,50}, k∈{5,10,15} and for the tricky set n∈{50,100}, k=20. In addition we varied the percentage of contracted edges (in a single tree in the first run and in both trees in the second run).

In Table 2 we show running times and solution quality of our algorithm when one of the input trees is binary. For the simple set of instances (regardless of the percentage of edge-contractions) 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 1. 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 1. 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 1.75 for the 4-approximation and 3 for the 6-approximation.

Table 3 shows our results on instances with two nonbinary trees. The exact algorithm for MAF is in this case much slower and this affects the running times even for the simple set. While the 4-approximation version has an average approximation ratio very close to 1 again, the running time is in this case worse than that of TERMINUSEST. For the tricky set the situation is even more significant; the exact MAF algorithm cannot deal with reticulation numbers above 15, while TERMINUSEST can get slightly further. On the other hand, the 7-approximation still runs much faster than TERMINUSEST, both for simple and tricky instances, while having an average approximation ratio of less than 2.6. The practical worst-case approximation ratio observed in these experiments was 1.5 for the 4-approximation and 4 for the 7-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 [18]. It is plausible that with 50%-contraction the instances suffer from the disadvantage of relatively high hybridization number without the branching advantages associated with (relatively) 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:

• 1000 leaves, 25%-contraction, on average 995 reticulations in 63 sec.

• 1000 leaves, 50%-contraction, on average 989 reticulations in 82 sec.

• 1000 leaves, 75%-contraction, on average 840 reticulations, in 656 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 [30] – version 3, 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 9903 rooted phylogenetic trees associated to the gene families contained in the database (the gene trees), along with the rooted phylogenetic tree describing the history of the twenty-two species contained in GreenPhylDB (the species tree). Note that the species tree for these species is not completely resolved, i.e. it is nonbinary. Among the gene trees, 2769 contain less than 3 species and they were discarded. Of the remaining 7134 trees, only 204 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 (MUL trees). We thus transformed the MUL trees into trees containing single copies of labels, applying the tools described in [31, 32] to the forest F of 7134 trees. As in Section 4.1 of [31], we obtained four data sets: F₁, F₂, $F_3^p$ and $F_3^s$, respectively containing 204, 1003, 5924 and 5789 trees. Note that only $F_3^s$ contains nonbinary trees. Finally, for each single labeled tree $G\in (F_1\cup F_2\cup F_3^p\cup F_3^s)$, we restricted the species tree S (containing 22 taxa) to the leaves of G and we applied our methods to all so obtained pairs (restricted S,G). The results are presented in Tables 4, 5, 6 and 7. 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.

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 2 from the optimum. In practice, the average-case approximation ratio that we observed was 1.003. 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 4 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-20. For other instances, NONBINARYCYCLEKILLER in its fastest mode is much faster than TERMINUSEST and produces solutions that are at most a factor 4 from the optimum (less than 2.6 on average).

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 1.75 from the optimum (less than 1.011 on average).

^a In [15], it has been shown that this number can go up to 40 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 [33] - to solve DFVS, rather than the polynomial-size formulation used by CYCLEKILLER. ILP is in neither case a bottleneck for the running time.

Leo van Iersel was supported by a Veni grant of The Netherlands Organisation for Scientific Research (NWO). Nela Lekić was supported by a Vrije Competitie grant of of The Netherlands Organisation for Scientific Research (NWO).

A preliminary version of this paper (restricted to the binary case) appeared in the proceedings of the 12th Workshop on Algorithms in Bioinformatics (WABI 2012) [28]. We thank Simone Linz and Leen Stougie for fruitful discussions. This publication is the contribution no. 2014-040 of the Institut des Sciences de l’Evolution de Montpellier (ISE-M, UMR 5554). This work has been partially funded by the French Agence Nationale de la Recherche, Investissements d’avenir/Bioinformatique (ANR-10-BINF-01-02, Ancestrome), and it has benefited from the ISE-M computing facilities.

Authors and Affiliations

1. Centrum Wiskunde & Informatica (CWI), P.O. Box 94079, 1090, GB, Amsterdam, The Netherlands

   Leo van Iersel

2. Department of Knowledge Engineering (DKE), Maastricht University, P.O. Box 616, 6200, MD, Maastricht, The Netherlands

   Steven Kelk & Nela Lekić

3. ISEM, CNRS – Université Montpellier II, Place Eugène Bataillon, 34095, Montpellier, France

   Celine Scornavacca

Corresponding author

Correspondence to Leo van Iersel.

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.

This article is published under license to BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly credited. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.

Reprints and permissions