Comparing the topology of phylogenetic network generators

Remie Janssen¹ Pengyu Liu²^*^**To whom correspondence should be addressed; e-mail: [email protected]. ¹Delft Institute of Applied Mathematics, Delft University of Technology, Delft, The Netherlands ²Department of Mathematics, Simon Fraser University, Burnaby, Canada

Abstract

Phylogenetic networks represent evolutionary history of species and can record natural reticulate evolutionary processes such as horizontal gene transfer and gene recombination. This makes phylogenetic networks a more comprehensive representation of evolutionary history compared to phylogenetic trees. Stochastic processes for generating random trees or networks are important tools in evolutionary analysis, especially in phylogeny reconstruction where they can be utilized for validation or serve as priors for Bayesian methods. However, as more network generators are developed, there is a lack of discussion or comparison for different generators. To bridge this gap, we compare a set of phylogenetic network generators by profiling topological summary statistics of the generated networks over the number of reticulations and comparing the topological profiles.

1 Introduction

A phylogenetic network is a directed acyclic graph that represents the evolutionary history of a set of species [26]. With the recent advancements in the theory of phylogenetic networks and the accumulation of evolutionary data, it is becoming increasingly viable and common to represent the evolutionary hierarchical relations with phylogenetic networks instead of phylogenetic trees, which cannot record reticulate evolutionary processes such as horizontal gene transfer (also known as lateral gene transfer), hybrid speciation, hybrid introgression, and recombination [17, 4, 5].

In evolutionary analysis, the stochastic process that generates random trees and random networks play an important role, especially in the reconstruction of phylogenies. Similar to phylogenetic tree reconstruction, there are two main types of methods for rebuilding evolutionary histories as phylogenetic networks: combinatorial methods and model based methods. Combinatorial methods are based on, for example, distances between taxa [7, 9, 31] or gene trees [2, 38, 50, 29, 28], while model based methods involve Bayesian statistics [53, 55] or likelihood calculations [52, 47] in models of sequence and lineage evolution. For both types of methods, random phylogenetic networks are needed in validating the reconstruction methods. Moreover, for Bayesian methods, the stochastic processes used to generate random networks also serve as priors [53]. In this paper, we call the stochastic processes that generate random phylogenetic networks the network generators. Note that this is different from the generator of a network defined in [27], which encodes the underlying reticulate structure of a phylogenetic network.

There are a few approaches that existing network generators usually take to constructing random networks. One approach is to alter a tree generator and add hybrid edges, for example, network generators [29, 43] that are based on the Yule-Harding model [23] or on a continuous-time Markov model of speciation to a birth-hybridization Markov model [53, 42]. The hybrid edges can either be added during the process of building a tree or after a tree is generated. For these methods, there is also an option to add the hybrid edges in a way that favours short-distance (local) hybrid edges [42, 43].

Other approaches that do not rely on tree generating processes include assembling a network directly from degree constraints [54] and sampling a network in the network space by random local surgeries, e.g., rearrangement moves which move one endpoint of an edge [21, 32]. The latter is used in Bayesian methods for sampling posteriors [53, 55]. In order to present a more comprehensive comparison of network generators, we include two new network generators in our analysis, in addition to the ones already available in the literature. The first is based on the beta-splitting model for trees [1], which we extend to networks by adding edges to a generated tree in various ways. The second is based on a birth-death process [24] which we extend with distance dependent hybrid speciation events.

For phylogenetic trees, it is well discussed which stochastic process generates phylogenetic trees that resemble trees reconstructed from data [6], and random tree generators are well compared from different points of view [48, 40, 12, 3]. However, as increasingly many network generators are developed, there is a lack of discussion about the differences and similarities of generated random networks. To bridge this gap, we examine and compare network generators by profiling topological summary statistics of the generated networks over the number of reticulations in the networks and comparing the topological profiles. The topological summary statistics used for the comparison include network balance, number of cherries, number of reticulated cherries, number, size and level of blobs as well as network class. Moreover, we generalize the recent polynomial representation for trees [39] to networks and leverage the high-resolution polynomial representation in comparing small phylogenetic networks.

We show that generators based on different types of approaches can generate networks with similar topological profiles. If a network generator is based on a tree generator, the topological characteristics are dominated by the topology of the underlying trees when the number of reticulations is small; and as the number of reticulations increases, the topological characteristics may diverge with respect to the integrated reticulation structures and the locality of added reticulations.

2 Methods

2.1 Topological profiles

A phylogenetic network or simply a network is a rooted directed acyclic graph whose out-degree zero nodes (leaf nodes) have in-degree one. We study binary networks, that is, every internal node in the network has in-degree one and out-degree two (tree nodes) or in-degree two and out-degree one (reticulation nodes). In a phylogenetic network, a reticulation event is represented by a reticulation node and the edges pointing to the reticulation node which are called hybrid edges. A reticulation node together with the hybrid edges incident to it is a reticulation structure.

The topological profile of a network generator is defined to be the topological summary statistics of generated networks as a function of number of reticulations. In other words, the topological profile of a network generator describes how the topological summary statistics change as the number of reticulations increases. We choose commonly used topological characteristics of phylogenetic networks for our comparison. These characteristics include network balance, number of cherries and reticulated cherries, size and level of blobs as well as proportions of generated networks that belong to various network classes.

2.1.1 Network balance

The balance of phylogenetic trees is a well studied topological characteristic. The balance for phylogenetic networks is recently revisited in [3], where the $B_{2}$ balance for phylogenetic networks [46] is reintroduced. The $B_{2}$ balance is an index based on the entropy of probability distribution on leaf nodes. Specifically, let $N$ be a network with the set of leaf nodes denoted by $L$ and $p_{l}$ be the probability of a directed random walk starting from the root and ending at the leaf node $l$ , then the $B_{2}$ balance of $N$ can be computed by $B_{2}(N)=-\sum_{l\in L}p_{l}\log_{2}(p_{l})$ .

The maximal value of $B_{2}$ balance for a network with $n$ leaf nodes is $\log_{2}(n)$ , which corresponds to the most balanced case where every leaf node has probability $1/n$ . For our comparison, we use the normalized $B_{2}$ balance, which is computed by dividing the $B_{2}$ balance a of network with $n$ leaf nodes by $\log_{2}(n)$ .

2.1.2 Cherries and reticulated cherries

The number of small subgraphs such as cherries and pitchforks is another well-studied topological characteristics for phylogenetic trees [12]. A cherry of a network is a subgraph consisting of two leaf nodes $x$ and $y$ that share a common parent node. For our comparison of phylogentic networks, in addition to the number of cherries, we also compute the number of reticulated cherries which are subgraphs consisting of two leaf nodes $x$ and $y$ , a reticulation node $w$ as the parent node of $x$ (or $y$ ), and the common parent of $y$ and $w$ (or $x$ and $w$ ) [7].

2.1.3 Blobs

A blob in a network is a maximal connected subgraph without cut edges and consisting of at least two nodes [22]. The size of a blob in a network is the number of nodes in the blob; the level of a blob is the number of reticulation nodes in the blob; and the level of a network is the maximal level over all its blobs.

2.1.4 Network classes

To overcome the obstacles faced in extending the mathematical theory of phylogenetic trees to phylogenetic networks, some subclasses of phylogenetic networks were introduced. These network classes include tree-child networks, stack-free networks, orchard networks, and tree-based networks. We examine the proportion of generated networks that belong to each of the classes as another topological characteristic in the profile of a generator.

A network is tree-based if it has an embedded spanning tree whose leaf nodes are exactly the leaf nodes of the network [20]. A network is orchard if it can be reduced to a network consisting of a single vertex by a sequence of cherry reduction operations [33, 18]. A stack in a phylogenetic network is a pair of adjacent reticulation nodes. A network is stack-free if there are no stacks in the network [45]. A network is tree-child if each internal node of the network has at least one child that is not a reticulation node [11].

For binary networks, a tree-child network is also a stack-free network as well as an orchard network [33]. Moreover, every stack-free network or orchard network is also tree-based [54, 25, 30].

