PhyloGen: Language Model-Enhanced Phylogenetic Inference via Graph Structure Generation

Authors: chenrui  duan
zelin  zang
siyuan  li
yongjie  xu
stan z li
Venue: NeurIPS_cc_2024_Conference
Published: 2024-04-29 10:29:33

Abstract: Phylogenetic trees elucidate evolutionary relationships among species, but phylogenetic inference remains challenging due to the complexity of combining continuous (branch lengths) and discrete parameters (tree topology). 
    Traditional Markov Chain Monte Carlo methods face slow convergence and computational burdens. Existing Variational Inference methods, which require pre-generated topologies and typically treat tree structures and branch lengths independently, may overlook critical sequence features, limiting their accuracy and flexibility.
    We propose PhyloGen, a novel method leveraging a pre-trained genomic language model to generate and optimize phylogenetic trees without dependence on evolutionary models or aligned sequence constraints. PhyloGen views phylogenetic inference as a conditionally constrained tree structure generation problem, jointly optimizing tree topology and branch lengths through three core modules: (i) Feature Extraction, (ii) PhyloTree Construction, and (iii) PhyloTree Structure Modeling. 
    Meanwhile, we introduce a Scoring Function to guide the model towards a more stable gradient descent.
    We demonstrate the effectiveness and robustness of PhyloGen on eight real-world benchmark datasets. Visualization results confirm PhyloGen provides deeper insights into phylogenetic relationships.

Main paper: 
 introduction: 
 Phylogenetic trees [42] (or evolutionary trees) are tree-structured graphs representing kinship relationships between species [11] , where each leaf node represents a distinct species and internal nodes represent evolutionary bifurcations. The tree's topology reflects evolutionary relationships based on their genetic characteristics, and branch lengths indicate evolutionary distances. Phylogenetics study is the foundation of evolutionary synthetic biology and is critical for tracking the evolutionary trajectories of species, analyzing the transmission pathways of newly discovered viruses, and providing meaningful insights for clinical applications [25, 13, 33] .Background. Despite DNA sequencing technologies [32] allowing us to construct evolutionary relationships based on molecular properties (DNA, RNA, and proteins), phylogenetic inference remains challenging due to its complex parameter space, which encompasses both continuous (branch lengths) and discrete (tree topology) components. Traditional MCMC-based methods like MrBayes [29] and RevBayes [9] , despite their ability to extensively explore the huge tree space, are hindered by slow convergence rates. Additionally, the number of possible tree topologies for n species grows factorially as (2n-5)!! for n ≥ 3, posing huge computational challenges. In contrast, VI-based methods leverage approximate inference to provide more efficient estimations. Depending on the data type and research objectives, these methods can be further categorized into tree representation learning and tree structure generation. Tree Representation Learning, as shown in Fig. 1(a) , intends to learn topological representations from existing tree structures. For instance, Subsplit The inputs are aligned sequences, and then tree structures and branch lengths are directly inferred by variational inference and biological modules. These methods optimize tree topology and branch lengths separately. (c) The inputs are raw sequences processed by a pre-trained language model to generate species representations. Then, an initial topology is generated through a tree construction module, and the topology and branch lengths are co-optimized by the tree structure modeling module.Bayesian Networks (SBNs) [45] focus on probabilistic representations from given tree structures without considering branch lengths. Variational Bayesian Phylogenetic Inference (VBPI) [46] and its extensions VBPI-GNN [44] utilize the topological probabilities provided by SBNs and jointly model tree structures in latent space by deriving branch lengths through variational approximations. However, these methods require pre-generated topology, based on which better topological representations are then learned. Tree Structure Generation, as shown in Fig. 1 (b), aim to infer tree structures directly from sequences. PhyloGFN [48] integrates VI and reinforcement learning to construct tree topologies, simplifying branch lengths into discrete intervals, thus necessitating posterior data for inference. VaiPhy [16] introduces the SLANTIS sampling strategy and basic biological models (e.g., Jukes-Cantor(JC) model [24] ) for topology and branch length estimation. ARTree [36] builds tree topologies recursively using a graph autoregressive model. GeoPhy [22] models tree topologies in a continuous geometric space. However, these methods require input sequences to be aligned to equal lengths, and tree topology and branch lengths are optimized separately. They also ignore sequence features that are critical for accurate phylogenetic analysis.Our Method. As depicted in Fig. 1(c ), we propose a novel approach based on a pre-trained genome language model. Our model does not rely on evolutionary models or the requirement to align input sequences to equal lengths and fully exploits the prior knowledge embedded in biological sequences. PhyloGen models phylogenetic tree inference as a conditional-constrained tree structure generation problem, aiming to generate and optimize the tree topology and branch lengths jointly. We map species sequences into a continuous geometric space and perform end-to-end variational inference without restricting topological candidates. To ensure the topology-invariance of phylogenetic trees, we incorporate distance constraints in the latent space to maintain translational rotation invariance. Our approach demonstrates effectiveness and efficiency on the eight real-world benchmark datasets and verifies its robustness through data augmentation and noise addition [41] . In addition, we propose a new scoring function to guide the model towards a more stable and faster gradient descent. 
