Finding a most parsimonious or likely tree in a network with respect to an alignment.

Phylogenetic networks are often constructed by merging multiple conflicting phylogenetic signals into a directed acyclic graph. It is interesting to explore whether a network constructed in this way induces biologically-relevant phylogenetic signals that were not present in the input. Here we show that, given a multiple alignment A for a set of taxa X and a rooted phylogenetic network N whose leaves are labelled by X, it is NP-hard to locate a most parsimonious phylogenetic tree displayed by N (with respect to A) even when the level of N-the maximum number of reticulation nodes within a biconnected component-is 1 and A contains only 2 distinct states. (If, additionally, gaps are allowed the problem becomes APX-hard.) We also show that under the same conditions, and assuming a simple binary symmetric model of character evolution, finding a most likely tree displayed by the network is NP-hard. These negative results contrast with earlier work on parsimony in which it is shown that if A consists of a single column the problem is fixed parameter tractable in the level. We conclude with a discussion of why, despite the NP-hardness, both the parsimony and likelihood problem can likely be well-solved in practice.

Introduction

Rooted phylogenetic networks are generalizations of rooted phylogenetic trees that allow horizontal evolutionary events such as horizontal gene transfer, recombination and hybridization to be modelled (Huson et al. 2010; Morrison 2011; Gusfield 2014). This is achieved by allowing nodes with indegree 2 or higher, known as reticulation nodes. Recent years have seen an explosion of interest in constructing rooted phylogenetic networks, fuelled by the growing awareness that incongruence in phylogenetic and phylogenomic datasets is not simply a question of evolutionary “noise”, but sometimes the result of evolutionary phenomena more complex than speciation and mutation (e.g. Zhaxybayeva and Doolittle 2011; Abbott et al. 2013; Vuilleumier and Bonhoeffer 2015).

Although many modelling questions surrounding the construction of phylogenetic networks are still to be answered, it is commonplace to associate a rooted phylogenetic network with the set of rooted phylogenetic trees that it contains (“displays”). Informally speaking, a rooted phylogenetic network displays a rooted phylogenetic tree if the tree can be topologically embedded inside the network. A network is not necessarily defined by the set of trees it displays, but the notion of display is nevertheless a recurring theme in the literature, since networks themselves are often constructed by merging phylogenetic trees subject to some optimality criterion. It is well-known that it is NP-hard to determine whether a network displays a given tree (Kanj et al. 2008), although on many restricted classes of phylogenetic networks the problem is polynomial-time solvable (Van Iersel et al. 2010; Fakcharoenphol et al. 2015; Gambette et al. 2015).

Although it is important to be able to determine whether a network displays a given tree, we may also wish to ask what the “best” tree is within the network, subject to some optimality criterion. For example, given a network N and a multiple alignment A, we may wish to ask for a tree T displayed by N with lowest parsimony score with respect to A. Similarly, if the network is decorated by edge lengths or probabilities, we may wish to identify a most likely tree displayed by the network, that is, a tree that maximizes the probability of generating A, under a given model of evolution. Such questions are natural, as the two following examples show. First, phylogenetic networks are often constructed by topologically merging incongruent phylogenetic signals (e.g. Kelk and Scornavacca 2014), and it is insightful to ask whether the network thus constructed displays trees which have interesting properties (such as a low parsimony score or a high likelihood) which were not in the input. Second, we may wish to perform classical phylogenetic tree construction under criteria such as maximum parsimony or maximum likelihood (e.g. Jin et al. 2006, 2007), but within the restricted space of trees displayed by a given network.

Problems of the above kind are already known to be NP-hard, since it is NP-hard to determine a most parsimonious tree T displayed by a given network N even when the alignment A consists of a single column and the network is binary (Fischer et al. 2015). However, the gadgets used in that hardness reduction produce networks with very high level, where level is the maximum number of reticulation nodes in a biconnected component of the network. On the positive side, the same article shows that the problem on an alignment consisting of a single column is FPT in the level of the network. This means that, on a network with level k, the problem can be solved in time where f is a function that depends only on k and n is the size of the network. Such results are useful in practice, when (as is often the case) k is small.