2.2 Polynomial comparison

A polynomial representation for phylogenetic trees was introduced in [40]. The polynomial representation is an isomorphic invariant for trees, that is, two trees are isomorphic if and only if they have the same polynomial representation. Moreover, the polynomials for rooted trees are irreducible, that is, they cannot be factored as a multiplication of smaller polynomials [39]. We leverage these results to extend the polynomial representation for further comparing phylogenetic network generators.

A natural way to connect phylogenetic networks with trees is by considering their sets of spanning trees. Let $N$ be a phylogenetic network and $\mathcal{T}_{N}$ be the deck of rooted spanning trees of $N$ , that is, $\mathcal{T}_{N}$ may contain identical elements. It is natural to represent $N$ with the polynomial $P(N,x,y)=\prod_{T\in\mathcal{T}_{N}}P(T,x,y)$ . It is unknown in general whether graphs can be reconstructed from their decks of spanning trees. However, it has been proved that two tree-child networks are isomorphic if and only if they have the same set of embedded spanning trees [19]. If we compute the polynomial for tree-child networks based on the set of embedded spanning trees, it is immediate that two tree-child networks are isomorphic if and only if they have the same polynomial. However, although not all networks are tree-child networks, it is still useful to compute their polynomials for comparison.

To compare the networks, We use the polynomial distance defined in [40], that is, the Canberra distance [37] between the polynomial coefficient vectors of the networks. Then we use multi-dimensional scaling [13] to visualize the distances, and we deploy k-medoids clustering [35] on the distances to determine if networks generated by different generators can be distinguished by the polynomial distance.

2.3 Network generators

Phylogenetic network generators typically alter tree generators by introducing reticulation event to the process. An important phylogenetic tree generator is the beta-splitting model which is a generalization of the Yule model and the PDA model [1, 6].

For our comparison, we select the LGT network generator [43] and the ZODS network generator [53] as representations of network generators based on the Yule model (i.e., beta-splitting model with $\beta=0$ ). We compare these network generators to a new network generator which is an extension of Heath’s tree generator [24]. Another network generator we choose is the NTK network generator [54] which is described not in terms of a branching process, but it utilizes the degree constraints of networks for direct sampling. Moreover, for comparison, we also introduce two additional network generators as references. The first is an extension of the beta-splitting tree generator. Instead of adding reticulations in the process, the beta-splitting network generator adds reticulations to generated trees in various ways. The second is the MCMC network sampler which is commonly used in Bayesian methods for sampling posteriors, but we employ it as a standalone generator.

For network generators that are based on a phylogenetic tree generator, we will discuss three types of reticulation structure integrated in the branching process. See Figure 1. For a pair of leaf nodes $x$ and $y$ of the network under generation at the current state, the n-type reticulation adds an edge from $x$ to a new leaf node $x^{\prime}$ , an edge from $y$ to a new leaf node $y^{\prime}$ and a horizontal edge from $x$ to $y$ (or $y$ to $x$ ); the y-type reticulation glues the leaf node $x$ and $y$ and produces a new leaf node $z$ ; the m-type reticulation branches at $x$ and $y$ and glues a pair of new leaf nodes producing three new leaf nodes $x^{\prime}$ , $y^{\prime}$ and $z$ . The n-type reticulation structure has a natural interpretation as an HGT or an hybrid introgression event, and the m-type can be interpreted as a hybrid speciation event.

Figure 1: Different types of reticulation structures, where reticulation nodes and hybrid edges are in red.

2.3.1 LGT network generator

The LGT network generator is a network generator designed to model the lateral gene transfer (LGT) events. The generator is introduced in [43], and, from its description, the speciation events are designed as in the Yule model, that is, a leaf node at the current state is chosen uniformly at random for branching. The LGT network generator uses n-type reticulations and has two parameters governing the introduction of reticulation events: a parameter $\alpha\in[0,1]$ controls the probability of the next event being a reticulation event instead of a speciation event, and a parameter $\gamma\geq 0$ regulates the probability of a reticulation event being within an existing blob, where smaller $\gamma$ corresponds to higher probability of a local reticulation event. Hence, if the number of events is fixed, the parameter $\alpha$ influences the number of reticulations in the network, and the parameter $\gamma$ influences the level of the network. For our comparison, we modified the code of the generator so that we have networks with the same number of reticulations and number of leaf nodes, making the $\alpha$ parameter irrelevant for our comparison.

2.3.2 ZODS network generator

The ZODS network generator is a birth-hybridization process introduced in [53] as a prior in the Bayesian method. The generator uses y-type reticulations and has a constant speciation rate $\lambda$ and a constant hybridization rate $\nu$ , that is, at a step when there are $k$ extant leaf nodes in a network, the speciation rate is $\lambda k$ and hybridization rate $k(k-1)\nu/2$ . Besides the reticulation type, another difference between the ZODS network generator and the LGT network generator is that the hybridization rate. The hybridization rate in the ZODS network generator depends quadratically on the number of leaf nodes, whereas in the LGT network generator, the hybridization rate depends linearly on the number of leaf nodes.

2.3.3 NTK network generator

A binary phylogenetic network has the degree constraints that a tree node has in-degree one and out-degree two, a leaf node has in-degree one and out-degree zero, and a reticulation node has in-degree two and out-degree one. The NTK network generator introduced in [54] does not rely on a branching process but utilizes these degree constraints and generates networks as follows. To generate a random network of a certain size, the generator selects uniformly at random $n$ nodes as leaf nodes and $r$ nodes as reticulation nodes from $2(n+r)-1$ nodes, and randomly inserting edges between a pair of nodes so that the degree constraints are met.

2.3.4 Beta-splitting network generator

The beta-splitting network generator is based on the beta-splitting model for trees [1, 6, 44]. The beta-splitting model for trees is a generalization of the Yule model ( $\beta=0$ ) and the PDA model ( $\beta=-1.5$ ). Here, we introduce three methods to create a network from a generated tree by repeatedly adding edges to the tree. The first method picks two edges uniformly at random, subdivides these two edges with two vertices respectively and adds another new edge between these two vertices with a random direction that does not introduce cycles (Algorithm 1); the second chooses edges more locally governed by a single parameter $P_{\text{stop}}\in[0,1]$ that stochastically determines how local the cycles are, a high $P_{\text{stop}}$ leads to more local reticulations (Algorithm 2); the third picks two leaf nodes, subdivides their incoming edges with two vertices and adds a horizontal edge between the two vertices with a random direction (Algorithm 3). We call the process of generating networks based on the PDA trees in these ways the PDA network generator, and the generated networks the PDA networks.

2.3.5 MCMC network sampler

The Markov-chain-Monte-Carlo (MCMC) network sampler performs a random walk through the space of phylogenetic networks by applying rSPR moves, small changes to the network that move one endpoint of an edge to another edge [21]. By selecting networks after sufficiently many random rSPR moves, the MCMC network sampler is able to sample (almost) uniformly from the set of (internally labelled) networks, as this space is connected [32]. The MCMC network sampler lies at the basis of a common method in Bayesian analysis for sampling posteriors, as in PhyloNet [49] and SpeciesNetwork [53].

2.3.6 Heath network generator

Finally, we present a new network generator based on the branching process introduced in [24], which uses auto-correlated speciation and extinction rates. For our analysis, we set the extinction rate to zero. The branching process is extended by adding hybrid speciation events, i.e., we integrate the m-type reticulation in this model. We introduce a new distance dependent hybridization rate $f_{hyb}(d)$ between two leaf nodes at the current state, which depends on the weighted distance $d$ between the leaf nodes, where $d$ is the sum of the lengths of all up-down paths between the two leaves, weighted by the probability of the up-down path, i.e., the product of the edge-probabilities on the path. The function $f_{hyb}(d)$ is piecewise linear, and it is linearly decreasing in $d$ (for $h_{l}\leq d\leq h_{r}$ ) with a minimum and maximum rate $h_{lr}$ and $h_{rr}$ .

