Orientability of Undirected Phylogenetic Networks to a Desired Class: Practical Algorithms and Application to Tree-Child Orientation

Phylogenetic networks are a powerful tool for representing complex evolutionary relationships between species that cannot be adequately modelled by trees. These networks are particularly useful in the presence of reticulate events, such as hybridisation, horizontal gene transfer (HGT) and recombination. Hybridisation refers to the interbreeding of individuals from different species, which can lead to the formation of a new hybrid species that shares genetic material from both parent species. HGT is the transmission of copied genetic material to another organism without being its offspring, a process particularly common in bacteria and archaea. Recombination is the exchange of genetic material between different genomes and is a common occurrence not only in viruses but also in bacteria, eukaryotes, and other organisms. While phylogenetic trees depict the hierarchical branching of evolutionary lineages, networks can represent both the hierarchical and non-hierarchical connections resulting from reticulate events, providing a more nuanced view of evolutionary history.

However, constructing directed phylogenetic networks from biological data remains a challenging task. Distance-based methods, such as Neighbor-Net, are widely used because they are scalable and helpful in visualising the data. However, the resulting networks are inevitably undirected, often making it difficult to interpret the evolutionary history. To provide a more informative representation of the data, it is meaningful to develop a method for transforming undirected graphs into directed phylogenetic networks in a way that ensures the resulting network has a desired property.

Recently, Huber et al. [1] introduced two different orientation problems, each considering orientation under a different constraint. The Constrained Orientation problem asks whether a given undirected phylogenetic network can be oriented to a directed network under the constraint of a given root edge and desired in-degrees of all vertices, and asks to find a feasible orientation under that constraint if one exists. In [1], it was shown that such a feasible orientation is unique if one exists, and a linear time algorithm for solving this problem was provided. However, when an undirected network has been created from data, it is not usually the case that there is complete knowledge of where to insert the root and which vertices are reticulations. The $\mathcal{C}$-Orientation problem does not constrain the position of the root or the in-degree of the vertices. Instead, it asks whether a given binary network can be oriented to a directed phylogenetic network belonging to a desired class $\mathcal{C}$.

The complexity of $\mathcal{C}$-Orientation is not fully understood for many classes $\mathcal{C}$, and no study has discussed practically useful methods. In [1], Tree-Based Orientation was shown to be NP-hard. Maeda et al. [2] conjectured that Tree-Child Orientation is NP-hard. Bulteau et al. [3] studied a similar problem and showed that for a graph with maximum degree five, determining whether the graph can be oriented to a tree-child network with a designated root vertex is NP-hard (Corollary 5 in [3]). Without the assumption of maximum degree five, the complexity of this problem is still unclear. Indeed, [3] was concluded by noting that even for graphs of maximum degree three (i.e. binary networks), finding an orientation to a desired class, such as tree-child, tree-based or reticulation-visible networks, may not be easy. Although [1] provided FPT algorithms for a special case of $\mathcal{C}$-Orientation where $\mathcal{C}$ satisfies several conditions, which are theoretically applicable to various classes including tree-child networks but are not easy to implement, there is no FPT algorithm to solve $\mathcal{C}$-Orientation in its general form. Also, no studies have pursued practically useful heuristics.

In this paper, we provide a practically efficient exponential time algorithm for $\mathcal{C}$-Orientation (Algorithm 1), which is FPT in both the reticulation number and the maximum size of minimal basic cycles selected by the algorithm. We also present a heuristic method for Tree-Child Orientation (Algorithm 2) which, although still exponential, runs very fast in practice because it only considers reticulation placements to maximise the sum of their pairwise distances. Using artificially generated networks, we compare the accuracy and execution time of the proposed methods for solving Tree-Child Orientation with those of the existing exponential time algorithm for $\mathcal{C}$-Orientation (Algorithm 2 in [1]). Our theoretical and empirical results demonstrate the usefulness of Algorithm 1, especially for relatively large input graphs with $5$ or more reticulations, where the exponential time method in [1] becomes computationally infeasible. Our Algorithm 2, while much faster than Algorithm 1, tends to decrease in accuracy as the reticulation number increases.