related works: 
 MCMC-based methods, such as MrBayes [29] and RevBayes [9] , have been widely used for phylogenetic inference due to their ability to explore vast tree spaces.VI-based methods offer a more efficient alternative to MCMC by leveraging approximate inference techniques. These methods can be categorized into two main approaches: structure representation and structure generation. Tree Representation Learning. This approach focuses on extracting information from existing tree structures. SBNs [45] capture the relationships between existing subsplits without addressing branch lengths. VBPI [46] and its variants VBPI-NF [43] and VBPI-A. Feature Extraction We utilize DNABERT2 [49] , a genome language model, to transform molecular (DNA) sequences Y = {y i } N i=1 into genomic embeddings E = {e i } N i=1 . These embeddings discern complex patterns and long-range dependencies, serving as the basis for generating the latent variables z * , which dynamically inform both the topology and the branch lengths. By introducing DNABERT2, we redefine the construction of phylogenetic trees as a continuous optimization problem within a biologically meaningful embedding space.To construct the phylogenetic tree, genomic embeddings E are input into an MLP network to derive the parameters µ and σ of the latent space, representing topological embeddings z * . This latent space effectively captures the evolutionary relationships that guide the subsequent tree structure generation process. The latent variable is then generated using the reparameterization trick [15] :z * = µ + σ ⊙ ε, ε ∼ N (0, I). Then, distance matrix D is computed using z * : D(i, j) = N i,j=1 z * i ⊕ z * j ,where ⊕ means an XOR operation reflecting the nucleotide mismatches. This distance matrix is fed into the Neighbor-Joining (NJ) algorithm [31] to generate the initial tree topology τ (z * ). The tree's topology is structured using a parent-child relationships list[L[i], R[i], P [i]],where P [i] denotes the parent of node i, and L[i] and R[i] denote the left and right child node. This representation highlights the tree's hierarchical nature, optimizing its structure based on the evolutionary patterns in the data.C. PhyloTree Structure Modeling The purpose of the tree structure modeling module is to jointly optimize both the phylogenetic tree's discrete topology and continuous branch lengths.C.1. Topology Learning. The first component involves a TreeEncoder R(z|τ (z * )) and a TreeDecoder Q(τ (z)|z) to learn the tree topology. Concretely, the encoder R receives an initial tree topology τ (z * ) and genomic embeddings E from DNABERT2 as inputs, conditions the latent state z to refine topological embeddings in a continuous space. We introduce variational inference to enhance the quality of the tree structure and adapt to data uncertainty and complexity by minimizing the Kullback-Leibler (KL) divergence. Additionally, the encoder R acts as a regularization mechanism, reducing overfitting and enhancing gradient stability. To refine the topology, the decoder Q samples from the probability distribution parameterized by the encoder's latent state. This process optimizes model parameters by minimizing the KL divergence between the true distribution P (τ (z)) = P (y|τ (z), B τ ) and the variational distribution Q(τ (z)) = q(τ (z), B τ ). To facilitate backpropagation through stochastic nodes, the latent variable z is sampled using the reparameterization trick: z = µ + σ ⊙ ε, where ε ∼ N (0, I) introduces controlled stochasticity while maintaining differentiability, ensuring that z effectively captures refined topological embeddings.C.2. Branch Length (Blens) Learning. The second component utilizes the inferred topology τ and genomic embeddings E from DNABERT2 to generate and adjust branch lengths.Step1: Node Feature via Linear-Time Dual-Pass Traversal. We utilize a linear time O(n) dual-pass traversal method, combining postorder (bottom-up) and preorder (top-down) strategies, guaranteeing each node is processed only once.Postorder Traversal (bottom-up) aggregates information from leaf nodes toward the root node.EQUATIONwhere c i is the scaling factor for node i, influenced by the node's connectivity, and K = 3 due to binary tree properties.EQUATIONwhere f i incorporates contributions from children nodes f j and its own initial feature e i from a pre-trained genome model. Initially,c i = 0, f i = e i .Preorder Traversal (top-down) propagates information from the root node toward the leaf nodes, enhancing each node's feature x i :EQUATIONwhere e p[i] is the feature of node i's parent.Step2: Feature Enhancement with Dynamic Graph Convolution (DGCNN) : We use the edge convolutional layer of DGCNN [35] to enhance node features. This approach captures both local interactions (between neighboring nodes) and global structural dependencies (reflecting the entire tree structure) within the phylogenetic tree. For layer L, inputs{x L i ∈ R F } N i=1 and outputs {x L+1 i ∈ R F ′ } Ni=1 have feature dimensions F = 768 and F ′ = 100. The transformations are:x L+1 i = DGCNN(x L i ), m ij = h(x L i , x L j ),where h : R F × R F → R F ′ is an MLP-shared asymmetric edge function.Topology-Invariance Property. Despite the non-uniqueness of 2D coordinates due to infinite possible translational rotations, the distances between species remain constant. We define this as topological invariance:d 2 ij = ∥z i -z j ∥ 2, where z i and z j are coordinates in the hidden space, maintaining accurate topological relationships. The function h(x i , x j ) incorporates global and local features through: h(x i , x j ) = h(x i , x i -x j , d 2 ij ). Aggregation of Features. The node features are aggregated using a MAX operation across all edges:EQUATIONwhere ⊕ denotes the MAX aggregation function, focusing on the most significant features.Step3: Reparametrization of Branch Length. Node features x L i obtained from edge convolution layers are further fed into the M network, parameterizing mean and log-variance for branch lengths:EQUATION(1)EQUATIONwhere h i are node features processed through an additional MLP 1 to capture inter-node interactions. Branch lengths are updated by reparametrization b = exp(µ + exp(σ) * rvs), where rvs ∼ N (0, I), ensuring differentiability and capturing the probabilistic nature of estimates. Throughout, the prior P (B τ ) is assumed exponential, while the posterior Q(B τ ), learned via a graph neural network, reflects inferred tree topologies and node characteristics.  
learning objectives: 
 As discussed in Appendix Background A, our primary goal is to maximize the expected marginal likelihood of the observed species sequence Y via max log p(Y |(τ (z), B τ )). The posterior distribution as p(τ (z), B τ |Y ) is difficult to infer directly, so we utilize variational inference to approximate it as q(τ (z), B τ |Y ). The detailed deviation is in Appendix D.2.To minimize the KL divergence between the true prior and the approximate posterior distributions, we start from the joint probability distribution:p(Y, τ (z), B τ ) = p(Y |τ (z), B τ )p(B τ |τ (z))p(τ (z)),where p(Y |τ (z), B τ ) represents the conditional probability of the observed data Y . We assume the tree topology τ (z) and the branch lengths B τ are conditionally independent.We introduce a variational distribution q(τ (z), B τ ) = q(B τ |τ (z))q(τ (z)) to approximate the true posterior p(τ (z), B τ |Y ). The ELBO loss can initially be written as:EQUATIONTo improve training variance and gradient stability, we introduce a regularization term R(z|τ (z * )) and reformulate the ELBO loss as:EQUATIONFor better performance and reduced variance, we adopt a multi-sample approach [23] :EQUATIONwhere z * k and B k τ are samples from the variational distributions Q(z * ) and Q(B τ |τ (z)), respectively. K represents the number of Monte Carlo samples from the variational distribution to compute the empirical average, providing a more robust approximation of the ELBO. Typically, K is chosen to be moderate to balance between computational efficiency and estimation accuracy [18] . And the multi-sample ELBO reduces to the standard ELBO formula when K=1.The total loss L total is then defined as:EQUATIONwhereL KL = -KL{q ϕ (τ (z)|y)||p θ (τ (z))} -KL{q ϕ (B τ |y)||p θ (B τ )}, minimizing the KL divergence, and L(S) = -N i=1 f (z i ) accounts for the scoring function loss of the embeddings.We derive the gradients of model parameters θ as follows:EQUATIONwhereH[Q θ (z * )] is the entropy of Q θ (z * ).The derivation of the model parameters ϕ:EQUATIONThe derivation of the model parameters ψ:EQUATION 