f_{hyb}(d)=\begin{cases}h_{lr}&\mbox{if }d\leq h_{l}\\ h_{lr}+\frac{(h_{rr}-h_{lr})(d-h_{l})}{h_{r}-h_{l}}&\mbox{if }h_{l}<d<h_{r}\\ h_{rr}&\mbox{if }d\geq h_{r}\end{cases}.

See A.2.1 for more details of the generator, and for an efficient method to update the weighted distance $d$ after several types of evolutionary events. Using these, our distance dependent extension is also viable for other tree generators than the Heath tree generator.

2.4 Data

To compare the network generators, we generate a large number of networks with each of them, varying the parameters if there are any. To keep the comparison straightforward, we generate networks with 100 leaf nodes for this purpose, and we vary the number of reticulations. When using the polynomial to compare networks, we use smaller networks with 10 leaves.

LGT networks

For the LGT network generator, we use three different parameter values ( $0.01$ , $0.1$ , and $1.0$ ) for the parameter $\gamma$ that regulated the locality of the reticulations. The parameter $\alpha$ becomes irrelevant as we fix the number of reticulation nodes and leaf nodes as mentioned before. For each of these parameter values, we generate 500 networks with $r=1$ , 2, 3, 4, 5, 10, 20, 50, 100, 200, 500 reticulations. This results in a set of $16500$ networks.

ZODS networks

Setting the number of leaf nodes and reticulation nodes, the parameters become irrelevant as we have argued before. Hence, we generate networks with a set number of reticulations using a reticulation rate that gives a large probability to generate such a network. We used the parameters $\nu=0.0005,0.001,0.002$ to create 500 networks for each reticulation number $r=1,\ldots,5$ ; we used $\nu=0.005$ to create 500 networks with 10 and 20 reticulations, $\nu=0.01$ for 500 networks with 20 and 50 reticulations, $\nu=0.015$ for 500 networks with 50 and 100 reticulations, and, finally, $\nu=0.02$ for 500 networks with 100 and 200 reticulations. Hence, the full set of networks totalled $23\cdot 500=11500$ networks.

NTK networks

For $r=1,2,3,4,5,10,20,50,100,200$ reticulations, we downloaded 500 networks with 100 leaf nodes from the webpage applet. This generator has no further parameters, so we have a total of $5000$ networks for this generator.

Beta-splitting networks

We used the beta values $0.0$ , $-0.5$ , $-1.0$ , and $-1.5$ . For each of these values, we created 50 trees and added $r=1,2,3,4,5,10,20,50,100,200,500$ reticulations to these trees to create 10 networks with each of the following methods: Uniformly , Horizontally , and locally with stopping probability $P_{\text{stop}}\in\{0.2,0.02,0.002\}$ . This results in a total of $4\cdot 50\cdot 11\cdot 10\cdot 5=110000$ networks.

MCMC sampled networks

We sample networks with $r=1,2,3,4,5,10,20,50,100,200,500$ reticulations, starting with a network sampled by the beta-splitting network generator with $\beta=-1.0,-1.5$ adding edges uniformly. For each of these parameter combinations, we have sampled $500$ networks. This results in a set of $1000$ networks for each reticulation number, and $11000$ networks in total.

Each sample is taken after $5000$ proposals. Because about half of the proposals were rejected as invalid rSPR moves, we took about $2500$ steps through the space of networks between each sample. As the diameter of the space of networks with $n$ leaf nodes and $r$ reticulation nodes is at most $2n+3r-2$ ; see [34]. This is reasonable to make the samples independent for $n=100$ and $r\leq 500$ .

Heath networks

For each of these options, we create $500$ networks where we add $r=1$ , 2, 3, 4, 5, 10, 15, 20, 30, 40, 50 reticulations locally, and similarly for less local reticulations. This results in $3\cdot 500\cdot 11\cdot 2=33000$ networks. For parameter sets used, see Table 3.

Note that we do not generate networks with 100, 200, and 500 reticulations with this generator. This is because the Heath network generator uses m-type reticulation events, which add a leaf to the network. Hence, to generate a network with $r$ reticulations, the network must contain at least $r+2$ leaves, so it is impossible to generate networks with 100 leaves and 100 reticulations.

3 Results

3.1 Topological profiles

3.1.1 Balance

It is displayed in Figure 2 that the balance of the generated networks depends primarily on the balance of the underlying tree model. For the generators based on the Yule model, that is, the LGT network generator and the ZODS network generator, the $B_{2}$ balance of networks with few reticulations are similar to the beta-splitting networks at $\beta=0$ (Yule model).

Similarly, for beta-splitting networks, the values of the parameter $\beta$ determining the balance of the underlying trees also determine the $B_{2}$ balance of the network when the number of reticulations is small. As the number of reticulations increases, the balance differs with respect to the method of adding reticulations. Adding reticulations horizontally at leaf nodes makes the network more balanced, while adding reticulations locally leads to less balanced networks.

Regarding the MCMC network sampler, networks with more reticulations are more balanced on average than networks with fewer reticulations. This is because, in the MCMC network sampler, the networks with few reticulations are relatively unbalanced compared to the other network generators, as unbalanced networks are over-represented among all distinguishable arrangements. In fact, when samples are taken after sufficiently many steps, the MCMC network sampler theoretically behaves exactly the same as the PDA model (beta-splitting model with $\beta=-1.5$ ) on trees, which is known to be relatively unbalanced.

Refer to caption — Figure 2: Dependence of the $B_{2}$ balance on the reticulation number in network generators.

3.1.2 Cherries and reticulated cherries

The relations between the number of cherries, the number of reticulated cherries and the number of reticulations for different network generators are displayed in Figure 3. For all network generators examined, the number of cherries decreases as the number of reticulations increases. This is because when adding hybrid edges, existing cherries can be removed while new cherries can never be created. It is clear that for network generators based on a branching processes, the initial number of cherries depends on the underlying tree shapes. In particular, the network generators based on the Yule model, namely the LGT and ZODS network generators, have similar initial number of cherries to the beta-splitting network with $\beta=0$ (Yule model). It can also be observed that the initial number of cherries of the NTK network generator is similar to the beta-splitting network with $\beta=0$ as well, and the number of cherries of the MCMC network sample resembles the PDA networks (the beta-splitting network with $\beta=0$ ).

For the LGT, ZODS, and Heath network generators, the number of reticulated cherries increases with the number of reticulations. The rate with which the number of reticulated cherries increases depends on the reticulation structure used in the process. Specifically, the ZODS network generator uses y-type reticulation structure, and the number of reticulated cherries increases the slowest; The LGT network generator uses n-type reticulation structure, and the number of reticulated cherries increases faster; the Heath network generator uses m-type reticulation structure, and the number of reticulated cherries increases the fastest. Also note that adding reticulations locally makes the number of reticulated cherries go up faster in the LGT network generator.

3.1.3 Blobs

We examine the relations between the number of blobs, the level, the blob size of networks and the number of reticulations in the networks generated by different generators. Figure 4 displays the results. For all the generators examined, The blob size increases slower than the number of reticulations, but eventually the ratio between the blob size and the number of reticulations becomes constant as the number of blobs in the networks becomes one. More precisely, in a network with $n$ leaf nodes and $r$ reticulations, as $r$ increases and the number of blobs in the networks becomes one, the number of nodes in the largest blobs converges to $n+2r\approx 2r$ and the ratio between the level of the network and the number of reticulations $r$ converges to $1$ .

Trivially, adding edges locally results in more and smaller blobs. The type of locality does affect the profile for the number of blobs subtly as well. The LGT network generator explicitly forces local reticulation events to stay within a blob. Hence, blobs cannot merge in a local event and networks often have multiple blobs even when they have a large number of reticulations. The beta-splitting network generator and the Heath network generator enforce locality differently, where blobs are allowed to merge, which nearly always leads to networks with one blob when a large number of reticulations are added.