The rest of the paper is structured as follows. Section 2 provides the necessary mathematical definitions and notation, including the definition of phylogenetic networks. Section 3 briefly reviews relevant results from [1] and formally states the problems of interest. Section 4 gives the theoretical background of our proposed methods, including the concept of ‘cycle basis’ and a theorem that allows us to reduce the search space (Theorem 4). Section 5 describes the proposed exact method (Algorithm 1) and a heuristic method (Algorithm 2). Theorem 7 ensures that the heuristic is correct when $r\leq 2$. We analyse the time complexity of Algorithm 1. Section 6 explains the experimental setup, including the method used to create undirected graphs (details can be found in the appendix), and presents the results of three experiments. Section 7 discusses the limitations and biological application of the heuristic method. Finally, Section 8 concludes the paper and outlines future research directions.

For an undirected phylogenetic network $N$ on $X$, orienting $N$ is a procedure that inserts the root $\rho$ into a unique edge $e_{\rho}=\{u,v\}$ of $N$ by replacing the edge $\{u,v\}$ with the arcs $(\rho,u),(\rho,v)$ and then orienting the other edges so that the resulting graph $\vec{N}$ is a directed phylogenetic network on $X$. In other words, a directed phylogenetic network $\vec{N}$ on $X$ is an orientation of $N$ if the underlying graph $U(\vec{N})$ becomes isomorphic to $N$ after suppressing the root $\rho$ of $U(\vec{N})$ (i.e. replacing the undirected path $u-\rho-v$ with the edge $\{u,v\}$).

Determining the orientability of undirected graphs to directed graphs with desired properties is a classic topic in graph theory and its application (e.g. [6]), but orientability problems for phylogenetic networks have only recently been studied [1, 3]. Huber et al. [1] discussed two types of problems, Constrained Orientation and $\mathcal{C}$-Orientation.

Before discussing Problem 3, we consider a general setting where we want to orient any undirected phylogenetic network $N=(V,E)$ on $X$ to a rooted directed phylogenetic network $\vec{N}$ on $X$. Then $\vec{N}$ must contain $r=|E|-|V|+1$ reticulation vertices (which can be easily verified by induction, and can also be derived from the equations in Lemma 2.1 in [9]), and so we need to decide which $r$ vertices among $|V|-|X|$ internal vertices of $N$ will have in-degree 2 in $\vec{N}$. The number of possible ways to select $r$ vertices from non-leaf vertices of $N$ is $\binom{|V|-|X|}{r}$, which is exponential. We will now consider how we can reduce the number of candidates to examine.

Lemma 2 allows us to exclude some of the inappropriate reticulation placements (specifically, cases where there is a cycle having no reticulation) that make orientation impossible. However, checking all cycles in a graph is computationally inefficient. To reduce the search space, we will now introduce some relevant concepts.

For a connected undirected graph $N=(V,E)$, the cycle rank of $N$ is defined to be the number $r:=|E|-|V|+1$ (e.g. p.24 in [10]), which is also known as circuit rank, cyclomatic number and the (first-order) Betti number of $N$. Note that $r$ is zero if $N$ is a tree and that $r$ is the number of desired reticulation vertices if $N$ is a binary phylogenetic network that is an instance of Problem 1. The cycle rank $r$ of $N$ can also be interpreted as the rank of a vector space called the ‘cycle space’, where each of the $r$ basis vectors, called a basic cycle, corresponds to the edge-set of a simple cycle in $N$. If we define the summation of cycles $C$ and $C^{\prime}$ as the cycle induced by the symmetric difference of their edge-sets, then the cycle space of $N$ can be identified with the set of even-degree subgraphs of $N$, so any cycle in $N$ can be expressed as a sum of basic cycles in a cycle basis $\mathcal{S}$ (for details, see Section 1.9 of [11], Section 4.3 of [12] and Section 6.4.2 of [10]).

Proposition 3 will be used in the proof of the main theorem (Theorem 4).

Theorem 4 provides a necessary condition for the feasibility of a reticulation placement $V_{R}$ in Problem 1 (binary version), where feasibility means that $(N,e_{\rho},V_{R})$ admits an acyclic orientation for some $e_{\rho}$. This implies that we need not consider all $\binom{|V|-|X|}{r}$ reticulation placements in solving Problem 2, regardless of the class $C$ we are interested in. Therefore, we may use an arbitrary cycle basis $\mathcal{S}$ of $N$ and choose exactly one reticulation vertex from each of the $r$ basic cycles, thereby reducing the search space for feasible reticulation placements in both Problems 2 and 3.

