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.

nodes. Recent years have seen an explosion of interest in constructing rooted phyloge-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 network N 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 character f is a surjective mapping f : X → S where S is a set of discrete states. When S contains two states we say that f is a binary character. Given a tree T = (V , E) and a character f , both on X , we say thatf : V → S is an extension of f to T iff (x) = f (x) for all x ∈ X . The number of mutations induced byf (on T ), denoted lf (T ), is the number of edges {u, v} ∈ E such thatf (u) =f (v). The parsimony score of f with respect to T , denoted l f (T ), is the minimum number of mutations induced ranging over all extensionsf 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 N displays T if there exists a subtree T of N such that T is a subdivision of T . An equivalent definition of "displays" relies on the notion of a switching, where a switching is a subtree N 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 N of N and a subdivision T of T such that T is a subgraph of N . In both definitions we say that T is an image of T inside N .
The softwired parsimony score 1 of a network N with respect to f is the minimum, ranging over all trees T displayed by N , of l f (T ).
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 f ∈ A 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 l A (T ), is simply f ∈A l f (T ).
When extending this concept to networks, two definitions have been proposed: the parsimony score of a network with respect to an alignment A, denoted l A (N ), can be defined as 1.
where T (N ) 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). 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 l f (T ) when a character f : X → S maps some of its taxa to the gap symbol, we can run Fitch's algorithm with a slight modification to the bottomup phase: for each taxon x such that f (x) ="-", 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 {0, 1, "-"}. Let k denote the total number of gap symbols in A. In polynomial time we can construct an alignment A on 2|X | taxa, which uses only states {0, 1}, and a network N on 2|X | taxa, such that there is a polynomial-time computable bjiection g mapping trees displayed by N to trees displayed by N . This bijection g has the property that, for each tree T displayed by N , l A (g(T )) = l A (T ) + k. 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 N (wrt A ).
Proof To obtain N from N we split each taxon x i into a cherry {x 1 i , x 2 i }. If, for a given character, x i had state 0 (respectively, 1), we give both x 1 i and x 2 i the state 0 (respectively, 1). If x i had state "-" we give x 1 i state 0 and x 2 i state 1. The idea is that by encoding a gap symbol as a {0, 1} cherry a single mutation is unavoidably incurred (on one of the two edges leading into x 1 i and x 2 i ) 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 {0, 1} 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 G obtained by replacing each edge {u, v} of G with exactly one of the two arcs (u, v) or (v, u). Given an orientation G 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 G of G, of the sum of the number of sources and sinks in G . 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.

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 {0, 1} and does not contain gap symbols. 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 {0, 1, "-"} consisting of 6|E| sequences, each sequence of length |V |. (We will remove the "-" symbols later). One can thus view A as a {0, 1, −} 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 e ∈ E replace the taxon x e of the caterpillar with a copy N e of the network shown in Fig. 1. The 6 taxa within N e are denoted x e,i , i ∈ {1, . . . , 6}. We use r e to refer to the root of N e .
Given that each N e contains 2 reticulations, there are 2 2 = 4 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 e = {u, v} ∈ E, switching 1 or 3 is used. Let T be the tree induced by this switching. Fix any optimal extension of A to T . Let T e be the subtree of T rooted at r e ; at least 5 mutations will be incurred on the edges of T e (with respect to the extension; see property (i) of Fitch's algorithm). Consider now the states allocated to r e 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 T e 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 r e ) 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 Switching 4 4 mutations

Fig. 2
The four switchings possible for N e . 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 overall number of mutations does not increase. Summarizing, whichever combination 00, 01, 10, 11 occurs at r e , 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 e ∈ E either switching 2 or 4 is used.
Observe that if, for a given e = {u, v}, the network N e uses switching 2, the bottom-up phase of Fitch's algorithm will allocate 01 (in columns u and v) to r e . 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 N e induced by the switching). See Fig. 2 for elucidation.
The central idea is that, since, for an edge e = {u, v}, 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 e 1 , e 2 , e 3 becomes a sink, the states at the roots of N e 1 , N e 2 , N e 3 (in column v) will all be 0, and for each e / ∈ {e 1 , e 2 , e 3 } the states at the root of N e (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 N e 1 , N e 2 , N e 3 (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 N e 1 , N e 2 , N e 3 (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 have Each edge in the graph will induce 4 mutations (within the N e part), and |E| = 3|V |/2, 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 (|V | − msso(G)).
Clearly msso(G) can easily be calculated from l A (N ). Finally, we can apply Observation 1 to obtain a network N and A without "-" symbols such that The 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: APXhardness 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
Proof 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 (α, β) Lreduction (Papadimitriou and Yannakakis 1991), where α, β ≥ 0, is defined as follows.
Definition 1 Let A, B be two optimization problems and c A and c B 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 O PT A is the optimal solution value of problem A and similarly for B.
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 mp(N , A) ≤ α · msso(G). Now, we know that msso(G) = maxcut(G)−|V |/2 (see appendix) and that maxcut(G) ≥ 2/3|E| = |V | (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, Hence taking α = 14 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 |msso(G) − s| ≤ β| p − mp(N , A)|. The back-mapping function g first ensures that all the N e gadgets are using type 2 or type 4 switchings, which might reduce the number of mutations to p ≤ p, and then extracts an orientation of G (thus establishing condition 2). Now, s = 7|V | − p and mp(N , So taking β = 1 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 T = (V , E) be a rooted binary phylogenetic tree on X . We associate probabilities p = ( p e ) e∈E ∈ [0, 1 /2] |E| to the edges of T and denote this (T , p). 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 1 /2, and then, for each vertex v below the root, either copy the state of the parent of v or flip it, with probabilities 1 − p e and p e , respectively. The restriction p e ≤ 1 /2 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 (T , p) with respect to f , denoted L f (T , p), and can be calculated as follows: Heref 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 (T , p) with respect to A, denoted L A (T , p), 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 A of an alignment A to a tree T = (V , E) is a set of functionsf : V → {0, 1} obtained by taking exactly one extension of each character in A. In practice, A 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 e = (u, v) ∈ E, we denote by h e ( A) the number of differences (that is, the Hamming distance) between the sequences that A associates to u and v. Finally, let l A (T ) denote e∈E h e ( A) = f ∈ A lf (T ). Note that the parsimony score l A (T ) is the minimum of l A (T ) over all extensions of A. Given these notations, we can express the likelihood of (T , p) as follows, where m = |A| = | A|, and A 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 A 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 p e ∈ [0, 1 /2] 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 γ v ∈ (0, 1) (and 1 − γ v 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  . 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) includes taking into account the probability of change p e along the edge connecting the selected parent to v. The inheritance probability parameters γ v (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 (N , p ). 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 p e 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 p e > 0. 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 v 1 and v 2 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 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) process (Zhang et al. 2017) where reticulate edges of zero length are practically impossible.
Trees displayed by a network with edge probabilities. We say that a network (N , p ) displays a tree (T , p), if N displays T in the usual topological sense (i.e. some subdivision T of T is a subtree of N ) and (T , p) can be obtained from T by repeatedly suppressing vertices with indegree-1 and outdegree-1, where here suppression also updates the probabilities. Specifically, if T contains two edges e 1 = (u, v) and e 2 = (v, w), where v has indegree-1 and outdegree-1, the suppression operation replaces these two edges with a single edge e = (u, w) and assigns it the probability p e = p e 1 (1 − p e 2 ) + (1 − p e 1 ) p e 2 . (2) This 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 (N , p ) 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 (T , p) displayed by (N , p ). Figure 4 shows, as an example, the four trees displayed by the network in Fig. 3.
Note that, in general, a tree T can have multiple distinct images T in the network, so it can occur that (N , p ) displays (T , p) for multiple different p. Also note that because p e ≤ 1 /2 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 p e > 1 /2 from edge probabilities that are at most 1 /2. It is also easy to see that if 0 < p e 1 , p e 2 < 1 /2, then max{ p e 1 , p e 2 } < p e < p e 1 + p e 2 . These observations lead to the following one, which will be useful later on: Observation 2 Let (N , p ) be such that for every edge of N , 0 < p e < 1 /2. Let (T , p) be a tree displayed by (N , p ) and e an edge of T . Finally, let E (e) be the subset of the edges of N whose probabilities contribute to p e . Then, p e < 1 /2 and max e ∈E (e) We say that (T * , p * ) is a most likely (ML) tree displayed by (N , p ) (with respect to A) if it maximizes L A (T , p), ranging over all (T , p) displayed by (N , p ). 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 c ∈ (0, 1), so the form c → 0 is to be understood as c approaches 0 to the right. Also, cp simply denotes the product between the scalar c and vector p.

Lemma 2 The function f (c) = L A (T , cp) is (c l A (T ) ) as c → 0.
Proof Write L A (T , cp) using Eqn. (1): where we have used e∈E h e ( A) = l A (T ). Note that the products in the second expression tend to a constant as c → 0. As a consequence, the term for A in the sum has order (c l A (T ) ) as c → 0. Since the lowest degree dominates, their sum is (c l A (T ) ).

Corollary 3 Let A be an alignment and T 1 and T 2 two trees such that l A (T 1 ) < l A (T 2 ).
Then, for any p 1 and p 2 , 2 ) for c sufficiently close to 0. cp 1 ) converges to 0 at a lower rate than L A (T 2 , cp 2 ). Thus there exists a neighborhood of 0 in which L A (T 1 , cp 1 ) > L A (T 2 , cp 2 ) 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 |X | = n 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 . (N , c) be a network of height d N , where all the edges are assigned a constant probability c e = c, with 0 < c < 1 /2. Let (T , p) be a tree displayed by (N , c). Then,

Lemma 3 Let
Proof Using Observation 2, we note that, for any edge e ∈ E of T = (V , E), c < p e < min{cd N , 1 /2}.
We begin by proving the upper bound in the statement. From Eqn.
(1) and the fact that (1 − p e ) m−h e ( A) < 1, we get the first inequality in the following: where the last inequality is obtained by noting that the sum has 2 m(n−1) terms (there are n − 1 internal nodes in a rooted binary tree, and thus 2 m(n−1) different extensions of A), and that l A (T ) ≤ l A (T ) ≤ m(2n − 2) (there are 2n − 2 branches in a rooted binary tree, and thus we cannot have more than 2n − 2 changes per character). The upper bound in the statement is larger than the one above.
As for the lower bound, if we use c < p e < 1 /2 in Eqn. (1): The lemma above shows the order of convergence to 0 of the likelihood L A (T , p) of a tree displayed by (N , c) as c → 0. 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. (N , c) be a network of height d N , where all the edges are assigned a constant probability c e = c, with 0 < c < d −2mn N 2 −3mn . If (T * , p * ) is a most likely tree displayed by (N , c), then T * is a most parsimonious tree displayed by N .

Proposition 1 Let
Proof Suppose that (T * , p * ) is a most likely tree displayed by (N , c), but not most parsimonious. That is, there exists (T , p) displayed by (N , c) with l A (T ) ≤ l A (T * )−1. But then, by using the lower bound in Lemma 3: Now apply the upper bound in Lemma 3 to T * , and combine it with c < d −2mn N 2 −3mn : The last terms of the two chains of inequalities above are equal, thus proving L A (T , p) > L A (T * , p * ). Since this contradicts the assumption that (T * , p * ) 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 (N , c) with respect to A, where c = (c e ) is such that c e = c, and 0 < c < d −2mn N 2 −3mn . 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 (T , p) displayed by a rooted phylogenetic network (N , p ) 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 {0, 1} 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 2 r 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 2 r (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?
A Appendix: Fitch's algorithm Fitch's algorithm (Fitch 1971) has two phases. In the first phase, known as the bottomup phase, we start by assigning the singleton subset of states { f (x)} to each taxon x. The internal nodes of T are assigned subsets of states recursively, as follows. Suppose a node p has two children u and v, and the bottom-up phase has already assigned subsets F(u) and F(v) to the two children, respectively. If F(u) ∩ F(v) = ∅ then set F( p) = F(u) ∩ F(v) (in which case we say that p is an intersection node). If F(u) ∩ F(v) = ∅ then set F( p) = F(u) ∪ F(v) (in which case we say that p is a union node). The number of union nodes in the bottom-up phase is equal to l f (T ). To actually create an optimal extensionf , we require the top-down phase of Fitch's algorithm. Start at the root r and letf (r ) be any element in F(r ). For an internal node u with parent p, we setf (u) =f ( p) (iff ( p) ∈ F(u)) and otherwise (i.e. f ( p) / ∈ F(u)) setf (u) to be an arbitrary element of F(u). For each node u of the tree, let ∪(u) be the number of union events in the subtree rooted at u. The following well-known properties of Fitch's algorithm are used repeatedly in the main hardness proof of this article: (i) every extension (optimal or otherwise) must incur at least ∪(u) mutations on the edges of the subtree rooted at u; (ii) an extension created by Fitch's algorithm induces exactly ∪(u) mutations on the edges of the subtree rooted at u (and u is assigned a state from F(u) in this extension).

B Appendix: NP-hardness and APX-hardness of MAX-SOURCE-SINKS-ORIENTATION
The following result is based on a sketch proof by Colin McQuillan 3 . We have been unable to find an original reference and hence have reconstructed the proof in detail. The APX-hardness proof is original.

Lemma 1. MAX-SOURCE-SINKS-ORIENTATION is NP-hard on cubic graphs.
Proof Recall that the classical MAX-CUT problem asks us to bipartition the vertex set of an undirected graph G, such that the number of edges that cross the bipartition is maximized. We reduce from the NP-hard problem CUBIC-MAX-CUT which is the restriction of the MAX-CUT problem to cubic graphs. Given an undirected cubic graph G, we simply write maxcut(G) to denote the number of edges in the maximum-size cut.
We reduce CUBIC-MAX-CUT to MAX-SOURCE-SINKS-ORIENTATION. Specifically, given an undirected cubic graph G = (V , E) (i.e. an instance of CUBIC-MAX-CUT) we will show that msso(G) = maxcut(G) − |V |/2, from which the hardness will follow.
We start by proving that msso(G) ≥ maxcut(G) − |V |/2. Fix an arbitrary cut C of G and let (U , W ) be the corresponding bipartition. If some vertex of U or W has more neighbours on the other side of the partition than its own, move it to the other side of the partition: this will increase the size of the cut. We repeat this until it is no longer possible and let C and (U , W ) refer to the cut and its induced partition at the end of this process. Note that now each vertex in U (respectively, W ) will have at most one neighbour in U (respectively, W ). We proceed by orienting the edges in the cut from U to W . Now, the remaining edges are either internal to U or internal to W . These edges must form a matching (i.e. they are node disjoint). For each such edge in U (respectively, W ), exactly one endpoint will become a source (respectively, sink). Nodes in U (respectively, W ) that are not adjacent to internal edges will already be sources (respectively, sinks) due to the orientation of the cut edges from U to W . Hence, if we write |C| to denote the number of edges in the cut C, we obtain an orientation of G with at least (|E| − |C|) + (|V | − 2(|E| − |C|)) sources and sinks. Hence, For the other direction, fix an arbitrary orientation of G and let s be the number of sources and sinks created by the orientation. We write V i (i ∈ {0, 1, 2, 3}) to denote those vertices of G which have indegree i. Let U = V 0 ∪ V 1 and let W = V 2 ∪ V 3 . Whenever an edge of G has been oriented from W to U , reverse its orientation: this only decreases the indegrees of the vertices in U and increases the indegrees of vertices in W so it cannot destroy any sources or sinks and it cannot cause a node to be on the "wrong" side of the bipartition. (In fact, it will cause the number of sources and sinks to increase, so this situation can only occur if the orientation was not optimal). Let s now refer to the number of sources and sinks once all arcs have been oriented from U to W . The edges (u, w) such that u ∈ U and w ∈ W form a cut; it remains only to count how many of these edges there are. We first count from the perspective of the vertices in U . The nodes in V 0 each generate 3 outgoing arcs. Let k be the total number of edges of the form (u 0 , u 1 ) where u 0 ∈ V 0 and u 1 ∈ V 1 . Note that each node in V 1 that does not receive any of these k arcs, must receive an arc which is outgoing from some other node in V 1 . It follows that the number of edges in the cut is (3|V 0 | − k) + (2|V 1 | − (|V 1 | − k)) = 3|V 0 | + |V 1 |.
If we count in a symmetrical fashion from the perspective of W , and let be the number of arcs of the form (u 2 , u 3 ) where u 2 ∈ V 2 and u 3 ∈ V 3 , it follows that the number of edges in the cut is (3|V 3 | − ) + (2|V 2 | − (|V 2 | − )) = 3|V 3 | + |V 2 |.
The hardness of msso can be strengthened to the following inapproximability result. Note that one consequence of APX-hardness is that msso does not permit a PTAS, unless P = N P.

Corollary 1 MAX-SOURCE-SINKS-ORIENTATION is APX-hard on cubic graphs.
Proof Note that the constructions and transformations used in the proof of Lemma 1 are all constructive and can easily be conducted in polynomial time. Moreover, they apply to arbitrary cuts/orientations, and not just optimal ones. This allows us to easily strengthen the described reduction to obtain a (1, 1) L-reduction from CUBIC-MAX-CUT to MAX-SOURCE-SINKS-ORIENTATION (see the main text for the definition of L-reduction). From this APX-hardness will follow, since CUBIC-MAX-CUT is APX-hard (Alimonti and Kann 1997;Berman and Karpinski 1999) and L-reductions are APX-hardness preserving. The (1, 1) means that the inapproximability threshhold for MAX-SOURCE-SINKS-ORIENTATION is at least as strong as that for CUBIC-MAX-CUT.