Despite the qualitative similarity in the profiles, the Heath network generator with m-type reticulations has distinct profiles with a more significant change in the number of blobs. Similarly, the beta-splitting network generator with different parameter $\beta$ also generates qualitatively similar profiles but with differences in the quantities; see Figure 2

3.1.4 Network classes

Lastly, We study the proportions of generated networks that belong to each of the network classes. Figure 5 displays the proportions with respect to the number of reticulations. As the number of reticulations increases, most of the generated networks fall out of any of the four classes.

For the beta-splitting networks, the beta parameter had negligible influence, so we combined data from all values of betas; see the left panel in Figure 5. The curves for all different methods of adding reticulations are very similar. One exception is that adding edges horizontally always leads to networks that are orchard, and thus also tree-based. In Lemma B.2, we show that it is in fact impossible for this method to produce non-orchard networks.

For the other network generators, the general results are similar. The exceptions are LGT networks and Heath networks. The LGT network generator, like the beta-splitting model that adds edges horizontally, only produces orchard networks. The Heath network generator is the most restricted, because these only produce tree-child networks as a result of its m-type reticulations.

3.2 Polynomial comparison of small networks

We deploy the polynomial as an alternative comparison between generated networks. We compare the network generators with similar profiles, namely the ZODS network generator versus the NTK network generator and the MCMC network sampler versus PDA network generator (with reticulations added uniformly). As a reference, we also compare the LGT network generator (with reticulations added uniformly, i.e. $\gamma=1.0$ ) vs the NTK network generator. Note that other than the network classes, the LGT network generator and the NTK network generator have similar topological profiles. Figure 6 displays the misclassification rates from k-medoids clustering sets of networks, the experiment is performed on networks with 10 leaves, and repeated 50 times for each number of reticulations. A misclassification rate of $0.5$ means that the two sets of generated networks cannot be distinguished, while a misclassification rate of $0$ means the two sets of generated networks can be completely distinguished. When the number of reticulations in the networks reaches 10, the polynomial distance can distinguish the LGT networks and NTK networks. However, the polynomial distance cannot successfully distinguish the ZODS networks from the NTK netowrks or MCMC networks from the PDA networks.

4 Discussion

We have compared the topological profiles of several network generators. In general, the topological characteristics are related to the underlying tree shapes when the number of reticulations is small. For generators based on a branching process, as the number of reticulations increases, the topological characteristics change with respect to the type of integrated reticulation structures and the locality of added reticulations. For example, y-type reticulations (ZODS network generator) can be used to create any network in any network class, and they lead to a relatively large number of cherries which seems to persist even for very large numbers of reticulations. This large number of cherries can be explained by the relatively large number of speciation events needed to generate networks with larger numbers of reticulations. Using n-type reticulations (LGT network generator) always leads to orchard networks, and a relatively high balance compared to the other generators. While using m-type reticulations (Heath network generator) always leads to tree-child networks, and to a sharp increase in the number of reticulated cherries, because each m-type reticulation involves two reticulated cherries.

For generators that are not based on a branching process, we have found resemblances between the NTK network generator and the ZODS network generator, as well as between MCMC network sampler and the PDA network generator (with reticulations added uniformly). The former resemblance could be explained by seeing the ZODS network generator as having a predefined set of nodes with degree constraints, which are connected in a top-down fashion. This makes it quite similar to the NTK network generator, but whose description does not explicitly include a top-down order for connecting the nodes. The resemblance between the MCMC network sampler and the PDA network generator is less obvious, as there only a clear relation between these methods when there are no reticulations. It remains an open problem whether the theoretical equivalence for these models on trees can be extended formally to an equivalence on networks.

These resemblances have been further analyzed by comparing small networks generated by the set of generators using the polynomial distance. We found that among the pairs of network generators with similar topological profiles, the polynomial distance can distinguish the LGT networks from the NTK networks, but cannot successfully distinguish ZODS networks from the NTK networks or MCMC networks from PDA networks. We limited our analysis to small networks due to the computation complexity of the polynomial. For larger phylogenetic networks, we can compute the average of coefficients over the set of spanning-tree polynomials to avoid too many multiplications of polynomials. Consider the polynomial coefficients as vectors, we can also compute the diversity [10] of the set of spanning-tree polynomials of a network. This would provide additional information about the topology of a network.

Perhaps surprisingly, we observed that adding reticulations locally does not only affect the number of blobs and their sizes, but also affects the balance of the networks. This effect is more prominent in the beta-splitting and Heath network generators than in the LGT network generator. It would be interesting to see whether it can be explained mathematically why local reticulations lead to networks with lower balance, and why this effect is not as strong as in the LGT network generator.

As we mentioned in the introduction, there are two main applications for network generators in phylogeny reconstruction: as a tool to validate network reconstruction methods, and as prior for Bayesian methods. The profiles and the analysis of the network topology presented in this paper can help select suitable generators in these tasks. It is an important question in validation which generator produces realistic networks that resemble networks reconstructed from data. This question is discussed for phylogenetic trees [6], but is not investigated for phylogenetic networks. One of the gaps in this investigation is the lack of realistic networks reconstructed from data. Once sufficient realistic networks become available, for example through reconstructions [41, 51] or through simulated evolution [14, 15, 16], the topological profiles and the polynomial can be utilized to compare the natural evolutionary processes with simulated evolution models or the random network generators, and networks reconstructed from data with the generated random networks.

To provide a more flexible framework for adding hybridization events to tree generators, we also developed a new network generator based on the Heath branching process. This network generator uses the m-type reticulation structure, and hybrid events are more likely for closely related species, which is modelled by a distance that we update efficiently after each event. Consequently, the Heath network generator has different topological profiles from the other network generators. For all these network generators based on a branching process, we did not include extinction events, which would change the topological profiles; see C.1. Moreover, extinction would allow network generators using m-type reticulations to have more reticulations than leaves, and network generators with m-type and n-type reticulations to generate networks outside the class of tree-based networks.

In the perspective of Bayesian methods, our results suggest that the priors introduced in [49] and [53] can behave quite differently in terms of balance and number of (reticulated) cherries, for in the former an explicit prior is absent so the actual prior is the MCMC network sampler, and in the latter, the prior is the ZODS network generator. Moreover, if a generator that is unable to generate all networks (e.g. the LGT network generator) is used as prior, then there exist networks that cannot be traversed by the random walk in sampling posterior distribution. Using such a prior thus requires new investigations on the connectedness of the spaces of networks that can be generated by such a generator, like for tree-child networks in [8, 36].

Acknowledgments

R.J. was supported by the Netherlands Organization for Scientific Research (NWO) Vidi grant 639.072.602. P.L. was supported by the grant of the Federal Government of Canada’s Canada 150 Research Chair program to Prof. C. Colijn.

References