The question emerges whether the positive FPT result goes through when A does not consist of a single column, but potentially many columns – a problem introduced more than one decade ago (Nakhleh et al. 2005). Here we show that this is not the case. We prove the rather negative result that locating a most parsimonious tree in a rooted binary network N is NP-hard, even under the following restricted circumstances: (1) each biconnected component of the network contains exactly one reticulation node (i.e. is “level-1”); (2) each biconnected component of the network has exactly three outgoing arcs; (3) the alignment A consists of two states. If indel symbols are permitted then the problem is not only NP-hard, but also difficult to approximate well (APX-hard). If any of the conditions (1)-(3) are further strengthened (respectively: the network becomes a tree; the reticulation nodes become redundant; the alignment becomes uninformative), the problem becomes trivially solvable, so in some sense this is a “best possible” (or the “worst possible”, depending on your perspective) hardness result. Next, we consider the question of identifying a most likely tree in the network. We obtain NP-hardness under the same restrictions (1)-(3), subject to the simple binary symmetric model of character evolution. It is no coincidence that restrictions (1)-(3) again apply, since the hardness of the likelihood question is established by reducing the parsimony variant of the problem to it. Specifically, we show that a most likely tree displayed by a network with sufficiently short branches is necessarily also a most parsimonious tree.

Although the main results in this paper are negative, some reasons for hope are given in the conclusion.

Preliminaries

A rooted binary phylogenetic networkN on a set X of taxa is a directed acyclic graph where the leaves (nodes of indegree-1 and outdegree-0) are bijectively labelled by X, there is a unique root (a node of indegree-0 and outdegree-2) and all other nodes are either tree nodes (indegree-1 and outdegree-2) or reticulation nodes (indegree-2 and outdegree-1). For brevity we henceforth simply use the term network. A rooted binary phylogenetic tree (henceforth tree) is a phylogenetic network without any reticulation nodes. A cherry is a pair of taxa that share a common parent. A rooted binary caterpillar is a tree with exactly one cherry.

The level of a network N is the maximum number of reticulation nodes in a biconnected component of the undirected graph underpinning N. In this article we will focus exclusively on level-1 networks. In level-1 networks, maximal biconnected components that are not single edges are simple cycles that contain exactly one reticulation node; such biconnected components are called galls. An arc whose tail (but not head) is a node of a gall is called an outgoing arc.

A characterf is a surjective mapping where S is a set of discrete states. When S contains two states we say that f is a binary character. Given a tree and a character f, both on X, we say that is an extension of f to T if for all . The number of mutations induced by (on T), denoted , is the number of edges such that . The parsimony score of f with respect to T, denoted , is the minimum number of mutations induced ranging over all extensions of f. Any extension that achieves this minimum is called an optimal extension. An optimal extension can be computed in polynomial time using Fitch’s algorithm (Fitch 1971), which for completeness we describe in the appendix along with some of its relevant mathematical properties. (Note that there potentially exist optimal extensions that cannot be generated by Fitch’s algorithm.)

For a network N and a tree T, both on X, we say that NdisplaysT if there exists a subtree of N such that is a subdivision of T. An equivalent definition of “displays” relies on the notion of a switching, where a switching is a subtree of N obtained by, for each reticulation node u, deleting exactly one of u’s incoming edges. N displays T if and only if there exists some switching of N and a subdivision of T such that is a subgraph of . In both definitions we say that is an image of T inside N.

The softwired parsimony score1 of a network N with respect to f is the minimum, ranging over all trees T displayed by N, of .

We now extend the above concepts to alignments. An alignment A is simply a linear ordering of characters. In this paper the linear ordering is irrelevant so we can arbitrarily impose an ordering and write without ambiguity. An alignment can naturally be represented as a matrix with |X| rows and |A| columns; we therefore use the terms characters and columns interchangeably (and, following the use of alignments in practice, we sometimes refer to the rows of the matrix as sequences). The parsimony score of a tree T with respect to A, denoted , is simply .

