Which Phylogenetic Networks Are Level-$k$ Networks with Additional Arcs? Structure and Algorithms

Suppose $X$ represents a non-empty finite set of present-day species. A rooted almost-binary phylogenetic network (on a leaf-set $X$) is defined to be a finite simple directed acyclic graph $N$ with the following properties (P1)–(P3). The vertex $\rho$ is called the root of $N$, and any vertex $v$ with ${\mathrm{indeg}}_N(v)>1$ is called a reticulation of $N$.

(P1)

There exists a unique vertex $\rho$ of $N$ with ${\mathrm{indeg}}_N(\rho )=0$ and ${\mathrm{outdeg}}_N(\rho )\in \{1,2\}$;

(P2)

The set of leaves of $N$ is identical to $X$;

(P3)

For any $v\in V(N)\setminus (X\cup \{\rho \})$, ${\mathrm{indeg}}_N(v)\in \{1,2\}$ and ${\mathrm{outdeg}}_N(v)\in \{1,2\}$;

The above definition allows $N$ to contain a vertex $v$ with ${\mathrm{indeg}}_N(v)={\mathrm{outdeg}}_N(v)=2$ or ${\mathrm{indeg}}_N(v)={\mathrm{outdeg}}_N(v)=1$. In this sense, it slightly generalizes the notion of rooted binary phylogenetic networks $N$ on $X$, which is defined by (P1), (P2) and (P4). When $N$ has the properties (P1), (P2) and (P5), $N$ is particularly called a rooted binary phylogenetic tree on $X$.

(P4)

For any $v\in V(N)\setminus (X\cup \{\rho \})$, $\{{\mathrm{indeg}}_N(v),{\mathrm{outdeg}}_N(v)\}=\{1,2\}$.

(P5)

For any $v\in V(N)\setminus (X\cup \{\rho \})$, $({\mathrm{indeg}}_N(v),{\mathrm{outdeg}}_N(v))=(1,2)$.

While the concept of support trees was originally defined for rooted binary phylogenetic networks in [8], the theoretical framework developed in [9], which includes the structure theorem and associated algorithms and forms the foundation of the present work, applies equally to rooted almost-binary networks. Therefore, in this paper, following the approach taken in [15], we will consider support trees of almost-binary (including binary) networks, unless stated otherwise. For brevity, we omit “rooted almost-binary” when no confusion is likely to arise.

As a support tree $T$ is a spanning tree of $N=(V,E)$, each $T$ is specified by its edge-set. This leads to the natural question: Which subsets $S$ of $E$ yield a support tree $N[S]=(V,S)$ of $N$? Francis and Steel [8] proved that such “admissible” subsets $S$ of $E$ are characterized by the following conditions (C1)–(C3). As in [9], we slightly generalize the original definition in [8] so that we can consider admissible edge selections for any subgraph of $N$.

Theorem 3 and Remark 4 allow us to identify each support tree $T$ of $N$ with its edge-set $E(T)$. In other words, counting or listing support trees of $N$ refers to counting or listing admissible subsets of $E(N)$.

We recall the relevant materials from [9, 15]. Given a rooted almost-binary phylogenetic network $N$, a connected subgraph $Z$ of $N$ with $m:=|E(Z)|\ge 1$ is a zig-zag trail (in $N$) if the edges of $Z$ can be permuted as $(e_1,\mathrm{\dots },e_m)$ such that either $\mathrm{head}(e_i)=\mathrm{head}(e_{i+1})$ or $\mathrm{tail}(e_i)=\mathrm{tail}(e_{i+1})$ holds for each $i\in [1,m-1]$. A zig-zag trail is represented by an alternating sequence of (not necessarily distinct) vertices and distinct edges, e.g., $(v_0,(v_0,v_1),v_1,(v_2,v_1),v_2,\mathrm{\dots },(v_m,v_{m-1}),v_m)$, but this can be more concisely written as or its reverse, where each edge $(v_i,v_{i+1})$ is represented by $v_i>v_{i+1}$. A zig-zag trail $Z$ in $N$ is maximal if $N$ contains no zig-zag trail $Z^{\prime }$ such that $Z$ is a proper subgraph of $Z^{\prime }$. Each maximal zig-zag trail in $N$ falls into one of the following four types. A crown is a maximal zig-zag trail $Z$ that has even $m:=|E(Z)|\ge 4$ and can be written in the cyclic form . A fence is any maximal zig-zag trail that is not a crown. An N-fence is a fence $Z$ with odd $m:=|E(Z)|\ge 1$, expressed as . A fence $Z$ with even $m:=|E(Z)|\ge 2$ is an M-fence if expressed as , and it is a W-fence if expressed as .