Aldous [1996] D. Aldous. Probability distributions on cladograms. In Random discrete structures, pages 1–18. Springer, 1996.
Baroni et al. [2005] M. Baroni, S. Grünewald, V. Moulton, and C. Semple. Bounding the number of hybridisation events for a consistent evolutionary history. Journal of mathematical biology, 51(2):171–182, 2005.
Bienvenu et al. [2020] F. Bienvenu, G. Cardona, and C. Scornavacca. Revisiting shao and sokal’s $b\_2$ index of phylogenetic balance. arXiv preprint arXiv:2010.08079, 2020.
Blair and Ané [2020] C. Blair and C. Ané. Phylogenetic trees and networks can serve as powerful and complementary approaches for analysis of genomic data. Systematic biology, 69(3):593–601, 2020.
Blais and Archibald [2021] C. Blais and J. M. Archibald. The past, present and future of the tree of life. Current Biology, 31(7):R314–R321, 2021.
Blum and François [2006] M. G. Blum and O. François. Which random processes describe the tree of life? a large-scale study of phylogenetic tree imbalance. Systematic Biology, 55(4):685–691, 2006.
Bordewich and Semple [2016] M. Bordewich and C. Semple. Determining phylogenetic networks from inter-taxa distances. Journal of mathematical biology, 73(2):283–303, 2016.
Bordewich et al. [2017] M. Bordewich, S. Linz, and C. Semple. Lost in space? generalising subtree prune and regraft to spaces of phylogenetic networks. Journal of theoretical biology, 423:1–12, 2017.
Bordewich et al. [2018] M. Bordewich, K. T. Huber, V. Moulton, and C. Semple. Recovering normal networks from shortest inter-taxa distance information. Journal of mathematical biology, 77(3):571–594, 2018.
Bryant and Tupper [2012] D. Bryant and P. F. Tupper. Hyperconvexity and tight-span theory for diversities. Advances in Mathematics, 231(6):3172–3198, 2012.
Cardona et al. [2008] G. Cardona, F. Rosselló, and G. Valiente. Comparison of tree-child phylogenetic networks. IEEE/ACM Transactions on Computational Biology and Bioinformatics, 6(4):552–569, 2008.
Choi et al. [2020] K. P. Choi, A. Thompson, and T. Wu. On cherry and pitchfork distributions of random rooted and unrooted phylogenetic trees. Theoretical population biology, 132:92–104, 2020.
Cox and Cox [2001] M. Cox and M. Cox. Multidimensional Scaling. Chapman and Hall, 2001.
Cuypers and Hogeweg [2015] T. Cuypers and P. Hogeweg. Endless evolutionary paths to virtual microbes. In Workshop, First EvoEvo Workshop, satellite workshop of ECAL2015, 2015.
van Dijk et al. [2016] B. van Dijk, T. Cuypers, and P. Hogeweg. Evolution of r-and k-selected species of virtual microbes: a case study in a simple fluctuating 2-resource environment. In 2nd EvoEvo Workshop, satellite workshop of CCS2016, 2016.
van Dijk et al. [2020] B. van Dijk, P. Hogeweg, H. M. Doekes, and N. Takeuchi. Slightly beneficial genes are retained by bacteria evolving dna uptake despite selfish elements. Elife, 9:e56801, 2020.
Elworth et al. [2019] R. L. Elworth, H. A. Ogilvie, J. Zhu, and L. Nakhleh. Advances in computational methods for phylogenetic networks in the presence of hybridization. In Bioinformatics and Phylogenetics, pages 317–360. Springer, 2019.
Erdős et al. [2019] P. L. Erdős, C. Semple, and M. Steel. A class of phylogenetic networks reconstructable from ancestral profiles. Mathematical biosciences, 313:33–40, 2019.
Francis and Moulton [2018] A. Francis and V. Moulton. Identifiability of tree-child phylogenetic networks under a probabilistic recombination-mutation model of evolution. Journal of theoretical biology, 446:160–167, 2018.
Francis and Steel [2015] A. R. Francis and M. Steel. Which phylogenetic networks are merely trees with additional arcs? Systematic biology, 64(5):768–777, 2015.
Gambette et al. [2017] P. Gambette, L. van Iersel, M. Jones, M. Lafond, F. Pardi, and C. Scornavacca. Rearrangement moves on rooted phylogenetic networks. PLoS computational biology, 13(8):e1005611, 2017.
Gusfield and Bansal [2005] D. Gusfield and V. Bansal. A fundamental decomposition theory for phylogenetic networks and incompatible characters. In Annual International Conference on Research in Computational Molecular Biology, pages 217–232. Springer, 2005.
Harding [1971] E. Harding. The probabilities of rooted tree-shapes generated by random bifurcation. Advances in Applied Probability, pages 44–77, 1971.
Heath et al. [2008] T. A. Heath, D. J. Zwickl, J. Kim, and D. M. Hillis. Taxon sampling affects inferences of macroevolutionary processes from phylogenetic trees. Systematic biology, 57(1):160–166, 2008.
Huber et al. [2019] K. T. Huber, L. van Iersel, R. Janssen, M. Jones, V. Moulton, Y. Murakami, and C. Semple. Rooting for phylogenetic networks. arXiv preprint arXiv:1906.07430, 2019.
Huson et al. [2010] D. H. Huson, R. Rupp, and C. Scornavacca. Phylogenetic networks: concepts, algorithms and applications. Cambridge University Press, 2010.
van Iersel et al. [2009] L. van Iersel, J. Keijsper, S. Kelk, L. Stougie, F. Hagen, and T. Boekhout. Constructing level-2 phylogenetic networks from triplets. IEEE/ACM Transactions on Computational Biology and Bioinformatics, 6(4):667–681, 2009.
van Iersel et al. [2019a] L. van Iersel, R. Janssen, M. Jones, Y. Murakami, and N. Zeh. Polynomial-time algorithms for phylogenetic inference problems involving duplication and reticulation. IEEE/ACM transactions on computational biology and bioinformatics, 17(1):14–26, 2019a.
van Iersel et al. [2019b] L. van Iersel, R. Janssen, M. Jones, Y. Murakami, and N. Zeh. A practical fixed-parameter algorithm for constructing tree-child networks from multiple binary trees. arXiv preprint arXiv:1907.08474, 2019b.
van Iersel et al. [2020a] L. van Iersel, R. Janssen, M. Jones, Y. Murakami, and N. Zeh. A unifying characterization of tree-based networks and orchard networks using cherry covers. arXiv preprint arXiv:2004.07677, 2020a.
van Iersel et al. [2020b] L. van Iersel, V. Moulton, and Y. Murakami. Reconstructibility of unrooted level-k phylogenetic networks from distances. Advances in Applied Mathematics, 120:102075, 2020b.
Janssen [2021] R. Janssen. Rearranging Phylogenetic Networks. PhD thesis, PhD thesis, Delft University of Technology, 2021.
Janssen and Murakami [2021] R. Janssen and Y. Murakami. On cherry-picking and network containment. Theoretical Computer Science, 856:121–150, 2021.
Janssen et al. [2018] R. Janssen, M. Jones, P. L. Erdős, L. van Iersel, and C. Scornavacca. Exploring the tiers of rooted phylogenetic network space using tail moves. Bulletin of mathematical biology, 80(8):2177–2208, 2018.
Kaufman and Rousseeuw [2009] L. Kaufman and P. J. Rousseeuw. Finding Groups in Data: An Introduction to Cluster Analysis. Hoboken, New Jersey: John Wiley & Sons, Inc., 2009.
Klawitter [2020] J. Klawitter. Spaces of phylogenetic networks. PhD thesis, PhD thesis, University of Auckland, in preparation, 2020.
Lance and Williams [1966] G. N. Lance and W. T. Williams. Computer Programs for Hierarchical Polythetic Classification (“Similarity Analyses”). The Computer Journal, 9(1):60–64, 1966.
Linz and Semple [2019] S. Linz and C. Semple. Attaching leaves and picking cherries to characterise the hybridisation number for a set of phylogenies. Advances in Applied Mathematics, 105:102–129, 2019.
Liu [2021] P. Liu. A tree distinguishing polynomial. Discrete Applied Mathematics, 288:1–8, 2021.
Liu et al. [2020] P. Liu, P. Biller, M. Gould, and C. Colijn. Polynomial phylogenetic analysis of tree shapes. bioRxiv, 2020.
Morando et al. [2020] M. Morando, M. Olave, L. J. Avila, J. W. Sites Jr, and A. D. Leaché. Phylogenomic data resolve higher-level relationships within south american liolaemus lizards. Molecular phylogenetics and evolution, 147:106781, 2020.
Morin and Moret [2006] M. M. Morin and B. M. Moret. Netgen: generating phylogenetic networks with diploid hybrids. Bioinformatics, 22(15):1921–1923, 2006.
Pons et al. [2019] J. C. Pons, C. Scornavacca, and G. Cardona. Generation of level-k lgt networks. IEEE/ACM Transactions on Computational Biology and Bioinformatics, 17(1):158–164, 2019.
Sainudiin and Véber [2016] R. Sainudiin and A. Véber. A beta-splitting model for evolutionary trees. Royal Society open science, 3(5):160016, 2016.
Semple and Simpson [2018] C. Semple and J. Simpson. When is a phylogenetic network simply an amalgamation of two trees? Bulletin of mathematical biology, 80(9):2338–2348, 2018.
Shao [1990] K.-T. Shao. Tree balance. Systematic Zoology, 39(3):266–276, 1990.
Solís-Lemus and Ané [2016] C. Solís-Lemus and C. Ané. Inferring phylogenetic networks with maximum pseudolikelihood under incomplete lineage sorting. PLoS genetics, 12(3):e1005896, 2016.
Steel [2016] M. Steel. Phylogeny: discrete and random processes in evolution. SIAM, 2016.
Wen et al. [2016] D. Wen, Y. Yu, and L. Nakhleh. Bayesian inference of reticulate phylogenies under the multispecies network coalescent. PLoS genetics, 12(5):e1006006, 2016.
Whidden et al. [2013] C. Whidden, R. G. Beiko, and N. Zeh. Fixed-parameter algorithms for maximum agreement forests. SIAM Journal on Computing, 42(4):1431–1466, 2013.
Yang et al. [2021] W. Yang, N. Feiner, C. Pinho, G. M. While, A. Kaliontzopoulou, D. J. Harris, D. Salvi, and T. Uller. Extensive introgression and mosaic genomes of mediterranean endemic lizards. Nature communications, 12(1):1–8, 2021.
Yu and Nakhleh [2015] Y. Yu and L. Nakhleh. A maximum pseudo-likelihood approach for phylogenetic networks. BMC genomics, 16(10):1–10, 2015.
Zhang et al. [2018] C. Zhang, H. A. Ogilvie, A. J. Drummond, and T. Stadler. Bayesian inference of species networks from multilocus sequence data. Molecular biology and evolution, 35(2):504–517, 2018.
Zhang [2016] L. Zhang. On tree-based phylogenetic networks. Journal of Computational Biology, 23(7):553–565, 2016.
Zhu et al. [2018] J. Zhu, D. Wen, Y. Yu, H. M. Meudt, and L. Nakhleh. Bayesian inference of phylogenetic networks from bi-allelic genetic markers. PLoS computational biology, 14(1):e1005932, 2018.