When extending this concept to networks, two definitions have been proposed: the parsimony score of a network with respect to an alignment A, denoted , can be defined aswhere is the set of trees displayed by the network. According to the first definition (introduced in Hein (1990)), each character can follow a different tree displayed by the network, while in the second one (introduced in Nakhleh et al. (2005)) all characters of the alignment follow the same tree. In this paper, we will adopt the latter definition, and a tree T that is the minimizer of this sum is called a most parsimonious (MP) tree displayed by N (with respect to A).

or

Note that in applied phylogenetics alignments often contain indels, encoded using a gap symbol “-”. From the parsimony perspective it is not uncommon to treat these symbols as wildcards that do not induce mutations; the taxon “does not care” what state it is assigned. (Note however that extensions are not allowed to contain gap symbols). To compute when a character maps some of its taxa to the gap symbol, we can run Fitch’s algorithm with a slight modification to the bottom-up phase: for each taxon x such that “-”, we assign the entire set of states S to x. Moreover, as the following observation shows, the use of “-” symbols does not make the problem of identifying a most parsimonious tree displayed by a network significantly harder.

Observation 1

Let A be an alignment for a set of taxa X and let N be a phylogenetic network on X. Suppose A uses the states . Let k denote the total number of gap symbols in A. In polynomial time we can construct an alignment on 2|X| taxa, which uses only states , and a network on 2|X| taxa, such that there is a polynomial-time computable bjiection g mapping trees displayed by N to trees displayed by . This bijection g has the property that, for each tree T displayed by N, . Consequently, T is a most parsimonious tree displayed by N (wrt A) if and only if g(T) is a most parsimonious tree displayed by (wrt ).

Proof

To obtain from N we split each taxon into a cherry , . If, for a given character, had state 0 (respectively, 1), we give both and the state 0 (respectively, 1). If had state “-” we give state 0 and state 1. The idea is that by encoding a gap symbol as a cherry a single mutation is unavoidably incurred (on one of the two edges leading into and ) and thus the state of the parent of the cherry in any (optimal) extension is irrelevant. The parent thus simulates the original gap symbol: the bottom-up phase of Fitch’s algorithm will always allocate the subset of states to the parent. (The bijection g, and its inverse, are trivially computable in polynomial time by splitting each taxon into a cherry, or collapsing cherries, respectively).

We defer preliminaries relating to likelihood until Sect. 4.

Let G be an undirected graph. An orientation of G is a directed graph obtained by replacing each edge of G with exactly one of the two arcs (u, v) or (v, u). Given an orientation of G, a source is a node that has only outgoing arcs, and a sink is a node that has only incoming arcs. Let msso(G) denote the maximum, ranging over all possible orientations of G, of the sum of the number of sources and sinks in . MAX-SOURCE-SINKS-ORIENTATION is the problem of computing msso(G). A cubic graph is a graph where every node has degree 3.

The proofs of the following are deferred to the appendix. These two results form the foundation of the hardness results given in the next section.

Lemma 1

MAX-SOURCE-SINKS-ORIENTATION is NP-hard on cubic graphs.

Corollary 1

MAX-SOURCE-SINKS-ORIENTATION is APX-hard on cubic graphs.

Hardness of finding a most parsimonious tree displayed by a network

In this section we will build on Lemma 1 and Corollary 1 to prove that computing a most parsimonious tree displayed by a rooted phylogenetic network N with respect to an alignment A is NP-hard and APX-hard already for highly restricted instances.

Theorem 1

It is NP-hard to compute a most parsimonious tree displayed by a rooted phylogenetic network N with respect to an alignment A, even when N is a binary level-1 network with at most 3 outgoing arcs per gall and A consists only of two states and does not contain gap symbols.

