Average-Tree Phylogenetic Diversity of Networks

Human-driven habitat degradation and environmental change have led to accelerated rates of species extinction, posing a significant threat to global biodiversity [18]. In response, conservation efforts face both logistical constraints – such as limited financial resources – and conceptual challenges, notably in identifying effective prioritization strategies. This has motivated the development of formal criteria and quantitative indicators to guide decision-making and assess biodiversity in a systematic and reproducible manner. Here, we study formal measures for quantifying the biodiversity of a given set of species and algorithms for computing and maximizing these measures.

Phylogenetic Diversity (PD) is a well-studied measure used to quantify the biodiversity of a set of species using their evolutionary relationships [22]. When these relationships are described by a rooted phylogenetic tree (with edge lengths) on a set $X$ of species, the phylogenetic diversity of a subset of the species $A\subseteq X$ is simply the total length of all edges that lie on at least one path from the root of the tree to a leaf in $A$ [8]. The main intuition behind this is that species that are evolutionarily further apart together contribute more diversity than species that are evolutionarily closely related. Phylogenetic diversity can more formally be interpreted as follows in terms of features (e.g., genes or morphological features). Suppose that (1) the number of features introduced on each edge of the tree is proportional to the length of that edge, (2) features are never lost, and (3) no feature is introduced independently on multiple branches. Then, the number of features in $A$ is proportional to its phylogenetic diversity. While the assumptions on the model may seem somewhat unrealistic, it has been used prominently in the literature [31, 3, 5, 24].

In many cases, however, describing evolutionary relationships by trees is too simplistic [7]. Rooted phylogenetic networks generalize rooted phylogenetic trees and are more suited to represent complex evolutionary scenarios involving, for example, hybridization events [12, 1, 26]. Roughly speaking, a rooted phylogenetic network is a directed acyclic graph with a single root and labeled leaves. Vertices with two (or more, in the nonbinary case) incoming arcs are called reticulations and can be used to model hybridization events, lateral gene transfer, or other types of horizontal evolutionary events. Generalizing the notion of phylogenetic diversity to networks is not obvious, and multiple variants have recently been proposed for rooted [32, 3, 5] and unrooted [20, 2, 4] networks.

The currently best-studied variants are AllPathsPD and NetworkPD [3, 14, 27, 5, 24]. Both definitions generalize PD on trees, in the following ways. $\mathrm{AllPathsPD}(A)$ is the total length of all edges that lie on at least one path from the root of the network to a leaf in $A$. NetworkPD generalizes AllPathsPD by allowing arcs to have inheritance probabilities, representing the probability that a given feature in the parent is inherited by the child. The idea of $\mathrm{NetworkPD}(A)$ is, roughly speaking, to compute the expected number of distinct features in $A$ by using the inheritance probabilities to compute, for each arc $uv$, the probability ${\gamma }_{uv}$ that at least one copy of a feature introduced on the arc $uv$ is inherited by a leaf in $A$. More precisely, ${\gamma }_{uv}$ is defined as ${\gamma }_{vw}$ times the inheritence probability of $uv$ if $v$ is a reticulation with child $w$, as $1-{\prod }_i(1-{\gamma }_{v{w}_i})$ if $v$ is a tree node with children $w_i$, and as $1$ or $0$ if $v$ is a leaf, depending on whether or not $v\in A$.

Note that, in the NetworkPD setting, features are inherited independently from different parents. For example, in the network in Figure 1, a feature introduced on the arc $\rho u$ is present in all parents of the reticulation above species $B$ but it has only a $0.75$ probability of being inherited by $B$. In addition, in NetworkPD it is possible to inherit the same feature from both parents, resulting in the presence of multiple copies of that feature. While this can be realistic in certain scenarios (e.g. allopolyploidy), in other cases it is more realistic to assume that each feature is inherited from at most one parent.

A less studied alternative averages the phylogenetic diversity score over all trees displayed by the network [31]. Each tree has an assigned probability equal to the product of the inheritance probabilities of the reticulation edges used by its embedding, and the PD-score of each tree is weighted by this probability. Thus, the resulting measure can be interpreted as the expected diversity of any tree displayed by the network. In this model, features inherited from different edges incoming to the same reticulation vertex are not considered to be independent. The assumption is basically that each feature is inherited from exactly one parent (or introduced on that edge). Note however that Wicke and Fischer [31] use an uncommon definition of displayed trees, in which trees are only considered to be displayed if they contain all vertices of the network and all their leaves are in $X$. In particular, a network that is not tree-based [13] would have no displayed trees and an undefined phylogenetic diversity score. Since we see no reason to ignore certain displayed trees, we will take the (weighted) average over all displayed trees in this paper. We call the resulting phylogenetic diversity measure on phylogenetic networks the average-tree phylogenetic diversity, APD for short, see Figure 1 for an example.

