From single-sequences to evolutionary trajectories: protein language models capture the evolutionary potential of SARS-CoV-2

Domain: chemistry
Published: 2026-02-19
Editorial outcome: major_revision

Abstract
Protein language models (PLMs) capture features of protein three-dimensional structure from amino acid sequences alone, without requiring multiple sequence alignments (MSA). The concepts of grammar and semantics from natural language have been suggested to have the potential to capture functional properties of proteins. Here, we investigate how these representations enable assessment of variation due to mutation. Applied to the SARS-CoV-2 spike protein via in silico deep mutational scanning (DMS), the PLM ESM-2 captures evolutionary constraints directly from sequence context, recapitulating what normally requires MSA data. Unlike other state-of-the-art methods which require protein structures or multiple sequences for training, we show what can be accomplished using an unmodified pretrained PLM. Applied to SARS-CoV-2 variants across the pandemic, we demonstrate that ESM-2 representations encode the evolutionary history between variants, as well as the distinct nature of variants of concern upon their emergence, associated with shifts in receptor binding and antigenicity. ESM-2 likelihoods can also identify epistatic interactions among sites in the protein. Our results here affirm that PLMs like ESM-2 are broadly useful for variant-effect prediction, including unobserved changes, and can be applied to understand novel viral pathogens with the potential to be applied to any protein sequence, pathogen or otherwise.

Introduction
The conventional study of protein sequence variation requires the alignment of amino acid sequence data. Alignments help identify where variation accumulates within sequences, revealing the evolutionary constraints that dictate observed variation. In human genome studies, impactful missense mutations, non-synonymous substitutions resulting in amino acid changes, are primarily linked to disease, while in viruses, amino acid changes are more often linked to adaptive evolution1,2. Alternatively, these changes (and in fact the majority of possible missense mutations) are deleterious resulting in, e.g., misfolding and a decrease in fitness or are “neutral” and have little to no consequence.

To date approaches for assessing the impact of mutations are premised on the availability of comparative genome sequence data to assess features like the relative proportions of non-synonymous to synonymous substitutions (dN/dS) or evolutionary conservation at individual sites. Similarly, genome sequence-based surveillance methods usually require aligned real-time sequence data to assess the relative growth rate of a pathogen lineage3,4. Experiments on novel variants can be initiated when first detected, but they take considerable time to complete.

Natural language processing methods can be applied to biological sequences like DNA and proteins5–7. Biological sequences are physical molecules represented by arbitrary characters. This is convenient since these characters, representing nucleotides or amino acids, can be processed in the same way as words in natural language. Hie et al.6 show how, much like natural language, biological sequence context can be used to create informative embeddings that reflect complex properties of the sequence. Specifically, concepts like grammar (the structural rules of a language) and semantics (the meaning of the language) can purportedly be applied to changes in protein amino acid sequences to capture ‘fitness’ and mutation effects, respectively6. The ‘grammaticality’ of a sequence is captured by its probability, indicating how well it conforms to the rules of protein formation and composition. The semantic difference between protein sequences can be measured by computing distances, a ‘semantic score’, between sequences within the embedding space.

Evolutionary scale model 2 (ESM-2)5,7 is a protein language model (PLM) trained on approximately 65 million unique protein sequences5,7. Here, we demonstrate that ESM-2 can infer fundamental properties that classically require a multiple sequence alignment. Applied to SARS-CoV-2 we show that the grammaticality and semantic score metrics reveal characteristic properties of the SARS-CoV-2 spike protein sequence (Fig. 1). We demonstrate these metrics have a broader meaning when computed using ESM-2 compared with previous models, likely due to differences in training data6, before assessing how to use these metrics for variant and sequence horizon scanning. Crucially, these metrics, produced from a single protein sequence using pre-trained models, bypass the need for alignment data. We make a case for their future application in characterising and assessing the constraints on proteins from emerging pathogens, before assessing how to use these metrics for variant and sequence horizon scanning.

Results
ESM-2 captures evolutionary potential
We first assess how well ESM-2 captures constraint on sequence variation from protein structure. Specifically, does the model assign appropriate values for relative grammaticality (the log-likelihood ratio between the mutation and the reference amino acid) and semantic score (the distance between sequence embeddings, e.g., mutant and wild-type) at different sites in the protein sequence. To do this analysis, we utilised an in silico deep mutational scanning (DMS) approach8 using embeddings for every possible single amino acid replacement in the SARS-CoV-2 spike protein. In vitro deep mutational scans involve experimentally substituting each site in a protein with every other amino acid and measuring the substitutions’ effect on phenotype. We do this computationally by embedding each of these mutated sequences using ESM-2 and then calculating grammaticality and semantic scores for each sequence.

The ESM-2 scores correlate well with the structural subunits of the spike (Fig. 2). Using our DMS approach to produce relative grammaticality and semantic scores, we observe that the protein has two main regions that broadly correspond to the two protein subunits S1 and S2 (Fig. 2B, C). Relative grammaticalities decrease when substitutions are introduced in the S2 indicating that the region is less tolerant to change compared to S1. We observe a statistically significant difference between the mean relative grammaticalities of the S1 (n=685 positions) (Supplementary Fig. 1) and S2 (n=588 positions) subunits (p=2.60e-160, two sided Mann–Whitney U test). Spike is a trimer with the core of the S2 subunit being formed by three closely interacting alpha helices from each monomer. As such, changes in this region could disrupt the formation of a stable and functional spike trimer by disrupting the inter-monomer interactions. On the other hand, in the S1 subunit, mutations in the N-terminal Domain (NTD) and receptor-binding domain (RBD) have higher relative grammaticalities. Despite being important to protein function, these regions need to be sufficiently flexible to facilitate receptor binding interactions and accommodate mutations to evade host immunity. This mirrors the accumulation of sequence variation in alignments. Specifically, using entropy, a measure of site-specific variation computed from multiple sequence alignments, there is significantly higher entropy observed at sites in S1 compared to S2 (Fig. 2A), with each region found to be significantly different (p=1.121e-10, two sided Mann–Whitney U test) (Supplementary Fig. 1). Thus, ESM-2 has correctly identified S1 as a more variable region than S2 (an analysis based on permutations of a single sequence) and indicates that ESM-2 has learnt this property about spike. Available full spike experimental DMS results from Dadonaite et al.9 show binding is the only measure with a significant difference between subunits (p=1.70e-06, two sided Mann–Whitney U test), which might be expected given that S1 contains the spike binding domain.

The semantic rank also increases in the NTD region of the S1 (Fig. 2A, B), indicating that changes here may produce large structural changes and are more likely to be accommodated due to higher relative grammaticalities across the region; presumably reflecting properties of the NTD, i.e., it is a region on the surface of the protein and contains many antigenic epitopes10. These results are consistent with prior results that show PLMs like ESM-2 can be used as an effective predictor of variant effect11. If PLMs can model regional constraints, they may also be able to identify epistatic interactions that are linked with protein structural constraints.

Our results show ESM-2 models spike protein constraints, mapping low relative grammaticalities to the S2 and high relative grammaticalities to the S1, which is heavily targeted by host antibodies (NTD, RBD)10,12–15 and predominantly on the protein surface (Fig. 2). Hie et al.6 use the semantic score as a proxy for antigenic change, however, this is not necessarily relevant for protein sites that are not involved in antibody binding. Since ESM-2 is trained on a wide diversity of proteins, semantic scores produced by this model represent shifts in structure/function rather than specifically antigenic change. Many high semantically ranked changes occur in the NTD (Fig. 2B) which contains a “supersite” with several epitopes targeted by neutralising antibodies10. The NTD slow relative grammaticality and high semantic score suggests it allows large structural changes to alter its antigenic properties without greatly affecting function.

ESM-2 identifies epistatic interaction sites
A consequence of protein structural constraint is that compensatory mutations are often required before a change can be accommodated. Such epistatic effects can be due to direct or indirect mechanisms linked to residues in close contact in the protein structure conformation and/or protein stability16,17. Here, we present an approach for assessing epistasis, where we find putative epistatic interactions among sites by identifying where changes in the amino acid at one site cause significant changes in amino acid likelihoods at another. Language models compute likelihoods for each amino acid in a sequence based on the context of the rest of the sequence. This means that the likelihoods depend on the other amino acids in the sequence and will change if other positions mutate. Positions that are unaffected by the mutated site would be expected to change minimally, whereas those affected may change their probabilities, i.e., lower or higher depending on whether an epistatic interaction is antagonistic or synergistic, respectively. A method by Zhang et al.18 uses a similar approach to predict contact maps, however, here we want to predict mutation specific effects instead of site specific contacting residues.

The BA.1 Omicron, variant of concern (VOC) is defined by over 30 unique substitutions in its spike protein sequence. By reverting these defining mutations to the Wuhan-Hu-1 reference sequence amino acid, we measure the effect this reversion has on the probabilities of all other BA.1 amino acids computed by the model. To illustrate these results we highlight the interactions of three key substitutions E484A, D614G and N969K19–22; selected due to their effects on antigenicity, conformational changes and presence in the S2 domain. For ease of interpretation, we inversed these values so that positive scores can be interpreted as a synergistic and negative scores antagonistic when the mutation is acquired. To illustrate these results we highlight the key substitutions E484A, D614G and N969K19–22; selected due to their effects on antigenicity, conformational changes and presence in the S2 domain.

The E484A substitution in the RBD of spike is involved in modulating binding to the ACE-2 receptor and is a good example of a substitution with localised effect (Fig. 3A). The VOCs Beta and Gamma contained the E484K substitution and this contributed to escape from several RBD-neutralising antibodies23. Interestingly, A, D, G and K substitutions at position 484 all confer resistance to human convalescent sera14,24. The mutation primarily changes amino acid likelihoods within the RBD (Fig. 3A).

Interestingly, RBD positions 482, 485, 486, 488 and 500 have the largest absolute changes in probability following the E484A mutation’s occurrence: 485 and 482 decrease their likelihood of changing indicating incompatibility with the E484A. Positions 486, 488, 500 increase their likelihoods of changing following the E484A mutation, indicating that this mutation is synergistic with these positions, possibly improving the structure or at least minimising disruption. Omicron mutation Q498R also increases its likelihood in response to E484A. E484A shows decreased immune evasiveness relative to E484K, however interactions between N501Y and Q498R mutations recover this20. While 501Y is not identified, position 500 increases in likelihood suggesting synergistic interactions with this region.

Our results show S1 mutations like E484A have shorter distances to their affected positions in sequence space and in three-dimensional space (Supplementary Fig. 3) likely due to it being the subunit comprising most of the accessible surfaces of the protein. Changes more ‘central’ to the protein structure like those in the S2 domain have a greater chance to affect more amino acids due to being surrounded by other sites. This can produce more impactful changes that have knock-on effects at distant locations.

The D614G substitution emerged early in the pandemic25, is now present in all circulating lineages, and is an example of a substitution with broad effect (Fig. 3A). Specifically, D614G has two main functional consequences: (1) to remove a hydrogen bond between 614 and 859 on the adjacent spike monomer, and (2) at 647 to contribute to the structure of the C-terminal domain. 614G increases the stability of the spike protein19 as well as the infectivity of the virus at the cost of increased susceptibility to neutralising antibodies19,26. Its reversion affected 46 other positions (Supplementary Fig. 4), nine of which are in the RBD, with a further 30 in S1 and seven in S2. While 859 and 647 were not affected positions, 615, 617 and 649 did change their likelihoods. 649 is an internal amino acid located behind 614 while 615 is in direct contact with the site and 617 is slightly downstream. 617 and 615 both increase their likelihoods, indicating that the D614G acquisition was synergistic with those sites. 649 decrease its likelihood by 3.7% on reversion suggesting an antagonistic effect with 614G.

The N969K reversion affected 46 positions on the protein (Fig. 3, Supplementary Figs. 3 and 4) and has been linked to the change in Omicron entry route21. N969K forms electrostatic contacts with Q755 on the adjacent protomer in the pre-fusion state27, as well as interacting with and displacing the HR2 backbone in the post-fusion spike structure22. It has long-range interactions both in linear sequence and in the 3D protein structure context (Supplementary Fig. 3A) and also has the greatest number of affected sites compared with the other BA.1 mutations. Position 976 had a likelihood increase of 21% while position 964 had a 15.9% increase, indicating a strong synergistic effect with the N969K substitution. N969K is positioned in the middle of a loop connecting two alpha helices, with 976 and 964 positioned at the end of the first helix and the start of the second. The large likelihood shift indicates strong synergy with the N969K mutation and could be as yet undiscovered epistatic positions. Both D614G and N969K have interacting residues on other monomers, however ESM was only given a single monomer sequence. The lack of other monomers in the embedding may prevent the detection of these interactions.

To further validate that ESM-2 was detecting epistatic effects with other sites, we compared the RBD positions identified by our BA.1 reversions with experimental data from Moulana et al.28, which characterised the epistatic interactions involved in BA.1’s increased ACE-2 binding affinity. We used two cutoffs, a relaxed threshold using a single standard deviation and a stringent two standard deviation threshold. We found that the identified RBD mutational interactions appear to be enriched with high ranked epistatic changes (Fig. 3B–D) in both thresholds. Using a hypergeometric test, the stringent threshold was found to be 6.43 fold over enriched (n=267, p-value 1.46e-10) with RBD-RBD BA.1 mutation interactions. This was also true for strong RBD-RBD mutation epistatic interactions (n=267, an over enrichment of 8.21 fold, p-value 0.00562) as well as for all BA.1 mutation-mutation interactions (n=653, 1.6 fold, p-value 0.01687). Strong epistatic interactions were defined as 1 standard deviation from the mean absolute experimental epistatic interaction. These remained significant after applying a Benjamini-Hochberg multiple test correction.

The mutation sets with the largest pairwise effects 496–501 and 501–505 are identified in the most stringent cutoff, with the next biggest impacting pairs 496–498 and 498–501 identified in the less stringent set (Fig. 3B–D). While all mutational pairs measure some difference, the ESM-identified changes appear to be predominantly impactful changes. Interestingly, the strongest epistatic effects that ESM-2 identifies between positions 371, 373 and 375 (Fig. 3B) are only found to be weakly involved in ACE-2 binding epistatic effects (Fig. 3C, D). However, recent experimental data from Furnon et al.29 has shown that these sites as well as 339 appear to be critical determinants of BA.1’s enhanced replication in nasal epithelial cells, a major phenotypic feature of Omicron SARS-CoV-2 viruses. These mutations are also highlighted by Moulana et al.28 due to their significant interactions with each other, despite their smaller effect on ACE-2 binding. As our measure of epistasis is not directly linked to ACE-2 binding and is simply a change in likelihood of the amino acid at the site, it is probable that ESM-2 is picking up on multiple epistatic effects rather than just binding affinity. The lack of direct relatedness to binding affinity may also explain why predicted effects in our analysis (Fig. 3B) do not always match the direction (positive or negative) to the results reported by Moulana et al.28. Several factors are at play in determining the suitability of mutations within the RBD including immune escape potential as well as binding affinity. Encouragingly, we also identify a larger set of the Moulana et al. epistatic interactions compared to an epistatic model30 utilising the pandemic sequence data to identify interactions (Supplementary Fig. 12).