Let be a cubic instance of MAX-SOURCE-SINKS-ORIENTATION. We will start by building a binary level-1 network N with 6|E| taxa and 2|E| reticulations, and an alignment A on states consisting of 6|E| sequences, each sequence of length |V|. (We will remove the “-” symbols later). One can thus view A as a matrix with 6|E| rows and |V| columns, or equivalently as a set of |V| characters for the 6|E| taxa of N.

To construct N, we start by taking a rooted binary caterpillar on |E| taxa. For each replace the taxon of the caterpillar with a copy of the network shown in Fig. 1. The 6 taxa within are denoted , . We use to refer to the root of .

Fig. 1

Although each sequence has length |V|, only columns u and v are shown. For this edge, the other symbols are “-”

To construct the alignment, we write to refer to the sequences, and write to refer to the state in its vth column. These states are assigned as follows. For each edge , we set the states of the 6 taxa () to be 1, 1, 0, 0, 0, 1, the states of the 6 taxa () to be 0, 1, 1, 1, 0, 0, and for each , we set the states of the 6 taxa () to all be “-”. Given that each edge is incident to exactly 3 edges, there are exactly “-” symbols in A.

Given that each contains 2 reticulations, there are different switchings of these reticulations possible, shown in Fig. 2. Note that switchings 1 and 3 both induce 5 mutations, while switchings 2 and 4 both induce 4 mutations. (Here by “induce mutations” we are referring to properties (i) and (ii) of Fitch’s algorithm, described in the appendix). We now claim that there exists an optimum solution in which only switchings 2 and 4 are used. Suppose, for some , switching 1 or 3 is used. Let T be the tree induced by this switching. Fix any optimal extension of A to T. Let be the subtree of T rooted at ; at least 5 mutations will be incurred on the edges of (with respect to the extension; see property (i) of Fitch’s algorithm). Consider now the states allocated to in columns u and v. There are four such uv combinations: 00, 01, 10, 11. If it is combination 01 or 10, we could replace with the subtree corresponding to switching 2 or 4 (respectively). This replacement subtree incurs only 4 mutations on its edges, so the total number of mutations in T decreases. If it is combination 00 or 11 we can use switching 2. This might induce a new mutation (on the edge incoming to ) but we again save at least one mutation on the edges of the subtree (because at most 4, rather than at least 5 mutations are incurred there), so the overall number of mutations does not increase. Summarizing, whichever combination 00, 01, 10, 11 occurs at , we can replace it with switching 2 or 4 without increasing the total number of mutations. Iterating this procedure proves the claim. Henceforth we can thus assume that for each either switching 2 or 4 is used.

Fig. 2

The four switchings possible for . The interior nodes are labelled by the output of the bottom-up phase of Fitch’s algorithm, for the two characters concerned. The symbol denotes where union events occur (i.e. mutations are incurred). The critical point is that both switching 2 and 4 incur the fewest number of mutations, and these select for 01 and 10 at the root, respectively, representing the choice of which way to orient edge e

Observe that if, for a given , the network uses switching 2, the bottom-up phase of Fitch’s algorithm will allocate 01 (in columns u and v) to . If, on the other hand, switching 4 is used, Fitch’s algorithm will allocate 10. In both cases, exactly 4 union events are generated on the nodes (of the subtree of induced by the switching). See Fig. 2 for elucidation.

The central idea is that, since, for an edge , a state 0 (resp. 1) in v implies a state 1 (resp. 0) in u and vice versa, we can use the choice of whether to use switching 2 or 4 (for each of the |E| reticulation pairs) to encode a choice as to which way to orient the corresponding edge. Without loss of generality we use state 0 to denote incoming edges, and state 1 to denote outgoing edges. Consider the bottom-up phase of Fitch’s algorithm. Observe that, if a vertex v incident to three edges becomes a sink, the states at the roots of (in column v) will all be 0, and for each the states at the root of (in column v) will be “-” i.e. “don’t care”2. Continuing Fitch’s algorithm along the backbone of the caterpillar shows that no mutations will be incurred on the edges of the caterpillar in column v. A completely symmetrical situation holds if a vertex becomes a source: the states at the roots of (in column v) will all be 1, and again no mutations are incurred on the edges of the caterpillar. On the other hand, if a vertex v is neither a source nor a sink, then the states assigned by the bottom-up phase of Fitch’s algorithm to the roots of (in column v) will consist of 0 (twice) and 1 (once) or 1 (twice) and 0 (once). Either way exactly 1 mutation is then incurred on the edges of the caterpillar (as can be observed by running the top-down phase of Fitch’s algorithm).