experiments: 
 In this section, we conduct experiments to demonstrate the effectiveness of our proposed PhyloGen.We aim to answer seven research questions as follows:RQ1: How effective is PhyloGen in generating tree structures under the benchmark datasets? RQ2: How diverse are the tree topologies generated by PhyloGen? RQ3: How consistent is the PhyloGen-generated tree structure compared to the MrBayes method? RQ4: How robust is PhyloGen to species sequences? RQ5: How does each PhyloGen's module affect its performance? RQ6: What evolutionary relationships between species does PhyloGen learn? RQ7: How do key hyper-parameters affect PhyloGen's performance? The closer the two curves are, the better, which suggests that our method is highly consistent with the gold standard MrBayes approach.performance improves in the later stages but exhibits large fluctuation. GeoPhy performs the worst, with consistently the lowest and more fluctuating ELBO values. The right plot of Fig. 4 demonstrates our model's advantages in the MLL metric, where it rapidly achieves and maintains high-performance levels. In contrast, ARTree and GeoPhy have lower MLL values, especially GeoPhy, which has the weakest performance throughout. To evaluate the topological diversity of trees generated by PhyloGen on DS1, we use three metrics: Simpson's Diversity Index [5] , Top Frequency, and Top 95% Frequency, as detailed in Tab. 3. A higher Diversity Index, which approaches 1, suggests broad diversity among generated tree topologies. The lower Top Frequency suggests a balanced distribution, preventing single tree structures from being overly dominant. Furthermore, the presence of 149 distinct topologies within the Top 95% Frequency underscores PhyloGen's ability to generate a diverse range of topologies.Fig. 5 shows the bipartition frequency distributions of trees inferred by PhyloGen for datasets DS1, DS2, and DS3. The horizontal axis indicates the ranking of the bipartitions in the tree topology, and the vertical axis indicates the normalized frequency of occurrence of the corresponding bipartitions. The similarity of our method's curves to those of MrBayes underscores its accuracy, demonstrating that PhyloGen consistently captures evolutionary patterns with reliability comparable to the gold standard. This indicates a robust validation of PhyloGen's phylogenetic inference capabilities. More detailed information is provided in Appendix E.3. Tab. 5 shows the performance of our model compared with the removal of the KL loss and the Scoring Function S, respectively. Ours performs best on ELBO and MLL metrics, with a significant decrease in performance after the removal of the KL loss, suggesting that the KL loss plays a key role in regularizing the model and avoiding overfitting. While removing the S module had a large impact on MLL, the ELBO impact was relatively small, indicating that the impact of the S module is more complex and may be related to specific feature extraction functions. Our future work will explore further enhancements to the S module and investigate other regularization techniques to refine the model's performance. Fig. 7 shows a phylogenetic tree constructed from DS1 dataset, where each leaf node represents a specific species, and the text next to the node indicates the species name. The branch lengths reflect the genetic distances, with shorter branches indicating recent evolutionary history and longer branches indicating greater genetic differences. The phylogenetic tree shown in Fig. 7 places Siren intermedia and Trachemys scripta, both aquatic organisms, on adjacent branches, reflecting our model could capture their adaptive evolutionary information to aquatic environments. Meanwhile, the reptilian species Heterodon platyrhinos and Trachemys scripta are also on neighbouring branches, suggesting a relatively recent common ancestor compared to other amphibian and reptilian species. Notably, our method does not require padding the sequences to a unified length, which effectively reflects actual sequence variation.For a detailed species sequence and topology comparison, please refer to Appendix E.6.Question: Does the paper fully disclose all the information needed to reproduce the main experimental results of the paper to the extent that it affects the main claims and/or conclusions of the paper (regardless of whether the code and data are provided or not)? Answer: [Yes] Justification: We detail dataset information, network architecture, various hyperparameters, and algorithms in the main text and appendices. Guidelines:• The answer NA means that the paper does not include experiments.• If the paper includes experiments, a No answer to this question will not be perceived well by the reviewers: Making the paper reproducible is important, regardless of whether the code and data are provided or not. • If the contribution is a dataset and/or model, the authors should describe the steps taken to make their results reproducible or verifiable. • Depending on the contribution, reproducibility can be accomplished in various ways.For example, if the contribution is a novel architecture, describing the architecture fully might suffice, or if the contribution is a specific model and empirical evaluation, it may be necessary to either make it possible for others to replicate the model with the same dataset, or provide access to the model. In general. releasing code and data is often one good way to accomplish this, but reproducibility can also be provided via detailed instructions for how to replicate the results, access to a hosted model (e.g., in the case of a large language model), releasing of a model checkpoint, or other means that are appropriate to the research performed. • While NeurIPS does not require releasing code, the conference does require all submissions to provide some reasonable avenue for reproducibility, which may depend on the nature of the contribution. For example (a) If the contribution is primarily a new algorithm, the paper should make it clear how to reproduce that algorithm. (b) If the contribution is primarily a new model architecture, the paper should describe the architecture clearly and fully. (c) If the contribution is a new model (e.g., a large language model), then there should either be a way to access this model for reproducing the results or a way to reproduce the model (e.g., with an open-source dataset or instructions for how to construct the dataset). (d) We recognize that reproducibility may be tricky in some cases, in which case authors are welcome to describe the particular way they provide for reproducibility.In the case of closed-source models, it may be that access to the model is limited in some way (e.g., to registered users), but it should be possible for other researchers to have some path to reproducing or verifying the results. Justification: We detail the network architecture of the model, the various hyperparameter settings and the algorithmic pseudocode in the main text and in the appendix.Guidelines:• The answer NA means that the paper does not include experiments.• The experimental setting should be presented in the core of the paper to a level of detail that is necessary to appreciate the results and make sense of them. • The full details can be provided either with the code, in appendix, or as supplemental material. 
conclusion: 
 Contributions In this study, we introduced PhyloGen, a novel approach leveraging pre-trained genomic language models to enhance phylogenetic tree inference through graph structure generation. By addressing the limitations of traditional MCMC and existing VI methods, PhyloGen jointly optimizes tree topology and branch lengths without relying on evolutionary models or equal-length sequence constraints. PhyloGen views phylogenetic tree inference as a conditionally constrained tree structure generation problem, jointly optimizing tree topology and branch lengths through three core modules: (i) Feature Extraction, (ii) PhyloTree Construction, and (iii) PhyloTree Structure Modeling. These modules map species sequences into a continuous geometric space, refine the tree topology and branch lengths, and maintain topological invariance. Our method demonstrated superior performance and robustness across multiple real-world datasets, providing deeper insights into phylogenetic relationships.Limitations and Future Works While our model demonstrates outstanding performance on standard benchmarks, it may benefit from using more expressive distributions or incorporating prior constraints to better capture complex dependencies and interactions in the latent space. Additionally, although the Neighbor-Joining (NJ) algorithm is effective for iterative tree construction, it is computationally intensive. We are exploring efficient data structures and parallel processing techniques to address this bottleneck. Furthermore, our model has primarily been applied to genomic data, and further research is needed to extend its applicability to diverse biological data, such as protein and single-cell data. 
experiment statistical significance: 
 Question: Does the paper report error bars suitably and correctly defined or other appropriate information about the statistical significance of the experiments?Answer: [Yes]Justification: In the experimental section, we provide standard deviations for the ELBO and MLL metrics and, in the appendix, we mention the choice of parameters and describe how we ran the various methods.Guidelines:• The answer NA means that the paper does not include experiments.• The authors should answer "Yes" if the results are accompanied by error bars, confidence intervals, or statistical significance tests, at least for the experiments that support the main claims of the paper. • The factors of variability that the error bars are capturing should be clearly stated (for example, train/test split, initialization, random drawing of some parameter, or overall run with given experimental conditions). • The method for calculating the error bars should be explained (closed form formula, call to a library function, bootstrap, etc.) • The assumptions made should be given (e.g., Normally distributed errors). • It should be clear whether the error bar is the standard deviation or the standard error of the mean. Guidelines:• The answer NA means that the authors have not reviewed the NeurIPS Code of Ethics.• If the authors answer No, they should explain the special circumstances that require a deviation from the Code of Ethics. • The authors should make sure to preserve anonymity (e.g., if there is a special consideration due to laws or regulations in their jurisdiction). 10. Broader Impacts Question: Does the paper discuss both potential positive societal impacts and negative societal impacts of the work performed? Answer: [Yes] Justification: We illustrate in the Appendix, in the Broader Impacts section, about the societal impact of our work, which has the potential to greatly improve our understanding of the evolution, origins, and transmission mechanisms of viruses and bacteria, particularly in the areas of public health and disease control. Guidelines:• The answer NA means that there is no societal impact of the work performed.• If the authors answer NA or No, they should explain why their work has no societal impact or why the paper does not address societal impact. • Examples of negative societal impacts include potential malicious or unintended uses (e.g., disinformation, generating fake profiles, surveillance), fairness considerations (e.g., deployment of technologies that could make decisions that unfairly impact specific groups), privacy considerations, and security considerations.• The conference expects that many papers will be foundational research and not tied to particular applications, let alone deployments. However, if there is a direct path to any negative applications, the authors should point it out. For example, it is legitimate to point out that an improvement in the quality of generative models could be used to generate deepfakes for disinformation. On the other hand, it is not needed to point out that a generic algorithm for optimizing neural networks could enable people to train models that generate Deepfakes faster. • The authors should consider possible harms that could arise when the technology is being used as intended and functioning correctly, harms that could arise when the technology is being used as intended but gives incorrect results, and harms following from (intentional or unintentional) misuse of the technology. • If there are negative societal impacts, the authors could also discuss possible mitigation strategies (e.g., gated release of models, providing defenses in addition to attacks, mechanisms for monitoring misuse, mechanisms to monitor how a system learns from feedback over time, improving the efficiency and accessibility of ML). 