We also tested the sensitivity and specificity of the model in comparison to the mutual information based approach from Innocenti et al.30 and found that when looking within the experimental data for epistatic sites, ESM-2 (stringent 2 standard deviation cut-off) had an equivalent sensitivity (0.120 vs 0.114) and improved specificity (0.842 vs 0.777). This was calculated using the strongly epistatic set of RBD-RBD mutation interactions. While the specificity is good for both approaches, suggesting both methods work well for ruling out false positive epistatic interactions (or in this instance, weak interactions), the low sensitivity suggests that the task of identifying true/strong interactions remains challenging. Considering ESM-2s lack of SARS-CoV-2 pandemic sequence data in training, future work could be directed towards leveraging pandemic data to improve the predictive performance, although the lack of extensive validation/training experimental pairwise datasets may prove to be a bottleneck here. Also, it should be again noted that ESM-2 is not trained explicitly to predict effects on ACE-2 binding, meaning the epistatic signal may be from another biologically meaningful measure, like the 371, 373 and 375 mutations.

In our epistasis analysis, position 330 changes its likelihood in the mutations mentioned above (E484A, D614G and N969K) and is also associated with several BA.1 mutations. This prompted an investigation into whether certain sites were more prone to fluctuations than others, and whether these sites were of functional importance. We identified positions that significantly changed their likelihood in each DMS mutant (putative ‘critical’ sites), before filtering for sites that occur across many different mutations and investigating their functional relevance to the protein. We also investigated which sites when mutated significantly affect a larger number of other sites (putative ‘impactful’ sites).

The critical sites in the NTD and RBD are predominantly prolines and cysteines, while elsewhere there is a larger range of amino acid types, seven of which contain aromatic side chains (F, Y and W) (Fig. 4A, B). These amino acids have important structural features. Proline is a rigid amino acid due to its side-chain which reduces its flexibility31,32 and prevents it from forming stable α-helices31. Prolines are often found where sharp turns are necessary for structure. Cysteine thiol side-chain allows them to form covalent di-sulphide bridges between adjacent cysteines33. Cysteines are well conserved throughout proteins33, and all of the identified cysteines were involved in di-sulphide bridges (Fig. 4A, B).

The impactful sites are also heavily enriched with cysteines throughout the protein, affirming their critical importance to protein structure (Fig. 4C, D). Two positions within the furin cleavage site (FCS) 682 and 685 are also highlighted. Both are arginine, a necessary polybasic component of the cleavage site that allows cleavage by furin34. The FCS is an important feature of SARS-CoV-2 that enhances transmission35. Critical and impactful sites play key structural and functional roles in the protein, with clear identification of structurally important amino acids such as cysteines and functionally relevant positions within the FCS.

ESM-2 provides orthogonal metrics for assessing mutational impact
An important question is whether the grammaticality and semantic score represent a different and useful measure of change compared with existing experimental, computational or evolutionary measures of mutation impact. We also wanted to determine whether these scores equate to fitness and antigenicity6 when computed using a general PLM like ESM-2. Extensive high-throughput experimental investigations have revealed various aspects of spike protein biology. We sought to correlate these measurements with PLM metrics, to better understand what characteristics are captured by the model. We calculated spearman’s rank correlations between the embedding metrics and experimental and computational metrics. We then used the full embeddings and logits as well as score combinations to fit linear and support vector regressions (SVR) to each metric and again scored using a spearman’s rank correlation.

We benchmarked our results against two other computational scores: EVEscape from Thadani et al.36 and mutability scores from Rodriguez-Rivas et al.37. EVEscape is a variational autoencoder based on EVE38 that was trained using pre-pandemic spike protein sequences. These prior sequences are given to EVEscape as a one-hot encoded alignment, which the model uses to learn distributions and dependencies across sites in the protein. The variational autoencoder produces a mutation score that when compared with the score of the reference sequence approximates the relative fitness of the mutant, termed EVEscape Fitness. This is combined with a structure derived contact accessibility and amino acid dissimilarity metric to produce the full EVEscape score. These additional metrics aim to weight the model towards detecting more immune-evasive mutations36.

There are two mutability score variants37, the direct coupling analysis (DCA) score and the Independent Site Model (IND) score. Both approaches use existing protein sequence alignments to calculate site-wise frequencies across the protein sequence. Much like the EVEscape fitness, the IND score is a likelihood ratio between the reference and mutated sequence. The DCA additionally incorporates coupling analysis to try to incorporate signals of mutational epistasis. This involves adding coupling terms between sites in order to approximate the epistatic interactions between them. We also tested other methods for computing grammaticality: masked grammaticality suggested by Allman et al.39; mutated grammaticality; and the relative grammaticality measure supplemented with the dissimilarity and accessibility as shown in Thadani et al.36.

Experimental scores were used from three in vitro DMS studies, two of which were performed on the RBD while the other was performed on a full spike protein. “Escape”, “Entry” and “Binding” were determined from a full spike DMS by Dadonaite et al.9 and correspond to: immune escape from human sera, cell entry, and binding affinity to soluble ACE-2, respectively. “RBD Escape” was produced by Yisimayi et al.40 using an RBD only DMS. The “Wuhan” and “Variant” binding and expression scores are from an RBD only DMS from Starr et al.41 where binding to ACE-2 and expression of the protein were measured across five spikes. The Variant score is made from an average of both measurements across all five spikes.

Broadly, all metrics weakly correlate with experimental datasets (Fig. 5A). Semantic scores and grammaticality metrics are typically negatively correlated with each other. The semantic score exhibits a higher correlation on both “RBD Binding” datasets, while also achieving stronger correlations compared to grammaticality metrics on “RBD Expression” data. Grammaticalities tend to have higher correlations with Escape and full spike Entry and Binding data. However, the grammaticalities are typically only competitive with EVEscape and mutability scores when the accessibility and dissimilarity components of EVEscape were added to the relative grammaticality.

Regression correlation coefficients were comparable to the direct score correlations for most metrics (Fig. 5B). However, many were found to be insignificant upon multiple test correction. Full spike Entry and Binding were exceptions, with most predictions significantly yet weakly correlating. Logits and embeddings do significantly improve upon all single score metrics and were found to be significant in every experimental dataset. Correlation coefficients are also much higher, with all but full spike Escape and Binding achieving coefficients greater than 0.5. Embeddings and logits are dense representations, and much more descriptive than any single numeric metric. Non-linear regressions like SVR appear to extract the most from these representations, with linear regression typically achieving weaker correlations likely due to their inability to exploit non-linear relationships in these high-dimensional representations (Supplementary Fig. 11).

We also calculated several computational metrics that are commonly used for the analysis of protein sequences and estimation of evolutionary constraints. These provide quantitative measures of substitution probability and structural characteristics at every site in the protein and are intended to be representative of the wide range of techniques currently used to analyse protein structural evolution. This allowed for comparison of the PLM metrics with existing metrics. The crystal structure of the spike protein was used to calculate “Accessibility”, “B-Factor”, “ΔΔG” and the environment-specific substitution table (ESST). Accessibility measures the antibody accessibility at every position. The B-Factor is the temperature factor derived from the protein crystallography experiment used to determine the 6VXX structure42 and is a measure of the local fit of the structure to experimental data; it is often increased if there is static or dynamic disorder. Protein stability change (ΔΔG) was assessed with the FoldX software. The ESST were taken from Mizuguchi et al.43 and provides a statistical estimate of substitution likelihood based on the observed frequency of substitutions at amino acids in similar local amino acid structural environments. We also calculated two sequence-based measures of substitution likelihoods. The position-specific scoring matrix (PSSM) was calculated using DeltaBlast and calculates the log likelihood of substitutions occurring based on their frequency in related sequences42. “Entropy” was used to provide an estimate of the variability present at a site.

Grammaticalities correlate better with computational metrics (Fig. 5A), particularly B-Factor and ESST. This indicates that grammaticalities (which align well with conservation) may also be a reasonable proxy for assessing protein flexibility. EVEscape and relative grammaticality combined with EVE components performed well on accessibility, but this is to be expected given that both use antibody accessibility as a component of the score. The semantic score performed better with PSSM and ESST but struggled to achieve high correlations with other metrics. The relative grammaticality with EVE components was often competitive with EVEscape on most metrics, although typically EVEscape achieved slightly greater correlations more often.

Regression improve correlations for most computational metrics (Fig. 5B), particularly for ESM-derived scores. Despite their accessibility components, EVEscape and relative grammaticality with EVE components both achieve worse correlations than the grammaticalities or the semantic score with SVR regression. Linear regression results are similar (Supplementary Fig. 11), and both EVEscape and mutability scores are found to be not significant. Grammaticality predictions are less correlated with B-Factor than the grammaticality scores however, are all found to be significant (Fig. 5B). EVEscape remains the most predictive score for ΔΔG but is only significantly predictive for one other computational metric ESST both with linear and SVR regressions. Embeddings and logits exhibit higher correlations with computational metrics when fitted with regression models, with all correlations above 0.65 except for Entropy, which remains challenging to predict. Larger representations like logits and embeddings therefore remain the most predictive but require fitting to data unlike single scores.

We also calculated a statistical distance metric, the Jensen Shannon distance (JSD), between these computational metrics and the semantic score and grammaticality at every site. Assessing the mean JSDs across different regions of the spike protein reveals a statistically significant difference between regions in S1 and S2 for the grammaticality score versus other metrics (Supplementary Fig. 6). This supports the conclusion that grammaticality carries implicit context-aware information about differences in evolutionary constraint between protein regions, as indicated by 2B; similar information not captured by existing, site-independent computational metrics. JSD between the semantic score and the other metrics does not show pronounced differences between S1 regions and S2 as consistently, although there were still statistically significant differences between some of the protein regions, again demonstrating contextual structural awareness in the PLM embeddings.

We observe that the ESM-derived scores correlate weakly with the experimentally tested or computationally derived metrics (Fig. 5). However, when testing which mutations appeared in real data, i.e., SARS-CoV-2 pandemic sequences, we see that computational metrics are often equivalent to or better predictors than experimental measures, despite being significantly quicker to produce (Supplementary Fig. 10). To compare approaches equally, the DMS data were filtered to include measurements only available in all datasets. For comparison purposes the top decile of scores (as previously used by Thadani et al.36) were then selected for each method, and the number of mutations observed in the in VOC sequences (Supplementary Fig. 10A) and in the population (Supplementary Fig. 10B) were assessed.

The relative grammaticality with EVE annotations, EVEscape, and EVEscape Fitness were found to identify the most mutations out of the computational approaches (Supplementary Fig. 10A). ESM metrics identify many fewer of the VOC mutations without the additional accessibility and dissimilarity components used by EVEscape. This would suggest that most of the VOC mutations were driven by their accessibility to antibodies. We know that for several of the VOC defining mutations this is the case, with sites like 484 and 501 located within epitopes bound by multiple classes of antibody24. Both RBD binding measurements identify the greatest number of mutations for almost all measured variants with the top 10% of scores (Supplementary Fig. 10A). Full spike Binding and Escape were much less predictive and identified very few mutations in the RBD or in the rest of the spike. Full spike Entry identifies almost all of the mutations; however, this is mostly within the top 25% of scores. Given that experimental DMS was produced later in the pandemic from when a viable set of computational scores could have been created, computational methods remain competitive at identifying VOC-defining mutations.

Considering the broader set of mutations (Supplementary Fig. 10B), we see that ESM metrics outperform other DMS methods when identifying mutations that occurred at higher frequencies (at least 100 times) in the pandemic data. The masked and mutated grammaticalities identify the greatest proportion of mutations in their top 10%. The masked score outperforms all experimental measures in the RBD and in the rest of spike. RBD experimental measures again identify more mutations than full spike DMS results. EVEscape fitness remains competitive with ESM-metrics, however EVEscape performs significantly worse. The relative grammaticality with EVE annotations also performs much worse here, although remains competitive with EVEscape. Mutations in this set don’t all persist or become fixed unlike the VOC mutations. This indicates that many are likely to be close to neutral. This could explain why EVEscape and the relative grammaticality with EVE annotations perform poorly here, since these scores are weighted towards antibody accessible changes. Escape mutations often come at a cost, e.g., with many reducing binding affinity44. VOCs like Omicron include several mutations that individually seem deleterious but together efficiently evade host immunity45.

While experimental results are considered the benchmark for comparisons, these computational metrics appear to be as predictive at identifying future mutations as the experimental measures while being significantly cheaper and quicker to produce. Experimental measures are more interpretable than ESM-derived scores since they assess a single property of a mutation. Sequence based language models derive their predictive signal from millions of protein sequences, spread across the tree of life. Hie et al.6 show that their scores are associated with viral fitness and antigenicity, however, the alignments used to train their model contain changes that are heavily associated with immune evasion. This mirrors SARS-CoV-2’s evolution during the pandemic, where most spike protein changes were the result of intense selective pressure to evade the host immune response and maintain viral fitness1,24,44,45. ESM-2 is trained on proteins that have a variety of selective pressures that are distinct from those faced by viral glycoproteins. This would suggest that an increase in the grammaticality or the semantic score is not simply a change in binding, escape, or any other single other protein property. Rather, they are a composite of these features, where associations to fitness or antigenicity can be drawn thanks to the over or under representations of similar sequences in the training data. Alignment entropy is predictive in much the same way, where the identification of changeable sites can help develop an understanding of protein properties, despite entropy not necessarily correlating with experimental results.

Languagemodels recapitulate evolutionary relationships between related sequences
Given that PLM approaches appear useful for interpreting single mutations within a protein sequence, we tested how ESM-2 would interpret differences between real SARS-CoV-2 sequences with multiple mutations. The earliest known SARS-CoV-2 sequence for each Pango lineage was extracted from the global data (retrieved from GISAID) and their spike protein sequences were embedded using ESM-25,7. The evo-velocity package was then used to produce an evo-velocity UMAP for the sequence embeddings. Evo-velocity assigns a putative directionality between the embeddings that describes the flow of evolution through the UMAP embedding space. Firstly, evo-velocity46, confirming these authors results, shows that PLM embeddings can distinguish between meaningfully different spike proteins (Fig. 6A, B). Secondly, it illustrates that embeddings can help to understand the evolutionary landscape of spike and the directionality of its evolution.