This means that the parsimony score is minimized by creating as many sources and sinks as possible. Specifically we haveEach edge in the graph will induce 4 mutations (within the part), and , which explains the term 6|V|. As argued above, sources and sinks to do not increase the parsimony score, and all other vertices increase the parsimony score by exactly 1, hence the term .

Clearly msso(G) can easily be calculated from . Finally, we can apply Observation 1 to obtain a network and without “-” symbols such thatThe transformation does not raise the level of the network or the number of arcs outgoing from any biconnected component. NP-hardness follows.

If we do allow “-” symbols then the following slightly stronger result is obtained: APX-hardness implies NP-hardness but additionally excludes the existence of a Polynomial Time Approximation Scheme (PTAS), unless P=NP. APX-hardness does not obviously hold if we encode the gap symbols using Observation 1 because the additive O(|V||E|) term thus created distorts the objective function.

Corollary 2

It is APX-hard to compute a most parsimonious tree displayed by a rooted phylogenetic network N with respect to an alignment A, even when N is a binary level-1 network with at most 3 outgoing arcs per gall and A consists only of states .

We give a (14, 1) L-reduction from msso, which is APX-hard, to the parsimony problem. L-reductions preserve APX-hardness so the result will follow. An L-reduction (Papadimitriou and Yannakakis 1991), where , is defined as follows.

Definition 1

Let A, B be two optimization problems and and their respective cost functions. A pair of functions f, g, both computable in polynomial time, constitute an L-reduction from A to B if the following conditions are true:where is the optimal solution value of problem A and similarly for B.

For every instance x of A, f(x) is an instance of B,

For every feasible solution y of f(x), g(y) is a feasible solution of x,

For every instance x of A, ,

For every feasible solution of f(x) we have

For brevity we refer to the optimum size of the parsimony problem as mp(N, A). We use the reduction described in the proof of Theorem 1 (before the gap symbols have been removed) with some slight modifications. The forward-mapping function f (condition 1 of the L-reduction) is the same mapping used in the proof of Theorem 1. The back-mapping function g (i.e. condition 2) will be described below. To establish condition 3 for a given we need to prove that . Now, we know that (see appendix) and that (because every cubic graph has a cut at least this large simply by moving nodes which have more neighbours on their side of the cut, to the other side). Hence, . We know that . Trivially therefore . Hence taking is sufficient. For the other direction, we need to show that for an arbitrary solution to the parsimony problem, which induces p mutations, the back-mapping function yields an orientation of G with s sources and sinks such that . The back-mapping function g first ensures that all the gadgets are using type 2 or type 4 switchings, which might reduce the number of mutations to , and then extracts an orientation of G (thus establishing condition 2). Now, and soSo taking is sufficient to establish condition 4.

Hardness of finding a most likely tree displayed by a network

The likelihood of a tree. We now introduce the basic concepts and notation that are necessary to define the likelihood of a tree with respect to an alignment. We largely follow the simple formulation by Roch (2006). First, we need a probabilistic model describing how sequences evolve along a tree. Here we assume the simplest model available, known as the Cavender-Farris model (Farris 1973; Cavender 1978), which can be described as follows. Let be a rooted binary phylogenetic tree on X. We associate probabilities to the edges of T and denote this . Under the Cavender-Farris model, each character evolves independently, as follows: at the root pick randomly a state between 0 and 1, each with probability , and then, for each vertex v below the root, either copy the state of the parent of v or flip it, with probabilities and , respectively. The restriction corresponds to the fact that, in a symmetric model, no amount of time can make a character more likely to change state than to remain in the same state.