safeguards: 
 Question: Does the paper describe safeguards that have been put in place for responsible release of data or models that have a high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)?Answer: [NA] Justification: Our work has no such risks.Guidelines:• The answer NA means that the paper poses no such risks.• Released models that have a high risk for misuse or dual-use should be released with necessary safeguards to allow for controlled use of the model, for example by requiring that users adhere to usage guidelines or restrictions to access the model or implementing safety filters. • Datasets that have been scraped from the Internet could pose safety risks. The authors should describe how they avoided releasing unsafe images. • We recognize that providing effective safeguards is challenging, and many papers do not require this, but we encourage authors to take this into account and make a best faith effort.12. Licenses for existing assets Question: Are the creators or original owners of assets (e.g., code, data, models), used in the paper, properly credited and are the license and terms of use explicitly mentioned and properly respected?Answer: [Yes] Justification: We cite all research works.Guidelines:• The answer NA means that the paper does not use existing assets.• The authors should cite the original paper that produced the code package or dataset.• The authors should state which version of the asset is used and, if possible, include a URL. • The name of the license (e.g., CC-BY 4.0) should be included for each asset.• For scraped data from a particular source (e.g., website), the copyright and terms of service of that source should be provided. • If assets are released, the license, copyright information, and terms of use in the package should be provided. For popular datasets, paperswithcode.com/datasets has curated licenses for some datasets. Their licensing guide can help determine the license of a dataset. • For existing datasets that are re-packaged, both the original license and the license of the derived asset (if it has changed) should be provided. 
 
 Appendix: A.1 Graph Structure Generation (GSG)Let G = (A, X) denote a graph, where A ∈ R N ×N is the adjacency matrix and X ∈ R N ×F is the node feature matrix with x i ∈ R F being the embedding of node v i . Given a feature matrix X, the target of GSG is to directly learn a graph structure A * by jointly optimizing the graph structure and Graph Neural Networks (GNN) [19, 1, 3] .Feature Extraction To better model the similarity of node pairs, a feature extractor is usually needed to map node features from a high-dimensional input space to a low-dimensional embedding space [40] .Graph Construction A similarity metric function is generally used to compute the similarity between embedding pairs as edge weights. There are several ways to construct a sparse adjacency matrix from a fully connected similarity matrix. For instance, we can create a graph that selects only connected node pairs whose similarity exceeds some predetermined threshold. In addition, we can connect the k nearest nodes to one node, thus constructing a k-nearest neighbor (kNN) [28] .Graph Structure modeling The core of GSG is the structure learner, which models edge connectivity to refine the preliminary graph [12] . Metric-based and neural network-based approaches are generally employed to learn edge connectivity through parameterized networks to receive node representations and regenerate the adjacency matrix that optimizes the graph structure, which is able to reveal the real connectivity relationship between nodes better and can be widely used in various downstream tasks [26] .Variational Autoencoders (VAE) [15] is a deep generative models that learn the distribution of input data by encoding it into a latent space. In this process, the encoder maps each input x to a latent space defined by parameters: mean µ and variance σ. Latent variables z are then sampled from this distribution for data generation.VI is employed within VAE to handle the computational challenges of estimating marginal likelihoods of observed data. This involves computing the log of the marginal likelihood:EQUATIONwhere p θ (X, Z) represents the joint distribution of the observable data x e.g. a Gaussian distribution, N (x|µ, σ) and its latent encoding Z under the model parameter θ.Since the direct estimation of marginal likelihoods is typically infeasible, VI introduces a variational distribution q ϕ (z|x) to approximate the true posterior. The goal of VI is to maximize the Evidence Lower Bound (ELBO), formulated as:EQUATIONThe first term is the reconstruction log-likelihood, log p θ (x|z) can be considered as a decoder, i.e., the log-likelihood between the reconstructed data and the original data given the potential representation.The second term, the KL divergence, quantifies the difference between the variational posterior q ϕ (z|x) and the latent prior p(z). Usually, VAE utilizes a reparameterization trick for gradient backpropagation through non-differentiable sampling operations. Once trained, VAEs can generate new data by directly sampling from the latent space and processing it through the decoder.MCMC-based methods, such as MrBayes [29] and RevBayes [9] , have been widely used for phylogenetic inference due to their ability to explore vast tree spaces. However, these methods are limited by the high-dimensional search space and the combinatorial explosion of tree topologies, which makes estimating posterior probabilities using simple sample relative frequencies (SRF) problematic [8] . The accuracy of these methods is compromised in unsampled tree spaces, and they tend to be unstable when estimating low-probability trees, often requiring impractically large sample sizes to achieve reliable results.VI-based methods offer a more efficient alternative to MCMC by leveraging approximate inference techniques. These methods can be categorized into two main approaches: structure representation and structure generation.Tree Representation Learning. This approach focuses on extracting information from existing tree structures. Subsplit Bayesian Networks (SBN) [45] : Given a collection of trees, SBNs capture the relationships between existing subsplits, providing probabilistic representations of various tree shapes under given data. However, SBNs do not directly address branch lengths.Variational Bayesian Phylogenetic Inference (VBPI) [46] and its variants VBPI-NF [43] and VBPI-GNN [44] : These methods introduce a two-pass approach to learn node representations, including branch lengths. VBPI employs variational approximations to handle branch lengths, allowing for joint modeling of the tree's latent space.Tree Structure Generation. This approach aims to infer tree structures directly from sequence data. PhyloGFN [48] utilizes GFlowNet [10, 20, 21] , a combination of VI and reinforcement learning. PhyloGFN constructs tree topologies by simplifying branch lengths into discrete intervals, requiring posterior data for accurate inference. VaiPhy [16] incorporates the SLANTIS sampling strategy [2] and basic biological models (e.g., Jukes-Cantor model) to estimate topology and branch lengths. ARTree [36] employs a graph autoregressive model to build detailed tree topologies. Branch length estimation is conducted independently through log-probabilities. Once the topology is determined, classical evolutionary models (e.g., Jukes-Cantor [24] or GTR (eneralized time reversible) [34] ) are used to estimate branch lengths via maximum likelihood or Bayesian methods. GeoPhy [22] models tree topologies within a continuous geometric space, offering a different approach to the distribution of tree topologies. ARTree and GeoPhy differ primarily in their approach to modeling tree topologies. While ARTree uses a discrete autoregressive model, GeoPhy leverages continuous geometric representations, providing a more flexible framework for capturing the evolutionary relationships among species.Our model PhyloGen performs phylogenetic inference on biological sequence datasets of 27 to 64 species compiled in [17] . Notably, our approach does not require sequences to be of equal length, overcoming a common limitation in traditional phylogenetic analysis. In Tab. 6, we summarize the statistics of benchmark datasets. To address the convergence challenges often associated with the ELBO in VAE models, we incorporate a scoring function S, implemented via an MLP network. This function assesses each leaf node in the latent space z and provides additional gradient information, facilitating more efficient learning and To further evaluate the performance of each S configuration, we computed their scaled cosine similarity to ELBO. Scaled cosine similarity is calculated as: Scaled Cosine Similarity = x•y ∥x∥∥y∥ , where x and y are vectors of normalized performance metrics for S and ELBO, respectively. As depicted in Figure 8 , the horizontal axis measures the similarity ranging from 0 to 1, with values closer to 1 indicating a closer alignment with ELBO. The vertical axis displays the normalized difference percentage, highlighting how each configuration deviates from ELBO. All configurations maintain high similarity values, generally above 0.8, suggesting robust performance across different architectures. Given these results, the Fully Connected (w/ FC) configuration is selected for its closest to the ELBO curve and has better performance consistency. What's more, the number of layers in MLP has less impact on performance, as verified in Fig. 8 and Fig. 3 .It is important to emphasise that we no longer distinguish between z * and z in the derivation for simplicity and uniformly use the latent variable z to denote the topology.In variational inference, we typically start from the joint probability distribution:EQUATIONwhere θ could be any set of parameters or latent variables, and Y is the observed data.For a tree model, assume that θ includes the tree topology τ (z) and the branch lengths B τ . Thus, we have:EQUATION) where p(Y |τ (z), B τ ) represents the conditional probability of the observed data Y , given the tree structure and branch lengths, while p(τ (z), B τ ) is the joint prior probability of the tree structure and branch lengths.Further assuming that the tree topology τ (z) and the branch lengths B τ are conditionally independent, the prior can be decomposed as:EQUATIONCombining the above formula 16 and 17, we can rewrite the joint probability as:EQUATIONIn variational inference, we introduce a variational distribution q(τ (z), B τ ) to approximate the true posterior distribution p(τ (z), B τ |Y ), and often assume:EQUATIONThus, the ELBO can then be written as:EQUATIONInserting the joint distribution and variational distribution formulas, we obtain:EQUATIONThis can be simplified by integrating the logarithmic terms:EQUATIONThis step introduces the key P and Q ratio. We are essentially comparing the variational distribution q and the conditional probability model p for discrepancies.The conditional distribution C(z|τ (z)) is introduced to enhance the model's capacity for capturing complex dependencies within the data. This distribution is aimed at refining the approximation of latent variables given specific conditions. The revised ELBO, incorporating C, becomes:EQUATIONThis formulation explicitly considers the influence of C(z|τ (z)), reinforcing the variational framework's ability to more accurately model the posterior distributions conditioned on complex hierarchical data structures.Since τ (z) and B τ are sampled from their respective distributions and B τ explicitly depends on the former:EQUATIONBegin the derivation of the model parameters θ:EQUATIONEQUATIONwhereH[Q θ (z)] = -E Q θ (z) [log Q θ (z)] is the differential entropy. The derivative of H is ∇ θ H[Q θ (z)] = -∇ θ E Q θ (z) [log Q θ (z)].Then we use the chain rule to apply the derivative operation to the desired logarithmic term yields: z) , we rewrite the ∇ θ L as:∇ θ H[Q θ (z)] = -E Q θ (z) [∇ θ log Q θ (z)] and ∇ θ log Q θ (z) = ∇ θ Q θ (z) Q θ (EQUATIONEQUATIONBegin the derivation of the model parameters ϕ:EQUATIONWe apply the reparameterization technique to express variables involving stochasticity as a combination of a deterministic transformation and a noise distribution that does not depend on the model parameters. With this transformation, the gradients of the model parameters can be passed directly through subsequent computations without being interrupted by the random sampling process. Next, we assume that z can be expressed as a result of a deterministic transformation h θ (ε z ).Accordingly, B τ is no longer sampled directly from a simple normal distribution but is instead obtained by transforming it through the function h ϕ (•), which combines the independent noise variable ε B (usually from a simple distribution such as the standard normal distribution) and τ (z) to generate B τ .EQUATIONSuch treatment allows the model to take into account the fact that complex dependency structures into account, allowing the computation of gradients to be carried out in this reparameterized way, which is important for reflecting the actual complexity in biological data.The expectation can be rewritten as:EQUATIONSubstituting this into the derivation formula gives:EQUATIONC ψ (Z|τ ) Begin the derivation of the model parameters ψ:EQUATIONPhyloGen's algorithm flow is as shown in Algorithm 1.Algorithm 1 Phylogenetic Tree Generation Initialize model parameters θ, ϕ, ψ 9:while not converged do 10:Sample z * ∼ R ψ (z * |z) using the reparameterization trick.Sample z ∼ Q θ (z|Y ) using the reparameterization trick.12:Compute the loss L total (Eq. 9)13:Compute gradients using backpropagation and Update parameters using Adam optimizer:14:θ ← θ -η∇ θ L (Eq. 10)15:ϕ ← ϕ -η∇ ϕ L (Eq. 11)16:ψ ← ψ -η∇ ψ L (Eq. 12)17:end while 18:return Updated θ, ϕ, ψ 19: end procedure 20: return the generated phylogenetic tree (τ, B τ )We focus on the most challenging aspect of the phylogenetic tree inference task: the joint learning of tree topologies and branch lengths. For this, we employ a uniform prior for the tree topology and an independent and identically distributed (i.i.d.) exponential prior (Exp( 10)) for the branch lengths. We evaluate all methods across eight real datasets (DS1-8) frequently used to benchmark phylogenetic tree inference methods. These datasets include sequences from 27 to 64 eukaryote species, each comprising 378 to 2520 sites. Notably, our approach does not require sequences to be of equal length, thus overcoming a common limitation in traditional phylogenetic analysis. For our Monte Carlo simulations, we select K = 2 samples and apply an annealed unnormalized posterior during each i-th iteration, where λ n = min{1.0, 0.001 + i/H} acts as the inverse temperature. This parameter starts at 0.001 and gradually increases to 1 over H iterations, effectively simulating a cooling schedule commonly used in annealing algorithms, similar to the approach in [45] , with an initial temperature of 0.001, which gradually decreases over 100,000 steps.During the model training process, we utilize stochastic gradient descent to process a total of one million Monte Carlo samples, employing K samples at each training step. The stepping-stone (SS) algorithm [37] in MrBayes is viewed as the gold-standard value. All models were implemented in Pytorch [27] with the Adam optimizer [14] . The MLL estimate is derived by sampling the importance of 1000 samples, with the larger mean value being better. The learning rate is initially set to 1e-4 and is reduced by a factor of 0.75 every 200,000 training steps. Momentum is set at 0.9 to prevent the optimization process from becoming trapped in local minima. Utilizing the StepLR scheduler, the current learning rate is multiplied by 0.75 every 200,000 steps to ensure steady progression, detailed in Tab. 7 and Tab. 8. We mainly compare PhyloGen with three types of methods, including two MCMC-based methods (i.e., MrBayes, SBN), two tree-structure representation learning methods (i.e., VBPI, VBPI-GNN), and five tree-structure generation methods. It should be noted that the results of all baseline methods are not included in the MLL tables, as some of the baseline methods are not provided with source code, and the results of the MLL metrics are not shown in the original paper.A bipartition frequency comparison plot is used in evolutionary analysis to show the bipartition frequencies observed in a sample of tree topologies. Each bipartition represents a node in the tree and defines the division of two sets of taxonomic units (e.g., species) on either side of that node. The frequency of these bipartitions in the posterior tree topology distribution reflects how often each topology appeared during the MCMC sampling process, which in turn reveals the confidence level of the different topological features.We sample 1000 Monte Carlo trees from the posterior distribution Q(τ ), calculate the frequency of each bipartition within these samples, and then compare these frequencies with the bipartition frequencies obtained using MrBayes. These frequencies are then compared with those obtained using the MrBayes method to visualize the consistency and differences between the two methods in terms of tree topology inference through graphs. Fig. 5 shows the bipartition frequency distributions of trees inferred by PhyloGen for the DS1, DS2, and DS3 datasets. The horizontal axis indicates the ranking of the bipartitions in the tree topology, and the vertical axis indicates the normalized frequency of occurrence of the corresponding bipartitions, reflecting the prevalence of various topological features observed during the MCMC sampling process. The results in the figure show that our method's variation curves are consistent with those of MrBayes, which indicates that our method has a high degree of consistency in tree topology inference and is able to capture underlying evolutionary patterns with a confidence level similar to the established gold standard MrBayes method.Tab. 5 shows the performance of our model compared with the removal of the KL loss and the Scoring Function S, respectively. ours performs best on both the ELBO and MLL metrics, with a significant decrease in performance after the removal of the KL loss, suggesting that the KL loss plays a key role in regularizing the model and avoiding overfitting. While removing the S module had a large impact on MLL, the ELBO impact was relatively small, indicating that the impact of the S module is more complex and may be related to specific feature extraction functions. Our future work will explore further enhancements to the S module and investigate other regularization techniques to refine the model's performance.Fig. 6 presents a comparative analysis of the MLL metric across six configurations of our model on the DS1 dataset. The standard configuration ('Ours') demonstrates superior performance with the highest MLL value, indicating optimal fit and stability throughout the training process. Variants without Layer Normalization or with a reduced hidden dimension show significant decreases in MLL metric. Furthermore, models using other distance functions for tree construction-specifically, Euclidean and Cosine distances-along with those using one-hot encoding for leaf node representation demonstrate lower MLL values. Notably, the one-hot approach works the worst, suggesting that both the method of computing the distance matrix and the choice of feature representation for leaf nodes critically influence the model's accuracy. Additionally, the distance matrix (calculated from the latent space z) that we specifically designed as an input to the NJ algorithm effectively captures key information about evolutionary relationships that a simple Euclidean or Cosine distance matrix would not be able to capture, resulting in MLL values lower than our standard configuration.E.5 Hyper-Parameter Analysis (RQ7) The first part of Tab. 9 demonstrates the effect of different output dimensions on the model performance at a fixed hidden layer dimension (hd=256). It can be seen that the ELBO and MLL metrics show different trends when the output dimension is increased from 2 to 16, especially when the output dimension is 8. The ELBO and MLL values are optimal, suggesting that higher output dimensions may contribute to the model performance, but there are some fluctuations. The second and third parts of Tab. 9 compare the effects of different hiding dimensions on model performance with and without using layer regularization. When using layer regularization, the model reaches optimal values for both ELBO and MLL at a hidden layer dimension of 256. The overall decrease in performance when layer regularization is not used shows that layer normalization plays an important role in model stability and performance. We recognise the importance of addressing the societal impact of our work. Phylogenetic inference methods such as the one we propose have the potential to greatly improve our understanding of the evolution, origin, and transmission mechanisms of viruses and bacteria. Such an understanding could have far-reaching societal implications, especially in the fields of public health and disease control. Justification: Our abstract and introduction do contain all the claims, innovations and contributions mentioned in the paper.• The proofs can either appear in the main paper or the supplemental material, but if they appear in the supplemental material, the authors are encouraged to provide a short proof sketch to provide intuition. • Inversely, any informal proof provided in the core of the paper should be complemented by formal proofs provided in appendix or supplemental material. • Theorems and Lemmas that the proof relies upon should be properly referenced. 


Decision: The reviewers find the novelty of the method and its performance as notable strengths. (Minor) weaknesses include: the code being not yet mature enough to be released (it was only made available to the Area Chairs), missing comparisons to recent methods and otherwise limited empirical experiments, and some need for improvement in the presentation. 

Overall, the recommendation is unambiguously in favor of accepting.

(As a final detailed comment: The terminology around "(pre-) aligned" vs "unaligned" sequences should be explicitly defined. To someone familiar with phylogenetic terminology, aligned sequences suggest multiple sequence alignment, which, correctly done, should in principle be more informative about phylogenetic relationships than unaligned sequences as one reviewer points out.)

Reviews: review: paper_summary: paper_summary: This paper presents PhyloGen, a novel approach for phylogenetic tree inference using pre-trained genomic language models and graph structure generation. PhyloGen aims to jointly optimize tree topology and branch lengths without relying on evolutionary models or equal-length sequence constraints. The method demonstrates superior performance and robustness across multiple real-world datasets compared to existing approaches.
strength_weakness: strength_weakness: Limited discussion of potential limitations or failure cases

Lack of comparison to some recent methods in the field

Source code not provided
questions: questions: How does PhyloGen perform on larger datasets with hundreds or thousands of species?

Have you explored using other types of pre-trained language models besides DNABERT2?

How sensitive is the method to hyperparameter choices?

What are the main computational bottlenecks of the approach?

Have you tested PhyloGen on simulated datasets where the true phylogeny is known?
limitations: limitations: While some robustness tests were performed (node additions and deletions), the paper doesn't explore all possible types of data noise or perturbations.

The paper doesn't deeply explore the interpretability of the model's decisions or the learned representations.

The paper doesn't extensively discuss how well the method generalizes to different types of genetic data or organisms not represented in the test datasets.

The method relies on pre-trained genomic language models (specifically DNABERT2), but doesn't explore how performance might vary with different pre-trained models.
paper_summary: paper_summary: The authors propose a new method, PhyloGEN, for phylogenetic inference. The method is able to perform end-to-end variational inference in order to jointly optimize the tree topology and the branch lengths. To achieve this, the authors propose using a pre-trained genomic language model to extract genome embeddings, and the embeddings are then used to form an initial topology which is refined iteratively using the proposed topology and branch length learning modules.
strength_weakness: strength_weakness: - The methods section (Section 3) is difficult to read in that variables are introduced without leading sentences, some notation is undefined, and design choices are introduced without clear motivation.
  - What is the motivation for $f_i$ and how is it defined in (2) if there are multiple children $f_j$?
  - The function $h$ is used before its definition.
  - Which function $h$ is the latter $h$ in the definition $h(x_i, x_j) = h(x_i, x_i - x_j, d_{ij}^2)$?
  - Why is max aggregation a good choice?
  - The scoring function $S$ is not described in enough detail. What does it output and how does it provide additional gradient information?

- Most experiments, in particular RQ2, RQ4, RQ5, RQ6, and RQ7 are performed on only one data set (DS1), making it difficult to draw meaningful conclusions. I would also consider the runtime an important metric, but it is measured only in conjunction with the robustness assessment on only one data set.

- The authors do not provide code for their method or experiments (although the code is promised to be released later). Without the code, considering the lack of detail in the methods section, it would be difficult to reproduce the proposed method.

- Minor:
  - The red and yellow colors are opposite to what they should be in the caption of Table 1.
  - There is a subsection numbered 3.1 in Section 4.
  - RQ7 is only referred to in the Appendix, and not in the main text.
  - Typos: "TreeEnocoder" (Figure 2 caption), "Metrtic" (Figure 3 title)
questions: questions: - Could the authors explain that if MrBayes is considered the gold standard, and the proposed method is validated by comparing it against MrBayes in RQ2 and RQ3, why is the proposed method able to achieve such a big increase in MLL and ELBO while the other methods yield relatively similar results? Relatedly, do the resulting topologies look significantly different in comparison to the other methods?

- The proposed method has low standard deviations in Table 2 for multiple datasets. Is this due to the optimization being robust with respect to the random seed or could it not indicate that $Q(z)$ does not support multiple topologies? Do the MrBayes posteriors support multiple topologies?

- How does the Bipartition Frequency Distribution look for other methods like GeoPhy and ARTree?

- How does the runtime of MrBayes compare to the runtime of the proposed method?

- Can the authors hypothesize why a simple fully connected $S$ would work so well (Figure 3)?

- In the construction of the distance matrix, could the authors explain why using XOR makes sense since aren't the $z_i^*$ continuous? Moreover, would it make sense to try out other similar algorithms to neighbor-joining to initialize the topology?

- I would suggest that the authors move section 3.3 (with the derived gradients) to the Appendix, leaving space for a more detailed description of the various design choices.
limitations: limitations: The authors do not explicitly discuss the limitations of their work. As I am not an expert in phylogenetic inference, it is difficult to judge the assumptions made in the paper. Are there conditions that are required to be met for using a pre-trained language model for extracting embeddings, and is the assumption of the tree topology and branch lengths being conditionally independent reasonable? It would be helpful for the authors to elaborate on what they think the possible limitations of their own work are.
paper_summary: paper_summary: The paper propose phylogenetic tree inference by modeling it as a problem of conditional-constrained tree structure generation. Its goal is to jointly generate and optimize the tree topology and branch lengths. By mapping species sequences into a continuous geometric space, PhyloGen performs end-to-end variational inference without limiting the possible topologies. To maintain the topology-invariance of phylogenetic trees, distance constraints are applied in the latent space to preserve translational and rotational invariance. They have shown the proposed model  is  effective  across eight real-world benchmark datasets.
strength_weakness: strength_weakness: I am not sure if the formulation equation 7, adding the term R mathematically make sense, as orignal ELBO in equation 6, is the lower bound, now adding the term log of R which is negative, and trying to maximize equation 7 would result in minimizing R?  Also in equation 7, q(\tao (z)) seems missing, or is there a kl between R and Q from module B?
I think the paper presentation needs to be improved. The figure 2 is very helpful but the text  and some notations are a bit confusing.
questions: questions: 1. Regarding the statement, "Existing Variational Inference methods, which require pre-generated topologies and typi6 cally treat tree structures and branch lengths independently, may overlook critical sequence features, limiting their accuracy and flexibility", can you elaborate on what kind of features do you refer to that the method could overlook?

2. I am not sure fromt he distance matix to the tree topology, the  Neighbor-Joining (NJ) algorithm can be differentiated to be able to learm end to end?

3. I think the section C.1. Topology Learning is written somehow confusing. Especially from line 121 to line 124. Why the distribution  p(\tao(z)) is represented as the distribution of y conditioned on tree structure? same as the q(\tao(z), B_\tao) is not clear what distribution it represents.  

4. I am not sure in this kind of tasks, if the ELBO is a good metric to reflect the inference power of the model.

5. The VAE model basically takes the embedding of the genomics sequences, and the initial topology of the tree \tao, then try to reconstruct the tree (without the leaf length)? I was thinking what happens if we remove the learning of the first component, make the phylo tree construction constant using directly the embeddings from LLM?

6. I am not sure if I understood how the scoring function is implemented, f(z), z is the embedding of the genomic sequence, f applied on the leaf nodes? or all the nodes, and why minimize it?
limitations: limitations: 1. While PhyloGen outperforms many existing methods, it still requires substantial computational resources, particularly during the training phase. The complexity of jointly optimizing tree topology and branch lengths, along with maintaining topological invariance, sounds computationally heavy, for instance:
1.1  Monte Carlo simulations to sample from the posterior distributions of tree topologies and branch lengths.
1.2 The tree construction module computes distance matrices from latent variables and uses algorithms like Neighbor-Joining (NJ) to generate initial tree structures. 
1.3 The utilization of  a pre-trained genome language model (DNABERT2) to transform DNA sequences into genomic embeddings, 

2, Looking at the table 5, I am wondering how it works in inference time, would peformance drops when we have less or more specious than the one in the training data the model was trained on? Also how it translates to other specious? moreover, did you train one model per dataset, or is there a way to train all together?
paper_summary: paper_summary: The authors introduce PhyloGen, a new method that uses a pre-trained genomic language model to generate phylogenetic trees without relying on evolutionary models or aligned sequences. PhyloGen treats phylogenetic inference as a conditionally constrained tree structure generation problem, jointly optimizing tree topology and branch lengths through three modules: Feature Extraction, PhyloTree Construction, and PhyloTree Structure Modeling. A Scoring Function guides the model towards stable gradient descent. The method achieves state-of-the-art performance when evaluated on benchmark datasets.
strength_weakness: strength_weakness: 1. Ablation Experiments Clarity

The ablation experiments are unclear. My main concern is the lack of clarity regarding why the proposed method significantly outperforms other methods. Specifically:

(1) Table 5: This table shows the average metrics across the eight datasets. However, Table 1 and Table 2 only shows the metrics for each dataset individually without any average metrics. It’s unclear what the performance gap is between the ablation model and the baseline methods when KL or S are removed.

(2) Figure 6: The ablation studies in this figure only show the metrics for the DS1 dataset. Please include the average metrics across all eight datasets.

To improve the clarity of the ablation experiments, I strongly suggest making the evaluation settings consistent in Table 1, Table 2,Table 5, and Figure 6. At the very least, include the average metrics across the eight datasets. Additionally, provide more explanations based on the complete ablation experiments to clarify why the proposed method significantly outperforms other methods.

2.  Claim on Aligned Sequences vs. Unaligned Sequences

The claim that unaligned sequences better reflect actual sequence variation than aligned sequences is questionable. Generally, aligned sequences with direct site-to-site comparisons are more effective in comparative sequence analysis. The authors need to provide more direct experimental evidence to support the advantage of using unaligned sequences.

3. Writing Issues

There are many writing issues. For instance,

a. line 141. "m_ij" -> "h"

b. In the caption of table1, "VBPI-GNNuse" -> "VBPI-GNN use"

c. line 213, "3.1 Experiment Setup"  ->  "4.1 Experiment Setup"
questions: questions: 1. Could you clarify whether the computation of the marginal log likelihood involves an evolutionary model like [1]? Does the evolutionary model utilize aligned sequences? If so, does this contradict the statement in the article that sequence alignment is not required?

2. I noticed that the author provides a comparison of run times. Could you provide details about the GPU used in the experiments, as well as the memory utilization?


[1] M. Zhou, Z. Yan, E. Layne, N. Malkin, D. Zhang, M. Jain, M. Blanchette, and Y. Bengio. Phy-logfn: Phylogenetic inference with generative flow networks. arXiv preprint arXiv:2310.08774, 2023.
limitations: limitations: I have no concerns on the potential societal impact of their work.

Comments: title: Reviewers comment
Thanks for Raising the Score
General response to all the reviewers
Have our response addressed your concerns?
Summarized Response to the Unaddressed Concerns
comment: Thank you for your thoughtful feedback and for considering our responses. 

We’re pleased you find the idea interesting and want to assure you that we’re committed to improving the paper’s presentation, methodology, and evaluation metrics, as you suggested.


Your **positive opinion is crucial**, and we hope you **might reconsider your score** given our commitment to these improvements.


Thank you again for your valuable insights.
Thank you for authors answer to clarify my questions. After reading the answer, I decide to keep my score, I think the idea is interesting but the paper presentations needs to be improved and parts of the methodology as well as the evaluation metrics needs to be clarified.
Dear Reviewer MKto,

We would like to express our sincere gratitude for your thorough review and for increasing your score after considering our response, which is very important to us. Based on your suggestions, we will keep on polishing our manuscript and add relevant tables to the camera-ready version. Thank you again for your valuable review and efforts in improving our work.

Warm regards,

Authors
Thank you for your response. I raised my score. Please keep in mind to fix the typos and grammar errors in the final version of the manuscript.
Dear Reviewers,


We greatly appreciate your effort and valuable feedback. In our response, we carefully illustrated and answered all your questions in detail. Additional experiments and analysis results are also provided as you had requested. In addition, we further clarified several unclear statements in the paper. We have incorporated all changes into the revised manuscript for your consideration. We hope your concerns have been addressed.


As you may know, unlike previous years, the discussion period this year can only last until August 13, and we are gradually approaching this deadline. We take it seriously and would like to discuss this with you during this time. We would be happy to provide more information based on your feedback or further questions.

If you are satisfied with our response, please consider updating your score. If you need any clarification, please feel free to contact us.

Best regards,
Authors
Dear reviewer XKPL,


We have provided the complete results of the ablation study for the requested dataset (DS1-DS8) based on your suggestions.


We sincerely appreciate all your valuable feedback to help us improve our work. **We believe that the issues raised by the reviewers have been addressed in the revision.** We hope that you will reconsider our submission of this response. If you are satisfied with our response and efforts, please consider updating your rating. If you need any clarification, please feel free to contact us. We are very pleased and look forward to hearing from you!


Best regards,

Authors
Dear Reviewer QAeZ,


We highly value your feedback and suggestions, which are crucial for us to review further and enhance our work. **We have revised our manuscript and addressed all your concerns.** We sincerely invite you to review our responses, hoping that it will help dispel any misconceptions about our work. Below is a summary of our response:


- Performance on Larger Datasets: We clarified PhyloGen's scalability and ensured a fair comparison by testing it on the same datasets as other state-of-the-art methods.
- Pre-trained Models and Hyperparameter Sensitivity: We justify the selection of DNABERT2 and demonstrate PhyloGen's robustness across different hyperparameter settings.
- Interpretability and Generalizability: We highlight several case studies and visualizations that demonstrate PhyloGen's interpretability and generalizability in various biological contexts. We also expanded the robustness tests and addressed concerns regarding computational bottlenecks.


In this revised version, we have sincerely and diligently responded to your valuable comments. **We trust that our responses have clearly addressed your questions and invite you to review them.**


We are pleased to inform you that several reviewers have given **high scores**, with one even **raising their score** after reviewing our responses. Notably, the reviewer with **high confidence** recognized the strength and robustness of our approach. In light of this, we hope you reconsider your assessment and increase your score. If you have any further questions, we would be happy to address them.


Best regards,

Authors
**Additional Experiments on Ablation Study (Part 2):**

Following our recent update, we have also completed additional ablation studies for datasets DS5 through DS8 to demonstrate our model's effectiveness further. The results support the conclusions drawn from the initial datasets. These results will be made available in the final manuscript for a comprehensive review.



**DS5 (50):**

| Methods | ELBO (↑) | MLL (↑) |
| --- | --- | --- |
| Ours | -8053.23 | -7928.4 |
| Ours w/o LN | -8133.16 | -7947.05 |
| Ours-Hid=64 | -8232.91 | -7953.48 |
| Ours-Euclidean | -8884.67 | -8674.23 |
| Ours-Cosine | -8550.97 | -8308.62 |
| Ours-One-Hot | -9915.63 | -9728.71 |


**DS6 (50):**

| Methods | ELBO (↑) | MLL (↑) |
| --- | --- | --- |
| Ours | -6324.90 | -6330.21 |
| Ours w/o LN | -6748.26 | -6341.30 |
| Ours-Hid=64 | -7630.79 | -7337.23 |
| Ours-Euclidean | -9039.58 | -8811.08 |
| Ours-Cosine | -7822.54 | -7545.93 |
| Ours-One-Hot | -9656.64 | -9437.25 |


**DS7 (59):**

| Methods | ELBO (↑) | MLL (↑) |
| --- | --- | --- |
| Ours | -36838.42 | -36838.42 |
| Ours w/o LN | -37579.46 | -37085.76 |
| Ours-Hid=64 | -37997.11 | -37571.84 |
| Ours-Euclidean | -43182.38 | -42889.76 |
| Ours-Cosine | -41446.14 | -41165.21 |
| Ours-One-Hot | -45627.14 | -45358.70 |


**DS8 (64):**
| Methods | ELBO (↑) | MLL (↑) |
| --- | --- | --- |
| Ours | -8409.06 | -8171.04 |
| Ours w/o LN | -8842.80 | -8358.32 |
| Ours-Hid=64 | -9319.05 | -8911.28 |
| Ours-Euclidean | -10128.19 | -9777.77 |
| Ours-Cosine | -9812.72 | -9436.38 |
| Ours-One-Hot | -13475.34 | -13187.75 |
Dear Reviewer,

We sincerely appreciate your efforts and valuable feedback. If you are satisfied with our responses and our improvements, please consider updating your score. 

If you need further clarification, please don't hesitate to contact us. We are grateful for your time and look forward to your response!
Dear Reviewer,

We sincerely appreciate your efforts and valuable feedback. If you are satisfied with our responses and our improvements, please consider updating your score. 

If you need further clarification, please don't hesitate to contact us. We are grateful for your time and look forward to your response!
Dear Reviewer,

We sincerely appreciate your efforts and valuable feedback. If you are satisfied with our responses and our improvements, please consider updating your score. 

If you need further clarification, please don't hesitate to contact us. We are grateful for your time and look forward to your response!
**Additional Experiments on Ablation Study:**

Thank you for your feedback and for raising your score. I understand your concerns regarding the experiments in Sec. 4.5 and have conducted additional ablation studies on multiple datasets, summarized below:

The rows represent different configurations: **Ours** (full model); **Ours w/o LN** (without layer normalization); **Ours-Hid=64** (reduced hidden dimensions); **Ours-Euclidean** and **Ours-Cosine** (replacing the custom distance matrix with Euclidean and cosine metrics); and **Ours-One-Hot** (using one-hot encoding for leaf nodes). The numbers in parentheses following each dataset (e.g., DS1 (27)) represent the number of species sequences.


**DS1 (27):**

| Methods | ELBO (↑) | MLL (↑) |
| --- | --- | --- |
| Ours | -7005.98 | -6910.02 |
| Ours w/o LN | -7111.31 | -7187.40 |
| Ours-Hid=64 | -8081.97 | -7728.90 |
| Ours-Euclidean | -9284.57 | -9091.84 |
| Ours-Cosine | -8943.73 | -8757.81 |
| Ours-One-Hot | -10433.23 | -9282.87 |

**DS2 (29):**

| Methods | ELBO (↑) | MLL (↑) |
| --- | --- | --- |
| Ours | -26362.75 | -26257.09 |
| Ours w/o LN | -26443.78 | -26301.85 |
| Ours-Hid=64 | -27511.50 | -27373.58 |
| Ours-Euclidean | -29355.51 | -29240.52 |
| Ours-Cosine | -28265.33 | -28147.96 |
| Ours-One-Hot | -30954.67 | -30848.88 |

**DS3 (36):**

| Methods | ELBO (↑) | MLL (↑) |
| --- | --- | --- |
| Ours | -33430.94 | -33481.57 |
| Ours w/o LN | -33737.22 | -33481.57 |
| Ours-Hid=64 | -34400.24 | -34209.57 |
| Ours-Euclidean | -38067.04 | -37925.02 |
| Ours-Cosine | -35976.66 | -35823.20 |
| Ours-One-Hot | -39120.53 | -38989.87 |

**DS4 (41):**

| Methods | ELBO (↑) | MLL (↑) |
| --- | --- | --- |
| Ours | -13113.03 | -13063.15 |
| Ours w/o LN | -13395.25 | -13111.01 |
| Ours-Hid=64 | -13410.90 | -13115.28 |
| Ours-Euclidean | -15438.36 | -15270.67 |
| Ours-Cosine | -13417.38 | -13101.95 |
| Ours-One-Hot | -16760.54 | -16606.26 |

**Analysis:**

Using Euclidean or cosine distance metrics significantly lowers MLL, indicating their limitations in capturing complex evolutionary patterns. The one-hot encoding yields the lowest MLL, highlighting the superiority of our feature extraction module and pre-trained language model. Additionally, reducing hidden dimensions shows a clear performance drop, underscoring the importance of model capacity. 

We are conducting further experiments, and I will upload the full results as soon as they are available.
Thank you for your valuable feedback and for raising your score. We are conducting the additional Ablation Study experiments and will share the results with you as soon as they are ready. 

We are also refining the manuscript to address the presentation issues raised by you and the other reviewers, aiming to enhance its clarity.
Thank you for answering my questions in detail and for providing a discussion on the limitations of your method. In light of this, I am slightly raising my score.

I cannot comment on the standard practice of using only DS1 for the experiments, but given the various somewhat ad-hoc design choices in your method,  I would expect to see at least the experiments of Section 4.5 performed using more data sets in a final revision. I also hope the authors pay careful attention to the multiple presentation issues raised by me and the other reviewers.