Ultrafast classical phylogenetic method beats large protein language models on variant effect prediction
Authors: Sebastian Prillo
Wilson Y. Wu
Yun S. Song
Venue: NeurIPS 2024 poster (2024)
Keywords: protein language models
evolution
phylogenetics
rate estimation
variant effect prediction
Primary area: machine_learning_for_other_sciences_and_fields

Abstract: Amino acid substitution rate matrices are fundamental to statistical phylogenetics and evolutionary biology. Estimating them typically requires reconstructed trees for massive amounts of aligned proteins, which poses a major computational bottleneck. In this paper, we develop a near-linear time method to estimate these rate matrices from multiple sequence alignments (MSAs) alone, thereby speeding up computation by orders of magnitude. Our method relies on a near-linear time cherry reconstruction algorithm which we call FastCherries and it can be easily applied to MSAs with millions of sequences. On both simulated and real data, we demonstrate the speed and accuracy of our method as applied to the classical model of protein evolution. By leveraging the unprecedented scalability of our method, we develop a new, rich phylogenetic model called SiteRM, which can estimate a general site-specific rate matrix for each column of an MSA. Remarkably, in variant effect prediction for both clinical and deep mutational scanning data in ProteinGym, we show that despite being an independent-sites model, our SiteRM model outperforms large protein language models that learn complex residue-residue interactions between different sites. We attribute our increased performance to conceptual advances in our probabilistic treatment of evolutionary data and our ability to handle extremely large MSAs. We anticipate that our work will have a lasting impact across both statistical phylogenetics and computational variant effect prediction. FastCherries and SiteRM are implemented in the CherryML package https://github.com/songlab-cal/CherryML.

Ultrafast classical phylogenetic method beats large protein language models on variant effect prediction
Sebastian Prillo* University of California, Berkeley spri pprillo@berkele, edu
Wilson Wu
University of California, Berkeley nosliv@berkeley. edu
Yun S. Song
University of California, Berkeley yssüberkeley. edu
Abstract
Amino acid substitution rate matrices are fundamental to statistical phylogenetics and evolutionary biology. Estimating them typically requires reconstructed trees for massive amounts aligned proteins, which poses a major computational bottleneck. In this paper, we develop near-linear time method to estimate these rate matrices from multiple sequence alignments (MSAs) alone, thereby speeding up computation by orders of magnitude. Our method relies on a near-linear time cherry reconstruction algorithm which call FastCherries and it can be easily applied to MSAs with millions of sequences. On both simulated and real data, we demonstrate the speed and accuracy of our method applied to the classical model of protein evolution. By leveraging the unprecedented scalability of our method, we develop new, rich phylogenetic model called SiteRM, which can estimate general site-specific rate matrix for each column of an MSA. Remarkably, variant effect prediction for both clinical and deep mutational scanning data in ProteinGym, we show that despite being an independent- sites model, our SiteRM model outperforms large protein language models that learn complex residue- residue interactions between different sites. We attribute increased performance to conceptual advances in our probabilistic treatment of evolutionary data and our ability handle extremely large MSAs. We anticipate that work will have lasting impact across both statistical phylogenetics and computational variant effect prediction. FastCherries and SiteRM are implemented in the CherryML package https: //tp///gtttt///tt com/ /oonglg cal/ hal/CherryML.
Introduction
Amino acid substitution rate matrices are crucial to statistical phylogenetics and evolutionary biology. Perhaps best known for their use in phylogenetic tree reconstruction [1. 2, 3, 4, 5], these transition-rate matrices allow us infer the evolutionary history of protein family. This has given insights into protein evolution, function, and fitness, which among other important applications, has informed the development of vaccines during the COVID-19 pandemic [6]. Statistically speaking, transition-rate matrices are parameters of continuous- time Markov chain models describing protein evolution. One of the most popular such models to date is the seminal model of Le and Gascuel [7] (LG for short), which posits that each site protein evolves independently following a continuous -ime Markoy chain parameterized by a global transition-rate matrix Q E 202020 up site-specific rates.
"Equal contribution; authors listed alphahetically
38th Conference on Neural Information Processing Systems (NeurIPS 2024).

---

Unfortunately, estimating transition-rate matrices under the LG and related classical models from a set of multiple sequence alignments (MSAs) is challenging task. Indeed, one has to first estimate phylogenetic trees and site-specific rates for each MSA, and then estimate the transition rate matrix (or matrices) Q which best explain the data. The procedures of tree plus site-rate estimation and of transition-rate matrix estimation are usually performed in a coordinate-ascenn fashion until convergence. Both of these steps are computationally challenging: maximum likelihood estimation (MLE) of trees is in general NP-hard [8], while MLE of a rate matrix Q given MSAs, trees and site-rates has historically proven computationally demanding task [9, 10].
Fortunately, recent work has made it possible to perform MLE of transition-rate matrices given MSAs, trees and site-rates in a scalable fashion. The method, called CherryML [11], replaces the full joint likelihood of the MSA data with composite likelihood over cherries in the trees. When applied to the LG model, CherryML is orders of magnitude faster than traditional MLE with the Expectation-Maxian Maximization (EM), all while sacrificing relatively small amount of statistical efficiency, estimated to be around 50% on simulation studies. The method was also shown to perform on-par with EM on real datasets when evaluated on held-out likelihood. While the work of CherryML represents a significant advance, the method still requires expensive phylogenetic tree reconstruction as prerequisite step. This hinders the scalability of the method when applied to many MSAs, some with millions of sequences. In fact, in the original CherryML work, MSAs were subsampled down to a thousand sequences each to alleviate the burden of tree reconstruction.
In this work, we introduce new methodology to significantly speed up the end-to-end estimation of rate matrices under the LG model, and extend this methodology to estimate site-specific rate matrices. Specifically, we speed up the phylogenetic tree plus site-rate estimation step required by the CherryML method, using a method which we call FasrCherries. The resulting end-to-end method consisting of CherryML with FastCherries can, given set of starting MSAs, estimate site-specific rates and an accurate transition-rate matrix in nearly linear time in the input dataset size. The method is thus essentially computationally optimal up to logarithmic factors and other constants. On a real MSA with more than 450,000 sequences of length 364, FastCherries took just 1,000 seconds 1 CPU core, and we can estimate site-specific rate matrices at an additional cost of one second per site.
Thanks to its unprecedented scalability, we applied CherryML with FastCherries to estimate site- specific transition-ratt matrices for MSAs with hundreds of thousands of sequences. This model, which we call SiteRM for short, posits that each site i of protein family f evolves independently according to family- and site- site-speciiic rate matrix Qf- Strikingly, we show that this independent-sing model excels at variant effect prediction, outperforming the seminal EVmutation [12] model with epistatic interactions as well as many large protein language models such as an ESM-1v ensemble [13]. We attribute the increased performance of SiteRM to the principled probabilistic treatment of the evolutionary process, which is absent in competing approaches. This opens up new avenue to variant effect prediction via the use of probabilistic models of protein evolution in time. FastCherries and SiteRM are available at https: : / htt//ggtttttt giggth..bbb... com/ hta///ongabb-oo cal/ val/CherryML.
2 Background and Related Work
2.1 Classical models of protein evolution
As time passes. proteins experience mutations, many of which are neutral, having no effect on fitness. On the other hand, some mutations are deleterious, leading to decreased fitness and hence negative selection, while other mutations confer increased fitness and thus are favored by natural selection. This complex evolutionary process of mutation and selection leads to the vast biological diversity observed today, and modeling it is the goal of statistical phylogenetics. Formally, model of protein evolution is a conditional distribution p(y Jx. t. f) describing the probability of sequence y. the evolved from sequence after time I, in protein family f.
Standard probabilistic models of protein evolution. Early work on protein evolution relied on counting-base based heuristics, such as the seminal works of Dayhoff et al. [14] and JTT [15]. The first probabilistic model of protein evolution was proposed by Whelan and Goldman [16]. It posits that each site of any protein evolves i.i.d. following a reve-reee continuous-timg time Markov chain parameterized by global transition-ran rate matrix Q 200 normalized so that the expected number of substitutions in a unit of time equals 1. Each protein family f may evolve at a different rate ay. Letting ly be the number of columns of an MSA for protein family f. we have pwAG (y) zz t, f)

---

I' {xp(a)1()2z vil, where exp( M) (a, bj denotes the (a,b) entry of the matrix exponential exp(M). The seminal extension of the above model by Le and Gascuel [7] still uses a global, normalized transition-rate matrix Q, but incorporates site-specific rates (a a2)..., af,,) giving
15 PLolul ff = II exp(a (Q) 31].
This is one of the most popular models of protein evolution to date, used extensively in phylogenetic tree reconstruction. The LG model has total of 400 - Σイ-ı4 I parameters, where m denotes the
number of protein families under study. The site-specific rates a are typically constrained with a model of site-rate variation, such as the I' model [17] or the babbablity-distributto----
model.
Phylogenetic models. Unfortunately. estimating models of protein evolution is a challenging task. Indeed, one does not have access to training data of the form (æt, y, f).. Instead, one has access to MSAs. Formally. let D = (Di.... Dm)) be MSAs for m protein families. The data D are modeled with phylogenetic model [11], which is parameterized by (1) a model of protein evolution p(y]ar, (, (' (as described above). (2) phylogenetic trees Ti, T..... T. (one tree per MSA), and (3) root state distribution Tome (r). The phylogenetic model posits that the MSA data were generated by sampling the root state from each tree following Froa (x)) and then running the evolutionary model p(y|a, f) down each tree. We denote the phylogenetic model's likelihood as Pphyto Pphylo(DiT). To estimate a model of protein evolution p(y| t, f) with MLE from data from a phylogenetic model Pohylo (D)T), one thus has to deal with the nuisance parameters T = (Ti, T2,,............. ....... Tm). Typically, this is done with coordinate ascent, wherein one alternately optimizes the trees 'T plus site-rates af given the transition- rate matrix Q and then optimizes the transition-ral rate matrix Q given the tree estimates T plus site-rates. Both of these steps have historically proven to be computationally demanding.
2.2 CherryML
Recently, the work of CherryML [11] introduced a scalable and accurate method for estimating the transition-rate matrix (or matrices) Q given the trees T and site-rates a: It is assumed, as is typical in statistical phylogeneties, that Q is time-reversible. For the LG model, this method was shown to be orders of magnitude faster than the traditional EM method. In fact, CherryML's runtime for the LG model is linear in the input data size up to logarithmic factors, making it essentially computationally optimal. From statistical point of view, CherryML is just ~ 50% less efficient than full MLE.
Composite likelihoood. The key to CherryML lies in the use of composite likelihood, which generally leads to consistent parameter estimates under weak assumptions [18]. Adapt- ing the notation of CherryML [11], for generic phylogenetic model Ppaylo the full joint like- lihood of the data is E = log Pphyto (D(T) = E Pphyti (D)TT). CherryML replaces the MSA's log-likelihood log Pphyio (TTTT)) (D))) by a composite likelihood over cherries in the trees, where cherries are iteratively picked until either O or leaf remains. Specifically, letting Hur u/) (r)hisss, be the ef cherries in tree Ty, CherryML considers the composite likelihood 6oomp = Σラ=1 Ejz log Pphyto (Ds[ug] Dy(u5] Tr), where D3 [u/] denotes the sequence in MSA D, cor- responding to leaf u/. For stationary time-reversible model, the term Pplyto (Dylue) Dyluf Tr) is exactly equal to p(Dr 1 where t.f is the distance between uz and uy in tree Ty. Therefore, the composite likelihood reduces to fcomp ΣØ= E log (D (ug) (1 This way, Cherry ML reduces the problem of learning the transition--rtt matrix (or matrices) Q given the trees T. site-rates, and MSAs D to supervised learning problem over the data (a,, y, t, f) given by (Dy[u/), Dy [or]] t2, f). Cherries are considered in both directions when forming the composite likelihood, so that if (ar, y. f) is training datapoint, then so is (y, x t, f). Thus, CherryML can be viewed as modern, principled version of the JTT method [15] where sequences are paired and then rate matrix estimation proceeds by MLE rather than via counting heuristics.
Time quantization. Cle further simplifies the composite likelihood by quantizing (or dis- cretizing) time into a finite number b of values TI < 72 <... < ть- Letting q(() be the quan- tized value of t (which is chosen to minimize relative error), in the LG model with site-rates rf for site S I ly of family f. CherryML leads to the quantized composite likelihood
L foomp, quant Σm= Σj E's (ttespplor//t/)//////
D 14)) The terms in this expres-
3

---