Appendix A Description and pseudo-code for the new network generators

A.1 Beta-splitting networks

Here, we describe how the beta-splitting network generator adds edges to a tree or existing network. Algorithms 1 and 3 are quite straightforward: they either choose two edges uniformly at random from the network, or two leaf-edges uniformly at random. In Algorithm 2, which adds the reticulations locally, the two edges are chosen as follows. First, one edge of the network is chosen uniformly at random, and then a random walk is performed through the network (with uniformly chosen steps). This random walk stops after each step with a stopping probability $P_{\text{stop}}$ . The higher $P_{\text{stop}}$ , the more local the reticulations are. If $P_{\text{stop}}$ is low,this random walk may take a very large number of steps. Hence, in our implementation, we also use a fixed maximal number of steps for this random walk.

Data: A network

N

Result: A network

N^{\prime}

obtained from

N

by adding one arc.

1 Set

N^{\prime}=N

;

2 Pick two edges

e_{1},e_{2}

N^{\prime}

uniformly at random;

3 Subdivide these edges with nodes

v_{1}

v_{2}

;

4 if $v_{2}$ is above $v_{1}$ then

5 Add the arc

(v_{2},v_{1})

N^{\prime}

;

7 else

8 Add the arc

(v_{1},v_{2})

N^{\prime}

;

10return

N^{\prime}

;

Algorithm 1 AddEdgeUniform

(N)

Data: A network

N

and a probability

p_{stop}

Result: A network

N^{\prime}

obtained from

N

by adding one arc.

1 Set

N^{\prime}=N

;

2 while True do

3 Pick an edges

e_{1}

N

uniformly at random;

4 Let

(u,v)

be a random orientation of

e_{1}

;

5 Let

\mathrm{stop\_value}=1

;

6 while $\mathrm{stop\_value}>p_{stop}$ do

7 Set

u=v

;

8 Let

v

be a random neighbour of

v

;

9 Sample

\mathrm{stop\_value}

uniformly at random from

[0,1]

;

11 Let

e_{2}

be the orientation of

(u,v)

as found in

N

;

12 if $e_{1}\neq e_{2}$ then

13 Quit the While-loop;

16Subdivide

e_{1}

and

e_{2}

with nodes

v_{1}

v_{2}

;

17 if $v_{2}$ is above $v_{1}$ then

18 Add the arc

(v_{2},v_{1})

N

;

20 else

21 Add the arc

(v_{1},v_{2})

N

;

23return

N

;

Algorithm 2 AddEdgeLocal

(T,p_{stop})

Data: A network

N

Result: A network

N^{\prime}

obtained from

N

by adding one arc.

1 Set

N^{\prime}=N

;

2 Pick two leaves

l_{1},l_{2}

N^{\prime}

uniformly at random;

3 Subdivide the incoming edges of these leaves with nodes

v_{1}

v_{2}

;

4 Add the arc

(v_{1},v_{2})

N^{\prime}

;

5 return

N^{\prime}

;

Algorithm 3 AddEdgeHorizontal

(N)

A.2 Heath networks

In this section, we give more details pertaining the implementation of our network generator based on the Heath branching process [24]. Algorithms 4 and 5 contain the pseudo-code for the model, and in the next subsection, we show that updating the distances as we do in Algorithm 5 is correct.

Data: A rate

r

, parameters

r_{\alpha},r_{\beta}

for the prior gamma distribution of

r

, and parameters

u_{\alpha},u_{\beta}

for the update gamma distribution.

Result: A new rate

r^{\prime}

1 Sample

u

from

\Gamma(u_{\alpha},u_{\beta})

;

2 Let

r^{\prime}=ur

be the proposed new rate;

3 With probability

\min(1,\operatorname{pdf}_{\Gamma(r_{\alpha},r_{\beta})}(r^{\prime})/\operatorname{pdf}_{\Gamma(r_{\alpha},r_{\beta})}(r^{\prime}))

, return

r^{\prime}

;

4 Otherwise, return

r

;

Algorithm 4 UpdateRate

(r,r_{\alpha},r_{\beta},u_{\alpha},u_{\beta})

A.2.1 Updating distances

As an example of an extension of a Markov model to a more flexible birth-hybridization model, we extend the model presented in [24], which uses auto-correlated speciation and extinction rates. To create networks, we add the options of including hybridization events. These events can be added to other Markov splitting models in the same way. So, if different or more realistic splitting models are desired (e.g., with speciation bursts) hybridization can easily be added to those as well.

Hybridization in this model should be thought of as hybrid speciation, and not hybrid introgression, that is, m-type reticulation model are used. Indeed, a hybridization event in this generator takes two extant taxa and combines them into a new taxon. The two parental taxa stay present in the network, so the number of taxa in the network increases by one.

In our model, the hybridization rate between two taxa depends on their distance. To keep it simple, we use the piece-wise linear function

f_{hyb}(d)=\begin{cases}h_{lr}&\mbox{if }d\leq h_{l}\\ h_{lr}+\frac{(h_{rr}-h_{lr})(d-h_{l})}{h_{r}-h_{l}}&\mbox{if }h_{l}<d<h_{r}\\ h_{rr}&\mbox{if }d\geq h_{r}\end{cases}.

Recomputing the total rates for speciation $\operatorname{sp}(N)=\sum_{x}\operatorname{sp}(x)$ and extinction $\operatorname{ext}(N)=\sum_{x}\operatorname{ext}(x)$ is quite simple, as it only requires one summation over the individual speciation and extinction rates $\operatorname{sp}(x)$ and $\operatorname{ext}(x)$ of all the extant taxa $x$ . Fully recomputing $\operatorname{hyb}(N)$ would entail finding all up-down paths between each pair of leaves. As there can be exponentially many up-down paths between a given pair of leaves, this full recomputation is not feasible.