In solving Problem 3, we will also use the following results. From the forbidden structures of tree-child networks (Fig. 1), Lemma 5 follows immediately, and this leads to Theorem 6.

We note that Theorem 6 provides a sufficient condition for tree-child orientability, not a necessary condition. In fact, a tree-child network can contain a pair of reticulations whose distance is less than 3 (see Fig. 2).

Based on Theorem 4, we propose an exact method (Algorithm 1) for $\mathcal{C}$-Orientation (Problem 2) and a heuristic method (Algorithm 2) for Tree-Child Orientation (Problem 3). While Theorem 4 holds for any cycle basis, our algorithms use a ‘minimal’ one for computational efficiency. A cycle basis $\mathcal{S}$ of $N$ is called minimal if the sum of the lengths of all cycles in $\mathcal{S}$ is not greater than that of any other cycle basis of $N$, i.e. minimising $\sum_{C\in\mathcal{S}}{|E(C)|}$. Intuitively, focusing on a minimal cycle basis allows us to use basic cycles that do not overlap too much. The problem of computing a minimal cycle basis has been extensively studied (e.g. [13, 14, 15]; see Chapter 7 of [16] for a brief literature review) and polynomial time algorithms exist (e.g. [16, 17, 15]). By using an algorithm given in [17], a minimal cycle basis of $N=(V,E)$ can be computed in $O(|V||E|^{2}/\log|V|)$.

For both Problems 2 and 3, our proposed methods first compute a minimal cycle basis $\mathcal{S}=\{C_{1},\dots,C_{r}\}$ of a given network $N$. Then, for Problem 2, our exact method (Algorithm 1) repeatedly picks exactly one reticulation vertex $v_{i}$ from each basic cycle $C_{i}$ to obtain a candidate set $V_{R}$ of $r$ reticulations, and solves Problem 1 for all $(e_{\rho},V_{R})$ until it finds a $\mathcal{C}$-orientation of $N$. While this algorithm still requires exponential time, it is more efficient than the exponential time method in [1] because the search space for appropriate $V_{R}$ is smaller than $\binom{O(|V|)}{r}$.

For Problem 3, we first note that a straightforward way to reduce the search space for $V_{R}$ is to exclude any candidate sets containing adjacent reticulations, as such configurations are invalid in tree-child networks (recall Fig. 1, right panel). Beyond this basic constraint, we introduce a novel approach to further reduce the search space based on Theorem 6, which suggests that reticulations should be placed as far apart as possible. This insight leads to our heuristic method (Algorithm 2) that considers only those reticulation sets $V_{R}$ that maximise the sum of pairwise distances between reticulations. While this significant reduction in the search space makes the algorithm faster than Algorithm 1, it may fail to find an existing tree-child orientation; more precisely, while false positives cannot occur, a ‘NO’ output from Algorithm 2 should be interpreted as ‘Probably NO’. However, Theorem 7 ensures that Algorithm 2 works correctly when $r$ is very small.

Theorem 8 states that Algorithm 1 is FPT both in the reticulation number $r$ and the size $c$ of longest basic cycles in $\mathcal{S}$ used in the computation.

Theorem 8 implies that the unparameterised worst-case complexity of Algorithm 1 is $O(|V|^{r+2})$. This complexity is equivalent to that of the exponential time method (Algorithm 2 in [1]) that performs $O(|V||E|)$ time calculations for all $\binom{|V|}{r}=O(|V|^{r})$ reticulation placements. Although two FPT algorithms for a special case of $\mathcal{C}$-Orientation were proposed in [1], Algorithm 1 differs in its parameterisation from them. Specifically, our Algorithm 1 is parameterised by $r$ and $c$, while the FPT algorithms in [1] are parameterised by $r$ and by the level of $N$, respectively.

