Beyond Level-1: Fast Inference of Generic Semi-directed Phylogenetic Networks

1. Problem Statement

Phylogenetic networks are essential for modeling reticulate evolutionary events such as hybridization, introgression, and horizontal gene transfer, which cannot be represented by standard phylogenetic trees. However, existing inference methods face two major limitations:

• Topological Restriction: Most widely used methods, including the popular SNaQ (Species Networks applying Quartets) algorithm, are restricted to Level-1 networks. In Level-1 networks, every biconnected component contains at most one hybrid node. This restriction prevents the modeling of complex evolutionary scenarios where reticulation events are not isolated (e.g., nested or overlapping hybridization).
• Computational Bottleneck: Inferring arbitrary networks under the Multispecies Network Coalescent (MSNC) model is computationally prohibitive. The likelihood equations for arbitrary networks are complex and require recursive computation of expected quartet concordance factors (qCFs), leading to a drastic increase in computational burden compared to tree estimation or Level-1 network estimation.

2. Methodology

The authors present an extended version of the SNaQ method (implemented in SNaQ.jl) that overcomes these limitations through algorithmic and optimization improvements.

A. Algorithmic Extensions

• Generalized Composite Likelihood: The method extends the composite likelihood framework to arbitrary binary, metric, semi-directed phylogenetic networks, removing the Level-1 constraint.
• Recursive qCF Computation: To compute expected quartet concordance factors ($E_N(X_{qi})$) for arbitrary networks, the authors implemented a modified version of a recursive algorithm (based on [10]).
  • Optimization: They introduced a computation graph approach. Instead of repeatedly copying and reducing the network for every quartet, they construct a directed, multifurcating graph that stores intermediate values. This allows for the calculation of expected qCFs and their gradients in a single walk of the graph.
  • Tree-like Quartets: An immediate speedup is achieved by first assuming quartets are tree-like; if this fails, the recursive algorithm is used as a fallback.

B. Gradient-Based Optimization

• Manual Forward Differentiation: The original SNaQ used gradient-free optimization (BOBYQA). The new implementation manually computes the gradients of the composite log-likelihood with respect to network parameters (edge lengths $t$ and inheritance proportions $\gamma$).
• L-BFGS Algorithm: By utilizing these efficient gradients, the authors switch to the L-BFGS (Limited-memory Broyden–Fletcher–Goldfarb–Shanno) algorithm. This results in substantial speedups in parameter optimization, enabling the inference of larger and more complex networks.

C. Flexible Search Spaces

The framework allows users to restrict the search space to specific topological classes. The authors focus on Tree-Child and Galled (TCG) networks:

• Tree-Child: Every hybrid node is a parent to at least one tree node in every rooting.
• Galled: Every cycle contains at most one hybrid node.
• Identifiability: Recent theoretical results indicate that TCG networks have identifiable topologies under the MSNC, making them a robust target for inference. The tool also supports unrestricted searches or user-defined topological constraints.

3. Key Contributions

1. Scalable Inference Beyond Level-1: The first implementation capable of inferring arbitrary semi-directed networks under a composite likelihood framework, moving beyond the Level-1 restriction.
2. Computational Efficiency: The combination of recursive qCF computation with a computation graph and gradient-based optimization (L-BFGS) yields significant runtime reductions, making genome-scale studies feasible.
3. Robustness Analysis: A comprehensive simulation study evaluating the method's performance when the true network lies inside vs. outside the TCG search space.
4. Empirical Application: Re-analysis of the Xiphophorus (Poeciliidae) phylogeny, demonstrating that relaxing the Level-1 constraint reveals a more complex and biologically accurate history of hybridization.

4. Results

A. Simulation Study

• Topological Accuracy:
  • When the true network is within the TCG space, SNaQ.jl reliably recovers the topology, with accuracy improving as the number of taxa, loci, and sequence length increases.
  • Data Limitation Warning: With limited data (e.g., 100 gene trees), networks with poor topological accuracy (high Unrooted Hardwired Cluster Distance, UHWCD) often achieve better Negative Composite Log-Likelihood (NCLL) scores than the true generating network. This suggests that with sparse data, the model may fit noise rather than the true topology.
  • Model Violations: When the true network is not TCG (e.g., TCNG or NTCNG), the method cannot perfectly recover the topology (since the search space is restricted to TCG). However, it still recovers meaningful information about hybridization events.
• Parameter Estimation:
  • Edge length ($t$) and inheritance proportion ($\gamma$) estimates are generally accurate when the topology is correct.
  • However, $\gamma$ estimates show high variance when the topology is incorrect, even if the median error is low.
• Hybrid Descendant Recovery: Even when the full topology is incorrect, the method accurately identifies the presence and direction of hybridization events (high F1 scores) for TCG networks. For non-TCG networks, precision remains high, but recall drops (the method identifies some hybrids but misses others).

B. Empirical Analysis: Xiphophorus (Poeciliidae)

• Search Strategies: The authors compared four search spaces: Level-1 (L1), TCG, Unrestricted (U), and TCG initialized from Unrestricted (TCGU).
• Model Fit: Networks inferred in TCG, U, and TCGU spaces yielded significantly better NCLL scores than previous Level-1 models, particularly when the number of hybrid nodes ($h$) was $\ge 2$.
• Biological Insights:
  • The best-fitting model (selected via DDSE model selection) was a Level-2 network with 4 hybrid nodes, found in the TCGU space.
  • This model revealed hybridization events previously undetected by Level-1 constraints, including ancient hybridization between X. milleri and X. xiphidium, and contemporary hybridization among X. malinche and X. birchmanni.
  • Interestingly, despite searching unrestricted spaces, the model selection criteria consistently favored Tree-Child and Galled topologies, suggesting the true evolutionary history of Xiphophorus may naturally fit this class.

5. Significance

This work represents a major leap forward in phylogenetic network inference. By breaking the Level-1 barrier, the authors enable the reconstruction of richer reticulate histories from genomic data. The computational optimizations make it possible to apply these complex models to genome-scale datasets, addressing critical biological questions regarding hybridization and introgression that were previously intractable. The study also provides crucial caveats regarding data requirements, highlighting that while the method is powerful, large datasets are necessary to distinguish between competing complex topologies.

Keywords: Phylogenetic networks, Composite likelihood, SNaQ.jl, Tree-child, Galled networks, Hybridization, MSNC.