As established in [9] (and will be restated in Theorem 8), any rooted almost-binary phylogenetic network $N$ admits a unique decomposition $𝒵=\{Z_1,\mathrm{\dots },Z_d\}$ into its maximal zig-zag trails ($d\ge 1$), called the zig-zag trail decomposition of $N$ (see Figure 1 for an illustration). To understand how this decomposition arises, observe that every edge $e$ of $N$ belongs to an obvious zig-zag trail $\mathrm{tail}(e)>\mathrm{head}(e)$ in $N$. If $\mathrm{tail}(e)$ has out-degree two in $N$, or $\mathrm{head}(e)$ has in-degree two in $N$, then this trail is not maximal. In such cases, by repeatedly extending the zig-zag trail, one can eventually obtain a maximal one. By construction, if each vertex of $N$ has in-degree and out-degree at most two, then a maximal zig-zag trail containing $e$ is unique for each $e$. In other words, if $N$ is almost-binary, then no two distinct maximal zig-zag trails in $N$ have a common edge. It follows that $\{E(Z_1),\mathrm{\dots },E(Z_d)\}$ is a partition of $E(N)$, and thus $𝒵$ is indeed a decomposition of $N$.

As in the approach taken in [9, 10], we often consider a maximal zig-zag trail $Z$ as a sequence $(e_1,\mathrm{\dots },e_{|E(Z)|})$ of edges, ordered according to their appearance in the trail, where $e_1$ and $e_{|E(Z)|}$ are called the terminal edges of $Z$ when $Z$ is a fence. For any maximal zig-zag trail $Z=(e_1,\mathrm{\dots },e_{|E(Z)|})$ in $N$, any subset $S$ of $E(Z)$ is specified by a $0$-$1$ sequence $⟨b_1{b}_2\mathrm{\dots }{b}_{|E(Z)|}⟩$, where $b_i=1$ if $e_i\in S$ and $b_i=0$ otherwise. For example, given a crown $Z=(e_1,e_2,e_3,e_4)$, a subset $\{e_1,e_3\}$ of $E(Z)$ can be encoded as $⟨\mathrm{1\hspace{0.33em}0\hspace{0.33em}1\hspace{0.33em}0}⟩$. We often write $⟨{(10)}^2⟩$ to avoid a repetition such as $⟨\mathrm{1\hspace{0.33em}0\hspace{0.33em}1\hspace{0.33em}0}⟩$.

For the reader’s benefit, we now illustrate how to find admissible subsets using small examples. Let $𝒵=\{Z_1,\mathrm{\dots },Z_d\}$ be the maximal zig-zag trail decomposition of $N$. Suppose $Z_1=(e_1,e_2,e_3,e_4,e_5)$ is an N-fence with $\mathrm{head}(e_4)=\mathrm{head}(e_5)$. Then, $\{e_1,e_3,e_5\}$ is the only admissible subset of $E(Z_1)$ because (C1) requires inclusion of both terminal edges $e_1$ and $e_5$, (C3) excludes selecting both $e_4$ and $e_5$, and (C2) then requires $e_3$, which in turn excludes $e_2$ due to (C3). Thus, the family $𝒮(Z_1)$ of admissible subsets of $E(Z_1)$ consists of a single element $\{e_1,e_3,e_5\}$, and so $𝒮(Z_1)=\{⟨\mathrm{1\hspace{0.33em}0\hspace{0.33em}1\hspace{0.33em}0\hspace{0.33em}1}⟩\}$. Next, suppose $Z_2=(e_1,e_2,e_3,e_4)$ is a crown. Then, $\{e_1,e_3\}$ and $\{e_2,e_4\}$ are the only two admissible subsets of $E(Z_2)$, because they are the only subsets that satisfy both (C2) and (C3), and crowns have no edge subject to (C1). Thus, $𝒮(Z_2)=\{⟨\mathrm{1\hspace{0.33em}0\hspace{0.33em}1\hspace{0.33em}0}⟩,⟨\mathrm{0\hspace{0.33em}1\hspace{0.33em}0\hspace{0.33em}1}⟩\}$.