sion can be grouped together by quantization time The leading to
foomp, quant
exxeexTT
k=1
where Cj a 20 x 20 count matrix for the number of transitions between two amino acids at a quantized distance of Tı, and (.--) denotes the matrix inner product. This function can be evaluated in time (bbs²) where
20 the number of states, which is remarkable since this no longer depends on the input dataset size. It can be optimized using automatic differentiation software such as PyTorch [19] and Tensorflow [20]. CherryML leverages differentiable implementation of the matrix exponential operator based on robust Taylor expansions [21] and utilizes the Adam optimizer [22].
3 Method
3.1
The SiteRM model
In this work, we propose a richer model of protein evolution, which we call the SiteRM model:
1 (y)e,M(y|z, = II exp((t)) vil.
where Qf....,Q! 0f, are site-specific transition-rate matrices. Unlike Q-matrix described above, Of are not normalized and they subsume site-specific rates. This model has 400 x Σ/=1' Ly parameters and strictly generalizes the LG model. This a particular case of the 'partition model' available in IQTree [5], in which each site forms its own group in the partition.
3.2 Near-linear time end-to-end estimation of model parameters
Here, we present an algorithm for speeding up end-to-end estimation of model parameters. Related work speeding up the tree reconstruction step required for rate matrix estimation is described in Appendix 1. Although we focus on the LG model [7] for the sake of analyzing statistical efficiency, our method will also enable fitting site-specific rate matrices under the richly parameterized SiteRM model, as we illustrate on the task of variant effect prediction. For the LG model, our work speeds up tree and site-rate estimation given the rate matrix Q and MSAs D. Together with CherryML's procedure for estimating a transition- rate matrix Q given trees, site-rates, and MSAs, this leads to a full coordinate-ascent procedure with mear-lineat time updates. making it exceptionally scalable.
Method overview. The key idea of our method is that CherryML's composite likelihood depends on the trees T only through the cherries (u/ uz) and their quantized distances 1/., In other words, there is point in estimating the full trees T if only their cherries will be used. Our method thus proceeds by first grouping the sequences within each MSA D, into disjoint pairs using a near-linear time divide-and-conquer approach which tries to minimize total Hamming distance across all pairs. Letting y be the number of sequences in MSA Dy (so that MSA Dj has size nj x 1y), we obtain Cf [ms/2] /2] pairs (uf, )..),... (wL, vfy)). These pairs are the putative cherries of the tree Ty- Next, the aste-specffic rates of af...., af, and the quantized times t!....,,'' ,th, separating each pair of sequences are estimated using coordinate ascent. We call our method FastCherries. A schematic of the FastCherrress/iieRRM method is provided in Figure 1.
Divide-and-conquer To pair up the set of sequences S, we use an almost-linear time distance-based tree topology reconstruction algorithm based on divide-and-conquer The general operation of the algorithm is as follows. Let d be dissimilarity function between protein sequences. We use the normalized Hamming distance d ignoring gaps, but other dissimilarities may be used. Assuming [S] 2 3. we first try to find the diameter of the set S (i.e., the two furthest away sequences) as follows: (1) take a random sequence 2 E S, (2) find the furthest sequence to 20 in S: ž arg m d(2, 20); (3) find the furthest sequences to 21 in S: 22 = arg max.es d(=,21). Our putative diametrically opposite sequences are given by 21 and 22. Note that this is analogous to the well- linear-time time algorithm used to find the diameter of tree. We now use 21 and 22 pivots
to split the set S into two subsets S and Sy based on distance to 2] and 22:
S = {z E S d(zz d(z,22)) and S2 = {= E S d(z, 21) d(2, 22)).

---

FastTree, etc. hours
XRATE, etc:: hours
FastCherries (this work]: minutes
CherryML: minutes
Multiple Sequence Multiple Sequence
Sequence Multiple Sequence
Multiple Sequence
Tree/Cherries and site rate estimation
Multiple Sequence
Multiple Sequence
Multiple Sequence
Rate
Alignment
Alignment
matrix estimation
-30
FastCherries (this work); near-linear runtime: For each MSA:
Repeat until convergence (-10 times)
Branch length (re-Jestimation
Multiple Sequence
Divide
(x2, yb)
Alignment
Site rate Initialization
On log b)
Conquer
Pairer
Site rate (re Jestimation O(nir)
Oinkloge)
Old + le elog!)
(Rey Ui.)
A schematic illustration of our FastCherrie//iteRM method. Rate matrix estimation from set of MSAs classically proceeds in two steps: estimation, followed by rate matrix estimation. The recently proposed CherryML method [11] significantly speeds up the rate matrix estimation step. Since CherryML. only requires the cherries in the trees, propose FastCherries, new near-linear method that estimates only the cherries in the tree (as well as the site rates) rather than the whole tree. FastCherries proceeds in steps: divide-and-conquer pairing step based Hamming distance, followed by site rate and branch length estimation. Site rate and branch length estimation alternate until convergence. CherryML's speed allows estimating not only single global rate matrix, but also one rate matrix per site, which we call the SiteRM model. In this schematic, computational complexities for each step of FastCherries is annotated at each step; sequences in the MSA, number of sites the MSA, $ number of states (e.g., 20 for amino acids), number of site rate categories of the LG model (e.g., 4 or 20 is typical), b number of quantization points used quantize time by CherryML. Precomputation of the matrix exponentials which is shared across all MSAs - is excluded from the schematic and costs Orrbs"). MSA illustrations adapted from [23].
Figure I:
number of
We recurse the above procedure S1 and S2. The base case consists of set S with either or sequences. has 2 sequences, we pair them. Otherwise, we return the lone sequence, which will get paired with another lone sequence up in the recursive calls if both |S1 and [8yl are odd. This way, all sequences except at most one (when |S| is odd) will be paired together. Assuming balanced splits of size 0 < < |/1/] < 1/c c for some universal constant the method has near-linear runtime of Olyn logn; ny). Further details can be found in Appendix A.2.
Categorical site-rate variation model. Having paired the sequences, we proceed to estimate site-rates af and quantized distances 1 between pairs. This is done with coordinate ascent until convergence (no change in site-rates) or a maximum of 50 iterations; we find that convergence usually happens after around 10 iterations. We adopt the CAT model of site-rate variation used by FastTree [2], which poses that site-rate of takes of (a hyperparameter) values from (1+2+
1,...,7), with a [[3, 1/3) prior. We use this model of site-rate variation as it is the one used by FastTree, which is one of the fastest phylogenetic tree reconstruction algorithms available and the used in the original CherryML work [1 I]; however, our method can adapted to other models of site-rate variation.
geometric grid between 1/r and r, i.e.,
Site-rate estimation. To estimate site-rates af given quantized divergence time t, for each cherry j and the transition-rale matrix Q, we perform 0-th order optimization on CherryML's composite likelihood. Each site-rate af can be optimized independently. To optimize site-rate af we simply need to find which of the rates in R maximizes that site's site's composite likelihood under the Gamma
5

---

prior. In other words, letting (xx) be the density of the r(3, 1/3) distribution, we optimize for:
of = arg max log ((x) + log(e (Dr Dyle Dy(e/1lD) + log(e" (DDyl Dyluf()) PER
2=1
in time O(bss²), the above objective can be evaluated in time ((cy) = and thus optimized by brute-force search in time (rny)). Since there are a total of ( sites, optimizing the site-rates this way takes time Ollyrny). This is equal to the size of the MSA Dy. multiplied by the number of rate categories r. We initialize the site-rates ag using simple heuristic described in Appendix A.3.
By precomputing all the matrix exponentials (cxp(xт< :
O(ny)
Branch length estimation. To estimate the quantized divergence time 1/ of each pair given the site-rates af and transition-rate matrix Q, we again perform 0-th order optimization on CherryML's composite likelihood. This means optimize for
1f arg max Σ log(e of af (DyLu// Dy Dyll/1l) + logle aQ ( D,
=E(Te:1kk0)
By using precomputed matrix exponentials, this objective can be evaluated in time O(ly) and optimized with brute-force search in time r We empirically find that the objective function above is (quasi-)conceee) which provides a significant speedup. (We could also have done this for optimizing site-rates, but since is already typically small for the CAT model, we found it an unnecessary optimization.) This way, optimizing all quantized divergence times t/ within given family takes O(ny log(b)lz) time. This is equal to the size of the MSA D,, multiplied by log(b).
in æ, that we can optimize it with binary search in time (log(b))y))
Runtime. As detailed Appendix A.4, the above FastCherries algorithm nearly linear the size of the MSA, and applying it to a real MSA with 453,8 19 sequences of length 364 took only about 1000 seconds on CPU core. Furthermore, combined with CherryML to estimate transition rate matrix, this leads to an end-to-end estimation procedure that is nearly linear in the dataset size. Also. FastCherries uses linear space, and thus memory-effieient (up a constant), detailed in Appendix A.S.
3.3 Regularized estimation of site-specific transition-rate matrices under the SiteRM model
We can employ the above algorithm to derive a regularized procedure for estimation under the more richly parameterized SiteRM model. Recall that the SiteRM model uses family- and site-specific transition- matrices Q'- To estimate these matrices, first apply FastCherries algorithm to estimate cherries, distances, and site-rates under the LG model. Given these, we simply apply CherryML's transition-rate matrix estimation procedure [ at each site separately to obtain Q/-
Regularization. Of course, the model is over- parameterized in the sense that there is usually not enough data to estimate the site-specific rate matrices accurately. Indeed, Of has 400 parameters and the only information to support its estimation are the [n,/2] pairs of amino acids at that site. Even with n in the millions, most pairs will not contain mutations at the site since pairs are by definition close sequences, and thus small entries in the rate matrix will be susceptible to large errors. We address this by mixing the empirical counts with pseudocounts from prior. Specifically, let Qo be the transition-rate matrix used for regularizing the model we use the LG transition-rate matrix in applications. We mix the empirical count matrices Cs. with pseudocounts from the LG model with transition-rad rate matrix Qo and site-rates a.. Mixing is performed with a regularization coefficient [0, 1]. Formally, we take
Č = (1- + X||Ck|
where Ps is the 20 x 20 pseudocount matrix with || = 1 given by Pa[z,y] TQ, [x] exp(a! ToQo) )xp)aaTTQQo))z, w). where &Q. is the stationary distribution of Qo This choice has the desirable property that when À we recover the LG model: of a Qo for all i, f. contrast, when À = 0, the prior model is ignored completely. We find that À = 0.5 works well in practice.
The runtime analysis of the above procedure is provided in Appendix A.4.3, Our method can thus be thought of as fine-tuning the LG model to each site of each protein family. This fine-tuning is achieved by mixing empirical counts with pseudocounts from the 'global' LG model.

---

4 Results
4.1 Scalable estimation under the LG model
We applied Ch with FastCherries to the benchmarks from the CherryML paper [11]. These are described in detail in Appendix A.6.
Figure 2a and Figure 2b respectively illustrate the total runtime and accuracy of the end-to-end method for data simulated under the LG model, using log-log plot. It is important to note that while in the original CherryML paper the ground truth trees were used and thus only runtime and accuracy of the rate matrix estimation step were assessed, we benchmarked the end-to-end procedure without access to the ground truth trees. particular, we iterated the process of tree estimation and rate matrix estimation four times as is typical on real data, starting from the uniform rate matrix. Accuracy measured via the median relative error of all the off-diagonal entries in the estimated rate matrix.
Simulations results. As Figure 2b shows, CherryML with FastCherries is one to two orders of magnitude faster than CherryML. with FastTree [2]. Runtimes of the tree estimation step of CherryML with FastCherries are so fast that total runtime for small dataset sizes is dominated by the PyTorch first-order optimizer (whose runtime (gbs²) is independent of input data size). Thus, only for larger dataset sizes does the runtime of CherryML with FastCherries see any noticeable increase. In terms of accuracy, Figure 2 (a) shows that CherryML with FastCherries shows a small asymptotic bias of around 2% median relative error. Otherwise, the relative statistical efficiency of CherryML with FastCherries compared to CherryML using ground truth trees or FastTree is around 50%, meaning that CherryML with FastCherries requires approximately twice as much data to achieve the same error. This is similar to the relative statistical efficiency of Cherry ML compared to MLE with EM.
Real data results. On real data benchmarks, we observed similar end-to-end speedups of one to two orders of magnitude for Cherry ML with FastCherries compared to CherryML with FastTree, while obtaining similar likelihoods on held-out data. This is both for the original LG paper's 's Pfam dataset, shown in Figure 2c,d, well for the QMaker [24] datasets, shown in Supplementary Figure $1.
4.2 Variant Effect Prediction with the SiteRM model
Here, we summarize the performance of our SiteRM model on the task of variant effect prediction. Prior work on applying phylogenetic models to variant effect prediction is described in Appendix A.I.
Data. A notable resource for this purpose is ProteinGym [25], which features dozens of models benchmarked across both deep mutational scanning (DMS) data and human clinical variants. The DMS substitutions benchmark comprises 2.4M variants across 217 DMS assays, while the human clin- ical substitutions benchmark contains 63k variants across 2,525 proteins. On the DMS substitutions benchmark, model pathogenicity scores are evaluated against experimental measurements, and evalu- ation is performed using diverse metrics including Spearman correlation, AUC, MCC, NDCG@ @ 10, and top 10 recall. The metrics are aggregated over all 217 assays, in such way to give equal representation to different kinds of DMS assays, which include protein activity, binding, expression, organismal fitness, and thermostability. The clinical benchmark contains binary pathogenic/benign labels, so that performance is evaluated by averaging AUC across all 2,525 proteins.
Traditional scoring. Probabilistic models of proteins such as Potts models [12] and protein language models (e.g.. ESM [13]) have been used for the task of variant effect prediction. These models provide either a likelihood function p(x)) or, in the case of masked language models, a conditional probability distribution p( (x- where M is the set of masked indices. These models can be used for variant effect prediction by scoring via log-likelihood ratios, as follows, Letting mul denote the mutant and
x the wild-type sequence, p(r) may be used to score mutants via score( ,i) = log p(amr In the case of masked language models, conditional likelihoods p( (o- may be used to score the substitution of FM by z'ar via sc. x11 = loog
p(a"")
p(Z'M x-M] p(zM x=M)
Our new approach. We take a conceptually different approach to variant effect prediction by leveraging probabilistic models of protein evolution p(y) 2x, We propose to score mutants via
p( |x=1, ,t)
score( == log
p( 10,,
7

---

b 10a
CherrMLL FastCherries CherryML Fastiree
CherryML with FastCherries
10a
18437
363004
74027
Median relative error (%g
23.220. 33.220.315.33
15.3
CherryML with Fastiree
CherryML with trees
8687
10.5
10
4381
101
11120.35 100.116.
6.4
Runtime
102
1135
222T
10.2
4.2
3.5
288
2.7
1.7
1.8
324
95
124
191
319
100
2.2
1.3
10
130
189
0.7
67
68
67
70
71
77
10"1
23
2
2
2
Number of families
Number of families
2500
improvement over (.In nats
0.40
Reproduced
Tree Estimetion
Matrix Estimetion
0.35
Average persite AIC:
0.30
FastCherries
0.25
Rurtime (s)
1500
0.20
1000
0.10
0.05
500
WAG
LG rate matrix
LG model CheryMLL with FastTree
LG model
CheryMLL
model
model CherryML
matrix
with FastCherries
CherryLL
with FastCherries
FastTree
Figure 2: CherryML with FastCherries applied to the LG model. (a) End-to-end runtime and (b) median estimation error as function of sample size for CherryML with FastCherries CherryML with FastTree (as well as an oracle with ground truth trees and site rates). Practically, the loss of statistical efficiency for CherryML. with FastCherries relative to FastTree (which perform similarly) is = 50% with small asymptotic bias of around 2%, yet Cherry ML with FastCherries is orders of magnitude faster when applied 1,024 families. The bulk of the end-to-end runtime is taken by rate matrix estimation. The simulation setup is the same as in the CherryML paper [11]. (c) On the benchmark from the LG paper [7]. CherryML with FastCherries yields similar likelihood on held-out families compared to CherryML with FastTree, while (d) shows that CherryML with FastCherries is approximately 20 times faster end-to-end than CherryLM with FastTree, with the bottleneck now being rate matrix estimation.
ground truth trees
where t is hyperparameter that controls the size of the evolutionary neighborhood around the wild-type. Our intuition for why this is good approach is described below.
Consider how the predictions change as varies. Whent 0, obtain Dirac delta at xx1 that all mutants get extremely bad scores. As starts increasing, proteins which are evolutionarily close to the wild-type will see an increase in probability. These are proteins that our model thinks would have arisen from the process of mutation and natural selection, making our conditional probability an excellent candidate for variant effect prediction. Finally, when f +0,, the stationary distribution of the evolutionary model is being used make predictions. In this case, the wild-type is being completely ignored, assuming we are comparing variants with respect to the same wild-type. For classical phylogenetic models, may be interpreted the expected number mutations per site. We use f for variant effect prediction, which means we explore a neighborhood mutation per site. Note that more conserved sites will have smaller site-rates, they may not mutate at all in time unit, while other less conserved sites may mutate multiple times.
one expected
A key insight is that traditional scoring approaches with p(z) may be thought of as a particular case of our approach with / oc. Indeed, given an arbitrary model of p(x) one may construet a model of protein evolution Pevo (ul=, (.,() whose stationary distribution is p(æ) via, for example, Metropolis-Hasting
discrete diffusion model. This insight has been observed before, most notably the works of Evo Velocity [26] and John Ingraham's 's Ph.D. thesis [27], which construct the energy landscape of protein language models and Potts models respectively evolutionary model. Unfortunately, likelihoods under this evolutionary model pevo (y/ x, t) will in general be intractable. Nonetheless, the point is that traditional scoring using p(x mut) is then exactly
continuous-tine

---

Table 1: Despite being an independent-sitee model, SiteRM matches or outperforms many notable
models, including the large protein language model ESM-Iv [13], the alignment- -basedd epistatic
models EVmutation [12] and DeepSequence [28], and the inverse-folding model ESM-IFI [29]. Best
performance among these models are shown in boldface.
Model name
Spearman
AUC
MCC
NDCG
Recall
SiteRM
0.423
0.734
0.330
0.776
0.207
EVmutation
0.395
0.716
0.305
0.777
0.222
DeepSequence (ensemble)
0,419
0.729
0.328
0.776
0.226
ESM-Iv Iv (ensemble)
0.407
0.723
0.320
0.749
0.211
ESM-IFI
0.422
0.730
0.331
0.748
0.223
equivalent to scoring with Pevo (x=""l//"" -1= oo). The analogy goes both ways: our approach can be thought of as a particular case of scoring with models of p(r), but where p(x) is constructed by conditioning an evolutionary model around the wild-type. Indeed, letting Pesa (y) 6) be model of protein evolution, one may define a model via p(x) (ar " f). In this case, scoring with Povo se|x=1 t) is equivalent to scoring with p(z). This way, the key conceptual difference between our approach and traditional approaches is the act of conditioning on the wild-type.
Example. We provide an illustrative example that underscores the power of conditioning on the wild-type, Suppose given site in the MSA for some protein family has four amino acids 1, L, D, E, each with frequency 25%. Note that I and L have similar biochemical properties (they are both non- polar, uncharged, and branched-chain), while and are similar (both being negatively charged). Suppose that mutations between / and L or between D and E at this site are tolerated, while other mutations among these residues are deleterious. In this case, using log log p(="mut) may not be able tell apart the benign mutations at this site from the two pathogenic ones. Indeed, in the extreme case of the site-independent EV mutation model [12], all variants this site will have exactly the same score. In contrast, our SiteRM model, despite being an independent-site
model, will learn and predict that only mutations that preserve the biochemical properties of the residue will be tolerated. Thus, despite being trained on exactly the same data, our SiteRM model has the built-in capacity to make sensible predictions of this kind whereas EVmutation may not. Although a highly expressive model might be able to succeed at this toy example by leveraging the correlations between sites, it is evident that an evolutionary model even if site-independent - is better placed to leverage insights like this one.
Variant effect prediction results. In the case of the SiteRM model with t = 1, we use
score( (v:
wi log p(
J: 1 1)
log(exp(Q/ |x (x.""""""".!!! (
to score variants with respect to the same wild-type sequence xt, Note that x may differ from "t at multiple sites, and (1) is an independent-sire sites model, which uses additive per-site effects.
We trained the SiteRM model on the MSAS provided by the ProteinGym benchmark. Duplicate sequences were removed. We used FastCherries with the LG rate matrix and 20 rate categories to estimate cherries and site rates. Then, we use CherryML. to estimate site-specfiic rate matrices, with the LG rate matrix as prior and with regularization strength X = 0.5 such that exactly half the data comes from pscudocounts and the other half from empirical counts. We found that modeling the gap as its own state improved performance, so that our best model uses 21 x 21 site-specific rate matrices, with the 21> 21 rate matrix learnt from the LG paper's Pfam dataset as the prior.
Table summarizes the DMS substitution benchmark results for some select notable models. Despite being an independent-sites model, SiteRM outperforms the large protein language model ESM- 1v and its inverse- prerssefoldiig version ESM-IF1 on Spearman correlation, AUC, MCC and NDCG. SiteRM also outperforms EVmutation [12], which is based on the seminal Potts model with pairwise interactions, as well as the DeepSequence ensemble model [28]. Only on Recall@ 10 do we observe less competitive results. We attribute the remarkable performance of our independent-sites model to conditioning on the wild-type, which is absent from competing approaches, and our ability to estimate family- and site-specific rate matrices without overfitting. Full results with all models reported in the ProteinGym paper [25] are provided in Supplementary Table S1. The table also shows ablations for the SiteRM model where we reduce the size of the training dataset by subsampling sequences in the MSA, by using FastTree instead of FastCherries, by excluding gaps. The

---

largest MSA processed by FastCherries had approximately 454K distinct sequences, each of length 364. FastCherries took approximately 1,000 seconds on a single CPU core to estimate cherries and site rates for this MSA. Of these, 30 seconds were spent on the divide-and-conge second was spent per site-specific
pairing step, and
970 seconds on distance and site rate estimation. Subsequently, rate matrix estimation, meaning a total time of around 1,500 seconds. Extrapolating our estimates shown in Supplementary Table $1, FastTree would have taken around 50 times longer.
Finally, on the human clinical benchmark, Supplementary Table S2 shows that SiteRM achieves an AUC of 0.911, which is ~ 0.02 higher than ESM-Ib [13], and less than 0.01 below the state-of-the-ar
5 Discussion
We have introduced an end-to-end method estimate amino acid substitution rate matrices from MSAs alone. By combining CherryML [11] with our FastCherries algorithm, it achieves near-linear runtime. Through simulations, we rigorously studied the computational and statistical efficiency of the method, finding it to be around 10 to 100 times faster than CherryML with FastTree [2], all while being only 50% less efficient and having small asymptotic bias of around 2% median relative error which we find negligible in practical applications. CherryML with FastCherries easily runs on MSAs with hundreds of thousands of sequences. We plan to contribute to the open-source CherryML package with our FastCherries method so that all researchers may use the end-to-end pipeline easily.
By leveraging FastCherries' scalability, we performed regularized inference under the SiteRM model, providing site-specific rate matrices for given MSA. We applied this independent-sire sites model of protein evolution to the task of variant effect prediction and found that it outperforms many well- established models such as deep protein language models. We believe that this seemingly paradoxical result that an independent-site model outperforms epistatic models is largely explained by our new conceptual approach, which conditions model of protein evolution around the wild-type.
Future directions. Our work leaves several intriguing avenues for future research. First, although in this work we focused on classical independent-sites models of protein evolution, the large amount of data generated quickly by our FastCherries method may be used to train more complex models of protein evolution that go beyond the independent-sites assumption, such as deep neural networks. These may combine the best of both worlds; the power of modeling the evolutionary process in time and the ability to take into account complex correlations between protein sites. Our work also opens up new avenue in statistical phylogenetics by enabling tree reconstruction with site-specific rate matrices. Software such as IQTree [5] already allows phylogenetic tree reconstruction under a 'partition model', where sites may be grouped such that all sites in the same group evolve under the same transition-rate matrix. Exploring extreme cases of the partition model where each site is its own partition, as in our SiteRM model, is an exciting and new avenue of research that promises to improve phylogenetic tree reconstruction, and find its way into other applications such as ancestral sequence reconstruction. To this end, we plan to release the site-specific rate matrices we have estimated for all 15,051 protein families in the TrRosetta dataset [30]. More broadly, our idea of bypassing tree estimation through the lens of composite likelihood may enable further methodological developments in other related models, such as mixture models [31]. All in all, we expect our work to have lasting impact in both statistical phylogenetics and computational variant effect prediction.
Limitations. One limitation of our current method is that it assumes- just like essentially all of phylogenetics - that the model of protein evolution is time-reversible. thereby constraining the set of possible rate matrices. This may be an important source of model misspecification in some cases, and there is an emerging body of work extending phylogenetic methods to iirevirreeersib models [32]. Although we find that our method has outstanding performance on VEP, even when compared to full tree reconstruction methods such as FastTree, our method does exhibit a small but non-zero amount of asymptotic bias which may matter in some other downstream applications. For example, the rate matrices estimated with CLerryML using FastCherries may yield different, less accurate tree topologies from those obtained with CherryML using FastTree or other tree reconstruction methods (such as PhyML [4]). The extent to which this is true requires further investigation.
Acknowledgments and Disclosure of Funding
This research is supported in part by NIH grants R56-HG013117 and R01-HG013117. S.P. would like to acknowledge Anastasia Ignatieva for early discussion about the project's idea.
10

---

References
[1] Alexandros Stamatakis. RAxML version tool for phylogenetic analysis and posl-analysis of large phylogenies. Bioinformatics, 3019):111-1313, 01 2014.
[2] Morgan N. Price, Paramvir S. Dehal, and Adam P. Arkin. Fasttree 2 approximately maximum- likelihood trees for large alignments. PLoS ONE, 5(3): Mar 2010.
[3] Fredrik Ronquist and John P. Huelsenbeck, MrBayes 3: Bayesian phylogenetic inference under mixed models. Bioinformatics, 19(12):1572-157) 08 2003.
[4] Stéphane Guindon, Jean- -rann-Françoi Dufayard, Vincent Lefort, Maria Anisimova, Wim Hordijk, and Olivier Gascuel. New Algorithms and Methods to Estimate logenies: Assessing the Performance of PhyML 3.0. Systematic Biology, 59(3):307-321, 05 2010.
Phy-
[5] Bui Quang Minh, Heiko A Schmidt, Olga Chernomor, Dominik Schrempf, Michael D Wood- hams, Arndt von Haeseler, and Robert Lanfear. IQ-TREE 2: New Models and Efficient Methods for Phylogenetic Inference in the Genomic Era. Molecular Biology and Evolution. 153::1550-153,, 02 2020.
[6] Emma B. Hodcroft, Nicola De Maio, Rob Lanfear, Duncan R. MacCannell, Bui Quang Minh, Heiko Schmidt, Alexandros Stamatakis, Nick Goldman, and Christophe Dessimoz. Want to track pandemic variants faster? Fix the bioinformatics bottleneck. Nature, 591(7848):30-33, March 2021.
[7] Si Quang Le and Olivier Gascuel. An Improved General Amino Acid Replacement Matrix. Molecular Biology and Evolution, 25(7):1307-1320 03 2008
[8] S. Roch. A short proof that phylogenetic tree reconstruction by maximum likelihood is hard. IEEE/ACM Transactions on Computational Biology and Bioinformatics, ((01):92-94, jan 2006.
[9] Peter S. Klosterman, Andrew Uzilov, Yuri R. Bendaña, Robert K. Bradley, Sharon Chao, Carolin Kosiol, Nick Goldman, and Ian Holmes. Xrate: a fast prototyping, training and annotation tool for phylo-grammars. BMC Bioinformatics, 7(1):428, 2006.
[10] Ian H Holmes. Historian: accurate reconstruction of ancestral sequences and evolutionary rates. Bioinformaticx, 3378): 01 2017.
[11] Sebastian Prillo, Yun Deng, Pierre Boyeau, Xingyu Li, Po-Yen Chen, and Yun S. Song. Cher- ryML: scalable maximum likelihood estimation of phylogenetic models. Nature Methods, 20( 1236, August 2023.
[12] Thomas A. Hopf, John B. Ingraham, Frank J. Poelwijk, Charlotta P. I. Schärfe, Michael Springer, Chris Sander, and Debora S. Marks. Mutation effects predicted from sequence co-variation. Nature Biotechnology. 35(2):128- -25, February 2017.
[13] Joshua Meier, Roshan Rao, Robert Verkuil, Jason Liu, Tom Sercu, and Alex Rives. Language models enable shot prediction of the effects of mutations protein function. In M. Ranzato, A. Beygelzimer, Y. Dauphin, P.S. Liang, and J. Wortman Vaughan, editors, Advances in Neural Information Processing Systems, volume 34, pages 29287-29300 Curran Associates, Inc., 2021.
[14] M. O. Dayhoff and R. M. Schwartz. Chapter 22: A model of evolutionary change in proteins. In in Atlas of Protein Sequence and Structure, 1978.
[15] David T. Jones, William R. Taylor, and Janet M. Thomton. The rapid generation of mutation data matrices from protein sequences. Comput. Appl. Biosci,, 8(3):275-222, 1992.
[16] Simon Whelan and Nick Goldman. A General Empirical Model of Protein Evolution Derived from Multiple Protein Families Using a Maximum-Likelihood Approach. Molecular Biology and Evolution, 18(5):691-69) 05 2001.
[17] Ziheng Yang. Maximum likelihood phylogenetic estimation from dna sequences with variable rates over sites: approximate methods. Journal of Molecular Evolution, 39(3):306-33) -14, 1994.
11

---

[18] Cristiano Varin, Nancy Reid, and David Firth. An overview of composite likelihood methods. Statistica Sinica, 21(1):5-42, 2011.
[19] Adam Paszke, Sam Gross, Soumith Chintala, Gregory Chanan, Edward Yang, Zachary Devito, Zeming Lin, Alban Desmaison, Luca Antiga, and Adam Lerer. Automatic differentiation in pytorch. In Advances in Neural Information Processing Systems 30, 2017.
[20]
Martín Abadi, Ashish Agarwal, Paul Barham, Eugene Brevdo, Zhifeng Chen, Craig Citro, Greg S. Corrado, Andy Davis, Jeffrey Dean, Matthieu Devin, Sanjay Ghemawat, Ian Goodfellow, Andrew Harp. Geoffrey Irving, Michael Isard, Yangqing Jía, Rafal Jozefowicz, Lukasz Kaiser, Manjunath Kudlur, Josh Levenberg, Dandelion Mané, Rajat Monga, Sherry Moore, Derek Murray, Chris Olah, Mike Schuster, Jonathon Shlens, Benoit Steiner, Ilya Sutskever, Kunal Talwar, Paul Tucker, Vincent Vanhoucke, Vijay Vasudevan, Fernanda Viégas, Oriol Vinyals, Pete Warden, Martin Wattenberg. Martin Wicke, Yuan Yu, and Xiaoqiang Zheng. TensorFlow: Large- scale machine learning on heterogeneous systems, 2015. Software available from tensorflow.org.
[21] Philipp Bader, Sergio Blanes, and Fernando Casas. Computing the matrix exponential with an optimized taylor polynomial approximation. Mathematics, 7(12), 2019.
[22] Diederik P Kingma and Jimmy Ba. Adam: A method for stochastic optimization, In Proceedings of the 3rd International Conference on Learning Representations (ICLR), 2015.
[23] Nick Bhattacharya, Neil Thomas, Roshan Rao, Justas Daupras, Peter K Koo, David Baker, Yun S. Song, and Sergey Ovchinnikov. Single layers of attention suffice predict protein contacts, 2021.
[24] Bui Quang Minh, Cuong Cao Dang, Le Sy Vinh, and Robert Lanfear. QMaker: Fast and Accurate Method to Estimate Empirical Models of Protein Evolution. Systematic Biology, 1065):004-1106, 02 2021.
[25] Pascal Notin, Aaron Kollasch, Daniel Ritter, Lood van Nickerk, Steffanic Paul, Han Spinner, Nathan Rollins, Ada Shaw, Rose Orenbuch, Ruben Weitzman, Jonathan Frazer, Mafalda Dias, Dinko Franceschi, Yarin Gal, and Debora Marks. Proteingym: Large-scale benchmarks for protein fitness prediction and design. In A. Oh, T. Neumann, A. Globerson, K. Saenko, M. Hardt, and S. Levine, editors, Advances in Neural Information Processing Systems, volume 36, pages 64331-64379. Curran Associates, Inc., 2023.
[26] Brian L. Hie, Kevin K. Yang, and Peter S. Kim. Evolutionary velocity with protein language models predicts evolutionary dynamics of diverse proteins. Cell Systems, 13(4):24--28..c6 2022.
[27] Probabilistic Models of Structure in Biological Sequences ProQuest.
[28] Adam J. Riesselman, John B. Ingraham, and Debora S. Marks. Deep generative models of genetic variation capture the effects of mutations. Nature Methods, 15(10):816-82) 822, October 2018.
[29] Chloe Hsu, Robert Verkuil, Jason Liu, Zeming Lin, Brian Hie, Tom Sercu, Adam Lerer, and Alexander Rives. Learning inverse folding from millions of predicted structures. In Kamalika Chaudhuri. Stefanie Jegelka, Le Song. Csaba Szepesvari, Gang Niu, and Sivan Sabato, editors, Proceedings of the 39th International Conference on Machine Learning, volume 162 of Proceedings of Muchine Learning Reseurch, pages 8946-8970. PMLR, 17-23 Jul 2022.
[30] Jianyi Yang, Ivan Anishchenko, Hahnbeom Park, Zhenling Peng, Sergey Ovchinnikov, and David Baker. Improved protein structure prediction using predicted interresidue orientations. Proceedings of the National Academy of Sciences, 1 (17()):1496-1))
2020.
[31] Si Quang Le. Cuong Cao Dang, and Olivier Gascuel. Modeling Protein Evolution with Several Amino Acid Replacement Matrices Depending on Site Rates. Molecular Biology and Evolution, 29(10):2911-2936 04 2012.
[32] Cuong Cao Dang. Bui Quang Minh, Hanon McShea, Joanna Masel, Jennifer Eleanor James, Le Sy Vinh, and Robert Lanfear. nQMaker: Estimating Time Nonreversible Amino Acid Substitution Models. Systematic Biology, 7105)) 1105)1110-1112,, 02 2022.
12

---

[33] Cuong Dang, Le Vinh. Olivier Gascuel, Bart Hazes, and Quang Le. Fastmg: a simple, fast, and accurate maximum likelihood procedure estimate amino acid replacement rate matrices from large data sets. BMC bioinformatics, 15:341, 10 2014.
[34] Naruya Saitou and Masatoshi Nei. The neighbor-joining method: a new method for reconstruct- ing phylogenetic trees. Molecular biology and evolution, 4(4):406-425, 1987.
[35] Elodie Laine, Yasaman Karami, and Alessandra Carbone. GEMME: Simple and Fast Global Epistatic Model Predicting Mutational Effects. Molecular Biology and Evolution, 36(1 1):2604 2619, 08 2019.
[36] Kazutaka Katoh and Hiroyuki Toh. PartTree: algorithm build approximate tree from large number of unaligned sequences. Bioinformatics, 23(3)377--774, 11 2006.
[37] Jung-Eun Shin, Adam J. Riesselman, Aaron W. Kollasch, Conor McMahon, Elana Simon, Chris Sander, Aashish Manglik, Andrew C. Kruse, and Debora S. Marks. Protein design and variant prediction using autoregressive generative models. Nature Communications, 12(1):2403, April 2021.
[38] Jonathan Frazer, Pascal Notin, Mafalda Dias, Aidan Gomez, Joseph K. Min, Kelly Brock, Yarin Gal, and Debora S. Marks. Disease variant prediction with deep generative models of evolutionary data. Nature, 599(7883):91-95, November 2021.
[39] Ethan C. Alley, Grigory Khimulya, Surojit Biswas, Mohammed AlQuraishi, and George M. Church. Unified rational protein engineering with sequence-basee deep representation learning. Nature Methods, 16(12):1315-3333 December 2019.
[40] Kevin K. Yang, Nicolo Fusi, and Alex X. Lu. Convolutions are competitive with transformers for protein sequence pretraining. Cell systems, 15(3):866-24.e2,, 2024.
[41] Daniel Hesslow, Niccoló Zanichelli, Pascal Notin, Iacopo Poli, and Debora Marks. Rita: a study on scaling up generative protein sequence models. arXiv preprint arXiv:2205.08 057889 2022.
[42] Erik Nijkamp. Jeffrey A. Ruffolo, Eli N. Weinstein, Nikhil Naik, and Ali Madani. Progen2: Exploring the boundaries of protein language models. Cell systems, 14((1):968-97) 9716689-78...... 2023.
[43] Céline Marquet, Michael Heinzinger, Tobias Olenyi, Christian Dallago, Kyra Erckert, Michael Bernhofer, Dmitrii Nechaev, and Burkhard Rost. Embeddings from protein language models predict conservation and variant effects. Human Genetics, 141(10):1699-144) October 2022.
[44] Roshan M Rao, Jason Liu, Robert Verkuil, Joshua Meier, John Canny, Pieter Abbeel, Tom Sereu, and Alexander Rives. Msa transformer. In Marina Meila and Tong Zhang, editors, Proceedings of the 38th International Conference on Machine Learning, volume 139 of Proceedings of Machine Learning Research, pages 8844-8856. PMLR, 18-24 Jul 2021.
[45] Pascal Notin, Mafalda Dias, Jonathan Frazer, Javier Marchena-Hurtadee Aidan Gomez, Debora S Marks, and Yarin Gal. Tranception: protein fitness prediction with autoregressive transformers and inference- time retrieval. arXiv (Cornell University), 2022.
[46] Pascal Notin, Lood Van Niekerk, Aaron W Kollasch, Daniel Ritter, Yarin Gal, and Debora Susan Marks. TranceptEVE: Combining family-specific and family-agnostic models of protein sequences for improved fitness prediction. In NeurIPS 2022 Workshop on Learning Meaningful Representations of Life, 2022.
[47] Kevin K Yang, Niccolò Zanichelli, and Hugh Yeh. Masked inverse folding with sequence transfer for protein representation learning. Protein Engineering, Design and Selection, 36:gzad015, 10 2022.
[48] J. Dauparas, I. Anishchenko, N. Bennett, H. Bai, R. J. Ragotte, L. F. Milles, B. I. M. Wicky, A. Courbet, R. J. de Haas, N. Bethel, P. J. Leung, T. F. Huddy, Pellock, D. Tischer, Chan, B. Koepnick, H. Nguyen, A. Kang, B. Sankaran, A. K. Bera, N. P. King, and D. Baker. Robust deep Bamming-based protein sequence design using proteinmpnn. Science, 378(6615) 2022.
13

---

A
Appendix
A.1 Related Work
Speeding up the tree reconstruction step required by rate matrix estimation has been of interest in recent years as genomic datasets have increased in size. In this direction, one work related ours is that of FastMG [33], which proposes naively splitting the input MSAs into sub-MSAs, and then running an end-to-end estimation procedure using these MSAs. The MSA splitting procedure is performed by using partial tree reconstruction with the Neighbor- to the Hamming distance. In this sense, their method has similar flavor ours, since decomposing each MSA into pairs of sequences. However, there are a number of key differences between approach and that of FastMG. Firstly, because FastMG proposes to naively split the MSAs into sub-MSAs and treat each independently, FastMG cannot split the MSA into sub-MSAs that are too small. Indeed, it is statistically infeasible to estimate site-rates reliably from too few sequences, since the site-rate at each position is only supported by the few observations at that site. In particular, it is impossible to estimate the site-rates reliably with FastMG by splitting each MSA into pairs of sequences, since only two amino acids would support the rate that site. Therefore, FastMG typically uses sub-MSAs of size 16. Because of this, FastMG still requires full phylogenetic tree reconstruction (albeit smaller trees of size 16). A major limitation that FastMG does not estimate site-rates for the original MSAs - it only does so for the sub-MSAs. This means that estimation under the SiteRM model cannot be naturally performed, since the site-rates required by the IGG prior we use as a regularizer. Our method can thus be thought of as a version of FastMG with parameter sharing, wherein sub-MSAs of size 2 share the site-rate parameters, thus allowing to capture rich site-rate variation reliably and estimate site-rates for the original MSAs. Finally, FastMG uses NJ to decompose the MSA, whose runtime is in general O(n²), while we propose to use a near- linear time On log n) algorithm. Together, these ingredients lead to our near-linear time algorithm with massive scalability.
(NJ) algorithm [34] applied
The application of phylogenetic methods to variant effect prediction is not new. A notable method in this space is GEMME [35], which uses phylogenetic trees and hand-crafted statistics for the purpose of variant effect prediction. This method performs on-par with the state-of-the-art method TranceptEVE for some metrics on the DMS and clinical substitutions benchmark. The main limitation of GEMME is that it is based on cand-crafted features and lacks a probabilistic basis, limiting its ability to scale larger datasets in seamless way. We expect that by training more sophisticated. epistatic models of protein evolution on the millions of training pairs (z, t) generated by our method, one may be able to scale up our approach and further improve on variant effect prediction.
A.2 Further details of the divide-and-conquee algorithm
Although the divide-and-conquer algorithm described in Section 3.2 does automatically guarantee balanced splits, they can be achieved by modifying the splitting procedure by first sorting the points based on the ratio d(z, 21 ///d(z 22) [0, +oo], assigning the (E)) ones with smallest ratio to S₁ and the O(e) ones with greatest ratio to S2; the remaining points are assigned with the original procedure. Moreover, this modified procedure only needs to be applied if the original split is unbalanced to begin with. Empirically, enforeing theoretically balanced splits completely unnecessary, and we use the original simple splitting criteria, which empirically achieves near-linear runtime on real data inputs. A similar algorithm called PartTree [36] has been proposed in the past for scalable phylogenetic tree reconstruction, but it has seen little adoption due to its inferior topological accuracy compared methods [2). It is important note that we care about rate matrix estimation, not tree estimation (since the trees T are nuisance parameters in our model), which is why such near-linear time algorithm turns out be profitable in setting. As we demonstrate in the Results section, the tree estimates (specifically, cherries) obtained this way are good enough for our purposes.
A.3 Site-rate initialization
To initialize the site-rates af in our FastCherrics method, we use the following heuristic that tries to align the initial site-rates with the I(3,, 1/3) prior. First, for each site i, we count the number of pairs
14

---

(x. supporting mutation at that site:
muts!
E
1D + D
K<jckens
which can be easily computed in time Onnyly 1ys) by just counting the number of occurrences of each state at each site. Gaps are treated as missing data and are ignored when computing the above quantity. Next, we sort the sites based on the number of such mutations in time Oll (gg()) Recall that the site-rates in the CAT model with r rate categories are geometrically spaced between 1/r and r. Furthermore, 13,3, 1/3) prior is used. Call these rates 7 Tr. Using our sorted order for sites based number of mutations, we assign number of sites rate category % proportional
to the measure of the [33, 1/3) distribution between 1i-1% and Vii++i,, where "o = O and 7++1 +00.
A.4 Runtime
A.4.1
Runtime of the FastCherries algorithm
Letting be the number of coordinate ascent steps used for optimizing site-rates and divergence
times, the runtime of our FastCherries algorithm described in Section 3.2 for family f is
Oly logn/
rbs
+ tnyly + lys log
(Lyrns + +ny log(b)())),
Pairing step
Matris exponcntial
Site initialization
Coondinate ascent
Site
Distance estimation
precompuration
Fate
Which is nearly linear in nyly, the size of the MSA Dy. Large applications encountered in practice have n = 450, 000, ly = 300, r = 20, = 100, 10, s = 20. For a real MSA with over 450,000 sequences of length 364, we find that our method takes around 1000 seconds on CPU core.
A.4.2 End-to-end runtime
Finally, combined with CherryML's procedure to estimate a transition-rate matrix Q given pairs (uf, un ), divergences tf and site-rates af, we obtain an end-to-end procedure for estimating a rate matrix Q under the LG model. Assuming for simplicity as in CherryML's analysis that families have the same number of sequences nt n and lengths l = 1, denoting by the number of first-order gradient descent steps of CherryMI. (implemented in PyTorchh using the Adam optimizer), q the number of iterations of the whole process (outer coordinate ascent), and recalling that m is the number of families, obtain an end-to-end runtime of
O
mnd log +
rbs exponential precompotation
+ mnn + mls + ml log
Outer coordinte ascent
Pairing stop
Site Fale initialization
mndr mond log(b) mnd log b
gbs² D)-
Inner coordinate
Site rate estimation
Distance estimation
CherryMI counting step
Cherry
fusst-order
escent
optimizer
(PyTorch)
In practice, in the CherryML work, q 4,g 2000. Thus, this end-to-end runtime is nearly
linear in the dataset size mni up to logarithmic factor q log n and a constant factor ga(r + log(b)). The procedure is also embarrassingly parallel over families f.
A.4.3 Runtime for estimating site-speeific transition-rate matrices
Empirically, on a large MSA with around 450,000 sequences length 364, it takes around 1000 seconds to estimate branch length and site-rates, and thereafter, it takes around 1 second per site to estimate site-specific transition-rate matrix Q! using the method described in Section 3.3. Thus, for an MSA with 450,000 sequences of length 364, it takes approximately 1500 seconds end-to-end to estimate site- site-specfii rate matrices. We should remark that unlike in the original CherryML procedure, we (1) initialize the optimizer with af Qo instead of their JTT-IPW strategy [11], and (2) use only 100 steps of the Adam optimizer [22] instead 2000. Indeed, since we are regularizing with a Qo this is reasonable initialization, and further, we find that 2000 steps is unnecessarily large, with
15

---

convergence already happening at around 100 steps. This observation is helpful to further boost scalability.
A.5 Note space complexity
Our method FastCherries uses linear space, Briefly, the logarithmic factors in the computational runtime come from divide-and-conquer and binary searches, which do not translate to the space complexity. In more detail, the space complexity when processing the m MSAs in parallel using p processes is brpps"b-+-plll pls). Here, (ps"br) is the cost of the precomputed matrix expomentials, which is shared across all MSAs (and can also be shared across all processes, in which case the space drops to O(s²br) The second term (pnl)) is the cost to load the p MSAs into memory. For each MSA, pairing takes O(n) space, since takes ((n) space to store the distance between the pivot and all other sequences. The space for computing the initial site rates for an MSA is Olss + r). Here, the (Is) term comes from the count matrix count (8, which tracks the number of times that character j occurs at position i. Computing site rates given branch lengths for an MSA takes ( space, and computing branch lengths given site rates takes (L)) space. In both cases, we only track the best rate for the current site or best branch length for the current site. This way, our method FastCherries is space efficient (up to a constant).
A.6 CherryML benchmark details
There are kinds of benchmarks in the CherryML work [11] which we leverage: a simulated data benchmark and a real data benchmark. In the simulated data benchmark, a ground truth rate matrix (the LG rate matrix) is used to simulate MSAs from ground truth trees and site rates. Runtime and accuracy of rate matrix estimation methods are evaluated as a function of the number of MSAs simulated, thus providing insight into the computational runtime and the statistical efficiency of each method. The real data benchmark is a compilation of benchmarks originally reported in the LG paper [7] and in the QMaker paper [24]. In these benchmarks, real MSAs are used to estimate rate matrix. The quality of the estimated rate matrix is judged by the lig-likelihood of MSAs from held-out families, thus probing the quality of the method in real-life settings. Importantly, these real data benchmarks contain MSAs from diverse areas of life (mammalian MSAs, bird MSAs, insect MSAs, and plant MSAs), thus probing the generalization capabilities of the method to different real biological settings.
A.7 QMaker datasets benchmark
The QMaker work [24] uses diverse datasets to highlight the importance of learning rate matrices which are specific to different domains of life. The CherryML work [11] uses these datasets to showcase the ability of the CherryML method to generalize to these diverse datasets. Here, we use these datasets in the same way to highlight the ability of CherryML with FastCherries to achieve CherryML. with FastTree;; in the CherryML. paper, it was shown that
comparable performance CherryML. with FastTreee matches the performance of EM with FastTTree Supplementary Figure S1 shows the results, with log- lillliikeooods the left and runtimes on the right. We can that indeed CherryML with FastCherries achieves comparable performance to CherryML with FastTree, all while being 10- 100 times faster end-to-end.
The details for this benchmark are as follows: starting from the LG rate matrix, tree estimation (with either FastTree or FastCherries) and rate matrix estimation with CherryML were iterated 4 times. The number of rate categories used is 4 as in the CherryML paper and the original LG paper [7]. The learnt rate matrices are then evaluated by computing the log- likeliood of the MLE trees fitted using FastTree with 4 rate categories.
should be noted that the 1000 MSAs in each of these datasets are relatively small, with less than 100 sequences each, which is why FastCherries is so fast and the runtime is barely noticeable. Most of the runtime is spent on CherryML's rate matrix estimator which is implemented in Py PyTorch [19) using 2000 epochs of the Adam [22] optimizer. In our experience, this is more than enough to achieve convergence, and why we lower it to 100 epochs when training the SiteRM model.
16

---

Avaraga per site op dibe lloood
Improvenaan over 1att
Stondard rate manrix
FasTTree
Parti Enttration
3020
0.05
2000
2020
-0.0
-C20
rate matrix
with
with
with
maliix
CherryMi
Cherryi
Cherryll
CherryAt
bird
bid
Average par paraata tog Witthehthddd
IF corooamtnnt over In nats
Fastiree
rete
Rote Varte Eabmation
Iee Esteoriinn
FastCheriies
6,00
80C0
Bunlmm
-030
8000
2000
WAG rate mabik
LG rate mabiix
LG with CharrVLL FasIThce
LG with ChanyML FastCherries
LG with
LG with CharryML FastCherrie
ChanyML FasITroe
irarvnal
mammal
Auaroie par site los Wiballlloooo 1ood
Improvemnnt Oour malts
0.1
36200
Extimation
51000
Bunkime
-04
Standard
4200
FestCherries
2000
matrix
CherrvML wet Fast lree
CherryMi with FastCherries
CherryML with Fastiree
CherryAt with Faetrierries
Averaga sandait log Bellhood provoaateett ouer JTT, mats
FastTice
Standard rate matrix
4320
Rute Marts balmatin
FastCherries
3520
022
3500
2000
0
with CheryMIL
with CherryMLL
with ChenyML
with CherryML Fasti harries
mabrik
FastTTee
Supplementary Figure S1: CherryML with FastCherries matches the performance of Cherry ML with FrstTree on diverse datasets. We reproduced the QMaker results from the CherryML paper [11] and added Ml. with FastCherries method. We observe similar held-out likelihood at fraction of end-to-end runtime, which is about 10- 100 times faster.
17

---

A.8 Full DMS substitutions results table
Here we provide the full table of results containing all the models reported in the ProteinGym benchmark [25], with our SiteRM model added.
Supplementary Table S1: ProteinGym DMS substitutions benchmark results, with our SiteRM model and ablations added under the 'Protein Evolution, Classical' model type. The ablations explore using less data by subsampling the MSAs, or using FastTree instead of FastCherries, gaps as a state. We also benchmark the LG model, which uses a global matrix and does quite well already. Despite being independent-siter model, SiteRM outperforms all protein language models expect VESPA, on all metrics except occasionally Recall @ 10. The approximate CPU time in seconds required for tree/cherry estimation on the largest MSA is shown in the last column 'slowest tree'. Empirically, FastCherries around 50 times faster than FastTree. Estimating site-speeific rate matrices takes around 1 additional second per site.
treating
Model type
Model name
Spearman
AUC
MCC
NDCG
Recall
Slowest tree
Protein
LG, FastCherries 1M
0.377
0.710
0.297
0.763
0.185
1,000 sec
Evolution,
SiteRM, FasTTree 1K
0.397
0.721
0.311
0.774
0.203
350 sec
Classical
SiteRM, FastCherrics 1K
0.398
0.721
0.312
0.772
0.205
sec
SiteRM, FasTTree 10K
0.404
0.724
0.315
0.776
0.204
3,100 sec
SiteRM, FastCherries 10K
0.408
0.726
0.319
0.777
0,206
60 sec
SiteRM, FastCherries IM
0.412
0.728
0.321
0.777
0.207
1,000
SiteRM, FastCherries IM, w/gaps
0.423
0.734
0.330
0.776
0.207
1,000 sec
Alignment-
Site-Independent [12]
0.359
0.696
0.286
0.747
0.201
based
WaveNet [37]
0.373
0.707
0.294
0.761
0,203
EVmutation [12]
0.395
0.716
0.305
0.777
0,222
DeepSequence (ensemble) [28]
0.419
0.729
0.328
0.776
0.226
EVE (ensemble) [38]
0.439
0.741
0.342
0.783
0.230
GEMME [35]
0.455
0.749
0.352
0.777
0.211
Protein
UniRep [39]
0.190
0.605
0.147
0.647
0.139
language
CARP (640M) [40]
0.368
0.701
0.285
0.748
0.208
ESM-Ib [13]
0.394
0.719
0.311
0.747
0,203
ESM-2 (15B) [13]
0.401
0.720
0.314
0.759
0,208
RITA [41]
0.372
0.707
0.293
0.751
0.193
ESM-1v (ensemble) [13]
0.407
0.723
0.320
0.749
0,211
ProGen2 XL. [42]
0.391
0.717
0.306
0.767
0.199
VESPA [43]
0.436
0.742
0.346
0.775
0.201
Hybrid
UniRep evotuned [39]
0.347
0.693
0.274
0.739
0.181
MSA Transformer (ensemble) [44]
0.434
0.738
0.340
0.779
0.224
Tranception L [45]
0.434
0.739
0.341
0.779
0.220
TranceptEVE [46]
0.456
0.751
0.356
0.786
0.230
Inverse Folding
ESM-IF1 [29]
0.422
0.730
0.331
0.748
0,223
MIF-ST [47]
0.401
0.718
0.311
0.766
0.226
ProteinMPNN [48]
0.258
0.639
0.196
0.713
0.186
18

---

A.9 Clinical substitutions results table
Our independent-sites SiteRM model outperforms the large protein language model ESM- lb at human climical substitutions variant effect prediction (VEP) by nearly 0.02 AUC points, and is less than 0.01 AUC points away from state-of-the-ar methods. We again find that modeling gaps as a state in SiteRM (thus learning 21 x 21 transition matrices) helps.
Supplementary Table S2:
Model name
Mean AUC
TranceptEVE [46]
0.920
GEMME [35]
0.919
EVE [38]
0.917
SiteRM
0.911
SiteRM w/o gaps
0.903
ESM-1b [13]
0.892
A.10 Compute Resources
The experiments in section 'Scalable estimation under the LG model" were performed on 'MacBook Pro' model 'MacBookPrro8.2 with 'Apple MI Max' chip with 10 CPU cores and 32 GB RAM. The ProteinGym variant effect prediction benchmarks were performed on a machine with 32 CPU cores, model 'Intel Xeon Skylake 6230 @ 2.1 GHz', and with 384 GB RAM; however, this is more RAM than necessary and the experiments can likely be performed comfortably with much less RAM. The total compute required by the project in the paper, which was dominated by variant effect prediction on human clinical substitutions, which took around day when parallelized with 32 cores and was done last. The large runtime is due the large number (2525) of clinical MSAs.
not much larger than that of the experiments reported
19

---

NeurIPS Paper Checklist
1. Claims
Question: Do the main claims made in the abstract and introduction accurately reflect the paper's contributions and scope?
Answer: [Yes]
Justification: Extensive benchmarking on simulated and real data from prior work (Cher- ryML) are used to demonstrate the speed and accuracy of the method, and benchmarks the extensive community resource ProteinGym support the variant effect prediction results. The theoretical complexity of the algorithm is also analyzed to provide a proof of the near-linear runtime.
Guidelines:
The answer NA means that the abstract and introduction do not include the claims made in the paper.
The abstract and/or introduction should clearly state the claims made, including the contributions made in the paper and important assumptions and limitations. A No or
NA answer to this question will not be perceived well by the reviewers.
The claims made should match theoretical and experimental results, and reflect how much the results can be expected to generalize to other settings.
It is fine to include aspirational goals as motivation long is clear that these goals are not attained by the paper.
2. Limitations
Question:
Does the paper discuss the limitations of the work performed by the authors?
Answer: [Yes]
Justification: Our paper contains separate 'Limitations' section dedicated to the limitations, such as the time-reversibility assumption that is standard in statistical phylogenetics.
Guidelines:
The answer NA means that the paper has no limitation while the answer No means that the paper has limitations, but those are not discussed in the paper.
The authors are encouraged to create a separate "Limitations" section their paper.
The paper should point out any strong assumptions and how robust the results are to violations of these assumptions (e.g., independence assumptions, noiseless settings, model well-specification asymptotic approximations only holding locally). The authors should reflect how these assumptions might be violated in practice and what the implications would be.
The authors should reflect on the scope of the claims made, e.g., if the approach was only tested few datasets with few runs. In general, empirical results often depend on implicit assumptions, which should be articulated.
The authors should reflect on the factors that influence the performance of the approach. For example, facial recognition algorithm may perform poorly when image resolution low or images are taken in low lighting. Or speech-to-text system might not be used reliably to provide closed captions for online lectures because it fails to handle technical jargon.
The authors should discuss the computational efficiency of the proposed algorithms and how they scale with dataset size.
If applicable, the authors should discuss possible limitations of their approach to address problems of privacy and faimess.
While the authors might fear that complete honesty about limitations might be used by reviewers as grounds for rejection, worse outcome might be that reviewers discover limitations that aren't 't acknowledged in the paper. The authors should use their best judgment and recognize that individual actions in favor of transparency play an impor- tant role in developing norms that preserve the integrity of the community. Reviewers will be specifically instructed not penalize honesty concerning limitations.
3. Theory Assumptions and Proofs
20

---

Question: For each theoretical result, does the paper provide the full set of assumptions and complete (and correct) proof?
Answer: [Yes]
Justification: However trivial, we performed analysis of the runtime of our method in the appendix and list the assumptions needed to guarantee near-linear runtime, such as balanced splits. Note that our work is of an applied nature so the practical runtimes and accuracies are more significant than the theoretical algorithm analysis itself, which is why the theoretical analysis is deferred to the appendix.
Guidelines:
The answer NA means that the paper does not include theoretical results. All the theorems, formulas, and proofs in the paper should be numbered and cross- referenced.
All assumptions should be clearly stated or referenced in the statement of any theorems.
The proofs can either appear in the main paper or the supplemental material, but they appear in the supplemental material, the authors are encouraged to provide 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.
4. Experimental Result Reproducibility
Does the paper fully disclose all the information needed to reproduce the main ex- perimental results the paper to the extent that affects the main claims and/or conclusions of the paper (regardless of whether the code and data are provided or not)?
Question:
Answer: [Yes]
We provide all methodological details for our FastCherries method and SiteRM model, including detailed time complexity analysis. For each experiment, details on how to reproduce them are provided in either the main text or appendix as appropriate. All our experimental results build benchmarks from prior works which are themselves well well-documented, such as the ProteinGym benchmark and the CherryML benchmarks. Guidelines:
Justification:
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, descríbing the architecture fully might suffice, or if the contribution is specific model and empirical evaluation, it may be necessary to either make 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 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 large language model), releasing of a model checkpoint, other means that appropriate to the research performed.
While NeurIPS does not require releasing code, the conference does require all submis- sions to provide some reasonable avenue for reproducibility, which may depend on the nature of the contribution. For example
(a) the contribution is primarily a new algorithm, the paper should make clear how to reproduce that algorithm.
(b) the contribution is primarily a new model architecture, the paper should describe the architecture clearly and fully.
(c) the contribution is a new model (c.g., a large language model), then there should either be way to access this model for reproducing the results or way to reproduce the model (c.g., with an open-source datasct or instructions for how to construct the dataset).
21

---

(d)
We recognize that reproducibility may be tricky in some cases, in which case authors are welcome 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., registered users), but it should be possible for other researchers to have some path to reproducing or verifying the results.
5. Open access to data and code
Question: Does the paper provide open access to the data and code, with sufficient instruc- tions to faithfully reproduce the main experimental results, described in supplemental material?
Answer: [Yes]
Justification: We have submitted all code to reproduce results alongside our submission. The provided folder contains README clear instructions to set up environment, install system requirements, and then reproduce all our results. Our work builds on the user-friendly benchmarks from the ProteinCiym and CherryML works, which we are very grateful for.
Guidelines:
The answer NA means that paper does not include experiments requiring code.
Please see the NeurlPS code and data submission guidelines (https://nips.cr ca/
httg//gddcccdsseeeee wwt.sumbbssiiooiiccoo..oo......i for more details.
While we encourage the release of code and data, we understand that this might not be possible, "No" an acceptable answer. Papers cannot be rejected simply for including code, unless this is central to the contribution (e.g., for new open-source benchmark).
The instructions should contain the exact command and environment needed to run to reproduce the results. See the NeurIPS code and data submission guidelines (https: /tipsscc/ppuleee wttpcccccccpcccccccccccccccccccc haabllccggiideseee /ooeeSSuu ssionPolicy) for more details.
The authors should provide instructions on data access and preparation, including how to access the raw data, preprocessed data, intermediate data, and generated data, etc.
The authors should provide scripts to reproduce all experimental results for the new proposed method and baselines. only subset of experiments are reproducible, they should state which ones are omitted from the script and why.
At submission time, to preserve anonymity, the authors should release anonymized versions (if applicable).
Providing as much information possible in supplemental material (appended to the paper) recommended, but including URLs to data and code permitted.
6. Experimental Setting/Details
Question: Does the paper specify all the training and test details (e.g., data splits, hyper- parameters, how they were chosen, type of optimizer, etc.) necessary to understand the results?
Answer: [Yes]
Justification; We built on top of the existing benchmarks from the ProteinGym and CherryML works. The choice of our hyperparameters À = 0.5 and f = 1.0 for variant effect prediction is motivated in the main text.
Guidelines:
The answer NA means that the paper does not include experiments.
The experimental setting should be presented the core the paper to 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, as supplemental material.
7. Experiment Statistical Significance
Does the paper report error bars suitably and correctly defined or other appropriate information about the statistical significance of the experiments?
Question:
Answer: [No]
22

---

Justification: Error bars are not reported because it would be too computationally expensive, particularly for the slower CherryML + FastTree method. Note that our experiments are all over large and diverse datasets, so although we do not provide error bars, statistical significance is conveyed through diversity of large benchmarks.
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, confi- dence intervals, or statistical significance tests, 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.
It is OK to report I-sigma error bars, but should state it. The authors should preferably report a 2-sigma error bar than state that they have 96% CI, if the hypothesis of Normality of errors is not verified.
For asymmetric distributions, the authors should be careful figures symmetric error bars that would yield results that are out of range (e.g. negative error rates).
show in tables or
If error bars are reported in tables or plots. The authors should explain in the text how they were calculated and reference the corresponding figures or tables in the text.
8. Experiments Compute Resources
For each experiment, does the paper provide sufficient information on the com- puter resources (type of compute workers, memory, time of execution) needed to reproduce the experiments?
Question:
Answer: [Yes]
Justification: Computational resources used for all experiments are provided in the appendix separate section.
Guidelines:
The answer NA means that the paper does not include experiments.
The paper should indicate the type of compute workers CPU or GPU, internal cluster, or cloud provider, including relevant memory and storage.
The paper should provide the amount compute required for each of the individual experimental runs as well as estimate the total compute.
The paper should disclose whether the full research project required more compute than the experiments reported the paper (e.g., preliminary or failed experiments that didn't 'ttntt make it into the paper).
9. Code Of Ethics
Question: Does the research conducted the paper conform, in every respect, with the NeurlPS Code of Ethics https://neuripg:: /t///nnuip..../ http///riss.cppullliieeepoo/
Answer: [Yes]
Justification: Our work does not involve human subjects nor participants; all datasets used are public and accessible, and all code used is open-source under suitable license. The primary application of our work variant effect prediction, which meant to aid doctors in the clinical setting.
Guidelines:
The answer NA means that the authors have 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.
23

---

The authors should make sure to preserve anonymity (e.g., if there is special consid- eration due to laws or regulations 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: Our work is applied to variant effect prediction, an important problem with positive societal impact. While our models could potentially be used in harmful way, for instance by computing and using polygenic risk scores against individuals, secure access to genetic data is beyond our domain.
Guidelines:
The answer NA means that there is no socictal impact of the work performed.
If the authors answer NA or No, they should explain why their work has no societal impact why the paper does not address societal impact.
Examples of negative societal impacts include potential malicious unintended uses (e.g., disinformation, generating fake profiles, surveillance), fairness considerations (c.g.. degloymemt 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 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, 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 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 attacks, mechanisms for monitoring misuse, mechanisms to monitor how system learns from feedback over time, improving the efficiency and accessibility of ML).
11. Safeguards
Question:
Does the paper describe safeguards that have been put in place for responsible release of data models that have high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)?
Answer: [NA]
Justification: Our work poses suck 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 encourage authors to take this into account and make best faith effort.
12. Licenses for existing assets
Question: Are the creators or original owners of assets (c.g., code, data, models), used in the paper. properly credited and are the license and terms of use explicitly mentioned and properly respected?
24

---

Answer: [Yes]
Justification: All benchmarks and models used are cited in our work, and all code and packages we used are open-source.
Guidelines:
The answer NA means that the paper does 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 used and, possible, include 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 (c.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 in the package should be provided. For popular datasets, paperswi thcode. com/ datasets has curated licenses for some datasets. Their licensing guide can help determine the license dataset.
For existing datasets that are re-packaged, both the original license and the license of the derived asset (if has changed) should be provided.
this information not available online, the authors are encouraged to reach out to the asset's s creators.
13. New Assets
Question: Are new assets introduced in the paper well documented and is the documentation provided alongside the assets?
Answer: [NA]
Justification: No new assets are released. Our code will be contributed the open-source ProteinGym and CherryML. packages upon publication.
Guidelines:
The answer NA means that the paper does not release new assets.
Researchers should communicate the details of the dataset/code/moddl as part of their submissions via structured templates. This includes details about training, license, limitations, etc.
The paper should discuss whether and how consent was obtained from people whose asset used.
At submission time, remember to anonymize your assets (if applicable). You can either create anonymized URL or include an anonymized zip file.
14. Crowdsourcing and Research with Human Subjects
For crowdsourcing experiments and research with human subjects, does the paper include full text of instructions given to participants and screenshots, if applicable, as well as details about compensation (if any)?
Question:
Answer: [NA]
Justification: Our work does not involve crowdsourcing nor research with human subjects.
Guidelines:
The answer NA means that the paper does not involve crowdsourcing nor research with human subjects.
Including this information the supplemental material is fine, but if the main contribu- tion of the paper involves human subjects, then much detail as possible should be included the main paper.
According to the NeurIPS Code of Ethics, workers involved in data collection, curation, other labor should be paid at least the minimum wage in the country of the data collector.
15. Institutional Review Board (IRB) Approvals or Equivalent for Research with Human Subjects
25

---

Question: Does the paper describe potential risks incurred by study participants, whether such risks were disclosed to the subjects, and whether Institutional Review Board (IRB) approvals (or an equivalent approval/review based on the requirements of your country or institution) were obtained?
Answer: [NA]
Justification:
The paper does not involve crowdsourcing nor research with human subjects.
Guidelines:
The answer NA means that the paper does not involve crowdsourcing research with human subjects.
Depending on the country in which research is conducted, IRB approval (or equivalent)
may be required for any human subjects research. If you obtained IRB approval, you should clearly state this in the paper.
We recognize that the procedures for this may vary significantly between institutions and locations, and expect authors adhere the NeurIPS Code of Ethics and the guidelines for their institution.
For initial submissions, do not include any information that would break anonymity (if applicable), such the institution conducting the review.
26


Reviews:
[{"review_id": "cQfjoI6MK7", "reviewer": "Reviewer_qr3c", "summary": "This paper introduces a new method for estimating amino acid substitution rate matrices from multiple sequence alignments, speeding up computation by orders of magnitude. The method, called SiteRM, outperforms traditional methods and large protein language models in variant effect prediction, showing its speed and accuracy in evolutionary biology.", "questions": "*  it would be better to provide a more vivid and more understandable method to explain the provided algorithm, such as drawing some diagrams or flowcharts of the pipeline.\n* although the provided method can improve the end-to-end runtime of the process, its performance drops a little comparing with the original method. Could you please list some possible reasons to explain the performance decrease, and some possible solutions.", "limitations": "* the performance of the method may vary across different protein families or evolutionary contexts. It would be better to assess the generalizability of the approach to diverse datasets and evolutionary scenarios.\n* the accuracy of the estimated rate matrices and predictions might depend on the quality of the input multiple sequence alignments. Are there any ways to address potential biases or errors in the MSAs to improve the robustness?\n* the benchmark used to evaluate the performance may have limitations or biases. It would be better to compare the performance on other benchmarks.", "rating": 7, "confidence": 4, "soundness": 3, "presentation": 2, "contribution": 2, "strengths": "* the proposed method to calculate amino acid substitution rate matrices is designed to be computationally efficient, with a near-linear runtime while maintaining comparable performance. This efficiency enables the analysis of extremely large MSAs, making it suitable for high-throughput applications.\n* SiteRM has shown superior performance in variant effect prediction compared to large protein language models that incorporate complex residue-residue interactions, which can be attributed to conceptual advances in the probabilistic treatment of evolutionary data.\n*  by estimating site-specific rate matrices for each column of a multiple sequence alignment, SiteRM captures the evolutionary dynamics at a finer resolution, which allows for a more accurate assessment of the impact of variants on protein function.\n* SiteRM can deal with large datasets with millions of sequences, which is particularly useful for handling the vast amount of data generated in clinical and deep mutational scanning studies, where comprehensive variant effect prediction is crucial.", "weaknesses": "* while the paper demonstrates the effectiveness of SiteRM in variant effect prediction, further exploration of its applicability to other evolutionary biology tasks or datasets could further understand its capabilities and limitations.\n* there lacks more detailed information on the benchmarking process, including datasets used, and potential biases in the evaluation, which could improve the transparency and reproducibility of the results.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for your detailed review. We respond below:\n\n> while the paper demonstrates the effectiveness of SiteRM in variant effect prediction, further exploration of its applicability to other evolutionary biology tasks or datasets could further understand its capabilities and limitations.\n\nWe agree and plan to explore other applications of SiteRM such as improving phylogenetic tree inference in future work. Software such as IQTree and its partition model puts this application within reach.\n\n> there lacks more detailed information on the benchmarking process, including datasets used, and potential biases in the evaluation, which could improve the transparency and reproducibility of the results.\n\nWe focussed on standard benchmarks from prior work (CherryML and ProteinGym) with the hope that it maximizes transparency and reproducibility. We also provided code with detailed instructions to reproduce all our results. Nonetheless, we agree that our work is not self-contained and more details on the benchmarks could be included. We plan to incorporate these more detailed descriptions of the benchmarks into the final version.\n\n> it would be better to provide a more vivid and more understandable method to explain the provided algorithm, such as drawing some diagrams or flowcharts of the pipeline.\n\nThis is a great suggestion. The closest we currently have to this is the runtime analysis in Appendix A.4.2. We will include a more graphical depiction of the algorithm in the final version.\n\n> although the provided method can improve the end-to-end runtime of the process, its performance drops a little comparing with the original method. Could you please list some possible reasons to explain the performance decrease, and some possible solutions.\n\nFastCherries shows a small asymptotic bias due to the use of Hamming Distance (HD) in the pairing step. We explored other alternatives during the project which we decided not to include in the paper, such as (1) pairing based on maximizing the composite likelihood of the pairing, (2) pairing based on minimizing the MLE distance between pairs, and (3) even a random pairer. The random pairer actually worked quite well for the WAG model of protein evolution on some benchmarks, showing asymptotic consistency with a relative statistical efficiency of ~1/8th, but worked poorly for the LG model due to the challenges of estimating site rates. We found that approaches (1) and (2) can indeed provide more accurate estimates in some cases. Since pairing based on HD worked very well already, we decided to make it the focus of the paper. We plan to explore other pairing methods in future work. It would be exciting to find a variant of FastCherries that retains the near-linear runtime while showing an even smaller error.\n\n> the performance of the method may vary across different protein families or evolutionary contexts. It would be better to assess the generalizability of the approach to diverse datasets and evolutionary scenarios.\n\nOur benchmark on the QMaker datasets (Supplementary Figure S1) shows that our method performs well on datasets from diverse parts of life. We agree that more in-depth error analysis may provide further insights and improvements to the method, but we leave this for future research.\n\n> the accuracy of the estimated rate matrices and predictions might depend on the quality of the input multiple sequence alignments. Are there any ways to address potential biases or errors in the MSAs to improve the robustness?\n\nThis is an important question which pertains to much work done in statistical phylogenetics. Although we did not explore it in our work, the speed of our method makes it possible (unlike prior work) to obtain bootstrap confidence intervals for the rate matrix estimation, which should enable users to understand the extent to which the estimates are stable to e.g. subsampling the MSA (whether rows or columns) or changing the MSA building algorithm.\n\n> the benchmark used to evaluate the performance may have limitations or biases. It would be better to compare the performance on other benchmarks.\n\nAs mentioned above, the QMaker benchmark provides evidence in this direction showing generalization to diverse domains of life.", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "BYMzljDShP", "reviewer": "Reviewer_cZPC", "summary": "The authors devise a fast and accurate phylogenetic inference algorithm based on CherryML. Their scalability allows them to fit flexible models to large protein families, prohibitively expensive for previous methods. As an example, they estimate site-specific substitution rates and use these rates to predict the effects of mutations. This method of mutation effect prediction is competitive with some state of the art generative modeling methods. This suggests scalable and accurate phylogenetic inference as a potentially important component of", "questions": "How does this work relate to the conclusions of Weinstein, Eli N., Alan N. Amin, Jonathan Frazer, and Debora S. Marks. 2022. “Non-Identifiability and the Blessings of Misspecification in Models of Molecular Fitness and Phylogeny.” Advances in Neural Information Processing Systems, December.?\n\nHow do you parameterize the substitution matrix?", "limitations": "Adressed.", "rating": 7, "confidence": 5, "soundness": 4, "presentation": 4, "contribution": 3, "strengths": "The writing is very clear. In particular, the authors describe CherryML very clearly.\n\nThe author's method is efficient and accurate (figure 1).", "weaknesses": "The authors state \"The first probabilistic model of protein evolution was proposed by Whelan and Goldman\". It would strengthen the paper to better cite the models that Whelan and Goldman were inspired by.\n\nI would appreciate the author comparing to more state-of-the-art mutation effect predictors.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for your thoughtful review.  Please see our response below:\n\n> The authors state \"The first probabilistic model of protein evolution was proposed by Whelan and Goldman\". It would strengthen the paper to better cite the models that Whelan and Goldman were inspired by.\n\nThank you. We will cite the earlier models in the final version, in particular Dayhoff's and JTT's.\n\n> I would appreciate the author comparing to more state-of-the-art mutation effect predictors.\n\nIn Supplementary Table S1, we show results including all variant effect predictors from the ProteinGym paper.\n\n> How does this work relate to the conclusions of Weinstein, Eli N., Alan N. Amin, Jonathan Frazer, and Debora S. Marks. 2022. “Non-Identifiability and the Blessings of Misspecification in Models of Molecular Fitness and Phylogeny.” Advances in Neural Information Processing Systems, December.?\n\nWhile we agree with the theoretical results presented in the above paper, our work challenges the assumption in their theorems that there exists a unique, \"global\" fitness function f. Instead, our intuition is that the fitness function is contextual, i.e. it is different in different parts of the evolutionary tree. For example, in some subtree, a positively charged amino acid may be required at position $i$, while in another subtree a negatively charged one may be required. This may be because of different contexts the protein is evolving in, due to e.g. protein-protein interactions. Our approach to variant effect prediction captures this intuition by using the \"local\" fitness function $P(y | x, t=1)$ at each point in the tree. We do not think that the theoretical results in the referenced paper necessarily explain our strong performance on variant effect prediction. Instead, we believe that it is the use of \"local\" fitness functions that is leading to our improved performance. There is certainly exciting theoretical and empirical work to be done in this direction.\n\n> How do you parameterize the substitution matrix?\n\nWe use the parameterization from the CherryML paper, namely the off-diagonal elements of $Q$ are given by $\\text{diag}(1/\\sqrt{\\pi}) \\cdot S \\cdot \\text{diag}(\\sqrt{\\pi})$ where $\\pi = \\text{softmax}(\\theta)$ for an unconstrained vector $\\theta$ (here $\\pi$ will be the stationary distribution of $Q$), and $S$ is a symmetric matrix given by $S = \\text{SoftPlus}(\\Theta + \\Theta^T)$, where $\\Theta$ is an unconstrained upper triangular matrix. This unconstrained parameterization allows CherryML to use unconstrained first-order optimizers to quickly estimate $Q$, as implemented by libraries such as PyTorch and Tensorflow.", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "2N5gdFzoTL", "reviewer": "Reviewer_f7u7", "summary": "The paper proposed a fast method for phylogenetic estimation from MSA called FastCherries, this method significantly speeds up the computational process with high accuracy. This method was demonstrated to be orders of magnitude faster than existing methods while achieving similar statistical efficiency. Further, the authors proposed the SiteRM model for variant effect prediction, also outperforming existing models, particularly large language models.", "questions": "See above.", "limitations": "Yes, the limitations are well described.", "rating": 6, "confidence": 3, "soundness": 4, "presentation": 3, "contribution": 3, "strengths": "1. The proposed methods are significantly faster than existing methods and are scalable for large sets of MSA. Meanwhile, the methods maintain a high level of accuracy, making a good contribution to the field of protein evolution.\n2. The methods made a superior performance in variant effect prediction than the large protein language model, which is effective in particularly GPU resources.\n3. The methods were tested rigorously on both simulated and real data.", "weaknesses": "1. The paper lacks a comparison to other methods in terms of computational resources, such as memory usage. Could the author provide some results?\n2. The model is still highly parameterized. Can the SiteRM model avoid overfitting with extensive datasets?\n3. How does the method perform with datasets that have a large number of gaps in MSAs (for example, the sequence has very low homology)?\n4. Can the methods be extended to incorporate other types of biological data, such as DNA/RNA sequence or structural alignment?", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for your thoughtful review. We respond below:\n\n> The paper lacks a comparison to other methods in terms of computational resources, such as memory usage. Could the author provide some results?\n\nOur method uses linear space; the logarithmic factors in the computational runtime come from divide-and-conquer and binary searches. This makes our method space efficient (up to a constant). We will discuss this and add practical results on memory usage in the final version. In more detail, the space complexity is $O(s^2 br + nl + ls).$  Here, $O(s^2 br)$ is the cost of the precomputed matrix exponentials, which is shared across all MSA’s. $O(nl)$ is the cost to load one MSA into memory. With a runtime of $n \\log n$, pairing takes $O(n)$ space, since it takes $O(n)$ to store the distance between x and all other sequences, and the space complexity equals the cost to store the stack of recursive calls. $O(ls + r)$ is the space of computing the initial site rates. The $O(ls)$ term comes from the count matrix $\\text{count}[i,j] =$   the number of times character $j$ occurs at position $i$. Computing site rates given branch lengths takes $O(r)$ space and computing branch lengths given site rates takes $O(l)$ space. In both cases, we only track the best rate for the current site or best branch length for the current site.\n\n> The model is still highly parameterized. Can the SiteRM model avoid overfitting with extensive datasets?\n\nAs described in Section 3.3 of our manuscript, we avoid overfitting by using pseudocounts drawn from the LG model, which serves as a prior. For the results described in the paper, we find that mixing the pseudocounts with the data at a 1:1 ratio works well in practice, but it should be possible to tune the regularization coefficient $\\lambda$ using cross-validation. Other options we did not explore that would benefit from large datasets would be some form of weight-sharing of the rate matrices across different positions.\n\n> How does the method perform with datasets that have a large number of gaps in MSAs (for example, the sequence has very low homology)?\n\nThe QMaker \"insect\" dataset is especially notorious for containing a large number of gaps. We observe a similar performance to FastTree in Supplementary Figure S1.\n\n> Can the methods be extended to incorporate other types of biological data, such as DNA/RNA sequence or structural alignment?\n\nIndeed, by using compound states one can model the joint evolution of amino acids and structural states. Our software is not limited to amino acid alphabets and can be applied to learn site-specific rate matrices for other kinds of molecular data such as DNA/RNA or structural states. We plan to pursue this in future work.", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "t9lKBxHHAr", "reviewer": "Reviewer_BDPQ", "summary": "The paper explores a method for estimating the transition matrix from multiple sequence alignments (MSAs). It utilizes a phylogenetic tree model, which is parameterized by the transition matrix and site rates. To estimate the transition matrix via maximum likelihood, an alternating optimization method is employed: first, the model parameters are fixed, and a set of phylogenetic trees is constructed from the MSAs. Then, given the MSAs, the trees, and the site rates, the transition matrix is estimated by maximizing the likelihood of the data given the tree.\n\nA recent method, CherryML, speeds up this process by replacing the likelihood of the data given the tree with the composite likelihood of the cherries (pairs of adjacent leaves) given the tree. This likelihood can be simplified using time-reversible models and can be learned from the tree, MSA, and leaf sequences.\n\nBuilding on the concepts of CherryML, this paper introduces a new method called FastCherry, which further simplifies the composite likelihood estimation by eliminating the need for tree construction. Instead, it focuses on creating disjoint pairs of similar sequences and estimating the composite likelihood from these pairs. Additionally, FastCherry extends the original WAG phylogenetic model by considering per-site transition matrices, with all parameters estimated using maximum likelihood.\n\nThe authors demonstrate the significance of their work through two different evaluation tests. The first test evaluates the speedup achieved by FastCherry compared to CherryML. The second test examines variant prediction performance in ProteinGym, comparing FastCherry to baseline approaches such as ESM-1v and ESM-IF1.", "questions": "1. How could the authors provide a deeper analysis of why the phylogenetic model performs better than language models in certain benchmarks?\n\n2. Could you please explain a bit on why on Figure 1.b when the number of family increases the error rate seem increase for CherryML?", "limitations": "NA", "rating": 5, "confidence": 4, "soundness": 2, "presentation": 3, "contribution": 2, "strengths": "The paper addresses an important application problem, though its relevance to the machine learning community may be limited. It is noteworthy that a phylogenetic model can achieve results comparable to, and sometimes slightly better than, pretrained PLMs in a specific use case. The paper is well-written, providing careful explanations of technical points that require a background in protein science.", "weaknesses": "The idea of simplifying the trees by partitioning the MSA into pairs of similar sequences is an approximate representation of the cherries in the tree. Since the partitioning algorithm is a greedy one, it may overlook the full structure of the MSA, reducing the relationships between sequences to binary relations between pairs. That also explains why CherryML with FastCherry lose the accuracy compared to CherryML with FastTree in Figure 1.b.\n\nThis work represents an incremental improvement over the original CherryML proposal, and therefore, its novelty is somewhat limited. Additionally, the results demonstrated in the variant prediction benchmark show only marginal improvements compared to language models. For a NeurIPS paper, I would expect more significant advancements, both in methodology and in the results presented.\n\nFor the machine learning community, it is not sufficient to demonstrate that one method is better than another through benchmark results alone. A deeper analysis explaining why the phylogenetic model performs better than language models in certain benchmarks and when it fails to do so would provide valuable insights for the community to learn from.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for your thoughtful review. Please find our response below:\n\n> The idea of simplifying the trees by partitioning the MSA into pairs of similar sequences is an approximate representation of the cherries in the tree. Since the partitioning algorithm is a greedy one, it may overlook the full structure of the MSA, reducing the relationships between sequences to binary relations between pairs. That also explains why CherryML with FastCherry lose the accuracy compared to CherryML with FastTree in Figure 1b.\n\nThe whole purpose of developing FastCherries is to avoid doing expensive computations while assuring little reduction in accuracy. It is a worthwhile tradeoff that can benefit applications requiring scalability, such as the one considered in our manuscript. Indeed, for the variant effect prediction task, it can be seen in Supplementary Figure S1 that FastTree is too slow to scale up to the full dataset size. In contrast, FastCherries provides no loss in performance and can crunch all the available data, yielding the best results.\n\n> This work represents an incremental improvement over the original CherryML proposal, and therefore, its novelty is somewhat limited.\n\nWe respectfully disagree with this view. As our results clearly demonstrate (Figure 1, Supplementary Figure S1, and Supplementary Table S1) and as the other reviewers have noted, the performance of CherryML with FastCherries is comparable to that of CherryML with FastTree while being one to two orders of magnitude faster.  Furthermore, we are proposing a novel method (SiteRM) to estimate site-specific rate matrices and this framework has the potential to advance phylogenetic inference significantly; software such as IQTree and its partition model puts this application well within reach.  When compared to the seminal LG model, SiteRM provides a substantial improvement in variant effect prediction, as seen in Supplementary Table S1, suggesting that it can model site-specific evolutionary constraints more accurately than previous approaches.\n\n> For the machine learning community, it is not sufficient to demonstrate that one method is better than another through benchmark results alone. A deeper analysis explaining why the phylogenetic model performs better than language models in certain benchmarks and when it fails to do so would provide valuable insights for the community to learn from.\n\n> How could the authors provide a deeper analysis of why the phylogenetic model performs better than language models in certain benchmarks?\n\nThe improved performance of our method comes from conditioning on the wildtype to obtain a local fitness function $P(y | x, t=1).$ Please see our response to Reviewer 2 for more discussion on the intuitions. We plan to make this more clear in the final version. We do agree with the reviewer that deeper error analysis of large language models would be a valuable research direction. Ultimately, we believe that the best variant effect predictors will combine our evolutionary modeling $P(y | x, t)$ with protein language modeling.\n\n> Could you please explain a bit on why on Figure 1.b when the number of family increases the error rate seem increase for CherryML?\n\nThe fact that the (rounded) error rate is 1.7 for $2^{11}$ and 1.8 for $2^{13}$ is just a small sampling effect. At these larger sample sizes, the statistical risk of the estimator has converged to its bias term.", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}]

Author rebuttal:
We thank all the reviewers for taking the time to carefully read our manuscript and provide thoughtful feedback. Different reviewers have asked different interesting questions, to which we have replied individually; we will also clarify them in the final version of our paper. The only major criticism is from reviewer #4, who states that our contributions are incremental over the original CherryML. As we explain in our rebuttal to the reviewer, we respectfully disagree with this view. In fact, the other reviewers have pointed out that our contributions are significant.  For clarity, we highlight our contributions below:

* Our new FastCherries algorithm significantly speeds up the tree estimation step of the CherryML framework, giving a near-linear time algorithm for end-to-end rate matrix estimation from MSAs. As we have pointed out in our response to reviewer #3 – whom we thank for bringing up space complexity – our method further has linear (and thus optimal up to constants) space complexity. While there is a small loss of accuracy, the whole purpose of developing FastCherries is to avoid doing expensive computations while assuring little reduction in accuracy. It is a worthwhile tradeoff that can benefit applications requiring scalability, such as the one considered in our manuscript.

* In applications, as our results clearly show (Figure 1 and Supplementary Figure S1 and Supplementary Table S1) and the other reviewers have highlighted, the performance of CherryML with FastCherries is comparable to that of CherryML with FastTree, while being one to two orders of magnitude faster. These results include diverse MSAs (Supplementary Figure S1 uses the QMaker MSAs which come from diverse areas of life). Furthermore, the performance of SiteRM over LG in variant effect prediction is far from incremental, as can be seen in Supplementary Table S1.

We agree that there is much interesting work to be followed up, such as applying our method to improve phylogenetic tree inference. We plan to pursue this research in future work.

We thank all reviewers again for their time and dedication.

Meta-review: All reviewers find this paper strong. Minor critical remarks/questions were made related to the significance of the improvement over prior SOTA, applicability in challenging cases (very large datasets, large number of gaps, ...), and theoretical analysis to complement empirical validation. Nearly all issues were resolved in the rebuttal and discussion.

Overall, the decision is unanimously to accept. The authors can hopefully make good use of the detailed reviews in subsequent work.
Decision: Accept (poster)