The process described above eventually associates a state to each element of X at the leaves of the tree, that is, it generates a random binary character. The probability of generating the binary character f is called the likelihood of with respect to f, denoted , and can be calculated as follows:Here ranges over all extensions of f to T. Because the model assumes that the characters in a sequence evolve independently, the probability of generating the binary sequences in an alignment A, named the likelihood of with respect to A, denoted , can be obtained as(Here, and in the rest of this section, we assume that alignments do not contain gap symbols.)

We now introduce some more notation that will be useful in the following. An extension of an alignment A to a tree is a set of functions obtained by taking exactly one extension of each character in A. In practice, can be represented as a matrix with |V| rows and |A| columns, in which the rows corresponding to the leaves of T are identical to the rows of A. For , we denote by the number of differences (that is, the Hamming distance) between the sequences that associates to u and v. Finally, let denote . Note that the parsimony score is the minimum of over all extensions of A. Given these notations, we can express the likelihood of as follows, where , and ranges over all extensions of A:Each term in the sum in Eqn. (1) expresses the product of the probabilities of transition between the sequences associated by to the endpoints of the edges, times the probability of the sequence at the root.

Networks with edge probabilities. It is possible to extend the Cavender-Farris model to describe the evolution of binary sequences on a binary rooted network N. For every edge e of N, we define a probability which represents again the probability of change between 0 and 1 along edge e. The evolution of a single character follows the same rules as in the case of a tree, except that when setting the state at a reticulation vertex v, one of its two parents is randomly selected with a probability (and for the other parent), also given as a parameter of the model. The state at v is generated as if the selected parent of v were the only parent of v, as in the tree case. That includes taking into account the probability of change along the edge connecting the selected parent to v. The inheritance probability parameters (e.g. Yu et al. 2012; Wen et al. 2016; Zhang et al. 2017), and mechanisms to correlate inheritance between neighboring characters in a sequence will not be discussed in the remainder of this paper. These aspects of the model are necessary to define the likelihood of the network, but they are irrelevant for the likelihood of the trees displayed by the network, which is all that concerns us here. In the following, we denote a network N and the probabilities of change along its edges as .

We note that in some cases, the edges entering a reticulation node (the reticulation edges) may represent an event of instantaneous combination between the sequences at the tails of the reticulation edges. The probabilities of change for these reticulation edges will necessarily equal 0, as they represent an immediate transfer of genetic information, and there is no time for sequence changes along these edges. The edges entering node w in the network of Fig. 3 are an example of this. Not all edges entering a reticulation, however, need be of this type. For the sake of generality, the networks we consider here may have “non-istantaneous” reticulation edges, that is reticulation edges with . For example, consider the edges entering node v in the network of Fig. 3. Edges like these simply mean that before the reticulate event happened, the sequence in the edge leading to the reticulation evolved independently of the rest of the tree, potentially accumulating changes. At the end of the edge, the sequence that underwent the reticulate event left no descendant leading to a leaf, other than the sequence at the reticulation itself. Figure 3 also displays a different representation of the same network, showing separately the nodes and that underwent the reticulation event. Neither of these nodes left any other descendant than v within the network. In some biological contexts (for example when reticulations represent homologous recombinations), reticulation edges representing lineages that have existed for a strictly positive amount of time are the norm, and not the exception. Examples of this are the phylogenetic networks generated by the coalescent with recombination model (e.g. Griffiths and Marjoram 1997; Nordborg 2001), or by the birth-hybridization process (Zhang et al. 2017) where reticulate edges of zero length are practically impossible.

Fig. 3