Remark 7 provides a key idea behind the direct product decomposition in Theorem 8. After obtaining $𝒵=\{Z_1,\mathrm{\dots }{Z}_d\}$ for $N$, one can construct an admissible subset $S$ of $E(N)$ simply by computing $S=S_1\cup \mathrm{\cdots }\cup S_d$, where $S_i$ is any admissible subset of $E(Z_i)$. The family $𝒯$ of support trees of $N$ – the family of admissible subsets of $E(N)$ – by listing all such combinations $(S_1,\mathrm{\dots },S_d)$ with each $S_i\in E(Z_i)$.

Using the definitions, notation, and key ideas outlined above, we now state the structure theorem for rooted almost-binary phylogenetic networks. It provides a way to canonically decompose such networks $N$ and characterizes the family of admissible subsets $S$ of $E(N)$, namely, the family of edge-sets of support trees of $N$.

Theorem 8 and Proposition 6 have furnished optimal algorithms for various computational problems on support trees [9, 10]. For example, the number $|𝒯|$ of support trees of $N$ is given by $|𝒯|={\prod }_{i=1}^d|𝒮(Z_i)|$, and using (2), it can be computed in $\mathrm{\Theta }(|E(N)|)$ time.

By Lemma 11, one can find an $𝒜$-admissible subset of $E(N)$ in essentially the same manner as for an admissible subset of $E(N)$: first, compute the zig-zag trail decomposition $𝒵=\{Z_1,\mathrm{\dots },Z_d\}$ of $N$, and then select edges of $Z_i$ according to the conditions (C1^⋆) and (C2^⋆) in Definition 10. As explained in Remark 7, when $Z_i=(e_1,\mathrm{\dots },e_{|E(Z_i)|})$ is a fence, (C1^⋆) requires that the terminal edges $e_1$ and $e_{|E(Z_i)|}$ be marked $1$. In contrast, new condition (C2^⋆) prohibits unselected consecutive edges of $Z_i$ – consecutive $0$’s in the $0$-$1$ sequence – regardless of whether $Z_i$ is a crown or a fence.

Thus, the family ${𝒮}_{𝒜}(Z_i)$ of $𝒜$-admissible subsets of $E(Z_i)$ is expressed by (3) for each $Z_i\in 𝒵$. We note that, if $Z_i$ is a crown, different choices of $e_1$ result in different sequences representing the same subset of $E(Z_i)$. For example, both $⟨\mathrm{1\hspace{0.33em}0\hspace{0.33em}0\hspace{0.33em}1}⟩$ and $⟨\mathrm{0\hspace{0.33em}1\hspace{0.33em}1\hspace{0.33em}0}⟩$ represent the same subset that does not satisfy condition (C2^⋆).

Moreover, by Lemma 11, ${𝒜}_N$ is characterized by ${𝒜}_N={\prod }_{i=1}^d{𝒮}_{𝒜}(Z_i)$, similarly to $𝒯$ in Theorem 8. Therefore, $|{𝒜}_N|={\prod }_{i=1}^d|{𝒮}_{𝒜}(Z_i)|$. As in (3), each number $|{𝒮}_{𝒜}(Z_i)|$ depends on whether $Z_i$ is a fence or not. This number can be more explicitly represented using Fibonacci and Lucas numbers (OEIS A000045 and OEIS A000032, respectively) as follows. The proof of Theorem 12 is given in Appendix A.2.

Theorem 12 yields an obvious algorithm for counting $|{𝒜}_N|$. It first computes the maximal zig-zag trail decomposition $𝒵=\{Z_1,\mathrm{\dots },Z_d\}$ of $N$, and then calculates $|{𝒜}_N|$ using (4).

We can see that $𝒜$-admissible subset $S$ of $E(Z_i)$ is $ℬ$-admissible if and only if $S$ contains no three consecutive edges in the edge sequence $(e_1,\mathrm{\dots },e_{|E(Z_i)|})$ of $Z_i$. Thus, the family ${𝒮}_{ℬ}(Z_i)$ of $ℬ$-admissible subsets of $E(Z_i)$ is expressed as in (5) for each $Z_i\in 𝒵$.