Theorem 8 also shows that the size of search space depends on the choice of $\mathcal{S}$ at Line 2 of Algorithm 1 although $c$ never exceeds the size of longest cycles in $N$ (the same applies to Line 2 of Algorithm 2). To illustrate this, consider two minimal cycle bases $\mathcal{S}=\{C_{1},C_{2},C_{3}\}$ with $|V(C_{1})|=|V(C_{2})|=|V(C_{3})|=4$ and $\mathcal{S}^{\prime}=\{C_{1}^{\prime},C_{2}^{\prime},C_{3}^{\prime}\}$ where $|V(C_{1}^{\prime})|=3$, $|V(C_{2}^{\prime})|=4$ and $|V(C_{3}^{\prime})|=5$ (note that they have the same total length $4+4+4=3+4+5$). When the former $\mathcal{S}$ is selected, the number of elements of $S$ at Line 3 of either algorithm is $|V(C_{1})|\times|V(C_{2})|\times|V(C_{3})|=4^{3}$, whereas $|S|=|V(C_{1}^{\prime})|\times|V(C_{2}^{\prime})|\times|V(C_{3}^{\prime})|=3\times 4\times 5$ for the latter $\mathcal{S}^{\prime}$.

We implemented our two proposed methods (Algorithms 1 and 2) and the existing exponential time algorithm described by Huber et al. (Algorithm 2 in [1]) using Python 3.11.6. In the implementation of Algorithms 1 and 2, we used the minimum_cycle_basis function from the Python package networkx to compute a minimal cycle basis. We note that the algorithm implemented in networkx is not the $O(|V||E|^{2}/\log|V|)$ algorithm given in [17] but the $O(|E|^{3}+|E||V|^{2}\log{|V|})$ algorithm given in Section 7.2 of [16]. Undirected binary phylogenetic networks were generated as test data using a method described in Appendix. The source code, test data, and the program used to generate the data are available at https://github.com/hayamizu-lab/tree-child-orienter. The full details of the results can be found at https://github.com/hayamizu-lab/tree-child-orienter/tree/main/results.

The $C$-Orientation problem, which asks whether a given undirected binary phylogenetic network can be oriented to a directed phylogenetic network of a desired class $C$, is an important computational problem in phylogenetics. The complexity of this problem remains unknown for many network classes $C$, including the class of binary tree-child networks. A simple exponential time algorithm for $\mathcal{C}$-Orientation was provided in [1]. FPT algorithms for a special case of the problem were also proposed in [1], but their practical application is limited due to the intricate nature of the procedures and the challenges in implementation, despite the constraints imposed on $\mathcal{C}$. Additionally, no study has explored heuristic approaches to solve $\mathcal{C}$-Orientation in practice, even for a particular class $\mathcal{C}$ such as tree-child networks.

In this paper, we have proposed a simple, easy to implement, practical exact FPT algorithm (Algorithm 1) for $\mathcal{C}$-Orientation and a heuristic algorithm for Tree-Child Orientation (Algorithm 2) based on Theorem 4. They improve the search space of the existing simple exponential time algorithm by using a cycle basis to reduce the number of possible reticulation placements. Our experiments showed that Algorithm 1 is significantly faster in practice than a state-of-the-art exponential time algorithm in [1]. Algorithm 2 is even faster, with a trade-off between the accuracy and the reticulation number. Further research using larger and more diverse datasets could provide more insight into the strengths and limitations of the proposed methods. Their usefulness and effectiveness should also be tested in real-world data analysis. As tree-child networks can describe evolutionary histories that do not involve frequent reticulation events, one can find various tree-child networks in the literature other than the one we used in Section 7.2, such as a hybridisation network of bread wheat (Fig. 3 in [19]), a hybridisation network of mosquitoes (Fig. 15 in [20]) and an admixture network of diverse populations (Fig. 3 in [21]).

Although we have used Tree-Child Orientation as a case study for performance evaluation, Algorithm 1 is applicable to orientation problems for classes $\mathcal{C}$ other than tree-child networks. For example, orientation for stack-free networks [22] or tree-based networks [23] is expected to be a good application, because one can quickly decide whether a given network belongs to such a class. Speeding up Algorithm 1 and extending it to non-binary networks are topics for future research. It would be also interesting to consider alternative FPT algorithms that do not use a cycle basis.

Improving the accuracy of Algorithm 2 is also an interesting direction for future research. As discussed in Section 7, there are clearly rooms for improvement in the current objective function. Introducing a more suitable objective function may lead to the development of faster and more useful tree-child orientation heuristics.