Hence, to efficiently compute $\operatorname{hyb}(N)$ in our implementation, we keep track of the weighted distances between all pairs of leaves. Updating these distances requires few calculations for each of the types of events, as we will show next.

1 Let

N

be the tree on one leaf

l

;

2 Sample speciation rate

\operatorname{sp}(l)=\operatorname{sp}(N)

l

from

\Gamma(s_{\alpha},s_{\beta})

;

3 Sample extinction rate

\operatorname{ext}(l)=\operatorname{ext}(N)

l

from

\Gamma(s_{\alpha},s_{\beta})

;

4 Sample hgt rate

\operatorname{hgt}(l)=\operatorname{hgt}(N)

l

from

\Gamma(s_{\alpha},s_{\beta})

;

5 Set

\operatorname{hyb}(N)=0

;

6 Sample the next event time

t^{+}=

from an exponential distribution with parameter

1/(\operatorname{sp}(l)+\operatorname{ext}(l)+\operatorname{hgt}(l))

;

7 while $|\operatorname{Taxa}(N)|<$ max_taxa and $|\operatorname{Taxa}_{extant}(N)|>0$ do

8 Pick a random event by their weights

\operatorname{sp}(N),\operatorname{ext}(N),\operatorname{hgt}(N),

and

\operatorname{hyb}(N)

;

9 if speciation then

10 Pick a random extant leaf

l

, weighted by by their speciation rates;

11 Attach two new edges with length

0

to this leaf;

12 Update the rates of l to get rates for the new leaves;

13 Replace

l

with the two new leaves in the set of extant leaves;

15 else if extinction then

16 Pick a random extant leaf

l

, weighted by by their extinction rates;

17 Remove

l

from the set of extant leaves;

19 else if hgt then

20 Pick a random extant leaf

l

, weighted by by their hgt rates;

21 Pick another leaf

m

uniformly at random;

22 Add length-

0

edges

(m,l)

(l,l^{\prime})

, and

(m,m^{\prime})

, where

l^{\prime}

and

m^{\prime}

are new leaves;

23 Pick an inheritance probability

p

along

(m,l)

uniformly at random from

[0,hgt_{inh}]

;

24 Copy the rates of

m

m^{\prime}

;

25 Set the rates of

l^{\prime}

to be

pr(m)+(1-p)r(l)

, and then update these;

26 Replace

l

and

m

with

l^{\prime}

and

m^{\prime}

in the set of extant taxa;

28 else if hybridization then

29 Pick a random pair of extant leaves

(l,m)

, weighted by their hybridization rates

f_{hyb}(d(l,m))

;

30 Add length-

0

edges

(m,x)

(l,x)

(l,l^{\prime})

(m,m^{\prime})

, and

(x,q^{\prime})

, where

l^{\prime}

q^{\prime}

, and

m^{\prime}

are new leaves;

31 Pick an inheritance probability

p

along

(m,x)

uniformly at random from

[0,1]

;

32 Copy the rates of

m

and

l

m^{\prime}

and

l^{\prime}

;

33 Set the rates of

q^{\prime}

to be

pr(m)+(1-p)r(l)

, and then update these;

34 Replace

l

and

m

with

l^{\prime}

q^{\prime}

, and

m^{\prime}

in the set of extant taxa;

36 Extend the length of all pendant arcs leading to extant taxa by

t^{+}

;

37 Recompute

\operatorname{sp}(N)

\operatorname{ext}(N)

\operatorname{hgt}(N)

, and

\operatorname{hyb}(N)

;

38 Sample

t^{+}=

from an exponential distribution with parameter

1/(\operatorname{sp}(N)+\operatorname{ext}(N)+\operatorname{hgt}(N)+\operatorname{hyb}(N))

;

40return

N

;

Algorithm 5 Heath

(s_{\alpha},s_{\beta},e_{\alpha},e_{\beta},hgt_{\alpha},hgt_{\beta},hgt_{inh},u_{\alpha},u_{\beta},h_{l},h_{r},h_{lr},h_{rr})

We will now show how to efficiently update the distances between all pairs of leaves efficiently after each of the following events: speciation, extinction, horizontal gene transfer/hybrid introgression, and hybrid speciation. In all these events, we assume all new edges have length 0 and are extended only later.

Speciation

Suppose the speciation event splits the extant taxon $l$ into two new taxa $m$ and $m^{\prime}$ . As the length of both the two new edges is zero, the distance between $m$ and $m^{\prime}$ is 0. For any other leaf $q$ , we get $d(m,q)=d(m^{\prime},q)=d(l,q)$ ; and for each pair of leaves not involving $m$ or $m^{\prime}$ , the distance is not changed by the speciation event.

Extinction

An extinction event does not affect the distances between pairs of extant taxa. The only change is that the set of extant taxa changes, i.e., in the set of pairs that we compute the hybridization rate for.

HGT

Suppose the HGT event uses $l$ as accepting taxon, and $m$ as donating taxon. Adding the arc $(m,l)$ only affects the distances between $l$ and other leaves. Indeed, for pairs of leaves where $l$ is not involved, the set of up-down paths is not affected by adding $(m,l)$ . If the inheritance probability along $(m,l)$ is $p$ , then the new distances involving $l$ are $d(l,q)=pd(m,q)+(1-p)d(l,q)$ .

Hybridization

Suppose this event adds a hybrid $q$ between $l$ and $m$ . The new arcs $(l,q)$ and $(m,q)$ can only be part of up-down paths between $q$ and another leaf. Hence, the only affected distances involve $q$ and another leaf. If the inheritance probability along $(m,q)$ is $p$ , then the new distances involving $q$ are $d(q,x)=pd(m,x)+(1-p)d(l,x)$ .

Merging

Suppose this event merges $l$ and $m$ into a reticulation $r$ and introduces the new edge $(r,q)$ (of length 0). Note that the only affected distances involve $q$ and another leaf. Let $p_{l}$ and $p_{m}$ be the parents of $l$ and $m$ , and let the inheritance probability along $(p_{l},r)$ be $p$ , then the new distances involving $q$ are $d(q,x)=pd(l,x)+(1-p)d(m,x)$ .

Extending pendant arcs

By extending the pendant arcs of the extant taxa by $t^{+}$ , each up-down path between a pair of extant taxa becomes $2t^{+}$ longer. Hence, the weighted sums of the distances along these up-down paths also become $2t^{+}$ longer.

Appendix B Network generators and classes

In the following lemmas, we prove that certain processes inherently lead to certain classes of networks, and other are able to produce any network.

Lemma B.1.

Let $N$ be an arbitrary network topology, then using merging (adding edges y-type) and speciation can generate $N$ .

Proof.

We prove this lemma by induction on the number of reticulations. First, for the basis of the induction, note that each network with 0 reticulations (i.e., each tree) can be generated by only using speciation events.

Now suppose that each network with at most $r$ reticulations can be generated, and let $N$ be an arbitrary network with $n$ leaves and $r+1$ reticulations. If $N$ has a leaf that is directly below a reticulation node. Removing this leaf, and splitting its parent into two nodes, one for each incoming arc, we get a network $N^{\prime}$ with $n+1$ leaves and $r$ reticulations from which $N$ can be obtained by one merging event. By the induction hypothesis, $N^{\prime}$ can be generated and, therefore, so can $N$ .

If $N$ contains no leaf directly below a reticulation, then there must be a reticulation in $N$ with a pendant subtree directly below it. Let $N^{\prime}$ be the network obtained from $N$ by removing this subtree. It is clear that $N$ can be obtained from $N^{\prime}$ using only speciation events, as each tree can be generated using speciation events. Moreover, by the arguments in the previous paragraph, $N^{\prime}$ can be generated as it has $r+1$ reticulations and a leaf directly below a reticulation. We conclude that any network with $r+1$ reticulations can be produced using merging and speciation events. By induction, this implies that all networks can be generated thusly. ∎

Lemma B.2.

Using HGT (adding reticulations n-type at the leaves) together with speciation generates only orchard networks.