By Lemma 11, ${ℬ}_N$ is characterized by ${ℬ}_N={\prod }_{i=1}^d{𝒮}_{ℬ}(Z_i)$. Therefore, $|{ℬ}_N|={\prod }_{i=1}^d|{𝒮}_{ℬ}(Z_i)|$. This number is expressed using Padovan numbers (OEIS A000931) and Perrin numbers (OEIS A001608) as follows. The proof of Theorem 14 is in Appendix A.3.

The construction of $𝒞$-admissible subsets of each $E(Z_i)$ is the same as that of admissible subsets for support trees, with the only difference that $𝒞$-admissible subsets allow W-fences.

By Lemma 11, ${𝒞}_N$ is characterized by ${𝒞}_N={\prod }_{i=1}^d{𝒮}_{𝒞}(Z_i)$. Therefore, $|{𝒞}_N|={\prod }_{i=1}^d|{𝒮}_{𝒞}(Z_i)|$. From (7), we have $|{𝒮}_{𝒞}(Z_i)|=2$ for any crown $Z_i$, $|{𝒮}_{𝒞}(Z_i)|=1$ for any N-fence $Z_i$, and $|{𝒮}_{𝒞}(Z_i)|=|E(Z_i)|/2$ for any other $Z_i$. Thus, we obtain Theorem 16 and Proposition 17.

As noted in Section 4.1, the elements of ${𝒞}_N$ can also be listed with $\mathrm{\Theta }(|E(N)|)$ delay.

One can solve any optimization problem on support networks by listing all elements of ${𝒜}_N$, but this approach is generally infeasible due to its enormous size. Although all of $|{𝒜}_N|$, $|{ℬ}_N|$, and $|{𝒞}_N|$ grow exponentially with the size of $N$, we present a comparison of their practical sizes to motivate our approach to Problem 2 described in Section 6. The full implementation including the code used to generate the networks, along with the generated samples and detailed count data, is available in our GitHub repository.

For this case study, we generated rooted binary phylogenetic networks with $n$ leaves and $r=2(n-1)$ reticulations. Such a network has $2n-2+3r=8(n-1)$ edges, and thus $50\%$ of its edges are reticulation edges, representing a high level of reticulation. For each $3\le n\le 10$, we created 100 such networks by the following procedure: starting from a rooted binary phylogenetic tree $N_0$, the procedure repeatedly adds a new edge between two arbitrary edges of $N_i$, directing it so as to preserve the acyclicity of $N_{i+1}$, and returns $N_r$. Note that networks generated in this manner may or may not be tree-based, since every new edge is added not between edges of the initial tree $N_0$ but between edges of the current network $N_i$.

For each generated network $N$, we computed $|{𝒜}_N|$, $|{ℬ}_N|$, and $|{𝒞}_N|$ using the linear-time algorithms described in Sections 4.1–4.3. The minimum, maximum, and median counts across the 100 samples of each size are summarized in Table 1.

Table 1 shows that $|{ℬ}_N|$ is consistently much smaller than $|{𝒜}_N|$, and $|{𝒞}_N|$ is even smaller. This indicates that, even though $|{𝒜}_N|$ can be prohibitively large, both $|{ℬ}_N|$ and $|{𝒞}_N|$ are typically small enough to allow exhaustive enumeration. We also see that the number of support networks varies widely even among networks of the same size, suggesting that the computational cost of exhaustive search within ${ℬ}_N$ or ${𝒞}_N$ may vary substantially depending on the structure of $N$. For example, the number $|{𝒞}_N|$ of minimum support networks of $N$ with $r=18$ was 1728 for a certain sample but was only 4 for another sample.

An obvious exact algorithm for Problem 2 computes $\mathrm{level}(G)$ of each element $G$ in ${𝒜}_N$ to obtain the minimum value, but this method is clearly impractical. Algorithm 1 performs an exhaustive search in ${ℬ}_N$ instead of ${𝒜}_N$. Recalling the numerical results in Table 1, we know that this search space reduction is significant. We can prove that ${ℬ}_N$ contains a correct $G^{\ast }$, yielding Theorem 19. The proof is given in Appendix A.4.