While APD is a natural and biologically relevant generalization of phylogenetic diversity from trees to networks, we show that it comes with great challenges. Although for most PD measures on networks, such as AllPathsPD and NetworkPD, the diversity score of a given set $A\subseteq X$ of species can be easily computed in polynomial time, we will show in this paper that for APD this is already computationally hard (more precisely, #P-hard; a polynomial-time algorithm for it would imply a collapse of the polynomial hierarchy.) Consequently, the maximization problem associated to APD is also computationally hard. This problem, called Max-APD, aims to find a set $A\subseteq X$ of $k$ species with maximum APD score.

On the positive side, we show that, although both these problems are #P-hard, they can both be solved exactly and rather efficiently in practice. We do this by presenting a parameterized algorithm finding a size-$k$ set of leaves that maximizes the APD score in binary networks with $r$ reticulations and $n$ leaves in $𝒪(2^{6r}\cdot nk)$ time. This algorithm can easily be adapted to compute the score of a given set $A$. The running time of our algorithm scales linearly with the number of species but exponentially with the number of reticulations in the network. Using practical experiments on simulated data we show that the algorithm can solve instances with up to $500$ species and $55$ reticulations within $5$ minutes on 16 GB RAM and 8 CPUs computer, with $k\le 5$. In addition, the algorithm efficiently solves instances with $k$ up to at least $40$ if the number of reticulations is fixed as $12$.

In the following, we prove that computing ${\text{APD}}_{𝒩}$ is #P-hard, even if all taxa are being saved, in contrast to other phylogenetic diversity measures on networks like NetworkPD [3, 27] or AllPathsPD [3, 14, 16], that can be computed (not optimized) in polynomial time.

We prove the theorem by reduction from Counting Perfect Matchings. A perfect matching of an undirected graph $G=(V,E)$ is a set of $|V|/2$ edges $E^{\prime }\subseteq E$ such that each vertex of $G$ is incident with exactly one edge of $E^{\prime }$. Counting the number of perfect matchings in a cubic graph is #P-hard [6, Theorem 6.2].

We first observe that the number of perfect matchings in a graph $G=(V,E)$ is the same as the number of mappings $\varphi :V\to E$ that are perfect assignments, as defined next. A mapping $\varphi :V\to E$ is an assignment if $\varphi (v)$ is an incident edge of $v$ for all $v\in V$. An assignment $\varphi$ is perfect if $|\varphi (V)|=|V|/2$. As $|{\varphi }^{-1}(e)|\le 2$ for each edge $e$, an assignment $\varphi$ is perfect if and only if $|{\varphi }^{-1}(e)|=2$ for each $e\in \varphi (V)$. That is, for each $e\in E$, a perfect assignment $\varphi$ assigns either both or none of the incident vertices to $e$. Thus, mapping each perfect assignment $\varphi$ to $\varphi (V)$ constitutes an injection from perfect assignments to perfect matchings in $G$. On the other hand, if $E^{\prime }\subseteq E$ is a perfect matching in $G$, then a perfect assignment $\varphi :V\to E$ can be defined from $E^{\prime }$ by setting $\varphi (v)$ to be the unique edge in $E^{\prime }$ incident to $v$ for all $v\in V$. It is easy to see that this mapping is an injection from the set of all perfect matchings in $G$ to the set of perfect assignments. Thus, the perfect matchings in $G$ are in bijection with its perfect assignments, so their number is the same. In the following, we let $M_G$ denote this number.

Given an arbitrary cubic graph $G=(V,E)$, we construct a network $𝒩$ on a set of taxa $X$ in such a way that $M_G$ can be determined from ${\text{APD}}_{𝒩}(X)$. See Figure 3 for an example.

Let $X:=\{x_v\mid v\in V\}$, that is, $X$ has one element for each vertex in $G$. Now, for each $v\in V$, add two vertices $r_v^1$ and $r_v^2$, where $r_v^2$ is a parent of $r_v^1$, and $r_v^1$ is a parent of $x_v$. For each $e\in E$, add a vertex $u_e$ to $𝒩$. For any vertex $v\in V$, choose one of its three incident edges $e$ arbitrarily (recall that $G$ is cubic), and add an arc $u_e{r}_v^1$. Then, for the remaining two edges $e^{\prime }$ and $e^{\prime \prime }$ incident with $v$, add the arcs $u_{e^{\prime }}{r}_v^2$ and $u_{e^{\prime \prime }}{r}_v^2$. Observe that, so far, $u_e$ has out-degree $2$ and in-degree $0$ for each $e\in E$, and $r_v^1$ and $r_v^2$ are reticulations for all $v\in V$. Furthermore, $u_e$ has a directed path to $x_v$ in $𝒩$ if and only if $e$ is incident to $v$ in $G$.

Create two new special vertices $u_1^{\ast }$ and $u_2^{\ast }$ and add an arc $u_1^{\ast }{u}_2^{\ast }$ (this will be the only arc that gets a large weight). Add a root $\rho$ as the parent of $u_1^{\ast }$. Now, for each $e\in E$, add a new vertex $r_e$ as the parent of $u_e$. For each vertex $r_e$, add two parents $p_e$ and $q_e$. Add arcs such that $\{u_1^{\ast }\}\cup \{q_e\mid e\in E\}$ and $\{u_2^{\ast }\}\cup \{p_e\mid e\in E\}$ form paths. We now have that each vertex $r_e$ is a reticulation, and exactly one of the incoming arcs of $r_e$ is on a directed path from the root that passes through $u_1^{\ast }{u}_2^{\ast }$. In Figure 4 an illustration is given.

We define the arc weights by setting $\omega (u_1^{\ast }{u}_2^{\ast }):=H$ for $H$ a large value to be defined later, and setting $\omega (a):=1$ for all other arcs $a\in A(𝒩)$.

It remains to define the inheritance probabilities on the reticulation arcs.

First, consider $r_e$ for some $e\in E$. Then, set $\lambda (p_e{r}_e):=1-\delta$ and $\lambda (q_e{r}_e):=\delta$, where $\delta \in (0,1)$ is a small constant to be defined later. Note that, when we choose ${𝒯}_{\sigma }$ with probability $P(\sigma )$ then, for each $e\in E$, the probability that $u_2^{\ast }$ has a path to some $r_e$ is exactly $1-\delta$.

Now, consider $r_v^1$ and $r_v^2$ for some $v\in V$. Then, set $\lambda (r_v^2{r}_v^1):=2/3$, $\lambda (a_1)=1/3$, and set $\lambda (a_2):=\lambda (a_3)=1/2$, where $a_1$ is the other arc incoming at $r_v^1$ and $a_2$ and $a_3$ are the two arcs incoming at $r_v^2$. The idea here is that when we choose a switching $\sigma$ and corresponding tree ${𝒯}_{\sigma }$ with probability $P(\sigma )$, each vertex $u_e$ for $e\in E$ adjacent to $v$ has probability $1/3$ of having a path to $x_v$ in ${𝒯}_{\sigma }$. As exactly one such $x_e$ will have a path to $x_v$, we can think of the choice of $\sigma$ as corresponding to a choice of assignment (where $v$ is mapped to $e$ in the assignment if $u_e$ has a path to $x_v$ in ${𝒯}_{\sigma }$). Our choice of inheritance probabilities ensures that each of the three edges incident to $v$ is chosen with probability $1/3$, independently for each $v$, and so each possible assignment is chosen with equal probability. This completes the reduction.

Given a switching $\sigma$, let ${\chi }_{\sigma }:=1$ if $u_2^{\ast }$ has a path to some leaf $x_v$ in ${𝒯}_{\sigma }$, and ${\chi }_{\sigma }:=0$ otherwise. Thus, the expected contribution of $u_1^{\ast }{u}_2^{\ast }$ to ${\text{APD}}_{𝒩}(X)$ is $H\cdot {\sum }_{\sigma \in 𝒮(R)}P(\sigma )\cdot {\chi }_{\sigma }$. We first show that this value depends entirely on the possible assignments of $G$.

In this section, we show that Max-APD is fixed-parameter tractable (FPT) with respect to the number of reticulations. For the sake of brevity, we restrict our theoretical proof to simple networks – that is, $𝒩$ has no cut-arcs except for those incident to degree-1 vertices. However, we see no reason why this result could not be extended to non-simple networks. Indeed, our implementation, see Section 5, can deal with arbitrary networks.

We make use of the notion of generators, as introduced in [28]. This gives us a certain useful partition of the leaves of a network into sides. At the high level, our algorithm “guesses” the optimal solution $Z$ by guessing which sides of the network $𝒩$ contain a leaf from $Z$. As the number of sides is $𝒪(|R(𝒩)|$, there are at most $2^{𝒪(|R(𝒩)|)}$ such guesses to consider. Once it has been decided which sides contain a leaf from the solution, the problem reduces to a problem on forests, which can be solved in polynomial time.

In what follows, let $𝒩$ be a network and let $r:=|R(𝒩)|$, $n:=|V(𝒩)|$ and $m:=|A(𝒩)|$. That is, $r$, $n$, and $m$ denote the number of reticulations, vertices, and arcs, respectively.