Proof.

Note that speciation and adding reticulations using the n-type are exactly the methods to “add a pair to a network” in the sense of [33]. The lemma now follows immediately, because these additions of pairs can be reversed, which gives a cherry-picking sequence for the network. ∎

Lemma B.3.

Using hybridization (adding reticulations m-type at the leaves) together with speciation generates only temporal tree-child networks, and each (binary) temporal tree-child network can be generated this way.

Proof.

The first part follows because we can label each internal node with the number of the step (event) it was introduced in. This gives a temporal labelling of the network. Moreover, the such a network must be tree-child because an m-type reticulation event can never lead to two neighbouring reticulations or a tree node with two reticulation children. For the second part, note that the temporal labelling can be modified so that for each time $t$ , the internal nodes with label $t$ either represent a speciation event, or an m-type hybridization event. ∎

Lemma B.4.

HGT (n-type) or hybridization (m-type) together with extinction and speciation can generate any network.

Proof.

This follows directly from Lemma B.1, because each merging event can be substituted for a HGT event and one extinction event, or a hybridization event and two extinction events. ∎

Appendix C Supplementary figures and tables

In this section, we provide additional data in the form of tables and figures. In Table 1, for each generator, we give an estimate for the number of reticulations at which half of the generated networks are in each network class. This value is found by finding the intersection of the horizontal line where the fraction is a half with the lines in Figure 5. In Table 2, we show the average balance for each of the generators and each number of reticulations. The stars show that adding reticulations significantly impacts the balance of the networks compared to the networks with 1 reticulation for the same network generator. Figure 1 shows the profiles for the number of (reticulated) cherries in the beta-splitting model. Note that the initial number of cherries depends primarily on the beta parameter, while all profiles, except when edges are added horizontally, are qualitatively the same. Figure 2 shows the dependence of the number of blobs and their sizes for the beta-splitting model. Note that the profiles are qualitatively the same, although the exact values differ slightly, as a result of the underlying tree shape.

Table 1: The number of reticulations at which half the networks are in a class for the different generators and parameters. Obtained by interpolation from the data with log-transformed number of reticulations, because in those data, the curves look more or less sigmoid (Figure 5), and so that we do not underestimate the effects at small numbers of reticulations.

			Class
Generator	Parameter	Locality	Tree-based	Orchard	Stack-free	Tree-child
Beta-splitting	Beta=0.0	Uniformly	41.2	24.0	10.1	6.6
		P_stop=0.002	41.6	21.4	11.0	7.0
		P_stop=0.2	33.9	9.4	8.9	6.9
		Horizontal			12.6	7.6
Zhang Ntk			32.7	25.2	9.1	6.2
LGT Generator		Locality_0.01			10.5	6.9
		Locality_0.1			10.4	6.6
		Locality_1.0			12.5	7.4
ZODS			33.5	28.2	9.8	6.2
Heath	keep extinct	less local
		local
	remove extinct	less local	?	?	26.8	13.9
		local	?	?	24.6	13.9
	no extinction	less local
		local
MCMC			33.1	19.7	8.9	6.5

Table 2: The average balance of networks for different network generators. Stars indicate whether the balance is significantly different from the networks with 1 reticulation of the same generator (

*:p<0.01

**:p<0.001

		Reticulations
Generator	Parameters	1	2	5	10	20	50	100	200	500
Beta-splitting beta=0.0	Uniform	0.771	0.769	0.762	0.760	0.745**	0.730**	0.710**	0.706**	0.712**
	P_stop=0.002	0.770	0.769	0.763	0.757	0.744**	0.722**	0.707**	0.699**	0.707**
	P_stop=0.2	0.770	0.768	0.761	0.752*	0.737**	0.705**	0.671**	0.635**	0.587**
	Horizontal	0.772	0.773	0.774	0.778	0.782	0.796**	0.815**	0.848**	0.894**
Zhang_Ntk		0.779	0.778	0.773	0.769	0.756**	0.744**	0.747**
LGT_Generator	Beta=0.01	0.780	0.779	0.777	0.769	0.760**	0.755**	0.767*	0.796*	0.848**
	Beta=0.1	0.782	0.784	0.779	0.778	0.771	0.782	0.807**	0.845**	0.876**
	Beta=1.0	0.781	0.779	0.777	0.771	0.781	0.792	0.806**	0.847**	0.884**
ZODS		0.779	0.779	0.776	0.778	0.776	0.775	0.774	0.787
Heath+extinction	less local	0.506	0.512	0.521	0.505	0.526	0.547**
	local	0.506	0.510	0.521	0.508	0.525	0.526
Heath+extinction+remove extinct	less local	0.732	0.727	0.713*	0.736	0.700**	0.693**
	local	0.722	0.722	0.716	0.715	0.687**	0.639**
Heath	less local	0.749	0.737	0.730**	0.742	0.716**	0.711**
	local	0.736	0.734	0.736	0.732	0.703**	0.646**
MCMC		0.442	0.444	0.449	0.461**	0.478**	0.510**	0.540**	0.569**	0.579**

C.1 Heath Parameters and Extinction

In the main text of this paper, we have not considered extinction. Even though our implementation of the Heath network generator has an option to add extinction. In Figures 3, 4, 5, and 6, we show the effect of extinction on the profiles, when using the parameters for the Heath network generator as detailed in Table 3. We have two methods of implementing extinction: simply generating a network where the desired number of taxa consists of both extinct and extant taxa; or generating a network with extinct taxa, and removing them to obtain a network with the desired number of extant taxa. When extinct networks are removed, Table 1 does not show reticulation numbers at which half of the generated networks are orchard or tree-based. This is because, for the reticulation numbers we could test for this model, none gave a fraction of less than half for these classes (Figure 3). Perhaps, this is because our networks have a small number of extinct taxa, and removing extinct taxa is necessary to produce networks that are not orchard or tree-based when adding reticulations as hybrid speciation.

Table 3: Parameters used in the Heath network generator

Extinction	No						Yes - remove						Yes - keep
Reticulations	$[1,2,3,4,5]$		$[10,15,20]$		$[30,40,50]$		$[1,2,3,4,5]$		$[10,15,20]$		$[30,40,50]$		$[1,2,3,4,5]$		$[10,15,20]$		$[30,40,50]$
Locality	local	less local	local	less local	local	less local	local	less local	local	less local	local	less local	local	less local	local	less local	local	less local
update shape	3.0	3.0	3.0	3.0	3.0	3.0	3.0	3.0	3.0	3.0	3.0	3.0	3.0	3.0	3.0	3.0	3.0	3.0
speciation rate mean	1.0	1.0	1.0	1.0	1.0	1.0	1.0	1.0	1.0	1.0	1.0	1.0	1.0	1.0	1.0	1.0	1.0	1.0
speciation rate shape	5.0	5.0	5.0	5.0	5.0	5.0	5.0	5.0	5.0	5.0	5.0	5.0	5.0	5.0	5.0	5.0	5.0	5.0
extinction rate mean	0.0	0.0	0.0	0.0	0.0	0.0	0.5	0.5	0.5	0.5	0.5	0.5	0.5	0.5	0.5	0.5	0.5	0.5
extinction rate shape	0.0	0.0	0.0	0.0	0.0	0.0	5.0	5.0	5.0	5.0	5.0	5.0	5.0	5.0	5.0	5.0	5.0	5.0
$h_{l}$	0.3	1.0	0.3	1.0	0.3	1.0	0.3	1.0	0.3	1.0	0.3	1.0	0.3	1.0	0.3	1.0	0.3	1.0
$h_{r}$	1.5	5.0	1.5	5.0	1.5	5.0	1.5	5.0	1.5	5.0	1.5	5.0	1.5	5.0	1.5	5.0	1.5	5.0
$h_{lr}$	0.08	0.015	0.3	0.05	0.6	0.125	0.06	0.01	0.25	0.05	0.6	0.11	0.08	0.015	0.3	0.075	0.75	0.15
$h_{rr}$	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0