Algorithm 2 performs an exhaustive search in an even smaller search space ${𝒞}_N$. It repeats the same procedure as Algorithm 1, except its iteration number is reduced from $|{ℬ}_N|$ to $|{𝒞}_N|$. By Theorem 19, Algorithm 2 therefore runs in $O(|{𝒞}_N|\cdot |E(N)|)$ time.

Algorithm 2 does not always output a correct solution of Problem 2, because minimizing the number of reticulations does not imply minimizing the level of support networks. Consider a level-$1$-based network $N$ in Figure 4(a), which has the property that any $G\in {𝒞}_N$ satisfies $r(G)=2$ and $\mathrm{level}(G)=2$. Then, Algorithm 2 outputs a level-$2$ support network as illustrated in Figure 4(b). An optimal support network $G^{\ast }$ with $\mathrm{level}(G^{\ast })=1$ exists in ${ℬ}_N\setminus {𝒞}_N$ while it is not optimal in terms of reticulation numbers, as shown in Figure 4(c).

We evaluated the performance of our exact algorithm (Algorithm 1) and heuristic algorithm (Algorithm 2) using rooted binary phylogenetic networks with $n=8$ leaves and various reticulation numbers $r$. Although our experiment focuses on only binary networks, this simplification does not undermine the validity of our evaluation, as the search space sizes $|{ℬ}_N|$ for Algorithm 1 and $|{𝒞}_N|$ for Algorithm 2 are not affected by whether $N$ is binary or almost-binary.

For each $r$ ranging from 1 to 50, we generated 25 instances using the procedure described in Section 4.4. The experiments were conducted on a MacBook Pro equipped with an Apple M3 Max CPU (4.05 GHz) and 36 GB of memory, with a timeout limit of 30 minutes per instance. The Python code and experimental results, including the outputs of both algorithms for all instances, are available in our GitHub repository.

We first assessed computational efficiency of Algorithms 1 and 2 by measuring runtime on each instance. Figure 5 shows average runtime on a log scale with minimum and maximum values across 25 instances for each $r$. Both algorithms exhibit expected exponential runtime growth as $r$ increases. Algorithm 1 completed up to $r=36$, beyond which timeouts occurred, while Algorithm 2 successfully handled all instances up to $r=49$. The performance gap is particularly pronounced for larger networks. For example, at $r=35$, Algorithm 1 approaches its limits while Algorithm 2 completes within seconds on average.

Notably, both algorithms exhibit substantial runtime variability even for the same $r$, as shown by the wide error bars in Figure 5. This variation reflects runtime dependence on network structure beyond just reticulation number, consistent with our observations in Section 4.4 regarding the significant variability of $|B_N|$ and $|C_N|$ among instances with identical $r$ values.

While Algorithm 1 guarantees optimal solutions (Theorem 19), Algorithm 2 provides no such theoretical guarantee. To evaluate the practical quality of the heuristic solutions, we compared the outputs of Algorithm 2 against the optimal values ${\text{level}}^{\ast }(N)$ computed by Algorithm 1 for all instances with $r\le 36$ where both algorithms completed successfully within the time limit. Figure 6 provides a breakdown of the output differences $\mathrm{\Delta }$ for each $r$, where $\mathrm{\Delta }$ is defined as the difference between the value returned by Algorithm 2 and the optimal value ${\mathrm{level}}^{\ast }(N)$ computed by Algorithm 1. The bars show the proportion of instances (out of 25) with $\mathrm{\Delta }=0$ (exact match), $\mathrm{\Delta }=1$, $\mathrm{\Delta }=2$, and $\mathrm{\Delta }\ge 3$.

Although the accuracy gradually declines as $r$ increases, Algorithm 2 maintains practical effectiveness across the entire range of tested networks. In fact, up to $r\le 17$, Algorithm 2 achieved exact solutions ($\mathrm{\Delta }=0$) for all instances. Even at $r=30$, it still returned optimal values for approximately 80% of instances. Large deviations ($\mathrm{\Delta }\ge 3$) remain rare throughout the tested range, occurring in less than 10% of instances even for the largest networks.