Alternative representations of a phylogenetic network having some reticulation edges with strictly positive lengths. A reticulation edge with positive length should be interpreted as ending in a node that undergoes a reticulation event, but leaves no descendant in the network, other than the reticulation node itself. Both edges entering node v in the network on the left are an example of this. The representation on the right, which strictly speaking is not a phylogenetic network, makes the biological interpretation of these edges explicit. In this representation, dashed edges denote an instantaneous event, and their length is necessarily 0 (not shown)

Trees displayed by a network with edge probabilities. We say that a network displays a tree , if N displays T in the usual topological sense (i.e. some subdivision of T is a subtree of N) and can be obtained from by repeatedly suppressing vertices with indegree-1 and outdegree-1, where here suppression also updates the probabilities. Specifically, if contains two edges and , where v has indegree-1 and outdegree-1, the suppression operation replaces these two edges with a single edge and assigns it the probabilityThis expresses the probability of having different states at the endpoints of a two-edge path, under the Cavender-Farris model.

It is important to observe that, as a result of the definitions above, if a network describes the evolution of a set of characters, then any of these characters taken separately evolves according to the Cavender-Farris model for one of the trees displayed by . Figure 4 shows, as an example, the four trees displayed by the network in Fig. 3.

Fig. 4

The four trees with edge probabilities displayed by the network in Fig. 3. Note that the edge probability 0.356 in the fourth network is obtained by two applications of Eqn. (2)

Note that, in general, a tree T can have multiple distinct images in the network, so it can occur that displays for multiple different . Also note that because for all edges in the network, the same will hold for the edges of the trees it displays, as no repeated application of Eqn. (2) can produce a probability from edge probabilities that are at most . It is also easy to see that if , then . These observations lead to the following one, which will be useful later on:

Observation 2

Let be such that for every edge of N, . Let be a tree displayed by and e an edge of T. Finally, let be the subset of the edges of N whose probabilities contribute to . Then, and

We say that is a most likely (ML) tree displayed by(with respect toA) if it maximizes , ranging over all displayed by . In the remainder of this section we consider the problem of finding such a most likely tree given a network with edge probabilities and an alignment.

A link between likelihood and parsimony. There are well-known relationships between the likelihood and the parsimony of a tree that imply that under some conditions a most likely tree is also a most parsimonious one (Tuffley and Steel 1997). We now illustrate one such relationship (Corollary 3 below), which is based on the observation that as we reduce the scale of a tree, its likelihood converges to zero at a rate that only depends on its parsimony score. Although it shares similarities with the results by Tuffley and Steel, we are not aware that it has been explicitly stated in the literature. This result is not necessary to obtain the other results in this section, but it provides the intuition behind them.

In the following statements, we assume that , so the form is to be understood as c approaches 0 to the right. Also, simply denotes the product between the scalar c and vector .

Lemma 2

The function is as .

Write using Eqn. (1):where we have used . Note that the products in the second expression tend to a constant as . As a consequence, the term for in the sum has order as . Since the lowest degree dominates, their sum is .

Corollary 3

Let A be an alignment and and two trees such that . Then, for any and ,

As , is , while is . That is, converges to 0 at a lower rate than . Thus there exists a neighborhood of 0 in which holds.

The corollary above can be extended to any collection of trees: irrespective of the edge probabilities assigned to them, if the trees are rescaled by a sufficiently small c, the most parsimonious trees will have likelihoods greater than all the other trees, meaning that a most likely tree in the collection of rescaled trees will necessarily also be most parsimonious.

Proving the NP-hardness of finding an ML tree in a network. In the remainder of this section, namely in the statements of the next two formal results, we are implicitly given a network N on X with and an alignment A with m characters on X. The height of a network N is the maximum number of edges in a directed path in N.

Lemma 3

Let be a network of height , where all the edges are assigned a constant probability , with . Let be a tree displayed by . Then,

Using Observation 2, we note that, for any edge of , .