The evo-velocity accurately describes the evolution of SARS-CoV-2 moving from the earlier ‘non-VOC’ clusters into VOC clusters, consistent with their evolution during the pandemic (Fig. 6A, B). Omicron and Delta form distinct clusters while Gamma forms a concentrated cluster on the fringes of the non-VOC sequences. Beta and Alpha are less homogenous, and recombinants tend to fall with their parental lineages (primarily Omicrons). Sarbecoviruses (which include bat and pangolin viruses related to SARS-CoV-2) form another distinct cluster close to the early SARS-CoV-2 sequences. These are also identified as the earliest sequences by the evo-velocity pseudotime algorithm (Fig. 6B and Supplementary Fig. 5). Evo-velocity first uses diffusion analysis to identify root and endpoint sequences before estimating the order of evolution using a pseudotime simulation (Fig. 6B and Supplementary Fig. 5). Pseudotime analysis achieved a significant (p=3.24e-296) spearman’s rank correlation of 0.86 against sampling time, confirming that the model has inferred the evolution of the sequences in the correct order. Evo-velocity recapitulates the topology of the representative sequence phylogeny in Fig. 6C, with early VOCs more closely related to non-VOC sequences than Omicrons and recombinants. This shows PLM embeddings capture meaningful representations that differentiate between distinct sequences. The congruence between the phylogenetic tree topology and the evo-velocity derived structure provides further confirmation that this method captures the evolutionary history of the spike protein.

PLM metrics also differentiate between distinct SARS-CoV-2 spike proteins. Using the semantic score or the relative grammaticality, a non-VOC, early VOC (Alpha, Beta, Gamma and Delta), and an Omicron cluster can be identified (Fig. 6, Supplementary Figs. 7 and 8). Unlike the UMAP (Fig. 6A and Supplementary Fig. 5) which is a dimensionality reduction and projection from high-dimensional space, grammaticality and the semantic score are much simpler to produce. Despite this, they still recapitulate much of the information shown in the UMAP. Evo-velocity, thus, captures differences between VOC sequences, recapitulates the topology of the nucleotide sequence phylogeny and reconstructs the direction of evolution accurately (Fig. 6). The PLM scores can effectively group lineages into their VOC categories (Fig. 6D, Supplementary Figs. 7 and 8).