We begin by proving the upper bound in the statement. From Eqn. (1) and the fact that , we get the first inequality in the following:where the last inequality is obtained by noting that the sum has terms (there are internal nodes in a rooted binary tree, and thus different extensions of A), and that (there are branches in a rooted binary tree, and thus we cannot have more than changes per character). The upper bound in the statement is larger than the one above.

As for the lower bound, if we use in Eqn. (1):where the last inequality is obtained by taking only one term in the sum. The lower bound in the statement is smaller than the one above.

The lemma above shows the order of convergence to 0 of the likelihood of a tree displayed by as . The higher the parsimony score, the faster the convergence. As a consequence, for c sufficiently close to 0, a tree with a lower parsimony score than another will have a higher likelihood. The following lemma shows how close is “sufficiently close”, by providing an explicit upper bound to c.

Proposition 1

Let be a network of height , where all the edges are assigned a constant probability , with . If is a most likely tree displayed by , then is a most parsimonious tree displayed by N.

Suppose that is a most likely tree displayed by , but not most parsimonious. That is, there exists displayed by with . But then, by using the lower bound in Lemma 3:Now apply the upper bound in Lemma 3 to , and combine it with :The last terms of the two chains of inequalities above are equal, thus proving . Since this contradicts the assumption that is a most likely tree, the statement follows.

The proposition above shows that the NP-hard problem of finding a most parsimonious tree in a network N with respect to an alignment A can be reduced to the problem of finding a most likely tree in with respect to A, where is such that , and . Since the reduction preserves the network and the alignment, the main result of this section follows from Theorem 1:

Theorem 2

It is NP-hard to compute a most likely tree displayed by a rooted phylogenetic network with respect to an alignment A, even when N is a binary level-1 network with at most 3 outgoing arcs per gall and A consists only of two states and does not contain gap symbols.

Conclusions and open problems

We have shown that, given a phylogenetic network with a sequence for each leaf, finding a most parsimonious or most likely tree displayed by the network is computationally intractable (NP-hard). Moreover, this is the case even when we restrict to binary sequences and level-1 networks; the simplest networks that are not trees. However, many computational problems that can be shown to be theoretically intractable can be solved reasonably efficiently in practice (see e.g. Cautionary Tales of Inapproximability by Budden and Jones (2017)). We end the paper by discussing whether we expect this to be the case for our problem.

There is a dynamic programming algorithm, described in Theorem 5.7 of Fischer et al. (2015), for finding a tree in a network that is most parsimonious with respect to a single character. The running time is fixed-parameter tractable, with as parameter the level of the network. Hence, this algorithm is practical as long as the level of the network is not too large. This algorithm can easily be extended to multiple characters (that all have to choose the same tree) when the number of characters is adopted as a second parameter. Indeed, for every root of a biconnected component, we introduce a dynamic programming entry not just for every possible state but for every possible sequence of states. However, the running time of this algorithm would be exponential in the number of characters, which makes it useless for almost all biological data. Similarly, the Integer Linear Programming (ILP) solution presented in the same paper can also be easily extended to multiple characters. However, there does not seem to be an easy way to do that without having the number of variables growing linearly in the number of characters. Hence, this approach is also unlikely to be useful in practice.

In contrast, consider the simple algorithm that loops through the at most trees displayed by the network, with r the number of reticulation nodes in the network, and computes the parsimony or likelihood of each tree (this naïve FPT algorithm was presented in Nakhleh et al. (2005), where it is named Net2Trees). Ironically, this simple algorithm would outperform the approaches mentioned above for any kind of data with a reasonably large number of characters. Hence, the main open question that remains is whether there exists an algorithm whose running time is linear (or at least polynomial) in the number of characters and whose dependency on r is better than  (for example recently an algorithm with exponential base smaller than 2 was discovered for the tree containment problem (Gunawan et al. 2016), although this algorithm does not obviously extend to generating all trees in the network). Another question of interest that remains open is the following: does the parsimony problem under restrictions (1)-(3) listed in the introduction permit good (i.e. constant factor) approximation algorithms, and possibly even a PTAS, when the alignment A does not contain any indels?