Detecting the distinct nature of the variants of concern on emergence
The emergence of SARS-CoV-2 and its subsequent variants of concern (VOCs) caused large waves of infections during the COVID-19 pandemic. Predicting their fitness advantage just from their initial sequences, before the viruses were in widespread circulation, has proven difficult. PLMs could potentially be used to meet this challenge of rapid characterisation of individual pathogen genomes as we have implemented for relative grammaticality and semantic scores (see Observable notebook: https://observablehq.com/@cvr-bioinfo/from-a-single-sequence-nature-communications).

The PLM metrics produced for each Pango lineage can be used to assess SARS-CoV-2 variants since they detect differences between sequences. We observed a characteristic “jump” in both semantic score and relative grammaticality between non-VOC, early VOCs and Omicron sequence clusters (Fig. 6D and Supplementary Fig. 8). Earlier VOCs required fewer changes to be successful since most of the human population was still naive to infection or vaccine-derived immunity. The Omicron lineage emerged during high levels of both vaccination and infection. It contained many more substitutions in the spike protein compared to previous variants, it changed its entry mechanism preference and managed to largely evade previous immunity which resulted in an extended vaccination regimen of three doses being recommended47,48. Unlike the other VOC groups, Omicron sequences initially decreased their semantic scores. The BA.2 variant had fewer mutations than the initial BA.1 variant, however later sequences that increased their number of mutations, continued to decrease their semantic scores. Clearly the semantic score is not a proxy for mutation count (Supplementary Fig. 9).

Once regular sequencing and surveillance are underway, a pipeline for analysing emerging sequences could identify outliers that may become future VOCs. Despite PLM scores distinguishing between VOCs well, a larger score does not always mean a variant will be successful. This could be due to measuring against Wuhan-Hu-1, which is no longer circulating and increasingly irrelevant SARS-CoV-2 variant. To tackle this, we created a ‘dynamic’ embedding that represents the average sequence circulating at a given time. This means we can check whether sequences are both divergent from the SARS-CoV-2 reference or from what is currently or recently circulating. By looking at a UK subset of sequences, we can assess this approach thanks to the UK’s high sequencing capacity paired with well-defined lineage waves. UK sequences were clustered to 99.9% similarity and a representative haplotype sequence for each cluster was embedded using the model. Waves of semantic change take place between all major VOC waves in the UK (Alpha to Delta to BA.1 to BA.2) (Fig. 7A). The non-VOC to Alpha transition resulted in sequences deviating by a dynamic semantic score of >2, although this quickly increased to almost three after the emergence of Delta. Following BA.2, the semantic changes in Omicron have been more incremental, with a gradual increase up to March 2023, followed by a decrease until July 2023. There have been several smaller waves after Omicron’s emergence, which suggests repeated dominance and replacement of consecutive distinct sub-lineages throughout this period. The haplotype embeddings also uncover the large diversity of semantic scores contained outside just the Pango representative sequences from Fig. 6D. Alpha and Delta have a larger range of dynamic semantic scores, equalling some of the most divergent Omicron sequences despite pre-dating them by months (Fig. 7A). This is also true to a lesser extent for non-VOC sequences. What is clear is that as a new VOC began to take over, sequences diverge very quickly from the average circulating sequence embedding.

The dynamic semantic score can highlight outlier sequences, although these might be related to factors beyond sequence properties such as epidemiology which make variant prediction difficult for both experimental and computational methods. However, outlier detection can filter out VOC-like sequences which are likely to be a subset of observed outliers. Using the first instances of each lineage (Fig. 7B), sequences above one and two standard deviations from the mean dynamic semantic score at each 3-month timepoint are classified as sequences of concern. We propose using these deviation thresholds as a sequence ‘warning’ indicators. When we filter for the first occurrences of known VOCs and important variants (Fig. 7C), we see that the majority (56%) are above one of the warning lines. These include all, at the time, WHO assigned VOC sequences (Alpha, Beta, Gamma, Delta and Omicron), with four of the five in the red ‘warning’ section. Three VUIs (Zeta, Iota and Mu) fall below the line, along with the BA.2.75 major variants and most of the XBB subvariants. VUIs were classified often due to their mutation composition or their early growth rates, yet many of these failed to maintain sustained spread due to competition from other variants. The Alpha variant was the dominant lineage in the UK when the VUI sequences arrived, which may explain their lack of spread in the UK, as well as their lack of risk classification. Alpha has a much larger dynamic semantic score than any of the early VOCs or VUIs and dominated in the UK until Delta’s emergence.

The initial XBB lineage was detected as a sequence of concern by the outlier detector, however subvariants were almost exclusively below the thresholds for detection (except XBB.1.5). This is likely due to the high degree of similarity between the XBB spike and descendant lineages as most of the sub-lineages only differ by one or two sites. While BA.2.75 is not selected by its first sequences, several of its sub-lineages (including the second sub-lineage to emerge) are classified into the red warning zone (Fig. 7B). The heavily mutated BA.2.86 lineage acquired an L455S substitution49 to produce JN.1, which subsequently became the dominant lineage circulating. Both BA.2.86 and JN.1 were identified by the outlier detection.

Using the haplotype sequences above the outlier warning lines, we counted the number of unique mutations in the spike protein that occur within those sequences (Fig. 7D). The NTD and RBD contain the sites with the greatest quantity and diversity of changes at a single position. Of the mutations identified, 46% of the positions were also present in the first sequence of an identified major variant sequence, with 39% of the exact mutations present.

Emergent antigenically distinct SARS-CoV-2 variants are likely to cause reinfection50. Consequences of this are quickly becoming apparent through the necessity for updated vaccines, with recent variants evading previously neutralising antibodies44,51. We assert that the dynamic embedding approach is a viable horizon scanning method that does not simply count mutations but accurately models them and the shifting immune landscapes these lineages exist within. By identifying sequences of concern and contextualising their mutations within circulating viral sequences and host immune pressures, dynamic embeddings offer another useful and easy to implement technique in the event of a future outbreak.

In conclusion, we have shown that the PLM metrics grammaticality and semantic score6 reveal characteristic properties of the SARS-CoV-2 spike protein sequence. Using just a single spike sequence we were able to determine the likely regions of variability in the spike protein, and identify regions where mutations were most likely to impact the structure and function of the protein. We describe how ESM-2 can map out evolutionary landscapes, identify epistatic effects and provide a method for horizon scanning of viral sequences of interest. Unlike other methods for predicting potential SARS-CoV-2 variant success, although representing improvements, for example, Ito et al.52 and Thadani et al.36, these need protein structures and multiple sequence alignments for training. In our implementation, the use of the pre-existing model ESM-2, negates the need for training and permits characterisation based on a single sequence. While there are clear areas for further improvement, this represents a paradigm shift for interpreting variation within protein sequences and could be applied to a novel pathogen when data is extremely sparse. A major strength of ESM-2 is its ability to learn broader evolutionary signals, which stems from its training data that spans the tree of life. While models trained from exclusively viral sequences may be produced in future, the sparsity of viral data included in training remains a major issue. Viral datasets currently contain only a small set of available and sufficiently distinct viral sequences. This is likely to limit the quality of any model inferences, however fine-tuning on such a set of sequences is more likely to yield improved performance. Models such as CovFit52 and LucaVirus53 have already demonstrated the power of fine-tuning general PLMs (LucaOne54 and ESM-25) on viral datasets, with more models likely to be released in future.

While ESM-2 is a good predictor of the computational and experimental metrics used to understand how mutations impact protein structures, these metrics do not correlate strongly with any one feature, indicating that PLM-derived scores are capturing something different about protein sequences. PLMs are now being used for DMS and variant effect prediction tasks across an impressive array of datasets, and often outcompete other methods55. They are helping improve the effectiveness of antibodies56 and even generate new proteins57. Our results here affirm that these models are useful and can be applied to great effect to understand novel viral pathogens. PLMs, thus, offer an exciting glimpse into the future use of language modelling within biology, and we should continue to press on with understanding the possibilities as well as the limitations of what these models can do.

Methods
ESM-2 and in silico DMS
The ESM-2 model was used for most of the analysis in the paper. ESM-2 was trained on the September 2021 version of UniRef50 meaning that SARS-CoV-2 spike is present in the training data5. UniRef clusters are made using MMseqs2 to group together sequences with at least N% identity overlap. UniRef90 is accordingly a set of clusters each represented by a sequence that has at least 90% identity to all other sequences in the cluster. UniRef50 is then a further clustering of the UniRef90 sequences where a sequence is selected for each cluster where the minimum identity is at least 50%. ESM-2 increases its diversity of sequences by using different UniRef90 sequences in certain minibatches during training. Since SARS-CoV-2 has a UniRef90 and UniRef50 cluster of which it is the representative sequence in both, it is very unlikely that anything other than the original SARS-CoV-2 Wuhan-Hu-1 spike protein would have been included in the training data. Further, the Omicron variant of the virus was not detected until November 20214,58 and at the time was the most divergent of the known circulating variants. The BA.1 lineage Omicron had an additional 39 differences (substitutions, insertions and deletions) relative to the Wuhan-Hu-1 spike protein, far below the number of differences required for the sequence to become its own cluster at either UniRef90 or UniRef50.

SARS-CoV-2 spike proteins were acquired by filtering the GISAID database (https://gisaid.org) for the earliest sequences from each Pango lineage with a fully intact spike protein sequence. The sequence was then embedded in the ESM-2 model to produce an embedding for the sequence. For the DMS data and for the embedding scores, the ESM-2 three billion parameter variant was used. The semantic score is equivalent to the L1 (Manhattan) distance between the embedding of the reference sequence (the spike protein from the original Wuhan-Hu-1 SARS-CoV-2 genome) and each of the Pango spike proteins. The grammaticality of a sequence is calculated as the product of the probabilities of each amino acid at each position in the spike protein. The probabilities come from a softmax of the last layer of the embedding and range between 0 and 1. However, many of the probabilities are small and for numerical stability the probabilities are represented in the log space. Relative grammaticality is the same as grammaticality, except the probability of a reference sequence is subtracted from the probability of the variant sequence so that the score is relative, in this case, to the SARS-CoV-2 reference sequence Wuhan-Hu-1. The sequence grammaticality represents the summed log-likelihoods of every reference position in the sequence, rather than just the mutated positions. For the DMS computations, a sequence was produced for every potential amino acid at every position in the spike protein. Each sequence was then embedded to calculate semantic scores and relative grammaticalities for every mutation.

Extended grammaticalities
During the training of ESM, the model is asked to predict the amino acid present at mask tokens that account for a random 15% subset of each training sequence. Allman et al.39 suggests that using the mask token is more appropriate for producing DMS scores, since this closely mirrors how the model is trained. This is different from Hie et al.6, where the probabilities are not taken from the mask token at the site of interest, but from the reference token. We checked how the masked approach performs relative to the standard grammaticality from Hie et al.6, while also investigating whether using the mutated sequence to calculate mutant probability was an appropriate measure, as this considers the rest of the sequence and the mutated site implicitly. Finally, we produced a grammaticality score that is weighted by the additional dissimilarity and accessibility components detailed in Thadani et al.36 to see if ESM metrics can also benefit from these additions. We included these alternative grammaticalities in the benchmarking (Fig. 5, Supplementary Fig. 10).

Evo-velocity
We used the evo-velocity package46 to embed the initial sequences using the ESM-2 650M parameter variant, and then performed velocity analysis and spearman’s rank for the sample dates against the pseudotime. The smaller model was chosen primarily due to hardware limitations of using the larger model, yet in several cases it has matched or even outperformed the larger model.

Epistasis experiments
The epistasis experiments used a SARS-CoV-2 BA.1 spike protein sequence embeddings using the ESM-2 three billion parameter variant. Due to the EPA insertion, the logits for this position were subsequently removed to map logits to spike’s three-dimensional protein structure. The BA.1 sequence contains several mutations relative to the Wuhan-Hu-1 SARS-CoV-2 reference spike sequence. Each of these mutations was reverted one by one back to the reference position, and the likelihood differences upon reversion were recorded. To eliminate noise, changes less than two standard deviations from the mean across all reverted mutations were removed and deemed not significant.

Dynamic embeddings and horizon scanning
Dynamic embeddings were computed by first gathering UK GISAID data and clustering sequences into variants with 99.9% similarity. These haplotype variant clusters produced 11,272 sample date labelled haplotypes with a sequence returned for each cluster. Each haplotype embedding was measured against a mean of the embeddings from the prior three-month period using the L1 distance, i.e., the semantic score. For sequence grammaticalities, mean sequence grammaticalities were calculated in a similar way and differences measured.

Assessing embeddings scores with known metrics
The language model’s metrics were first assessed using a Spearman’s rank against several known biological scores. Next, Support Vector Regression (SVR) was used as a simple model to fit the model scores as well as the embeddings and logits to the biologically relevant metrics. Models were fit to the data using five-fold cross validation, and a linear kernel for the SVR. Model results were reported as the average spearman’s rank between the folds, with the error bars as +/- 1 standard deviation from the mean. To assess statistical divergence between metrics at every site, scores were normalised to between 0 and 1, preserving rank order. The Jensen Shannon distance for each comparison was then calculated for every site independently using the philentropy (v.0.8.0) package in R. Amino acid positions used to define protein regions were as follows: S1 N-terminal (14-306), S1 RBD (331-528), S1 C-terminal (529-686), S2 (687-1273), all other positions were classed as indeterminate. Greater JSD implies greater divergence between two distributions. Significance was determined by a Mann–Whitney U test with Bonferronicorrection.

Selection analysis signals and entropy
Signals of ancestral evolutionary selection in the animal (bat and pangolin) sarbecovirus most closely related to SARS-CoV-2 referred to as the “nCoV” clade (see Lytras et al.59) were inferred on a set of 167 sarbecovirus genomes, accounting for recombination by inferring selection separately in each non-recombinant segment. These results are published in Martin et al.60 and presented in more detail in the following Observable notebook: https://observablehq.com/@spond/ncos-evolution-nov-2021. Sites under negative selection were inferred using the FEL and sites under positive selection using MEME61 by testing on internal branches of the nCoV clade. Sites denoted as conserved have the same amino-acid residue among all sarbecovirus sequences in the analysis. The variability of each site in the SARS-CoV-2 sequence was obtained from the entropy of the predicted distribution of credible evolutionary states.

Antibody accessibility and substitution probabilities
Structure-based epitope score, referred to as “accessibility”, which approximates antibody accessibility for each spike protein amino acid position, was calculated using BEpro software62, see Harvey et al.24. Scores relating to substitution probabilities, namely, ESST probability, Log PSSM and predicted ΔΔG, were obtained for every possible single amino acid substitution for the 6VXX SARS-CoV-2 spike structure (note that only values for Chain A are included in the results as data is generally identical across all three chains). ESST probability values were calculated using Environment Specific Substitution Table (ESST)63 after local structural environments were calculated by JOY43 Log Likelihood substitution values were calculated using Position Specific Scoring Matrices (PSSM) with the DELTA-BLAST64 algorithm in BlastX. It should be noted that this is a sequence-based method, so residue numbering does not match the numbering in 6VXX and values are available for residues not described by the 6VXX PDB file. ΔΔG values were predicted by FoldX65 software. FoldX uses empirical energy functions to predict the energetic effect of mutations to protein stability. The predicted ΔΔG quantifies the change in the free energy of unfolding between the wild-type and mutated structure. The 6VXX structure was first repaired with the RepairPDB function to fix residues with bad torsion angles, van der Walls’ clashes or total energy. Substitutions were then performed using the repaired structures on all three chains simultaneously using the BuildModel function, giving the change in free energy of unfolding, with negative values implying stabilising mutations.

Deep mutational scanning data
The receptor binding DMS data was taken from experimental DMS studies of the SARS-CoV-2 from Yisimayi et al. and Starr et al.40,41. The Wuhan-Hu-1 scores were taken as is, while the variant average score was calculated by averaging the scores for each position between each of the SARS-CoV-2 variant specific DMS results41. For the RBD mutational escape values we utilised the high-throughput mutation antibody escape profiling results presented in Yisimayi et al.40. This study used a panel of 1,350 monoclonal antibodies against all possible RBD substitutions. The backbone virus used was the SARS-CoV-2 BA.5 variant instead of Wuhan-Hu-1, however, this should still provide the most comprehensive dataset of unique substitutions’ effect on antibody escape. Our mutational escape metric is the average of the raw antibody escape values for each substitution on each site of the RBD across all tested monoclonal antibodies. The full spike DMS data was taken from Dadonaite et al. and DMS scores (“Binding”, “Escape” and “Entry”) are from a BA.2 and XBB.1.5 spike backbone and averaged together to give the presented score9.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "In the paper entitled “From a single sequence to evolutionary trajectories: protein language models capture the evolutionary potential of SARS-COV-2 protein sequences”, the authors use pLMs (ESM2) to understand the virus evolutionary trajectory, especially focus on the SARS-COV-2 spike protein. While some figures look good, the paper lacks sufficient novelty, and several figures and analyses appear to duplicate previously published work. Below, I outline several serious concerns.\nMinor Issue:\n1. Figure 1 Clarity: The upper panel of Figure 1 appears misleading, as it does not illustrate how the transition to the next position occurs. Please consider revising the figure to clearly demonstrate the progression between positions.\nMajor Issues:\n1. Lack of Novelty: The manuscript lacks sufficient novelty in its current form. Forecasting the evolutionary potential of viruses using language models—whether supervised fine-tuned models (https://doi.org/10.1038/s41467-023-39199-6 ), unsupervised models(https://www.science.org/doi/10.1126/scitranslmed.abk3445, https://www.cell.com/cell-systems/fulltext/S2405-4712(22)00038-2 ), or variational autoencoders(https://doi.org/10.1038/s41586-023-06617-0 )—has been previously validated and demonstrated through various analyses. This work does not offer new insights or solutions beyond what has already been established in the literature. Additionally, many figures seem to replicate data from published datasets with published models. For example, Hie et al. have shown that semantic and grammaticality scores are effective for forecasting novel mutations with protein language models.\n2. Overstatement of Contributions: In the Abstract, the authors claim, \"We demonstrate that the grammaticality and semantic scores represent novel metrics.\" However, these metrics were introduced by Hie et al. in 2021, who also used self-supervised language models to perform similar analyses. This overstatement diminishes the perceived originality of the work.\n3. Training Data Overlap with ESM-2: The ESM-2 model was trained on data up to April 2021, by which time the pandemic had been ongoing for over a year and numerous variants had emerged. Consequently, it's not surprising that the model captures evolutionary potential, as these variants are included in its training data. This raises concerns about the model's ability to predict truly novel mutations that were not part of its training set, particularly those emerging at the early stages of the pandemic.\n4. Epistasis Interaction Analysis: While the proposed method for identifying epistatic interactions is interesting, the manuscript lacks quantitative metrics or head-to-head comparisons with other methods that capture epistatic relationships. For instance, this paper (https://doi.org/10.1038/s41588-019-0432-9) identifies epistatic relationships based on protein 3D structures and contact maps. Benchmarking your method against such existing approaches as well as the GB1 epistasis dataset they used would strengthen your analysis and highlight its advantages.\n5. Justification for Machine Learning Models: The necessity of using additional features to train logistic regression or SVR classifier models is not clearly justified. Both models perform well even with simple sequence logits, suggesting that the prediction task may be relatively straightforward. A crucial aspect of employing machine learning models is assessing their out-of-distribution (OOD) generalization capabilities, the current analysis fails to do so. Without demonstrating improved performance or generalization, the added complexity may not be warranted.\n6. Repetition of Existing Work: Figure 5 seems to replicate analyses performed by Hie et al. and Han et al., without introducing any algorithmic or analytical novelty. Simply changing the model type does not constitute a substantial advancement in the field. The manuscript would benefit from novel approaches or insights that differentiate it from existing literature.\nConclusion:\nIn summary, the major concern is the lack of novelty, as much of the analysis presented has been explored extensively in prior literature. To enhance the contribution of this work, I recommend that the authors introduce novel insights or methodologies rather than compiling existing methods and analyses. Focusing on unique aspects or providing innovative solutions would significantly improve the manuscript's impact.", "reviewer_2": "Summary:\nLamb et al. used ESM2, applied to SARS-CoV-2 Spike protein, and analyzed how grammaticality and semantics can potentially capture (1) the mutability of a specific position, (2) epistatic interactions, (3) new metrics for mutational impacts, (4) outlier sequences in real time to aid in surveillance efforts. However, the authors simply report on the results of the model without adequate validation or benchmarking of its use to predict SARS-CoV-2 mutations or variants of concern (i.e., could the model have been used to forecast mutations that would later become dominant and how many predictions would have needed to be tested to accurate flag a true variant of concern). Additionally, the authors use “PLMs” consistently throughout the manuscript when describing their results when only ESM2 was evaluated. ESM2 was trained up until April 2021 and so a large proportion of SARS-CoV-2 sequences seen during the pandemic are included in training, yet this is not considered in the paper at all.\nMajor comments\nESM2 capture evolutionary potential\nThe authors claim that ESM scores can be used to capture evolutionary potential. However, no validation is provided for this. They simply report that a statistically significant difference was observed between the ESM2 scores between the S1 and S2 domains. While the authors provide some qualitative description of why this might be expected, no attempt is made to quantitatively assess this claim (for example, they could have investigated if this is also true for experimental DMSs).\nESM2 provide “new” metrics for assessing mutational impact\nThe authors found that neither grammaticality nor semantics correlate well with any previously used metrics (experimental or computational). While this is a significant shortcoming, the authors report it as a discovery of “new metrics” with little attempt to describe what these metrics might be.\nDetecting the distinct nature of the variants of concern on emergence\nThe authors claim that ESM can identify differences between VoCs and state that “Clearly the semantic score is not a proxy for mutation count”. However, this is not “clear” at all from the analysis they presented. One possibility is to show what the plots (e.g. those in Fig 6a) would be if they plotted the number of mutations instead of the ESM metrics.", "reviewer_3": "In Lamb et al., the authors explore the use of protein language models (PLMs), specifically ESM-2, to characterize emerging mutations in SARS-CoV-2 spike protein. The paper introduces two novel metrics for assessing the effect of amino acid mutations: grammaticality (the sum of the log probabilities for each amino acid) and semantic score (the distance in embedding space). These scores are used to:\n- Identify regions of the spike protein which are more prone to accumulating mutations\n- Detect potential epistatic interactions between residues\n- Compare computational predictions with DMS studies and other quantities measured from the protein structure\n- Reconstruct the virus evolutionary history from the Wuhan-Hu-1 strain and predict emerging variants\nThe authors acknowledge and discuss that the grammaticality or semantic scores exhibit weak correlations with experimental data, such as DMS studies or quantities derived from protein structures. The approach is intriguing, as the pretrained ESM-2 model can be applied immediately upon identifying the first sequence of a new pathogen. However, in my opinion, the paper lacks benchmarking against other computational methods designed to predict the impact of mutations or epistatic interactions. Overall, the paper falls short of convincingly demonstrating that the proposed metrics provide meaningful predictive or descriptive value, nor does it present clear evidence of improvements over existing methods.\nSpecific comments that should be addressed are:\n- Epistatic interactions. The authors identify a set of epistatic interactions between residues in the spike protein. This finding is interesting and supported by the fact that many identified epistatic interactions appear to be close in the 3d structure. However, it is difficult to assess its robustness without a comparison to other computational methods designed to identify epistasis. Benchmarking against tools trained on sequences from emerging SARS-CoV-2 strains such as those described in https://www.pnas.org/doi/full/10.1073/pnas.2104241118 or https://genomebiology.biomedcentral.com/articles/10.1186/s13059-024-03355-y could strengthen their claims.\n- Benchmarking Grammaticality and Semantic Scores: The authors argue that their method differs from existing approaches because it does not require an MSA or protein structure. While this claim is valid, it is important to note that ESM-2 is still trained on protein sequence evolutionary data. I believe it would be valuable to compare the grammaticality and semantic scores with other established metrics, such as the Evo score provided by https://evescape.org/emergingvariants or the mutability score described in https://www.pnas.org/doi/full/10.1073/pnas.2113118119. These approaches, which rely on MSAs of homologous sequences (available before the pandemic), could serve as meaningful benchmarks. Examining correlations between these metrics and the grammaticality or semantic scores could provide deeper insights into their predictive utility for assessing the impact of mutations.\n- Grammaticality and Semantic Scores vs. Experimental Measures: The correlation of grammaticality and semantic scores with deep mutational scanning (DMS) studies is limited. Moreover, these metrics only partially reflect quantitative features derived from protein structures, such as accessibility, B-factors, or |ΔΔG| values. The authors note that “these metrics do not correlate strongly with any one feature, indicating that PLM-derived scores are capturing something novel about protein sequences.” However, similarly low correlation values have been reported for models trained on evolutionary data, such as those described in https://www.pnas.org/doi/full/10.1073/pnas.2113118119 and https://www.nature.com/articles/s41586-023-06617-0. The benchmarking approach proposed in the previous point could help clarify which constraints are specifically captured by the ESM-2 model.\n- Not working website. The website providing the ESM-2 scores (https://sars2.cvr.gla.ac.uk/) is currently non-functional (502 Bad Gateway). This limits the reproducibility and validation of their findings.\n- Predicting emerging mutations and fitness of VOC strains. To provide evidence that the model can be used for predictive purposes (predict emerging variants) the authors a) analyse the congruence between the phylogenetic tree topology and the evo-velocity derived structure and b) use dynamic embedding to compute the deviation from a average sequence circulating a given time. While this is interesting, it is difficult to assess where this model could be used to predict emerging mutations. A more compelling demonstration would involve testing whether these scores can prospectively identify mutations that emerged during SARS-CoV-2 evolution (similar to what done in https://www.pnas.org/doi/full/10.1073/pnas.2113118119 and https://www.nature.com/articles/s41586-023-06617-0).\nMinors:\n- The paper lacks clarity on which score - grammaticality or semantic - should be prioritized for predictive purposes. Resolving this confusion and providing a clearer explanation of their respective roles in prediction would enhance the overall impact of the work.\n- There is no title for the “Conclusion” section.", "reviewer_4": "I co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is part of the Nature Communications initiative to facilitate training in peer review and to provide appropriate recognition for Early Career Researchers who co-review manuscripts."}, "author_response": "Response to reviewer comments, NCOMMS-24-59219-T\nReviewer #1 (Remarks to the Author):\nIn the paper entitled “From a single sequence to evolutionary trajectories: protein language\nmodels capture the evolutionary potential of SARS-CoV-2 protein sequences”, the authors\nuse pLMs (ESM2) to understand the virus evolutionary trajectory, especially focus on the\nSARS-COV-2 spike protein. While some figures look good, the paper lacks sufficient novelty,\nand several figures and analyses appear to duplicate previously published work. Below, I\noutline several serious concerns.\nWe understand the reviewer’s concern over novel methodology. However, our purpose was\nnot to make new models but rather to explore the utility of existing ones. Specifically, our\npaper is concerned with applying the concepts of grammaticality and semantic scores, as\nsuggested by Hie et al in their Science 2021 paper, and we show the frontier protein\nlanguage model, ESM-2, can make useful and novel predictions. Especially, when applied to\nsingle sequence-based epistatic interaction detection and horizon scanning via dynamic\nreferencing. To the best of our knowledge, these results and the analytical framework are\nnovel.\nMinor Issue:\n1. Figure 1 Clarity: The upper panel of Figure 1 appears misleading, as it does not illustrate\nhow the transition to the next position occurs. Please consider revising the figure to\ndemonstrate the progression between positions clearly.\nApologies, the figure has been amended to illustrate the progression between positions\nmore clearly.\nFigure 1. Schematic summary of our methodology. Deep mutational scanning involves taking a\nsequence and mutating every position to each of the possible alternative amino acids. These are then\npassed to the PLM to produce embeddings and logits which can be used for downstream tasks or to\nproduce the metrics relative grammaticality and semantic score. Identifying epistatic interactions\ninvolves reverting the mutations from a SARS-CoV-2 variant (here BA.1) and measuring the effect this\nhas on the other likelihoods. Embeddings and log-likelihoods (logits) can be used for surveillance by\nproducing metrics over timescales for newly emergent protein sequences of concern, or by looking at\nevolutionary trajectories such as with evo-velocity.\nMajor Issues:\n1. Lack of Novelty: The manuscript lacks sufficient novelty in its current form. Forecasting\nthe evolutionary potential of viruses using language models—whether supervised fine-tuned\nmodels (https://doi.org/10.1038/s41467-023-39199-6 ), unsupervised\nmodels(https://www.science.org/doi/10.1126/scitranslmed.abk3445,\nhttps://www.cell.com/cell-systems/fulltext/S2405-4712(22)00038-2 ), or variational\nautoencoders(https://doi.org/10.1038/s41586-023-06617-0 )—has been previously validated\nand demonstrated through various analyses. This work does not offer new insights or\nsolutions beyond what has already been established in the literature. Additionally, many\nfigures seem to replicate data from published datasets with published models. For example,\nHie et al. have shown that semantic and grammaticality scores are effective for forecasting\nnovel mutations with protein language models.\nOur novelty in comparison with the papers listed here are the following:\n1. We focus on what can be accomplished with a single sequence and use of an\navailable existing model, ESM-2, eg, following the spillover of a novel virus with\nlimited data. This scenario occurs frequently as we are typically constrained by\navailable sequences for multiple alignment-based methods and limited experimental\ndata. Almost every other tool requires more than just a single viral sequence to\nperform any type of analysis, which makes using PLMs an attractive proposition. This\nis clearly different from the premise in supervised fine-tuned models\n(https://doi.org/10.1038/s41467-023-39199-6)\n2. We proposed a novel single-sequenced epistatic interaction detection analysis and\nvalidated it with real experimental data.\n3. We proposed a novel horizon-scanning score that is agnostic to any experiment-\nderived features such as ACE binding affinities and epidemiological features\nconsidered in https://www.science.org/doi/10.1126/scitranslmed.abk3445.\n4. Hie et al (https://pubmed.ncbi.nlm.nih.gov/33446556/) proposed that grammar and\nsemantic scores are analogous to protein fitness and antigenicity. We show that\nthese metrics do not correspond to these properties when applied to a generic PLM,\nsuch as ESM-2. This is likely due to the training datasets not being constructed from\nsingle viral protein alignments dominated by antigenic mutations. This is not\nsomething that to our knowledge has been presented previously. We assess the\nutility of these metrics and identify what they capture given a more foundational PLM,\nESM-2. We then propose approaches that can be used when more sequence data is\navailable and show that sequence logits potentially identify epistatic sites within the\nprotein, a technique not mentioned in Hie et al.\nWe do agree with the reviewer that the Evescape model is a good comparison in the single\nsequence scenario. We have conducted a systematic benchmarking experiment with\nEvescape on predicting mutations that occurred during the pandemic.\nSupplementary Figure 10: Data from each of the mutational scans were filtered so that there is a\nmeasurement for every mutation in each dataset to allow for equivalence. The full spike mutations\nexclude the RBD Experimental mutations since the mutation set could only increase selection of RBD\nmutants given the lack of other regions in this data. (A) Barchart for each feature showing the number\nof mutations from variant of concern sequences present in the top 10% of shared feature predictions.\nThe top set of bars shows the shared mutations within the RBD, while the bottom shows mutations\nshared across the whole of the spike protein. (B) A cumulative sum of shared DMS mutations that\nhave appeared at least 100 times during the pandemic that appear in the top 10 % of all features.\n2. Overstatement of Contributions: In the Abstract, the authors claim, \"We demonstrate that\nthe grammaticality and semantic scores represent novel metrics.\" However, these metrics\nwere introduced by Hie et al. in 2021, who also used self-supervised language models to\nperform similar analyses. This overstatement diminishes the perceived originality of the\nwork.\nWe are not claiming that these are our new metrics and apologise for any confusion; we\nhave removed this sentence for clarity. Rather we show these are ‘new' metrics in that they\ncapture something new biologically. Our aim is to extend the work of Hie et al who showed\nsemantic score can be used to measure antigenicity. We show that this property is not\nrestricted to antigenicity when computed using another model such as ESM-2. This\nrelationship appears derived from the training data, which in Hie et al was virus protein\nspecific. Mutations in these virus datasets are mostly antigenic, hence we suspect this is\nwhere the semantic score to antigenicity relationship is derived.\nMore importantly, we demonstrate the base protein language model, ESM-2, can be used to\nrecapitulate information from a single sequence that would normally require alignment data,\nand we demonstrate this can be used to assess epistasis incredibly from data obtained from\na single sequence. We believe this demonstrates the power of protein language models and\nbuilds on Hie et al.’s suggestion to use these metrics.\n3. Training Data Overlap with ESM-2: The ESM-2 model was trained on data up to April\n2021, by which time the pandemic had been ongoing for over a year and numerous variants\nhad emerged. Consequently, it's not surprising that the model captures evolutionary\npotential, as these variants are included in its training data. This raises concerns about the\nmodel's ability to predict truly novel mutations that were not part of its training set,\nparticularly those emerging at the early stages of the pandemic.\nThe reviewer is correct that ESM-2 was trained after the pandemic started. However, the\ntraining data only contains the Wuhan-Hu-1 sequence, ie, just one sequence, based on\nUniRef selection criteria. Specifically, ESM-2 was trained on the Sep 2021 version of\nUniRef50 meaning that the SARS-CoV-2 spike was present at the time of training. However\nmore importantly, UniRef clusters are made using MMseqs to group together sequences\nwith at least N% identity overlap. UniRef90 is accordingly a set of clusters each represented\nby a sequence that has at least 90% identity to all other sequences in the cluster. UniRef50\nis then a further clustering of the UniRef90 sequences where a sequence is selected for\neach cluster where the minimum identity is at least 50%. ESM-2 increases its diversity of\nsequences by using different UniRef 90 sequences in certain mini batches during training.\nSince SARS-CoV-2 has a UniRef90 and 50 cluster of which it is the representative sequence\nin both, then only the original SARS-CoV-2 Wuhan-Hu-1 spike protein would have been\nincluded in the training data.\nFurther, the Omicron variant of the virus was not detected until November 2021 and at the\ntime was the most divergent of the known circulating variants. The BA.1 lineage Omicron\nhad an additional 39 differences (substitutions, insertions and deletions) relative to the\nWuhan-Hu-1 spike protein, far below the number of differences required for the sequence to\nbecome its own cluster at either UniRef90 or 50.\n4. Epistasis Interaction Analysis: While the proposed method for identifying epistatic\ninteractions is interesting, the manuscript lacks quantitative metrics or head-to-head\ncomparisons with other methods that capture epistatic relationships. For instance, this paper\n(https://doi.org/10.1038/s41588-019-0432-9) identifies epistatic relationships based on\nprotein 3D structures and contact maps. Benchmarking your method against such existing\napproaches as well as the GB1 epistasis dataset they used would strengthen your analysis\nand highlight its advantages.\nWe agree with the reviewer that additional benchmarking is required for predicting epistatic\ninteractions. However, the GB1 analysis and experiment was not suitable for this study\nbecause it focuses specifically on contact map prediction, which is a subset of our broader\ndefinition of epistatic interactions, and does not concern a viral protein. We did investigate\nproducing a contact map for GB1, since there are PLM based methods for predicting contact\nmaps (namely https://www.pnas.org/doi/10.1073/pnas.2406285121). We found that ESM-2\ncan be used to identify contacts, as shown below. However for inclusion in our paper, we\ndecided to focus on a more directly relevant SARS-CoV-2 dataset.\nGround Truth ESM-2 650M ESM-2 650M with 0.5 threshold\nInstead in this revision, we have conducted an analysis using a specific SARS-CoV-2\nexperimental dataset which investigated a very similar question to what we were attempting\nto answer (https://www.nature.com/articles/s41467-022-34506-z). Moulana et al. produce\nexperimental data of all possible intermediate sequences between Wuhan-Hu-1 and the\nBA.1 Omicron RBD. For each variant they measure ACE-2 binding affinity, and using a\nbiochemical epistasis model process their data into a sum of mutational effects. These begin\nwith single mutational effects, increasing to pairwise effects and up to 5th-order effects.\nSpecifically, we used the pairwise effects to compare to our BA.1 reversion data and show\nthat several of the top pairwise epistatic interactions between BA.1 mutations are shown by\nour method. We report these findings in the updated Figure 3. We also discuss the caveats\nof this analysis, which ties into our benchmarking results. This includes a discussion about\nthe measured epistatic effects being related to ACE-2 binding affinity, and how this is but\none of many possible epistatic effects that a mutation can have. Our scores of course cannot\nmeasure all of these effects at the same time, however, they can pick up a disruption that\nlikely affects at least one of these features of the protein. Often regions involved in binding\nare also relevant to escape for example. We believe our approach here is interesting given\nits lack of knowledge about sequences spreading in the pandemic, and its relative ease of\nuse to identify putative epistatic effects between sites.\nFigure 1. (A) Monomeric structures showing the changes in probabilities for three example\nsubstitutions in SARS-CoV-2's spike protein: E484A, D614G and N969K. This is the negative of the\nreversion likelihood. The mutation site is coloured yellow, red sites increase in probability while blue\nsites decrease. Mutation probabilities were only shown if they were outside two standard deviations of\nthe mean change from all BA.1 reversions. (B) Heatmap of BA.1 RBD mutations and the change in\nlikelihoods at other BA.1 mutation sites. Larger circles indicate the sites were two standard deviations\naway from the average change, smaller circles indicate one standard deviation away. (C)\nExperimental epistasis data from Moulana et al.21 showing the pairwise effect of mutations on ACE-2\nbinding affinity. Positive values indicate a pairwise increase in ACE-2 binding affinity. Larger circles\nindicate the sites were identified by ESM using two standard deviations away from the average\nchange, smaller circles indicate one standard deviation away. The smallest points were not identified\nby ESM. (D) Experimental epistasis data [2] coloured based on the ranking of the absolute epistatic\neffect values. This allows large negative values to be counted as high ranking and vice versa.\nWe have also shown favourable comparisons to an existing computational method Innocenti\net al. which relies on mutual information from a large number of pandemic sequences for\npredicting epistasis shown below. Our approach captures more experimentally validated\nevents.\nSupplementary Figure 12. Results from Moulana et al. filtered by the RBD epistatic interactions\nidentified by Innocenti et al.\n5. Justification for Machine Learning Models: The necessity of using additional features to\ntrain logistic regression or SVR classifier models is not clearly justified. Both models perform\nwell even with simple sequence logits, suggesting that the prediction task may be relatively\nstraightforward.\nThe “additional features”, i.e., the embedding and the sequence logits are not in fact\nadditional. The semantic score and grammaticalities are calculated using these features, i.e.,\nthe semantic score is just the L1 distance between embeddings and the grammaticality is\njust the logit at the mutated position. The purpose of including these features during the\nregression section was to reaffirm that the model representations are meaningful and\ncontain much of the information that is required for successful prediction. The reason for\nselecting a linear regression and an SVR regression was to show that two of the simplest\nmethods (one linear, the other non-linear) can produce meaningful results. It also highlights\nthat reducing complex representations (2048 dimensional vector embeddings and sequence\nlength vectors of probabilities) into single scores is often simplistic.\nA crucial aspect of employing machine learning models is assessing their out-of-distribution\n(OOD) generalization capabilities, the current analysis fails to do so. Without demonstrating\nimproved performance or generalization, the added complexity may not be warranted.\nAs mentioned above, the SARS-CoV-2 spike protein is not well represented in the ESM-2\ntraining dataset. Further, the maximum sequence length used to train ESM-2 was 1024\namino acids, which is smaller than the spike protein. The vast majority of the training dataset\n(UniRef50 or UniRef90) is not comprised of viral proteins, and an even smaller proportion of\nthose are viral glycoproteins. As such, most of what we have tested can be considered\nOOD already and we only show held-out test scores.\n6. Repetition of Existing Work: Figure 5 seems to replicate analyses performed by Hie et al.\nand Han et al., without introducing any algorithmic or analytical novelty. Simply changing the\nmodel type does not constitute a substantial advancement in the field. The manuscript would\nbenefit from novel approaches or insights that differentiate it from existing literature.\nEvo-velocity was used with the updated ESM-2 model to show that the model could\ndifferentiate between the circulating SARS-CoV-2 sequences in a meaningful way. We see\nthat the method closely approximates what we observe using a phylogenetic tree of the\nSARS-CoV-2 spikes, and the trajectories correlate well with the sequence sample dates\nduring the pandemic. While this in and of itself is not novel, we feel that it further\ndemonstrates that the models have some understanding of the relatedness of sequences.\nThis, in addition to showing the models can predict properties of sequence (correlation\nsection) helps to show that the models can be used for biological insights.\nConclusion:\nIn summary, the major concern is the lack of novelty, as much of the analysis presented has\nbeen explored extensively in prior literature. To enhance the contribution of this work, I\nrecommend that the authors introduce novel insights or methodologies rather than compiling\nexisting methods and analyses. Focusing on unique aspects or providing innovative\nsolutions would significantly improve the manuscript's impact.\nOverall, we hope the reviewer will agree with us about the novelty of our study: 1) The\nsingle-sequence premise. 2) Novel epistatic interaction detection 3) Novel horizon scanning\nanalysis. 4) Generalisation of existing metrics such as grammaticality and semantic scores.\nThese contributions have a broad impact on pandemic preparedness and response.\nReviewer #2 (Remarks to the Author):\nSummary:\nLamb et al. used ESM2, applied to SARS-CoV-2 Spike protein, and analyzed how\ngrammaticality and semantics can potentially capture (1) the mutability of a specific position,\n(2) epistatic interactions, (3) new metrics for mutational impacts, (4) outlier sequences in real\ntime to aid in surveillance efforts. However, the authors simply report on the results of the\nmodel without adequate validation or benchmarking of its use to predict SARS-CoV-2\nmutations or variants of concern (i.e., could the model have been used to forecast mutations\nthat would later become dominant and how many predictions would have needed to be\ntested to accurate flag a true variant of concern). Additionally, the authors use “PLMs”\nconsistently throughout the manuscript when describing their results when only ESM2 was\nevaluated. ESM2 was trained up until April 2021 and so a large proportion of SARS-CoV-2\nsequences seen during the pandemic are included in training, yet this is not considered in\nthe paper at all.\nThanks to the reviewer for their insightful comments. We have addressed all of them below,\nincluding adding a much expanded benchmarking section, replaced (where necessary)\nPLMs with the ESM-2 acroynm, and discussed the training data used for ESM-2.\nOn the ESM-2 training data point we have responded to this in Reviewer 1’s comments as\nwell but to reiterate here:\nThe reviewer is correct that ESM-2 was trained after the pandemic started. However, the\ntraining data only contains the Wuhan-Hu-1 sequence, ie, just one sequence, based on\nUniRef selection criteria. Specifically, ESM-2 was trained on the Sep 2021 version of\nUniRef50 meaning that the SARS-CoV-2 spike was present at the time of training. However\nmore importantly, UniRef clusters are made using MMseqs to group together sequences\nwith at least N% identity overlap. UniRef90 is accordingly a set of clusters each represented\nby a sequence that has at least 90% identity to all other sequences in the cluster. UniRef50\nis then a further clustering of the UniRef90 sequences where a sequence is selected for\neach cluster where the minimum identity is at least 50%. ESM-2 increases its diversity of\nsequences by using different UniRef 90 sequences in certain mini batches during training.\nSince SARS-CoV-2 has a UniRef90 and 50 cluster of which it is the representative sequence\nin both, then only the original SARS-CoV-2 Wuhan-Hu-1 spike protein would have been\nincluded in the training data.\nFurther, the Omicron variant of the virus was not detected until November 2021 and at the\ntime was the most divergent of the known circulating variants. The BA.1 lineage Omicron\nhad an additional 39 differences (substitutions, insertions and deletions) relative to the\nWuhan-Hu-1 spike protein, far below the number of differences required for the sequence to\nbecome its own cluster at either UniRef90 or 50.\nMajor comments\nESM2 capture evolutionary potential\nThe authors claim that ESM scores can be used to capture evolutionary potential. However,\nno validation is provided for this. They simply report that a statistically significant difference\nwas observed between the ESM2 scores between the S1 and S2 domains. While the\nauthors provide some qualitative description of why this might be expected, no attempt is\nmade to quantitatively assess this claim (for example, they could have investigated if this is\nalso true for experimental DMSs).\nWe thank the reviewer for the suggestion. We have now provided an explicit test for\nmutations captured by experimental DMS data between S1 and S2.\nThe DMS data comes from Dadonaite et al., which measures the mutational consequences\non ACE2 binding, antibody escape, and cell entry. Out of the three measurements, ACE2\nbinding affinity is significantly different between S1 and S2 regions, matching results from\nsingle sequence ESM2 analysis and entropy.\nCorrelations with experimental DMS (those that go outside the RBD), show only Binding has\na significant difference between subunits. Escape and Entry are not significant, while other\nexperimental datasets only measure within the RBD.\nThis is now mentioned in the paper:\n“Available full spike experimental DMS results from Dadonaite et al. show Binding is the only\nmeasure with a significant difference between subunits (p-value = 1.70e-06, two sided\nMann-Whitney U test), which might be expected given that S1 contains the spike binding\ndomain.”\nThe full details of the tests are the following:\nESM2 provide “new” metrics for assessing mutational impact\nThe authors found that neither grammaticality nor semantics correlate well with any\npreviously used metrics (experimental or computational). While this is a significant\nshortcoming, the authors report it as a discovery of “new metrics” with little attempt to\ndescribe what these metrics might be.\nApologies we meant new in the sense of a generically useful ‘new’ model for the community,\nnot a new model by us. We’ve also discussed the scores in more detail. The grammaticality\nscores for instance tend to align better with B-Factor and ESST scores. The Semantic score\ncorrelated better with PSSM and ESST. However, we make a point of saying that we expect\nthese correlations to be weak as nothing in the training data other than the bias of the\nsequences themselves explicitly enforce these metrics as a constraint. As such, it seems\nunreasonable to expect that the ESM-2 scores should strongly correlate with all features, let\nalone any single one. In fact, this would likely limit their usability and predictiveness, which in\nlater figures we show to be comparable with experimental measures. The majority of this\ndiscussion is included in the expanded section “PLMs provide orthogonal metrics for\nassessing mutational impact”.\nOur generic analysis has particular value compared to Hie et al, in which the scores were\noriginally proposed and the analysis was based on a viral-specific model. There,\ngrammaticality and semantic scores are equated to fitness and antigenicity. We show the\nmetrics are no longer specifically related to these two properties. Instead, they have the\npotential to be related to a boarder set of properties.\nDetecting the distinct nature of the variants of concern on emergence\nThe authors claim that ESM can identify differences between VoCs and state that “Clearly\nthe semantic score is not a proxy for mutation count”. However, this is not “clear” at all from\nthe analysis they presented. One possibility is to show what the plots (e.g. those in Fig 6a)\nwould be if they plotted the number of mutations instead of the ESM metrics.\nThis is a good suggestion, we have done this (see: Supplementary Figure 9).\nReviewer #3 (Remarks to the Author):\nIn Lamb et al., the authors explore the use of protein language models (PLMs), specifically\nESM-2, to characterize emerging mutations in SARS-CoV-2 spike protein. The paper\nintroduces two novel metrics for assessing the effect of amino acid mutations: grammaticality\n(the sum of the log probabilities for each amino acid) and semantic score (the distance in\nembedding space). These scores are used to:\n- Identify regions of the spike protein which are more prone to accumulating mutations\n- Detect potential epistatic interactions between residues\n- Compare computational predictions with DMS studies and other quantities measured from\nthe protein structure\n- Reconstruct the virus evolutionary history from the Wuhan-Hu-1 strain and predict\nemerging variants\nThe authors acknowledge and discuss that the grammaticality or semantic scores exhibit\nweak correlations with experimental data, such as DMS studies or quantities derived from\nprotein structures. The approach is intriguing, as the pretrained ESM-2 model can be applied\nimmediately upon identifying the first sequence of a new pathogen. However, in my opinion,\nthe paper lacks benchmarking against other computational methods designed to predict the\nimpact of mutations or epistatic interactions. Overall, the paper falls short of convincingly\ndemonstrating that the proposed metrics provide meaningful predictive or descriptive value,\nnor does it present clear evidence of improvements over existing methods.\nThanks for these helpful comments, we’ve addressed these concerns, see responses below.\nWe hope you agree these changes have improved the manuscript and make the paper more\ncomparable with alternative methods such as Evescape and epistasis predictors.\nSpecific comments that should be addressed are:\n- Epistatic interactions. The authors identify a set of epistatic interactions between residues\nin the spike protein. This finding is interesting and supported by the fact that many identified\nepistatic interactions appear to be close in the 3d structure. However, it is difficult to assess\nits robustness without a comparison to other computational methods designed to identify\nepistasis. Benchmarking against tools trained on sequences from emerging SARS-CoV-2\nstrains such as those described in https://www.pnas.org/doi/full/10.1073/pnas.2104241118\nor https://genomebiology.biomedcentral.com/articles/10.1186/s13059-024-03355-y could\nstrengthen their claims.\nAs described above we have added comparison to experimental data as a benchmark for\nour epistasis method, since this has the most direct one-to-one comparison. This data is also\nused by Innocenti et al. as part of their validation. We compared their graph of RBD\nmutations to our set of identified epistatic interactions between RBD mutations, and show\nthat our method identifies a larger number of those interactions. This figure is in the\nsupplementary.\nOur results: Innocenti et al.\n- Benchmarking Grammaticality and Semantic Scores: The authors argue that their method\ndiffers from existing approaches because it does not require an MSA or protein structure.\nWhile this claim is valid, it is important to note that ESM-2 is still trained on protein sequence\nevolutionary data. I believe it would be valuable to compare the grammaticality and semantic\nscores with other established metrics, such as the Evo score provided by\nhttps://evescape.org/emergingvariants or the mutability score described in\nhttps://www.pnas.org/doi/full/10.1073/pnas.2113118119. These approaches, which rely on\nMSAs of homologous sequences (available before the pandemic), could serve as\nmeaningful benchmarks. Examining correlations between these metrics and the\ngrammaticality or semantic scores could provide deeper insights into their predictive utility\nfor assessing the impact of mutations.\nIn this revision, we have added a systematic benchmarking study comparing the\nEvescape(https://evescape.org/emergingvariants) and mutability\nscores(https://www.pnas.org/doi/full/10.1073/pnas.2113118119), and have added this to the\npaper. Specifically, there is now a much-expanded correlation section.\nThe following has been added to the main text explaining the benchmarking:\n“We benchmarked our results against two other computational scores: Evescape from\nThadani et al. and mutability scores from Rodriguez-Rivas et al.. Evescape is a variational\nautoencoder that was trained using pre-pandemic spike protein sequences. It produces a\nfitness score (Evescape Fitness) that when combined via a normalised weighted sum with a\nstructure derived contact accessibility and amino acid dissimilarity metric produces the\nEvescape score. There are two mutability score variants35, the direct coupling analysis (DCA)\nscore and the Independent Site Model (IND) score. The IND score is a site independent\nmutability score, while the DCA score uses coupling analysis to try to incorporate signals of\nmutational epistasis. Similarly to Evescape, these scores are obtained using pre-pandemic\nviral protein sequences. We also tested other methods for computing grammaticality:\nmasked grammaticality suggested by Allman et al.; mutated grammaticality; the relative\ngrammaticality measure supplemented with the dissimilarity and accessibility as shown in\nThadani et al..”\nWe have expanded the Figure 5 correlation and regression analysis to include both Evescape\nand Mutability scores. Broadly, all scores correlate weakly with experimental and most\ncomputational metrics. ESM-2 scores with accessibility and dissimilarity are competitive with\nEvescape, the best competitor model. ESM-2 embeddings and logits outperform all individual\nscores significantly, suggesting these representations contain much richer information than\nsingle scores.\nFigure 5: Comparison using (A) Spearman's Rank correlations between the semantic score,\ngrammaticality and relative sequence grammaticality and traditional metrics. Bars with an asterix\nmeans the correlation was found to be significant after a Bonferroni correction (multiple test\ncorrection). And, (B) Spearman's Rank correlations between the language model metric and the\ntraditional metric. Each pair was fitted using a support vector regression model using an RBF kernel,\nwith 5-fold cross validation which was repeated on 3 randomised splits of the data. Bars represent the\nmean of the correlations, with the error bar +/- 1 standard deviation of the correlations. Bars with an\nasterix mean the correlation was found to be significant after a Bonferroni correction.\n- Grammaticality and Semantic Scores vs. Experimental Measures: The correlation of\ngrammaticality and semantic scores with deep mutational scanning (DMS) studies is limited.\nMoreover, these metrics only partially reflect quantitative features derived from protein\nstructures, such as accessibility, B-factors, or |ΔΔG| values. The authors note that “these\nmetrics do not correlate strongly with any one feature, indicating that PLM-derived scores\nare capturing something novel about protein sequences.” However, similarly low correlation\nvalues have been reported for models trained on evolutionary data, such as those described\nin https://www.pnas.org/doi/full/10.1073/pnas.2113118119 and\nhttps://www.nature.com/articles/s41586-023-06617-0. The benchmarking approach\nproposed in the previous point could help clarify which constraints are specifically captured\nby the ESM-2 model.\nAs above, further discussion and benchmarking on scores has been added.\n- Not working website. The website providing the ESM-2 scores (https://sars2.cvr.gla.ac.uk/)\nis currently non-functional (502 Bad Gateway). This limits the reproducibility and validation of\ntheir findings.\nApologies these online parts have now been relocated to an Observable dashboard where\nall of the plots are now available, along with the accompanying data used to produce the\nplots. See: https://observablehq.com/@cvr-bioinfo/from-a-single-sequence-nature-\ncommunications.\n- Predicting emerging mutations and fitness of VOC strains. To provide evidence that the\nmodel can be used for predictive purposes (predict emerging variants) the authors a)\nanalyse the congruence between the phylogenetic tree topology and the evo-velocity derived\nstructure and b) use dynamic embedding to compute the deviation from an average\nsequence circulating a given time. While this is interesting, it is difficult to assess where this\nmodel could be used to predict emerging mutations. A more compelling demonstration would\ninvolve testing whether these scores can prospectively identify mutations that emerged\nduring SARS-CoV-2 evolution (similar to what done in\nhttps://www.pnas.org/doi/full/10.1073/pnas.2113118119 and\nhttps://www.nature.com/articles/s41586-023-06617-0).\nWe have updated this section to include predictions based on our dynamic embedding\napproach and added further comparisons in a similar way to the suggested Evescape\n(https://evescape.org/emergingvariants) and Mutability\n(https://www.pnas.org/doi/full/10.1073/pnas.2113118119) methods. Mutation prediction\nresults are shown in the benchmarking section in Supplementary Figure 10.\nSupplementary Figure 10: Data from each of the mutational scans were filtered so that there is a\nmeasurement for every mutation in each dataset to allow for equivalence. The full spike mutations\nexclude the RBD Experimental mutations since the mutation set could only increase selection of RBD\nmutants given the lack of other regions in this data. (A) Barchart for each feature showing the number\nof mutations from variant of concern sequences present in the top 10% of shared feature predictions.\nThe top set of bars shows the shared mutations within the RBD, while the bottom shows mutations\nshared across the whole of the spike protein. (B) A cumulative sum of shared DMS mutations that\nhave appeared at least 100 times during the pandemic that appear in the top 10 % of all features.\nWe show that ESM-2 metrics typically perform worse than competitors when considering\nVOC mutations, with the exception of the ESM relative grammaticality when it is extended\nusing the accessibility and dissimilarity components from Evescape. When these\ncomponents are added, we see that comparable results are achieved relative to other\nexperimental and computational methods (Supplementary Figure 10A). When detecting\nmutations that appear during the pandemic as a whole, we show that standalone ESM\nmetrics outperform all other approaches (Supplementary Figure 10B).\nWe then show a new method to identify sequences of interest in an expanded “Detecting the\ndistinct nature of the variants of concern on emergence” section. The latest results are\npresented in Figure 7.\nWe have added the following text to describe the results:\n“Using the first instances of each lineage (Figure 7B), sequences above 1 and 2 standard\ndeviations from the mean dynamic semantic score at each 3-month time point are classified\nas sequences of concern. When we filter for the first occurrences of known VOCs and\nimportant variants (Figure 7C), we see that the majority (59%) are above one of the warning\nlines. These include all WHO-assigned VOC sequences (Alpha, Beta, Gamma, Delta and\nOmicron), with 4/5 in the red warning section.”\n“Using the sequences above the outlier warning lines, we counted the number of unique\nmutations in the spike protein that occur within those sequences (Figure 7D). The NTD and\nRBD contain the largest number of mutated sites, and often the sites with the greatest\nquantity and diversity of changes at a single position. Of the mutations identified, 43% of the\npositions were also present in the first sequence of an identified major variant sequence,\nwith 37% of the exact mutations present.”\nFigure 7. (A) UK SARS-CoV-2 spike sequences through the pandemic. Each point represents a\nsequence cluster with 99.9% sequence similarity. Dynamic semantic scores were calculated for each\nsequence cluster, with the black line showing the mean sliding score. (B) Haplotypes are filtered on first\noccurrence of a new Pango lineage with standard deviations computed for each 3-month window. The\norange warning area denotes sequences above 1 deviation from the mean, while the red warning area\ndenotes >2 standard deviations. (C) Haplotypes filtered on the first appearance of a VOC/major variant,\nwith deviation thresholds showing if they were detected on first appearance. (D) First haplotype\nsequence site mutation frequencies for sequences identified by warning thresholds mapped onto the\nspike monomer.\nMinors:\n- The paper lacks clarity on which score - grammaticality or semantic - should be prioritized\nfor predictive purposes. Resolving this confusion and providing a clearer explanation of their\nrespective roles in prediction would enhance the overall impact of the work.\nWe have added much more thorough benchmarking of each method to the paper to outline\nthe pros and cons of each metric. For predictive tasks for single mutations, we see that the\ngrammaticality scores are typically more useful. This is likely because large semantic scores\noften have big effects on the protein, and as such are not seen often. The grammaticalities\nalso allow for a measure of epistasis, not quite possible with the semantic score. However,\nthe semantic score is crucial for the dynamic embedding approach, since it allows for\ncomparisons between the embeddings. As such, one is not necessarily better than the other,\nbut there are situations where one might pick one score over the other.\n- There is no title for the “Conclusion” section.\nA title has now been added to the conclusion.\nReviewer #4 (Remarks to the Author):\nI co-reviewed this manuscript with one of the reviewers who provided the listed reports. This\nis part of the Nature Communications initiative to facilitate training in peer review and to\nprovide appropriate recognition for Early Career Researchers who co-review manuscripts.\n2nd response to reviewer comments, NCOMMS-24-59219-T\nWe are grateful to the reviewers for their time and effort in assessing our paper, and for their\ninsightful and constructive comments and suggestions. See our detailed responses\nhighlighted in blue below:\nReviewer #1 (Remarks to the Author):\nThank you for the revised manuscript. While some minor issues have been addressed,\nseveral fundamental concerns regarding the novelty of the proposed methods and the\nquantitative validation of the results remain. My specific comments are as follows:\n1. On the Use of Single Sequences and ESM-2:\na. Authors' statement: \"We focus on what can be accomplished with a single\nsequence and use of an available existing model, ESM-2, eg, following the\nspillover of a novel virus with limited data... Almost every other tool requires\nmore than just a single viral sequence... which makes using PLMs an\nattractive proposition.\nb. Comment: While the scenario described (limited data, single sequences for\ninference) is highly relevant in emerging pathogen contexts, the claim that the\nframework's effectiveness and novelty stem from the use of a pre-trained\nmodel like ESM-2 in this inference mode is problematic. The ability to perform\ninference on a single sequence after a model has been trained is a common\ncharacteristic of many models, including multi-sequence alignment-based\nmethods once trained (e.g., EVEscape performs inference on single\nsequences or new variants based on its training). The core contribution's\nnovelty must therefore lie specifically in the analytical methods\ndeveloped or applied using the single-sequence ESM-2 output, rather\nthan merely the choice of a pre-trained model or its application for\ninference.\nWhile the reviewer’s claim that models like EVEscape also do inference on a single\nsequence is true, that inference must be made from a model that has specifically been\ntrained on sequences relevant to the query sequence. This means that in the case of a novel\nvirus, EVEscape would have to be trained on related viral sequences in order to be used.\nYou cannot, for instance, use the influenza version of the EVE model on SARS-CoV-2\nsequences; these are 2 separate models.\nThis is not, however, the case for the use of ESM-2 as we demonstrate here, which is why\nusing a general-purpose PLM is attractive. Our paper details how you might use such a\nmodel in the context where it is of most value, i.e., where you do not have immediate access\nto an alignment or where you would have to wait for another model to be trained. This is\nwhat we mean by saying you can analyse a single viral sequence.\n2. On the Novelty of Epistatic Interaction Detection:\na. Authors' statement: \"We proposed a novel single-sequenced epistatic interaction\ndetection analysis and validated it with real experimental data.\"\nb. Comment: The proposed method for detecting epistatic interactions appears to be\nhighly similar, if not a direct adaptation, of the approach presented in Zhang et al. (PNAS\n2024, DOI: 10.1073/pnas.2406285121). Both methods utilize systematic perturbation\n(mutation) of the input sequence and compute changes in model output (likelihood\nchange/Jacobian) to infer interactions. While this paper applies the concept specifically to\nSARS-CoV-2, the underlying methodology does not seem conceptually distinct or novel\ncompared to the more general framework described by Zhang et al. To claim novelty, the\nauthors would need to clearly articulate and prove what is fundamentally different or\nimproved in their methodological approach compared to this prior work.\nSincere apologies for this oversight. We have now cited the work of Zhang et al, and added\nthe following text on page 6:\n“A method by Zhang et al. uses a similar approach to predict contact maps, however here\nwe want to predict mutation specific effects instead of site specific contacting residues.”\nWe had used their method in our previous response to reviewers' document, as we had\nused it to produce the SARS-CoV-2 contact map, as it shows how direct contacts can be\npredicted using ESM-2. Zhang’s method, while related to ours, differs significantly in intent:\nZhang et al. are looking specifically for contacts between sites. Their signal comes from\nmutating to every possible amino acid at a site, and doing this for every site in the sequence.\nThe Jacobian is then calculated from this set of measurements, giving you a site-based\ncontact map. Our approach simply mutates the site to the mutation of interest and identifies\npotential epistatic interactions from the perturbations from this single mutation. For\nidentifying direct contacts, it would be likely that Zhang et al. will perform better, as the\ncontact signal between sites will likely be better produced by observing effects from multiple\ndifferent mutations at a single site. However, as we mentioned, we are inferring mutation-\nspecific effects, not site effects.\n3. On the Novelty of the Horizon-Scanning Score:\na. Authors' statement: \"We proposed a novel horizon-scanning score that is agnostic\nto any experiment-derived features...\"\nb. Comment: The concept of developing scores to assess the potential risk or\nevolutionary trajectory of new variants (\"horizon scanning\") is not novel. Tools like\nEVEscape already provide regular, comprehensive screenings of novel variants on a weekly\nbasis. For the proposed score to be considered novel and a significant contribution, the\nauthors need to clearly define what specifically is novel and why the score is better\nabout the score's mathematical definition, the methodology used to derive it from the\nPLM outputs, or the biological/evolutionary property it uniquely captures, beyond\nsimply being derived from ESM-2 and being \"agnostic\" to specific experimental features\n(which is also a characteristic of many sequence-based or PLM-based scores).\nWhile the concept of horizon scanning is not novel, we believe the novelty arises from using\na sliding window with an embedding approach to determine risk in a way that moves away\nfrom referring back to the initial SARS-CoV-2 sequence. As mentioned, EVEscape do\nprovide screenings of new variants. This is not what our method does. Our method does not\nproduce a new deep mutational scan for each new emerging variant, although this could be\ndone with ESM-2. It creates a 3-month sliding window mean embedding to measure from,\nwhich appears to be a better proxy for local divergence than the original SARS-CoV-2\nsequence.\n4. On the Comparison to Hie et al. and Grammar/Semantic Scores:\na. Authors' statement: \"Hie et al... proposed that grammar and semantic scores are\nanalogous to protein fitness and antigenicity. We show that these metrics do not correspond\nto these properties when applied to a generic PLM, such as ESM-2... We assess the utility of\nthese metrics and identify what they capture given a more foundational PLM, ESM-2. We\nthen propose approaches that can be used when more sequence data is available and show\nthat sequence logits potentially identify epistatic sites within the protein, a technique not\nmentioned in Hie et al.\"\nb. Comment: The analysis applying concepts from Hie et al. (using PLM-derived scores\nas proxies for biological properties) to a different pre-trained model (ESM-2 instead of Hie's\nLSTMs) is presented. However, simply re-applying a previously proposed concept with a\ndifferent existing model does not inherently constitute significant scientific novelty, unless\nthis application reveals profound, previously unknown properties or limitations of\nfoundational models like ESM-2 regarding their representation of protein evolution,\nsupported by rigorous analysis. I did not see any more novelty from the current format.\nRespectfully we disagree. We have a significantly different scope of the utilities of the scores\ncompared to Hie et al. We focus on in silico mutational scan, epistatic interactions and\nhorizon scanning. Notably, we are interested in predicting variant effects more broadly (and\nusing the foundation model ESM-2) not just antigenicity, and we benchmarked our\npredictions against evolutionary, structural metrics and experimental DMS data.\nThe LSTM in Hie et al and ESM-2 differ significantly in model capacities. They also differ\nsignificantly in terms of the scope and scale of their training data. We assessed the\ncapabilities and limitations of ESM-2.\nImportantly, our analysis mimics the situation at the beginning of the pandemic, when only a\nsmall number of sequences were available, as we systematically tested the performance of\nESM-2 with single sequences.\nThese analyses and insights extend the work of Hie et al., confirming their utility applied to\nvirus data. These models are being increasingly used to analyse viral sequence data, and it\nis important to identify where they can be used, as well as to develop new methods that take\nadvantage of their unique benefits (alignment-free, single sequence inference etc).\n5. On Epistatic Interaction Benchmarking:\na. Authors' statement: \"We agree with the reviewer that additional benchmarking is\nrequired for predicting epistatic interactions. However, the GB1 analysis and experiment was\nnot suitable... We did investigate producing a contact map for GB1... We found that ESM-2\ncan be used to identify contacts, as shown below. However for inclusion in our paper, we\ndecided to focus on a more directly relevant SARS-CoV-2 dataset.\"\nb. Comment: I appreciate the authors' acknowledgment that additional benchmarking is\nneeded. The brief visual demonstration of GB1 that their method can identify contacts is not\na substitute for rigorous, quantitative evaluation. For a computational methods paper,\nvisualizations and anecdotal examples are useful for illustration but cannot objectively\ndemonstrate performance, quantify accuracy (e.g., precision, recall, correlation), or allow\ncomparison to other methods. It is imperative that the authors include a comprehensive\nquantitative evaluation of their epistatic prediction method on a relevant dataset (such\nas the SARS-CoV-2 data they mention) using appropriate statistical metrics to support any\nclaims about its efficacy.\nWe have addressed this with the inclusion of specificity and sensitivity comparisons on the\nBA.1 epistatic dataset with a competing mutual information method, as well as using a\nhypergeometric test to check for enrichment of mutations at important and experimentally\nconfirmed positions within the RBD. We have included the following text on pages 9 and 10.\n“Using a hypergeometric test, the stringent threshold was found to be 6.43 fold over enriched\n(p-value 1.462928853945546e-10) with RBD-RBD BA.1 mutation interactions. This was also\ntrue for strong RBD-RBD mutation epistatic interactions (an over enrichment of 8.21 fold, p-\nvalue 0.005625658517307826) as well as for all BA.1 mutation-mutation interactions (1.6\nfold, p-value 0.01687039075859053). Strong epistatic interactions were defined as 1\nstandard deviation from the mean absolute experimental epistatic interaction. These\nremained significant after applying a Benjamini-Hochberg multiple test correction. ”\n“We also tested the sensitivity and specificity of the model in comparison to the mutual\ninformation based approach from Innocenti et al.30 and found that when looking within the\nexperimental data for epistatic sites, ESM-2 (stringent 2 standard deviation cut-off) had an\nequivalent sensitivity (0.120 vs 0.114) and improved specificity (0.842 vs 0.777). This was\ncalculated using the strongly epistatic set of RBD-RBD mutation interactions. While the\nspecificity is good for both approaches, suggesting both methods work well for ruling out\nfalse positive epistatic interactions (or in this instance, weak interactions), the low sensitivity\nsuggests that the task of identifying true/strong interactions remains challenging.\nConsidering ESM-2s lack of pandemic sequence data, future work could be directed towards\nleveraging pandemic data to improve the predictive performance, although extensive\nvalidation/training experimental pairwise datasets may prove to be a bottleneck here. Also, it\nshould be again noted that ESM-2 is not trained explicitly to predict effects on ACE-2\nbinding, meaning the epistatic signal may be from another biologically meaningful measure,\nlike the 371, 373 and 375 mutations.”\n6. On the BA.1 Reversion Data Comparison:\na. Authors' statement: \"...we used the pairwise effects to compare to our BA.1\nreversion data and show that several of the top pairwise epistatic interactions between BA.1\nmutations are shown by our method. We report these findings in the updated Figure 3. We\nalso discuss the caveats... We believe our approach here is interesting given its lack of\nknowledge about sequences spreading in the pandemic, and its relative ease of use...\"\nb. Comment: Regarding the comparison to the BA.1 reversion data: While\ndemonstrating that \"several of the top pairwise epistatic interactions... are shown by our\nmethod\" provides some illustrative examples, this remains anecdotal evidence. Pointing to a\nfew matching interactions, even the top ones, does not constitute a robust quantitative\nvalidation of the method's overall performance. The heatmap visualization (Figure 3) itself\nappears to show a significant number of predicted interactions that do not correspond to the\nexperimental data (potential false positives). To support claims about the method's ability\nto identify epistatic effects, the authors must provide quantitative metrics (e.g.,\ncorrelation coefficient with experimental values, precision/recall if interactions are\ntreated as binary, or other relevant statistics) to objectively assess the agreement\nbetween their predictions and the experimental data. Claims about the method's \"interesting\"\nnature or \"ease of use\" cannot substitute for demonstrated quantitative accuracy.\nComment addressed above.\nReviewer #2:\n(This review was also done by reviewer #3).\nReviewer #2 had an overall comment regarding the lack of validation and benchmarking in\nthe previous draft, a point with which I agreed. I believe the new sections added by the\nauthors provide a satisfactory response to these concerns, as discussed above.\n1. ESM-2 training data: Reviewer #2 raised a valid point regarding the training set of\nESM. The inclusion of SARS-CoV-2 genomic data in the training set could undermine\nthe ability to assess the true predictive power of the proposed approach, as it\neffectively prevents a clear separation between training and test data. The authors’\nresponse is satisfactory; however, as mentioned above, I believe they should more\nclearly specify the differences between the training sets used across the various\ntools. ESM-2 is pretrained on the UniRef90 database, which includes a broad range\nof protein sequences from bacteria, eukaryotes, and other organisms - not just viral\nproteins. It would be interesting to explore how such a model would perform if trained\nexclusively on viral sequences. While I understand that this is beyond the scope\nof the current work, I would appreciate it if the authors could briefly comment\non this.\nWe have added to our Discussion on page 22:\n“A major strength of ESM-2 is its ability to learn broader evolutionary signals, which stems\nfrom its training data that spans the tree of life. While models trained from exclusively viral\nsequences may be produced in future, the sparsity of viral data include in training remains a\nmajor issue. Viral datasets currently contain only a small set of available and sufficiently\ndistinct viral sequences. This is likely to limit the quality of any model inferences, however\nfinetuning on such a set of sequences is more likely to yield improved performance. Models\nsuch as CovFit and LucaVirus have already demonstrated the power of finetuning general\nPLMs (LucaOne and ESM-2) on viral datasets, with more models likely to be released in\nfuture.”\n2. ESM2 capture evolutionary potential and “new” metrics: Authors now compare ESM-\n2 scores to DMS data from Dadonaite et al. and show ACE2 binding is significantly different\nbetween S1 and S2. I think the use of such DMS data plus additional experimental data\nstrengthens the authors’ claim.\nWe thank the reviewer for recognising this contribution.\nReviewer #3 (Remarks to the Author):\nI appreciate the effort made to address my previous comments. The revised manuscript\nshows significant improvement, particularly through the inclusion of comparisons to\nexperimental data, as well as benchmarking against other computational methods. Although\nthe overall correlation with experimental data remains low, and it is still uncertain whether\nsuch an approach could be applied in future pandemics, I believe the paper provides\nvaluable insights into how ESM2 can be used to detect epistasis and forecast mutations.\n1. Epistatic interactions: I appreciate that the authors incorporated experimental\nepistatic data from Moulana et al. to compare their computational predictions with\nmeasured ACE2-binding affinities, addressing Reviewer 1’s comments. I am, in fact,\nsomewhat surprised that the model is able to capture such effects at all, especially\nconsidering that, to my understanding, the ESM2 model is trained on individual\nprotein sequences rather than on interacting protein pairs.\nWe thank the reviewer for their comment. ESM-2 is indeed trained on individual proteins\nrather than pairs, however due to the structural and functional conservation that is learned\nby the model, it appears to detect interactions between sites on the spike protein that impact\nits interaction with ACE-2, despite having no direct knowledge of it.\n2. The authors rightly note that the ESM2 model is not directly linked to ACE2\nbinding. I think it is important to emphasize this point more explicitly: a priori, models\ntrained on evolutionary data from single proteins are not expected to capture effects related\nto protein-protein interactions - particularly those specific to the Spike protein - ACE2\ninterface. I would also like to thank the authors for including the comparison with Innocenti et\nal., which in my opinion provides a helpful benchmark against another computational method\nthat relies only on SARS-CoV-2-specific data.\nWe agree and have added a clarification to our main text on page 10 on this point:\n“Also, it should be again noted that ESM-2 is not trained explicitly to predict effects on ACE-\n2 binding, meaning the epistatic signal may be from another biologically meaningful\nmeasure, like the 371, 373 and 375 mutations.“\n3. Benchmarking with EVEscape, IND, and DCA scores: I think the inclusion of this\nbenchmarking analysis is valuable - it helps assess both the validity and limitations of the\nauthors' approach. That said, I have a few comments that could help clarify how these\ndifferent methods operate:\na. The IND (Independent) and DCA (Direct coupling analysis) models are both trained\nwith pre-pandemic MSAs of viral sequences (Coronaviridae). The IND is based on single-\nsite amino acid frequencies and identifies positions more likely to mutate based on\nconservation. DCA goes a step further by incorporating pairwise correlations between sites,\npotentially capturing signals of epistasis.\nb. A related approach was later introduced by EVEfitness, which also uses pre-\npandemic MSAs of viral sequences as input to train a variational autoencoder. In this sense,\nEVEfitness generalizes IND and DCA by replacing single-site amino acid frequencies and\npairwise correlations statistics with an autoencoder model. EVEscape builds upon\nEVEfitness by incorporating additional biological features - such as epitope accessibility and\nantibody binding propensity - to refine predictions of viral evolvability.\nc. ESM-2 is a transformer-based protein language model trained on millions of protein\nsequences—not limited to viral proteins.\nI think this section would benefit from a rewrite to more clearly explain each\napproach and the type of data each model is trained on.\nThank you for the constructive feedback, we have taken this on board and expanded on the\ndescriptions of both methods, and how they have been trained. The following text was added\non page 11:\n“We benchmarked our results against two other computational scores: Evescape from\nThadani et al and mutability scores from Rodriguez-Rivas et al. Evescape is a variational\nautoencoder based on EVE that was trained using pre-pandemic spike protein sequences.\nThese prior sequences are given to EVEscape as a one-hot encoded alignment, which the\nmodel uses to learn distributions and dependencies across sites in the protein. The\nvariational autoencoder produces a mutation score that when compared with the score of the\nreference sequence approximates the relative fitness of the mutant, termed Evescape\nFitness. This is combined with a structure derived contact accessibility and amino acid\ndissimilarity metric to produce the full Evescape score. These additional metrics aim to\nweight the model towards detecting more immune-evasive mutations.\nThere are two mutability score variants, the direct coupling analysis (DCA) score and the\nIndependent Site Model (IND) score. Both approaches use existing protein sequence\nalignments to calculate site-wise frequencies across the protein sequence. Much like the\nEVEscape fitness, the IND score is a likelihood ratio between the reference and mutated\nsequence. The DCA additionally incorporates coupling analysis to try to incorporate signals\nof mutational epistasis. This involves adding coupling terms between sites in order to\napproximate the epistatic interactions between them.\nWe also tested other methods for computing grammaticality: masked grammaticality\nsuggested by Allman et al.; mutated grammaticality; the relative grammaticality measure\nsupplemented with the dissimilarity and accessibility as shown in Thadani et al.”\n4. Moreover, the results presented in Supplementary Figure 10 are somewhat\nsurprising. Specifically:\na. EVEfitness appears to outperform EVEscape, and\nb. DCA and the IND model achieve equal performances.\nThese findings are unexpected. DCA - by incorporating epistatic interactions - has\ngenerally been shown to outperform the frequency-based IND model when it comes\nto predicting evolutionary trajectories. Likewise, In the original EVEscape paper,\nEVEscape consistently outperforms the simpler EVEfitness model. Could these\ndiscrepancies be due to the specific threshold used (e.g., selecting the top 10% most\nfrequently observed mutations)? Have the authors tested whether these rankings are\nconsistent across different thresholds or evaluation criteria? I would appreciate it if\nthe authors could comment on that.\nThe EVEscape fitness components improved performance in identifying mutations present in\nthe population is actually noted by the authors in the original paper. The authors mention\nthat the increase in performance comes from the identification of low-frequency mutations,\nand that while these mutations lack strong selective pressure to be maintained in the\npopulation, they don’t impair function.\nThe full EVEscape score is weighted towards sites that are accessible and amino acids that\nare more dissimilar. This weighting aims to skew EVEscapes preference towards immune\nevasive mutations, which means low frequency mutations with minimal impact tend rank\nlower when using this score. Conversely, this makes mutations that are impactful, but less\nlikely to occur rank more highly, making EVEscape worse at this task than its fitness\ncomponent. However, this may be why the EVEscape score performs better at identifying\nVOC mutations, as shown in Supplementary Figure 10B.\nDCA performs similarly, detecting more VOC mutations than IND while performing slightly\nworse on the frequency task. We did have a look at a 25% threshold for mutations, and DCA\ndoes then perform slightly but consistently better than IND in the frequency task. However,\nthe patterns for EVEscape and EVEscape fitness are maintained.\nIt should also be noted that in our benchmarks, we assessed each method on the set of\ncommon mutations between them. This was to ensure that methods were not doing better\ndue to numerical features i.e. a larger set of mutations comprising the top 10%, or in the\ncase of experimental measures, enrichment of RBD mutations when comparing across the\nwhole of the protein.\n5. Moreover, I think the authors should add a Methods section clarifying where the\nbenchmarking data were downloaded from, in order to improve reproducibility.\nA benchmarking section with links to the files used for each has been added and is as\nfollows:\n“The EVEscape benchmarking data (Supplementary Figure 10) were collected from the\nGitHub repository (commit 8238e4f,\nhttps://github.com/OATML-\nMarkslab/EVEscape/blob/main/results/summaries_with_scores/full_spike_evescape.csv,\nhttps://github.com/OATML-\nMarkslab/EVEscape/blob/main/results/summaries_with_gisaid/spike_dist_one_scores_gisai\nd.csv and https://github.com/OATML-\nMarkslab/EVEscape/blob/main/data/gisaid/single_mutant_count_by_month.csv ).\nData from DCA and IND mutability scores were collected from the GitHub repository (commit\naeffe23, https://github.com/GiancarloCroce/DCA_SARS-CoV-\n2/blob/main/data/data_dca_proteome.csv).\nSupplementary Figure 10A uses variants mentioned in the https://github.com/OATML-\nMarkslab/EVEscape/blob/main/results/summaries_with_gisaid/spike_dist_one_scores_gisai\nd.csv file, but calculates the mutations from a representative set of sequences from GISAID\n(https://doi.org/10.55876/gis8.240620pm) so that all mutations (not just single nucleotide\nchanges) could be identified.\nSupplementary Figure 10B uses the data from https://github.com/OATML-\nMarkslab/EVEscape/blob/main/data/gisaid/single_mutant_count_by_month.csv , since this\nincludes mutation frequency tracking for each month. This is restricted to just single\nnucleotide changes mutations.\nEpistasis benchmarking data from Innocenti et al.28 was downloaded from the paper\nsupplementary table 1 https://static-content.springer.com/esm/art%3A10.1186%2Fs13059-\n024-03355-y/MediaObjects/13059_2024_3355_MOESM1_ESM.xlsx .”\n6. Predicting emerging VOC strains: The approach presented in Figure 7 for identifying\nVOCs and Spike protein positions to monitor is compelling and suggests a potential use of\nESM for forecasting mutations - a point that, in my opinion, was not clearly conveyed in the\nprevious draft.\nWe thank the reviewer for recognising this contribution.\nReviewer #4 (Remarks to the Author):\nN/A\nReviewer #5 (Remarks to the Author):\nI have had the chance to carefully consider the authors' manuscript, first round of reviewer\nfeedback, the authors' rebuttal, and the second round of comments from reviewer 1. Overall,\nI disagree with reviewer 1 on points of novelty, but I do think some very basic statistical\nquantification should be provided, but this is very reasonable analysis to include in the\nmanuscript and would merit an additional round of revision before I can personally\nrecommend acceptance of the paper.\n1. Regarding the novelty of single-sequence models such as ESM-2, \"epistatic\"\ninteractions, horizon scanning, and grammaticality/semantic scores, while these have\nbeen previously applied in other contexts, the investigation of these scores in the\ncontext of SARS-CoV-2 after several years of evolution and extensive experimental\ndata characterization, is still interesting and worth communicating via a publication.\nAs long as the authors cite the appropriate papers, this should be fine. I do agree that\nthe categorical Jacobian method of Zhang et al., PNAS, 2024 needs a citation here.\nAlso, while both the reviewer and the authors have been referring to this as\n\"epistasis\", many of the contacts recovered from this analysis are of structural\nproximal residues and the authors should just be careful not to describe all of the\ninteractions predicted by this approach as \"epistatic.\" This can be addressed with\nchanges to the manuscript text.\nWe have added a reference for the Zhang et al. paper\n“A method by Zhang et al. uses a similar approach to predict contact maps, however here\nwe want to look at mutation specific effects instead of site specific contacting residues.”.\nFor the definition of epistasis, we do identify structurally proximal residues and include these\nin our definition of epistasis. To make it clear what definition we’re using we added a citation\nto a review paper from Harms and Thornton (https://www.nature.com/articles/nrg3540) which\ndiscusses several different types of epistasis: both structurally proximal interactions as well\nas longer-ranged indirect mechanisms included:\n“Such epistatic effects can be due to direct or indirect mechanisms linked to residues in\nclose contact in the protein structure conformation and/or protein stability (Harms and\nThornton, 2013.”\n2. Regarding the epistatic interaction benchmarking, I am not sure that the GB1\nanalysis suggested by the reviewer in the original referee report is that meaningful, and I\nagree with the authors that it is sufficient to simply cite prior work showing that ESM-2 can\nlearn 3d contacts for general proteins. However, I do agree with the reviewer that a\nquantitative analysis of this data should be added. This could look like a very simple and\nreasonable analysis showing statistical significance of the findings, as well as a\ncomparison to a non-neural and non-ML baseline -- for example, 3d distance correlates\nstrongly with many putative epistasis/allostery predictors.\nSee response above where this comment was addressed (Reviewer 1, Comment 5).\n3. Regarding the BA.1 reversion analysis, which is related to the above comment, I\nagree with the reviewer that the analysis could simply use a statistical test, where a null\ndistribution could potentially be generated with random permutations, alongside simple non-\nneural/non-ML baselines.\nSee response above where this comment was addressed (Reviewer 1, Comment 5)."}, {"round": 2, "reviewer_comments": {"reviewer_5": "The authors have adequately addressed my concerns and I have no further comments. I am happy to recommend publication of this manuscript."}, "author_response": ""}]

Decision letter