Predicting a Protein's Stability under a Million Mutations
Authors: Jeffrey Ouyang-Zhang
Daniel Jesus Diaz
Adam Klivans
Philipp Kraehenbuehl
Venue: NeurIPS 2023 poster (2023)
Keywords: stability
proteins
biology
physical
TL;DR: parallel decoding the change in thermodynamic stability under all single and double mutations

Abstract: Stabilizing proteins is a foundational step in protein engineering. However, the evolutionary pressure of all extant proteins makes identifying the scarce number of mutations that will improve thermodynamic stability challenging. 
Deep learning has recently emerged as a powerful tool for identifying promising mutations.
Existing approaches, however, are computationally expensive, as the number of model inferences scales with the number of mutations queried. 
Our main contribution is a simple, parallel decoding algorithm.
Mutate Everything is capable of predicting the effect of all single and double mutations in one forward pass. 
It is even versatile enough to predict higher-order mutations with minimal computational overhead.
We build Mutate Everything on top of ESM2 and AlphaFold, neither of which were trained to predict thermodynamic stability.
We trained on the Mega-Scale cDNA proteolysis dataset and achieved state-of-the-art performance on single and higher-order mutations on S669, ProTherm, and ProteinGym datasets.
Our code is available at https://github.com/jozhang97/MutateEverything.

Predicting a Protein's Stability under a Million Mutations
Jeffrey Ouyang-Zhang
Daniel Diaz UT Austin danny wwwzuutexas edu
UT Austin jozhangdutexas. edu
Adam R. Klivans
Philipp Krähenbühl
UT Austin
UT Austin philkr www.uueexeeeco edu
klivans@cs. wtw.us.eeeee...e edu
Abstract
Stabilizing proteins is foundational step in protein engineering. However, the evolutionary pressure of all extant proteins makes identifying the scarce number of mutations that will improve thermodynamic stability challenging. Deep learning has recently emerged as powerful tool for identifying promising mutations. Exist- ing approaches, however, are computationally expensive, as the number model inferences scales with the number of mutations queried. Our main contribution is a simple, parallel decoding algorithm. Our Mutate Everything is capable of pre- dicting the effect of all single and double mutations in one forward pass. It is even versatile enough to predict higher-order mutations with minimal computational overhead. We build Mutate Everything top of ESM2 and AlphaFold, neither which were trained to predict thermodynamic stability. We trained on the Mega- Scale cDNA proteolysis dataset and achieved sate-of-the-are performance on single and higher-order mutations on S669, ProTherm, and ProteinGym datasets. Our code is available https: //githut.b../ com/ wwwhangg//uutaeeeeeeeeeeeeeeee Butteetrrrtterrryttnnngggg
1 Introduction
Protein engineering is the discipline of mutating a natural protein sequence to improve properties for industrial [5, 78] and pharmaceutical applications [2, 25, 41]. However, evolution simultaneously optimizes several properties of a protein within its native environment, resulting in proteins with marginal thermodynamic stability (~ 5-15 kcal/mol) [39] which become non-functional in in- dustrial setting. Therefore, repurposing natural protein for biotechnological applications usually begins with identifying non-deleterious mutations that stabilize the structure. With stabilized structure, downstream engineering goals, which often require exploring destabilizing mutations, become tractable. Historically, this process has been bottlenecked by the requirements of extensive laboratory characterization of rational designs [33, 79] or directed evolution libraries [4, 12, 21, 22]. The recent explosion biological data [9, 29, 72, 73] has enabled deep learning frameworks to accelerate the identification of stabilizing mutations. A successfully stabilized protein often requires several mutations. However, current frameworks do not take into account epistatic interactions between multiple mutations. Thus, to significantly accelerate protein engincering, it is critical to efficiently navigate the epistatic landscape of higher-order mutations [10, 69). However, due to its combinatorial nature, thorough exploration quickly becomes computationally prohibitive.
In this paper, we introduce Mutate Everything to directly predict changes thermodynamic stability (AAG) for all single and higher-order mutations jointly. Mutate Everything is a parallel decoding algorithm that works in conjunction with a sequence model for predicting thermodynamic stability.
37th Conference on Neural Information Processing Systems (NcurIPS 2023).

---

K
P1K
A AKVMGH... -1.0 kcal/mol
P
E Figure 1: Mutate Everything efficiently predicts AAG, the change in thermodynamic stability of folding, for over a million mutations (e.g. all single, double mutations) in single inference step. This helps identify and prioritize stabilizing mutations (AAG < 0) in protein engineering efforts. The notation for a mutation AA from. pos, AAto (e.g. "PIK" mutates from K at position 1)
Mutate
K3E
PATVMGH... -0.5 kcal/mol
Everything
K
K W
M5K H7W
PAKVMGH...
PAKV t... +1.0 kcal/mol M H
Prior models for predicting thermodynamic stability work by producing an embedding or representa- tion of the input sequence. Our decoder takes this representation and uses linear network create further representations, one for every possible single mutation. These mutation-level representa- tions are then further aggregated form representations for higher-order mutations. We feed these higher-order representations into a lightweight multi-layer perception (MLP) head to output predicted AAG measurements. Since the mutation-level representations are computed only once, we are able to scale to millions of (higher-order) mutations, as the aggregation and MLP-heal computations are inexpensive. As such, given a fixed computational budget, our model evaluates millions more potential mutations than prior methods.
Mutate Everything estimates the effects of all single and double mutations for single protein in seconds on one GPU. This can efficiently compute the stability of all single and double mutants across all 20,000 proteins in the human proteome within short time frame. To the best of our knowledge, Mutate Everything is the first tool that renders the computational analysis of double mutants across the entire proteome tractable. Mutate Everything can be used in conjunction with any model that generates representations for protein sequences. In this paper, we use AlphaFold [31] as a backbone architecture and fine-tune it on a dataset labeled by experimentally derived AAG measurements. This results in the first accurate model of stability prediction based on AlphaFold.
We evaluate the performance of our Mutate Everything on well-established benchmarks for predict- ing experimentalll-vaaidaae ProteinGym [49]. On Pro'Therm high-order mutations (PTMul), our model achieves Spearman of 0.53, compared to 0.50 in the next best method. On S669 single mutations, our model achieves a Spearman of 0.56, outperforming the -fate-of-the-art at 0.53. On ProteinGym, the Mutate Every- thing outperforms state-of-the- thetth methods from 0.48 to 0.49. Our model makes finding the most stabilizing double mutations computationally tractable. We show on the double mutant subset of the cDNA- prote-pots dataset (cDNA2) [72], where only 3% of known mutations are stabilizing, that Mutate Everything effectively ranks stabilizing mutations ahead of destabilizing ones.
AAG measurements of mutations: ProTherm [48], S669 [54], and
2 Related Works
2.1 Protein Engineering
Protein engineering is the process of mutating a natural protein to improve particular phenotypes, such as thermodynamic stability and function. The field has historically relied on rational design and stochastic methods, such as error-prone PCR [1]. DNA shuffling [71], and directed evolution (DE) [4, 27], to identify gain-of-function function mutations. Rational design is limited to proteins with solved structures and requires an extensive understanding of biochemistry and specific knowledge of the particular protein to select mutations highly likely to improve the target phenotypes [33, 79J. Directed evolution requires little to no knowledge of the protein and instead generates a library of protein variants that are then screened for the target phenotypes (c.g. fluorescence brightness, antibiotic
2

---

resistance, stability, activity) [4, 27]. The library can be generated via site-saturated mutagenesis of handful of positions in the sequence [12] or DNA shuffling, in which the gene is digested into random fragments and reassembled into full-length sequences [71]. After screening the library, the most "fit" variant is then selected as the initial sequence for the next round of directed evolution. This iterative process repeats until it obtains a protein variant with the desired phenotypes.
Machine learning has demonstrated its ability augment rational design and accelerate the stabi- lization of variety of proteins [20, 28, 40, 52, 68]. Separately, machine learning-guided directed evolution (MLDE)[75] has been shown to improve the likelihood of obtaining the global fitness maximum by 81-fold compared to traditional DE [76]. MLDE has accelerated the engineering of several proteins, such as the enantioselectivity of enzymes for kinetic resolution of epoxides [23] and the activity and expression of a glutathione transferase [46]. Mutate Everything empowers the experimentalist to accelerate the stahilization of protein for both rational design and ML.DF.
2.2 Machine Learning for Protein Structure Prediction
Recent advances in machine learning have led to remarkable progress in protein structure prediction. AlphaFold [31] has demonstrated that deep learning is highly effective at predicting protein structures from a sequence by using evolutionary history via a multiple sequence alignment (MSA). AlphaFold passes the MSA and pairwise representation of the sequence into Evoformer to capture the co- evolutionary patterns between residues. The Evoformer output is then processed by the Structure Module, which predicts the protein's structure. We challenge prior works postulating that AlphaFold cannot be used for stability prediction and show that fine-tuning these rich co-evolutionary and structural features yield highly performant stability predictors [53].
Evolutionary Scale Modeling (ESM) [38, 62] has shown that protein structure prediction can be performed without MSAs and specialized architectures by leveraging large transformers pre-trained on masked token prediction. Other works extended this masked pre-training framework, including MSA-Transfon Transformer [60] which incorporates sequence's MSA as input, and Tranception [49], which develops a hybrid convolution and attention based autoregressive architecture. We show that fine- tuning these evolutionary representations excels at protein stability assessment without MSAs.
2.3 Protein Stability Assessment
Traditional approaches to protein stability assessment relied on a combination of physics-based methods, statistical analysis, and traditional machine learning techniques. Physics-based tools, such as FoldX, Rosetta, and SDM, utilize energy functions and statistical patterns to assess how mutations affect a protein's 's stability [32, 67, 77]. DDGun [44] directly inferred AAG from heuristics, including Blosum substitution scores, differences in interaction energy with neighboring residues, and changes in hydrophobicity. Many traditional machine learning approaches used support vector machines and decision trees with physics-based baseed feature engineering [11, 14, 15, 17, 37, 66). Others ensemble existing machine learning methods [35, 56, 59, 63, 64].
Recently, deep learning-basee based approaches have recently begun to outperform existing physics and traditional machine learning approaches [13, 19, 34]. ACDC-NN [6] trains an asymmetric convolu- tional neural network for predicting forward and reverse mutation effects. Thermonet [36] voxelizes and feeds both the wild-type and mutant protein structures into a 3D-CNN to regress a AAG value. PROSTATA feeds the wild-type and mutant protein sequence into a pre-trained ESM2 model and then regresses AAG value [74]. Stability Oracle takes the native and mutant amino acid type along with the local structure as input and regresses A AG [18]. Several deep learning frameworks [64, 80, 81] model multiple mutations in addition to single mutations. In this paper, we develop framework that models the protein to enable efficient enumeration of all mutation candidates.
3 Preliminary
Problem Setup. A protein w (wi...., wr) is sequence of amino acids wi E AA, where AA {A, C)..,,Y are the 20 different amino acid types encoded by genetic information. Let μ == (p, t) denote mutation substituting the amino acid wp at position p E [1,..., L] to amino acid type / E AA. Our goal is to determine the change in thermodynamic stability AAG E R for a protein under large number of mutation sets [M,..., MN) where mutation set M, Kk
3

---

MLP[V]
AAG
MLP[V]
for P1V
AAG for
MLP[E]
AAG for V4E
MLP[E]
+
MLP
P1V
V4E
AlphaFold
H7K
AAG for H7K
MLP[K]
MLP[K]
(a) Backbone
(b) Mutate Single
(c) Mutate Multiple
Figure 2: Mutate Everything decodes any mutation. The model takes a protein sequence of length input. (a) pre-trained feature extractor (e.g. ESM [38], AlphaFold [31]) computes per-position features for each amino acid in the sequence (yellow). Each token is reshaped by an MLP into an array of all possible mutations at that location. (b) A single mutation directly decoded its AAG. (c) For mutation set containing multiple mutations, selected mutations are aggregated and decoded obtain predicted AAG.
The mutation sets include higher-order mutations (K 1), which is beyond the scope of state-of-the art single mutation predictors (K = 1) [18, 36, 74]. These higher-order mutations are strictly more challenging than single mutations. The model must understand the effect of each single mutation and also the epistatic interactions that arise between individual mutations. We introduce a method that captures these interactions for any mutation and decodes a corresponding AAG.
Furthermore, existing paradigms cannot efficiently scale to a large number of mutations. Higher-order mutation models [16, 35, 44, 64] take protein and its mutant as input and predict the change in stability. These methods would need to enumerate and sort over a million model evaluations to find the most stabilizing mutations for small protein of 75 amino acids. Our model runs fine-tuned backbone only once and a lightweight decoder once for each mutation in parallel.
Feature Extraction. We fine-tune AlphaFold and other pre- pre-trained models to alleviate the scarcity of labeled data. AlphaFold's Evoformer learns co-evolutionary patterns, such residue interactions, and the Structure Module learns structural information, such as solvent accessibility or residue-residue distance. We demonstrate this is connected protein's 's stability and how mutations affect stability.
Each element of the protein sequence tokenized by its amino acid type wi and positionally encoded. This produces set of high-dimensional tokens for the transformer to take as input. The tokens are fed into a pre-trained model, which extracts per-position features (21)f=1 where æ IRd The transformer feature extractor is finetuned with our training objective. This is illustrated in Figure 2a.
Methodology
Our model uses per- per-posstion features (z1) determine the stability mutation set M. We distinguish between single point mutations K = M = 1, and higher-order mutations K > 1.
Both single and higher-order mutations share the same backhone architecture and encoding scheme, see Figure 2b. For each mutation, μ = (p, t) at position p to amino acid f, we compute a latent mutation representation z 2(4) E IRd, which captures the effect the mutation on the protein. We decompose the mutation representation =(x) f'(x,)) + ht' where sequence-dependder features f' {(f { project the per- per-position feature according to the mutated amino acid type, and eequuence-indpeedeeen embeddings unique to the amino acid types A, C, ...,, Y. Intuitively, f captures a contextualized feature at the position, and /t captures information about the mutated amino acid type. This mutation representation builds up feature representing the amino acid substitution.
features h' E {h "...... }
Decoding single mutations. Figure 2b illustrates single mutation decoding. The change in thermo- dynamic stability AAG for a single mutation is read out from its mutation representation 2(4). A lightweight linear head g² Rdd - I decodes this representation into the change in thermodynamic stability: P ==(2())).

---

A AG for P1C:V4D
P
A
Head
A
A
C
D
Y
K
A
Y
Amino Acid Expansion
AÁG for M5Y:G6A:H7C
V
A
Head
M
A
Y
AlphaFold
G
A
C
D
Y
A AG for H7Y
Head
H
A
...
Figure 3: Mutate Everything decodes many mutations. Each token output by the backbone zạ is expanded to 20 tokens 2(H), each corresponding to one unique amino acid type (yellow with amino acid types). For each mutation set M, we first aggregate the features corresponding each mutation's position and amino acid being mutated to (arrows). A lightweight head then decodes the aggregated features obtain the predicted AAG. These heads are illustrated in Figure 2b and Figure 2c. Parallel feature aggregation and decoding enable efficient stability prediction for millions of mutations.
Decoding higher-order mutations. Figure 2c illustrates higher-order mutation decoding. The AAG for higher-order mutation M {μr ) is read out from its K corresponding single point mutation representations ([()))] Mutate Everything aggregates the representations for each single point mutation within higher-order mutation. By aggregating the single mutations in the latent space, the model Iearns account for the interactions between them. We aggregate with summation for its simplicity. A lightweight head I > E then decodes the higher-order mutation seamlessly parameterize any higher-order
AAGPr g Ex 2(44)). This formulation allows mutation. For higher-order mutations (e.g. K > 1), we empirically find that is beneficial to predict the residual the sum of single mutation AAGs from our single mutation decoder (see Table 5d).
Decoding million mutations in one forward pass. This decoding scheme efficiently compute AAG for millions of mutations, whether single or higher-order (see Figure 3). Our insight is to precompute all mutation features and reuse them for all predictions. First, compute all Lx20 mutation representations Z = {=((p, t)) : E [1.... I.] E AA). This expands each token my to 20 representations that each capture the effect of mutating to specific amino acid (including the wild-type amino acid). The backbone and this expansion are executed once for a protein sequence.
A lightweight head decodes AAG for any higher-order mutation (a mutation set). For each mutation set, we index and add the relevant mutation representations in Z. Index and sum operations are parallelized with outer- summation for K 2 and the "gather" operation for K 2. Finally, the lightweight decoder computes the AAG for each mutation set in parallel.
Training. Our model is fine-tuned using supervised learning. We optimize the Huber loss.
L(w, (M, AAGE EED
|AAG (w, M.) AAGE
(1)
N
where is a protein sequence and (M, AAGE ) are the mutation sets for this protein with experimentally validated AAG. This approach efficiently backpropagates the loss through all labeled mutations for a given protein in a single pass, unlike prior work which backpropagates through only mutation in single pass. We train single and double mutants with average 2000 mutations for each protein. When training on double mutants, we learn the residual to the sum of the experimental (ground truth) AAG for the constituent single mutations. Note that we use the experimental A AG of the single mutants to train double mutants but during testing, we use the model's predicted AAG of the single mutants infer double mutant AAG. Thus, do not need experimental data for single effects to predict combined effects.
Only 3% of the double mutants in our training set are stabilizing (AAG < -0.5 kcal/mol [8]). To improve our performance on stabilizing mutations, we subsample the number of destabilizing double mutations to most & times the number of stabilizing double mutations.
5

---

Method
Stabilizing r,
nDCG
DetPr
AUC
MCC
Mean
0.00
0.02
0.03
0.50
0.00
MSA
0.04
0,06
0.02
0.60
0.01
ESM2 [38]
0.04
0.08
0.02
0.65
0.00
PROSTATA [74] (Additive)
0.08
0.15
0.05
0.75
0.05
Mutate Everything (Additive)
0.08(0.03)
0.25(0.01)
0.10(0.01)
0.81(0.01)
0.25(0.02)
Mutate Everything (Ours) 0.14(0.02)
0.43 (0.02) 0.16 (0.01) 0.84 (0.01) 0.27 (0.01)
Table 1: Stabilizing Double Mutation Results on cDNA2. The metrics evaluate the model's performance on stabilizing mutations. Our additive baseline naively adds predicted single mutation AAGs. Our model state-of-the-art at finding stabilizing mutations. The proteins in cDNA2 val have at most 36% sequence similarity to those in the training set. Parenthesis indicates standard error across 7 training runs.
5
Results
5.1
Implementation Details
Our primary feature extractor is AlphaFold. implemented in OpenFold 31]. The multiple sequence alignment (MSA) is computed using Colabfold [43]. Sequence features from the Evoformer and Structure Module are aggregated as input to the decoder. dimension to D 128. Our amino acid projections and single mutation decoders are implemented with a linear layer. The higher-order mutation decoder transforms the previous embedding with a 2-layer MLP. These representations are aggregated and fed into 3-layer MLP to predict AAG.
adapter maps the backbone hidden
We train the cDNA display proteolysis dataset [72], which leverages a high throughput proteolysis screen enabled by next generation sequencing (NGS) to extract noisy AAG values for 100 mini- proteins totaling over 100,000 single and double mutations. To evaluate generalization on unseen proteins, we remove proteins from our training set that have high sequence similarity using MMSeqs2 [70] to any protein in our evaluation benchmarks. The short protein lengths (< 100 amino acids) keep the memory requirement small for the higher-order decoder. See Section A. for more information about cDNA and other datasets.
We fine-tune a pre-trained backbone on single mutations for 20 epochs. Then, we finetune the model on both single and double mutations for 100 epochs using a cosine learning rate schedule with 10 warmup epochs. We use a batch size of 3 proteins due to the high memory requirements of AlphaFold. We learning rate of 3x 10-4 and weight decay of 0.5. Training takes 6 hours on 3 A100 GPUs.
5.2 Finding Stabilizing Double Mutations
We evaluate our model's ability to find the stabilizing mutations on cDNA2. cDNA2 is our validation split of the cDNA double mutations dataset, consisting of 18 proteins totaling 22,000 double mutations. The proteins in the validation set have at most 36% homology with those in the cDNA training set. Of these mutations, only 198 or 0.8% are stabilizing with AAG < -0.5 kcal/mol [8]. Like protein engineering, this scenario closely resembles the challenge of identifying small number of stabilizing mutations amidst a much larger pool of destabilizing or neutral mutations.
We evaluate different methods on stabilizing mutations in Table 1. Stabilizing r, is the Spearman coefficient on the experimentally stabilizing subset introduced in [65]. Normalized discounted cumulative gain (nDCG) [30] measures the quality of the ranked mutation list by taking into account its experimentally validated AAG and position in the list [59]. Detection Precision (DetPr) the proportion of experimentally validated stabilizing mutations among the top K = 30 predictions. We adapt previous state-of oftthe----heertttt art single mutation predictors by naively adding both single mutation predictions. A detailed description of the metrics is found in Section A.2.
Our model demonstrates exceptional performance in prioritizing stabilizing double mutations over destabilizing ones, achieving significantly higher normalized discounted cumulative gain of 0,43 compared 0.25, well as a superior detection precision of O. 16 compared to 0. 10. Our model additionally improves classification metrics Matthews Correlation Coefficient (MCC) and Area under Curve (AUC) by 0.02 and 0.03, respectively.

---

Method
I
AUC
MCC
RMSE
Stabilizing Es
Mcan (cDNA)
0.00
0.50
0.00
2.42
0.00
Mean (ProTherm)
0.00
0.50
0.00
2.26
0.00
MSA
0.07
0.51
-0.05
N/A
0.02
ESM2 [62]
0.05
0.51
0.04
N/A
-0.01
FoldX [67]
0.41
2.95
DDGun [44]
0.25
0.63
0.16
2.21
0.17
DDGun3D [44]
0.26
0.64
0.18
2.24
0.17
PROSTATA [74] (Additive)
0.21
0.60
0.05
2.25
0.00
Mutate Everything (Additive) 0.50(0.02)
0.76(0.01)
0.37(0.02) 0.78(0.01) 0.43(0.01) 2.04(0.01) 0.19(0.01)
2.02(0.03)
0.20(0.01)
Mutate Everything (Ours) 0.53(..01)
Table 2: Multiple Mutation Results on PTMul. Our additive baseline naively adds predicted single mutation AAGs. Our model presents strong mutation assessor. The proteins in PTMul have at most 35% sequence similarity to those in the training set. Parenthesis standard error across runs.
To the best of our knowledge, Mutate Everything is the first work that models all double mutations computationally tractable manner. Figure shows the runtime of several methods on protein of 317 amino acids on an A100 GPU, The dashed line indicates the transition from evaluating single mutants to double mutants. Mutate Everything predicts AAG for all single and double mutations in pass of the model. It runs in 0.6 seconds using an ESM2 backbone, and 12.1 seconds
Mutate Everything (AF)
Mutate Everything (ESM)
10a
PROSTATA
Runtime (Sec)
10
10
10a
10:
AlphaFold backbone. PROSTATA [74] also uses the ESM2 backbone but takes protein sequence and its mutated sequence input and outputs the change in stability. PROSTATA takes 306 hours evaluate all double mutations with batch size of 768. On an 8 GPU node, this will take 1.5 days.
10a
102
10f
10a
Number of Mutations
Figure 4:
Runtime analysis.
Our
model's runtime is constant.
5.3 Higher-order Mutation Results
We evaluate our model for predicting changes in thermodynamic stability for higher-order mutations in ProTherm. ProTherm is database of thermodynamic (AAG) and thermostability (ATm) ex- perimental characterization protein mutations curated from the literature. We consider a subset of this database that contains AAG values for higher-order mutations, named ProThhrm Multiple (PTMul). PTMul contains 858 mutations [44]. We keep 846 experimentally validated AAG for higher-order mutants, including 536 doubles and 183 triples, after removing 4 duplicate and 8 ambiguous mutations. PTMul proteins have at most 35% homology to the proteins used in training.
Table 2 shows our results PTMul. A description of baselines is provided in Section A.3. Our additive baselines naively add the sum single mutation AAGs. This baseline rigorously evaluates our model's s ability to learn epistasis, or the non-additive interactions when combining two mutations. Mean always predicts the global test dataset mean statistic. ESM and MSA perform poorly when comparing mutations across proteins. On ProTheem, Mutate Everything achieves the state-of-the-a performance 0.53 Spearman correlation T., compared to 0.50 of our additive baseline. A breakdown of PTMul performance double mutations and more (> 2) mutations is found in Table 6.
While other methods also handle multiple mutations, they usually adopt a mutation- or protein-level level train and test split. This results in data leakage since the same protein or homolog will have mutations both the training and testing sets, leading to inflated reported metrics [16, 18, 35, 64]. In our study, we use sequence similarity- splitting where our training proteins have at most 35% homology to those in MTMul. fairly evaluate generalization to new proteins, we exclude these inflated comparisons from our study.

---

Forward
Reverse
Method
AUC
MCC
RMSE
Ig
AUC
MCC
RMSE
Structure-based Methods
MCSM [57]
0.37
0.66
0.13
1.53
0.24
0.61
0.10
2.33
1-Mutant3.0 [11]
0.35
0.64
0.08
1.53
0.17
0.59
0.06
2.35
DUET [56]
0.42
0.68
0.19
1.52
0.26
0.63
0.12
2.16
FoldX [67]
0.27
0.62
0.14
2.35
0.32
0.67
0.20
2.54
MAESTRO [35]
0.46
0.69
0.26
1.45
0.22
0.61
0.11
2.12
PopMusic [17]
0.42
0.69
0.22
1.51
0.24
0.62
0.12
2.10
SDM [77]
0.39
0.67
0.21
1.67
0.14
0.59
0.12
2.15
INPS3D [66]
0.44
0.70
0.20
1.49
0.36
0.69
0.23
1.78
Dynamut [63]
0.38
0.68
0.20
1.59
0.37
0.67
0.17
1.69
ThermoNet [36]
0.38
0.69
0.21
1.62
0.35
0.66
0.18
1.66
PremPS [14]
0.42
0.66
0.20
1.50
0.43
0.66
0.22
1,49
DDGun3D [44]
0.43
0.71
0.27
1.60
0.41
0.71
0.23
1.61
ACDC-NN [7]
0.46
0.73
0.23
1.49
0.45
0.72
0.22
1,50
Stability Oracle [18]
0.53
0.75
0.34
1.44
0.53
0.75
0.32
1.43
Sequence-based Methods
MuPro [15]
0.27
0.61
0.03
1.60
0.22
0.64
0.08
2.41
I-Mutant3.0 00..0gS [11]
0.34
0.67
0.25
1.54
0.23
0.61
0.05
2.25
DDGun [44]
0.43
0.72
0.24
1.73
0.41
0.71
0.23
1.76
ACDC-NN-Son [55]
0.44
0.72
0.26
1.52
0.44
0.72
0.23
1.52
INPS-Seq [66]
0.44
0.72
0.25
1.52
0.44
0.72
0.21
1.53
PROSTATA [74]
0.50
0.73
0.28
1.44
0.50
0.73
0.29
1.44
Mutate Everything (Ours)
0.56
0.76
0.35
1.38
0.49
0.73
0.30
1.48
Table 3: Comparisons on Single Mutation Thermodynamic Stability prediction in S669. We report regression and classification metrics the forward and reverse datasets for structure-based approaches (top block) and sequence-based approaches (bottom block). r refers to Spearman corre- lation coefficient. Our training proteins have at most 30% sequence similarity with those in S669. Mutate Everything is state-ofthhe-ee on S669 forward.
5.4 Single Mutation Results
We additionally validate our model's ability to predict changes thermodynamic stability under single mutations. We compare Mutate Everything prior works the newly introduced S669 dataset, which was established to address the data Icakage issuc and enable fair comparisons between existing methods [54]. We also evaluate the commonly studied reverse dataset setting, in which the model takes the mutant sequence as input and predicts AAG for mutating back to the wild-type sequence. The experimental value is obtained by negating the original AAG value going from the wild-type sequence to the mutant sequence. S669 is well curated and extensive, containing 94 proteins totaling 669 mutations [54]. S669 proteins have most 30% sequence similarity with those in our training set to ensure separation between the training and validation sets.
Table 3 shows our comparisons with prior works. Mulate Everything is ttate-of-the-art S669, achieving a Spearman correlation of 0.56, where the prior art obtained 0.53. Mutate Everything outperforms the existing sequence model by large margin, with 6-point gain Spearman correlation and 3-point improvement on AUC. Our method even as powerful as the latest structure-based methods which have access atom coordinates.
Mutate Everything is competitive on the reverse dataset evaluation, in which the original and mutated sequences are flipped and AAG negated. Our performance drops on reverse mutations because reverse mutations are out of distribution for our model. We did not perform data augmentation to train on reversed mutations commonly done in the literature [18, 36, 74]. We found that it was beneficial bias our reverse predictions using AAG to the amino acid found the original sequence. We found the performance to be similar across 5 runs (<: 0.01 standard error).

---

TroteinGym-Stabiliy
ProteinGym
Model
ra
AUC
MCC
AUC
MCC
DeepSequence [61]
0.44
0.74
0.34
0.43
0.74
0.35
EVE [24]
0.47
0.75
0.36
0.46
0.76
0.36
ESMIv [42]
0.41
0.72
0.31
0.41
0.74
0.33
Progen2 [47]
0.45
0.74
0.35
0.42
0.74
0.34
MSA Transformer [60]
0.50
0.77
0.38
0.44
0.75
0.35
Tranception [49] Tranception+MSS
0.40
0.71
0.31
0.40
0.73
0.32
0.45
0.74
0.35
0.45
0.76
0.36
TranceptEVE [50]
0.50
0.77
0.38
0.48
0.77
0.38
MSA
0.41
0.72
0.32
0.39
0.72
0.31
Mutate Everything
0.52
0.78
0.39
0.38
0.72
0.30
Mutate Everything +MSA
0.52
0.78
0.41
0.44
0.75
0.36
Mutate Everything +
0.53
0.79
0.42
0.49
0.78
0.39
Table 4: Comparisons on mutation fitness prediction ProteinGym. We report ranking and classification metrics both subset of ProteinGym with stability-like phenotypes and the entire ProteinGym dataset. Our stability model can generalize to all phenotypes, outperforning existing methods stability-like phenotypes and performing similarly to existing methods otherwise. The maximum homology between our training set and ProteinGym proteins is 43%.
5.5 Generalization to other Phenotypes
We show that stability predictors generalize other phenotypes well. ProteinGym substitutions a dataset of 87 proteins with 1.6 million mutation sets labeled with a generic fitness score, such as activity, growth, or expression. The mutation sets have up to 20 substitutions and 14 proteins contain higher-order mutations. The ProteinGym-Stabie fitness labels that known to correlate with thermodynamic stability (thermostability, abundance, expression). ProteinGym proteins have at most 45% sequence similarity with those in our training set. Due memory constraints, we truncate proteins at most 2000 amino acids.
split contains 7 proteins and 26,000 mutations with
Table 4 compares our work with previous works. Our MSA baseline uses the empirical amino acid distribution in the MSA given position compute likelihood ratio. The MSA retrieval prior weights the scores using the per-position first-order statistics in a multiple sequence alignment. On ProteinGymm-Staillee
Mutate Everything obtains Spearman correlation of 0.52, compared to 0.50 of the next best method, TrancepLEVE. On the one protein with labeled stability. Mutate Everything obtains 0.51 1 vs 0.33 of Tranception with MSA retrieval [51]. On the full ProteinGym benchmark, we outperform the strongest methods by ensembling our predictions with Tranception [49]. Ensembling our method with evolutionary models by averaging scores gains 0.1 points on ProteinGym and 0.01 points on ProteinGym troteinGym-Stubilityy suggesting that our model learns complementary information to the evolutionary models. For comparison, TransEVE ensembles an evolutionary model (Tranception) and family-based model (EVE) for a 0.03 point gain [50]. We found the performance to be similar across 5 runs (< 0.01 standard error). An in-depth study on how protein stability related function can be found in [26].
5.6 Ablation Studies
We ablate several architectural design choices in Table 5. Core architectural designs are evaluated on S669. Architectural designs for only multiple mutations are evaluated on cDNA2.
Backbone. We analyze the effect of different feature extractors on single mutation stability assess- ment in Table 5a. ESM2 [38] takes only sequence as input whereas M Transformer [60] and AlphaFoold [31] take the sequence and its MSA as input. Both ESM2 and MSA-Transforme trained on masking residue modeling while AlphaFold is also trained on structure prediction. MSA input improves the stability prediction performance and AlphaFold's architecture and training further improve stability prediction.
are

---

backhone
AUC
aggregator
nDCG
DetPr
AUC
ESM2 [38]
0.47
0.72
outer prod. + flat.
0.39
0.15
0.82
MSA- Transformer [60]
0.53
0.75
product
0.42
0.14
0.85
AlphaFold (scratch)
0.36
0.65
0.43
0.16
0.84
AlphaFold [31]
0.56
0.76
whiten + sum
0.25
0.12
0.80
(a) Backbone on S669, MSA inputs improve performance. A pre-trained AlphaFold model performs the strongest.
(b) Mutation Aggregation cDNA2.
Higher- order mutations aggregate single mutation features. Adding these features competitive.
backbone opt
r
AUC
higher model
nDCG
DetPr
AUC
freeze
0.48
0.72
direct
0.40
0.13
0.85
finetune
0.56
0.76
multiply
0.17
0.08
0.77
add
0.43
0.16
0.84
(c) Backbone Optimization S669. Fine- tuning the AlphaFold backbone improves per- formance.
Higher-order Modeling cDNA2. We model AAG the double mutation as residual com- bined with the AAG its single mutation.
Table 5: Mutate Everything ablation experiments. Default Settings are marked in grey.
Backbone Optimization We try fine-tuning and freezing the AlphaFold backbone during stability fine- tuning in Table 5c. We find that fine-tuning the backbone improves stability prediction performance.
Aggregation The technique to aggregate mutation representations for higher-order mutations is ablated in Table 5b. To control for feature dimension size, we reduce the head dimension when aggregating with outer product and flattening. Aggregating with product and summation performs similarly on double mutations. We add the features for a natural extension to more mutations.
Higher Order Modeling In the direct prediction model, the multi-mutant head directly predicts AAG. In the multiply and add models, the multi-mutant head learns how interactions among single mutations affect the higher-order mutation's s AAG. In these models, the multi mutant head output learns additive bias multiplicative scaling to the sum of single mutation AAGs. In Table 5d, we show that learning a bias performs the strongest among these three options.
Limitations
First, our training dataset contains biases that may affect model performance. The training set contains only small proteins, which may limit performance on larger ones. Our model may exhibit biases towards certain types of mutations due to the data imbalance in training set. Second, the limited availability of experimental stability data poses a challenge for in-silico evaluation. Evaluation on larger and more diverse datasets is necessary to fully assess the generalizability of our model. In the future, we hope that tig-thrroughpt experimental assays will enable more rigorous evaluation and further improvements protein stability prediction.
7 Conclusion
present a method that efficiently scales thermodynamic stability prediction from single mutations to higher-order mutations. Our key insight is that the effects of mutations the same protein are correlated. Thus, for a target protein, it suffices to run a deep backbone once and decode the effect of all mutations simultancously using shallow decoder. With the AlphaFold model our method outperforms existing methods on a variety of single and multiple mutation benchmarks. Our method scales millions of mutations with minimal computational overhead and runs in fraction of the time it would take prior works.
backbone,
8 Acknowledgements
This work is supported by the NSF Al Institute for Foundations of Machine Learning (IFML).
10

---

References
[1] Randomization of genes by per mutagenesis. Genome research 2(1), 28-33 (1992) 2
[2] Adams, J.P., Brown, M.J., Rodrigudd A.. Lloyd. R.C., Roiban, G.D.: Biocataly- sis: A Pharma Perspective. Advanced Synthesis and Catalysis 361(11). 2421-2432 (2019).. http://tbb..ggg/1..00///aaaa.se
[3] Ahdritz, G., Bouatta, N., Kadyan, S., Xia, Q. Gerecke, W., O'Donnell, TJ., Berenberg, D., Fisk, I., M.M., Zhang, S.,
Zanichelli, N., Zhang, B., Nowaczynski, A., Wang. B., Stepniewska--ziubing
Ojewole, A., Guney, M.E., Biderman, S., Watkins, A.M., Ra, S., Lorenzo, P.R., Nivon, L.,, Weitzner, B.. Ban, Y.E.A., Sorger, P.K., Mostaque, E., Zhang, Z., Bonneau, R., AlQuraishi M.: Openfold: Retraining alphafold2 yields new insights into its learning mechanisms and capacity for generalization. bioRxiv (2022). https///di..or//0000 10.1 11110102..1111.. 21..0222.......5722111
[4] Amold, F.H.: Design by directed evolution. Accounts chemical research 31(3), 125-131 (1998) 1,
[5] Bell, E.L.,, Finnigan, W., France, S.P., Green, A.P., Hayes, M.A., Hepworth, L.J., Lovelock, S.L., Niikur,, H., Osuna, S., Romero, E., Ryan, K.S., Turner, N.J., Flitsch, S.L.:: Biocatalysis. Nature Reviews Methods Primers 1(1), 1-21 (2021). https://dii.rg//1000 https://bb.orgg11..0088844433333 1038/$3556-021 021-0004--
[6] Benevenuta, S., Pancotti, C., Fariselli, P., Birolo, G., Sanavia, T.: An antisymmetric neural network predict free energy changes protein variants. Journal Physics D: Applied Physies 54(24), 245403 (2021)3
Benevenuta, S., Pancotti, C.,, Fariselli, P., Birolo, G., Sanavia, T.: An antisymmetric neural network predict free energy changes protein variants. Journal Physics D: Applied Physies 54(24), 245403 (2021) 8
[8] Benevenuta, S., Birolo, G., Sanavia, T., Capriotti, E., Fariselli, P.: Challenges in predicting stabilizing variations: An exploration. Frontiers in Molecular Biosciences 9, 1075570 (2023) 5, 6, 17
191 Brückner, A., Polge, C., Lentze, N., Auerbach, D., Schlattner, U.: Yeast two-hybrid, a powerful tool for systems biology. International journal of molecular sciences 10(6), 2763-2788 2788 (2009)
[10] Cadet, F. Saavedra, E.. Syren, P.O., Gontero, B.: Machine learning, epistasis, and protein engineering: From sequence-structurg functionn relationships regulation of metabolic pathways. Frontiers in Molecular Biosciences 1098289 (2022)
[11] Capriotti, E., Fariselli, P., Casadio, R.: I-mutant2. O: predicting stability changes upon mutation from the protein sequence structure. Nucleic acids research 33(spppl_2) 2), W306-W110 (2005) 3,
[12] Chen, K., Amold, F.H.; Tuning the activity of an enzyme for unusual environments: sequential random mutagenesis of subtilisin for catalysis in dimethy/formmiide Procecdings of the National Academy of Sciences 90(12), 5618 062262222222 (1993) 1.
[13] Chen, T., Gong, C., Diaz, D.J., Chen, X., Wells, J.T., Wang. Z.. Ellingion, A., Dimakis, A., Klivans, A.. et al.: Hotprotein: A novel framework for protein thermostability prediction and editing. In: The Eleventh International Conference on Learning Representations (2022) 3
[14] Chen, Y., Lu, H., Zhang, N., 7hu, 7.,, Wang, S., Li, M.: Premps: Predicting the impact of missense mutations protein stability. PLoS computational biology 16(12), e1008543 (2020) 3.8
[15] Cheng. J.. Randall, A., Baldi, P.: Prediction of protein stability changes for single-site mutations using support vector machines. Proteins: Structure, Function, and Bioinformatics 62(4), 1125-1 132 (2006) 3,
[16] Dehghanpoor, R., Ricks, E., Hursh, K., Gunderson, S., Farhoodi, R., Haspel, N., Hutchinson, B., Jagodzinn ski, F. Predicting the effect of single and multiple mutations on protein structural stability. Molecules 23(2), 251 (2018) 7,7
[17] Dehouck, Y., Kwasigroch, J.M., Gilis, D., Rooman, M.: Popmusic 2.1: a web server for the estimation of protein stability changes upon mutation and sequence optimality. BMC bioinformatics 12(1). 1-12 (201 8
[18] Diaz, D.J., Gong, C., Ouyang-Zhang, J., Loy, J.M., Wells, J.T., Yang, D., Ellington, A.J. Dimakis, A., Klivans, A.R.: Stability oracle: A structure-based graph-transformat for identifying stabilizing mutations (2023). https://doi.org// 11100/023..0.1.544420...
3, 16, 17, 18, 19
[19] Diaz, D.J., Kulikova, A. A. Ellington, A.D., Wilke, C.O.: Using machine learning to predict the effects and consequences of mutations in proteins. Current Opinion Structural Biology 78, 102518 (2023) 3

---

[20] d'Oelsnitz, S., Diaz, D.J., Acosta, DJJ., Schechter, W.W. Minus. M.B., Howard, J.R., Do, H. Loy, J.. Alper. I.. Ellington, A.D.: Synthetic microbial sensing and biosynthesis of amaryllidaceae alkaloids. bioRxiv pp. 2023-04 (2023) 3
[21] d'Oelsnitz, S., Kim, W., Burkholder, N.T., Javanmardi, K., Thyer, R., Zhang, Y., Alper, H.S., Ellington, A.D.: Using fungible biosensors to evolve improved alkaloid biosyntheses. Nature Chemical Biology 18(9), 981 989 (2022)
[22] Ellefson, J.W., Gollihar, J., Shroff, R., Shivram, H. Iyer, V.R., Ellington, A.D.: Synthetic evolutionary origin of proofreading reverse transcriplase. Science 352(6293), 1590-1593 (2016) 1
[23] Feng, X., Sanchis, J., Reetz, M.T., Rabitz, H.: Enhancing the efficiency of directed evolution in focused enzyme libraries by the adaptive substituent reordering algorithm. Chemistry-A European Journal 18(18), 5646-5654 (2012)3 3
[24] Frazer, J., Notin, P., Dias, M., Gomez, A., Min, J.K., Brock, K., Gal, Y., Marks, D.S.: Disease variant prediction with deep generative models of evolutionary data. Nature 599(7883), 91-95 (2021) 9
[25] Gehauer, M., Skerra, A.: Engineered protein scaffolds as next-generation therapeutics. Annual Review of Pharmacology and Toxicology 60, 391-415 (2020). https://doi.orr/10.. http:///00...1////eeeeeeeaaaaaa 44 021118
[26] Gerasimavicius, L.,, Livesey, B.J., Marsh, J.A.: Correspondence between functional scores from deep
mutational scans and predicted effects on protein stability. bioRxiv pp. 2023-02 (2023) 9
[27] Giver, L., Gershenson, A., Freskgard, P.O., Amold, F.H.: Directed evolution of thermostable esterase. Proceedings of the National Academy of Sciences 95(22), 12809 12813 (1998) 2,
[28] Hie, B.L., Shanker, V.R., Xu, D., Bruun, T.U., Weidenbacher, P.A., Tang. S., Wu, W., Pak. J.E., Kim, P.S.: Efficient evolution of human antihodics from general protein language models. Nature Biotcchnology (2023)
[29] Jaroszewiez, W., Morcinck-Orlowska. J., Pierzynowska, K., Gaffke, L., Wegrzyn, G.: Phage display and other peptide display technologies. FEMS Microbiology Reviews 46(2), fuab052 (2022) 1
[30] Järvelin, K., Kekäläinen, J.: Cumulated gain-based evaluation of ir techniques. ACM Transactions on Information Systems (TOIS) 20(4), 422 446 (2002) 6, 17
[31] Jumper, J., Evans, R., Pritzel, A., Green, T., Figurnov, M., Ronneberger, O., Tunyasuvunakool, K., Bates, R., Zidek, A., Potapenko, A., el al.: Highly accurate protein structure prediction with alphafold. Nature 596(7873), 583-589 (2021) 2, 3, 4, 6, 9, 10
[32] Kellogg, E.H., Leaver-Fay, A., Baker, D.: Role of conformational sampling in computing mutation-induces changes in protein structure and stability. Proteins: Structure, Function, and Bioinformatics 79(3), 830-838 (2011)3
[33] Klabunde, T., Petrassi, H.M., Oza, V.B., Raman, P., Kelly, J.W., Sacchettini, J.C.: Rational design of potent human transthyretin amyloid discase inhibitors. Nature structural biology 7(4), 312-321 (2000) 1,
[34] Kouba, P., Kohout, P., Haddadi, F., Bushuiev, A., Samusevich, R., Sedlar, J., Damborsky, J., Pluskal, T., Sivie, J., Mazurenko, S.: Machine learning-guided protein engineering. ACS Catalysis 13, 13863-13895 (2023) 3
[35] Laimer, J., Hofer, H., Fritz, M., Wegenkittl, S., Lackner, P.: Maestro-mult agent stability prediction upon point mutations. BMC bioinformatics 16(1), 1-13 (2015) 3, 7,8
[36] Li, B., Yang, Y.T., Capra, J.A., Gerstein, M.B.: Predicting changes in protein thermodynamic stability upon point mutation with deep 3d convolutional neural networks. PLoS computational biology 16(11). el c1008291 (2020) 3, 4,
[37] Li. G., Panday. S.K., Alexov, E.: Saafec-seq: sequence-based method for predicting the effect of single point mutations on protein thermodynamic stability. International journal of molecular sciences 22(2), 606 (2021) 3
[38] Lin, Z., Akin, H., Rao, R., Hie, B., Zhu, Z., Lu, W., Smetanin, N., dos Santos Costa, A., Fazel- Zaraand, M., Sercu, T., Candido, S., et al.: Language models of protein sequences at the scale of evolution enable accurate structure prediction. bioRxiv (2022) 3, 4, 6, 9, 10, 18
12

---

[39] Loell, K., Nanda, V: Marginal protein stability drives subcellular proteome isoelectric point. Proceedings of the National Academy of Sciences 115(46). 11778-11783 (2018)
[40] Lu. H., Diaz, DJ.., Czarnecki, N.J., Zhu, C., Kim, W. Shroff. R., Acosta, DJJ., Alexander, B.R.. Cole. H.O., Zhang. Y. et al.: Machine learning-aide aided engineering of hydrolases for pet depolymerization. Nature 604(7907), 662-667 (2022) 3
[41] Meghwanshi, G.K., Kaur, N., Verma, S., Dabi, N.K., Vashishtha, A., Charan, P.D., Purohit, P., Bhandari, H.S., Bhojak, N., Kumar, R.: Enzymes for pharmaceutical and therapeutic applications. Biotechnology and Applied Biochemistry (4), 586-601 (2020). https://doi.orgg org/ 1002/bab.. 1919
[42] Meier, J.. Rao, R., Verkuil, R., Liu, J., Sercu, T., Rives, A.: Language models enable zero-shot prediction of the effects of mutations on protein function. Advances Neural Information Processing Systems 34, 29287-29303 (2021) 9, 18
[43] Mirdita, M., Schüitze, K., Moriwaki, Y., Heo, L., Ovchinnikov, S., Steinegger, M.: Colabfold: making protein folding accessible to all. Nature methods 19(6), 679-682 (2022) 6
[44] Montanucei, L.. Capriotti, E.,, Frank, Y., Ben- Tal, N., Fariselli, P.: Ddgun: an untrained method for the prediction of protein stability changes upon single and multiple point variations. BMC bioinformatics 20, 1-10 (2019) 3, 4, 7, 8, 16, 18
[45] Montanucci, L., Savojardo, C., Martelli, P.L.,, Casadio, R., Fariselli, P.: On the biases in predictions of protein stability changes upon variations: the inps test case. Bioinformatics 35(14), 2525-2527 (2019) 16
[46] Musdal, Y., Govindarajan, S., Mannervik, B.: Exploring sequence- space a poplar glutathione transferase using designed information-rich gene variants. Protein Engineering, Design and Selection 30()), 543-549 (2017)3 3
[47] Nijkamp, E., Ruffolo, J., Weinstein, E.N., Naik, N., Madani, A.: Progen2: exploring the boundaries of protein language models. arXiv preprint arXiv:2206. 13517 (2022) 9
[48] Nikam, R., Kulandaisamy, A., Harini, K., Sharma, D., Gromiha, M.M.: Prothermdb: thermodynamic database for proteins and mutants revisited after 15 years. Nucleic acids research 49(D1), D420-D442 (2021) 2, 16
[49] Notin, P., Dias, M., Frazer, J., Hurtado, J.M., Gomez, A.N., Marks, D., Gal, Y.: Tranception: protein fitness prediction with autoregressive transformers and inference-time retrieval. In: International Conference on Machine Learning- pp. 16990-17017. PMLR (2022) 2, 3, 16, 18
[50] Notin, P.M., Van Niekerk. L... Kollasch, A.W. Ritter, D., Gal, Y.. Marks, D.: Trancepteve: Combining samily-specific and family-agnostic models of protein sequences for improved fitness prediction. bioRxiv pp. 2022- 12 (2022)
[51] Nutschel, C., Fulton, A., Zimmermann, O., Schwaneberg, U., Jaeger, K.E., Gohlke, H.: Systematically scrutinizing the impact of substitution sites on thermostabilitty and detergent tolerance for bacillus subtilis lipase a. Journal of chemical information and modeling 60(3), 1568-1584 (2020)
[52] Paik, I., Ngo. P.H., Shroff, R., Diaz, DJJ., Maranhao, A.C., Walker, D.J., Bhadra, S., Ellington, A.D.: Improved dna polymerase variants derived via a machine learning approach. Biochemistry (2021)* 3
[53] Pak, M.A., Markhieva, K.A., Novikova, M.S., Petrov, D.S., Vorobyev. I.S., Maksimova, E.S., Kondrashov, F.A., Ivankov, D.N.: Using alphafold to predict the impact single mutations on protein stability and function. Plos one 18(3), e0282689 (2023) 3
[54] Pancotti, C., Benevenuta, S. Birolo, G., Alberini, V., Repetto,
V., Sanavia, T., Capriotti, E.,, Fariselli, P.: Predicting protein stability changes upon single-poin mutation: a thorough comparison of the available tools on new dataset. Briefings in Bioinformatics 23(2), bbab555 (2022) 2, 8, 16, 17
[55] Pancotti, C., Benevenuta, S., Repetto, V., Birolo, G., Capriotti, E., Sanavia, T., Fariselli, P.: A deep-Icarning sequence-based method to predict protein stability changes upon genetic variations. Genes 12(6), 911 (2021) 8
[56] Pires, D.E., Ascher, D.B., Blundell, T.L.: Duet: a server for predicting effects of mutations on protein stability using an integrated computational approach. Nucleic acids research 42(W1), W314-W319 (2014) 8
[57] Pires, D.E., Ascher, D.B., Blundell, T.L.: mesm: predicting the effects of mutations in proteins using brsed-based signatures. Bioinformatics 30(3), 335-342 (2014) 8
13

---

[58] Pucci, F., Schwersensky. M., Rooman, M.: Artificial intelligence challenges for predicting the impact of
mutations protein stability. Current opinion in structural biology 72, 161-168 (2022) 17
[59] Qiu, Y. Wei, G.W.; Persistent spectral theor--guided protein engineering. Nature Computational Science pp. 1-15 (2023) 3, 6, 17
[60] Rao, R.M., Liu, J., Verkuil, R., Meier, J., Canny, J., Abbeel, Sereu, T., Rives, A.: Msa transformer. In: International Conference on Machine Learning. pp.
PMLR (2021) 3, 9, 10
[61] Riesselman, A.J., Ingraham, J.B., Marks, D.S.: Deep generative models of genetic variation capture the effects of mutations, Nature methods 15(1)), 816-822 (2018) 9
[62] Rives, A., Meier, J., Sercu, T., Goyal, S., Lin. Z., Liu, J., Guo, D., Ou,, M., Zitnick, C.L., Ma, J., Fergus, R.: Biological structure and function emerge from scaling unsupervised learning to 250 million protein sequenees. PNAS (2019). httpss//doi.ory org/ 10.1100/622800 3,7
[63] Rodrigues, C.H., Pires, D.E., Ascher, D.B.: Dynamut:
predicting the impact of mutations protein conformation, flexibility and stability. Nucleic acids rescarch 1). 11350-W555 (2018) 8
[64] Rodrigues, C.H., Pires, D.E., Ascher, D.B.: Dynamut2: Assessing changes in stability and flexibility upon single and multiple point missense mutations. Protein Science 30(1), 60-69 (2021) 3,4,7
[65] Sanavia, T., Birolo, G.. Montanucci, L.. Turina, P. Capriotti, E., Fariselli, P.: Limitations and challenges in protein stability prediction upon genome variations: towards future applications in precision medicine. Computational and Structural Biotechnology Journal 18, 1968- 1979 (2020). https://dooi.org/10.. http://bbbboor///0000016....../p c0166f...bj..220.07....
https://dan orgg /0/ 1atpg//..11011//11111 csbjj 2020. 00.20.0..00 07.007.01111 011 6, 17
[66] Savojardo, C., Fariselli, P., Martelli, P.L., Casadio, R.: Inps-md: web server to predict stability of protein variants from sequence and structure. Bioinformatics 32((6), 2542-2544 (2016) 3,8
[67] Schymkowizz, J., Borg, Stricher, F., Nys, R., Rousscau, F., Serrano, I..: The foldx web server: an online force field. Nucleic acids research (3(supp 2), 9382-W388 (2005) 3, 7, 8, 18
[68] Shroff, R., Cole, A. W. Diaz, D.J. Morrow, B.R., Donnell, I,, Annapareddy, A., Gollihar, J., Ellington, A.D., Thyer, R.: Discovery of novel gain-of-function mutations guided by structure- deep learning.
ACS synthetic biology 1)1,, 2927-2935 (2020) 3
[69] Starr, T.N., Thorton, J.W.: Epistasis in protein evolution. Protein science 25(7), 1204-1218 (2016) 1
[70] Steinegger, M., Soding, J.: Mmseqs2 enables sensitive protein sequence searching for the analysis of massive data sets. Nature biotechnology 35(11), 1026-1028 (2017)
[71] Stemmer, W.P.: Rapid evolution of protein in vitro by dna shuffling. Nature 370(6488), 389-391 (1994)
3
[72] Tsuboyama, K., Dauparas, J., Chen, J., Laine, E., Mohseni Behbahani. Y., Weinstein, J.J. Mangan, N.M.,
Ovchinnikov, S., Rocklin, G.J.: Mega-scale experimental analysis of protein folding stability biology and protein design. bioRxiv pp. 2022- 12 (2022) 1, 2, 6, 16
[73] Tung, J. J.W., Heydari, K., Tirouvanziam, R., Sahaf, B., Parks, D.R., Herzenberg, L.A.,, Herzenberg, L.A.: Modern flow cytometry: a practical approach. Clinics in laboratory medicine 27(3), 453-468 (2007)
[74] Umerenkov, D., Shashkova, T.I., Strashnov, P.V., Nikolaev, F., Sindeeva, M., Ivanisenko, N.V., Kardymon, O.L.: Prostata: Protein stability assessment using transformers. bioRxiv pp. 2022 12 (2022) 3, 4, 6, 7, 8, 18, 19
[75] Wittmann, B.J., Johnston, K.E., Wu, Z., Arnold, F.H.: Advances in machine learning for directed evolution. Current opinion in structural biology 69, 11-18 (2021) 3
[76] Wittmann, B.J. Yue, Y., Arnold, F.H.: Informed training design enables efficient machine learning- assisted directed protein evolution. Cell Systems 12(11), 1026-1045 (2021) 3
[77] Worth, C.L., Preissner, R., Blundell, T.L.: Sdm server for predicting effects of mutations on protein stability and malfunction. Nucleic acids research 39(supl 2), W215-W222 (2011) 3, 8
[78] Wu, S., Snajdrova, R., Moore, J.C., Baldenius, K., Bornscheuer, U.T.: Biocatalysis: Enzymatic Syn- thesis for Industrial Applications. Angewandte Chemie Internatoona Edition 60(1). 88 119 (2021). https://doi.org/10.10000 wtwp:00000022aae2222200666666.
14

---

[79] Yeung, N.. Lin, Y.W., Gao. Y.G., Zhao, X., Russell. B.S., Lei, L... Miner, K.D., Robinson, H., Lu. Y.: Rational design of structural and functional nitric oxide reductase. Nature 462(7276), 1079-1082 (2009) 1.2
[80] Zhou, B.,, Lv. O., Yi, K.. Xiong. X., Tan, P. Hong. L.. Wang. Y.G.: Accurate and definite mutational effect prediction with lightweight equivariant graph neural networks. arXiv preprint arXiv:2304.06299 (2023)
[81] Zhou, Y., Pan, Q.. Pires, D.E., Rodrigues, C.H., Ascher, D.B.: Ddmut: predicting effects of mutations protein stability using deep leamnin.. Nucleic Acids Research p. gkad472 (2023)
15

---

A
Experimental Setting
We elaborate on our experimental protocol, including the datasets used for training and testing, the evaluation metries, and the baselines that calibrate our models.
A.1 Datasets
fundamental principle in machine learning is to partition data so that the model is not evaluated on the same data training data (in our case, to new proteins). Protein sequences are regarded as the same protein when they have large sequence similarity (e.g. 50%), even when they are not exactly the The impact of training on protein sequences similar to those in the test set has received considerable attention [18, 45, 54]. We list the datasets we use below and report the maximum sequence similarity of our training set with each test set. This ensures a clear separation between our training and test sets.
trained on. This is essential assess the model's ability generalize beyond its
string.
cDNA proteolysis [72] is large scale dataset containing mutant proteins with AAG measurements. The dataset is generated by cDNA display proteolysis, a massively parallel technique that measures thermodynamic folding stability for millions of protein variants [72]. The cDNA proteins average 56.1 amino acids in length with a maximum length of 72 and a minimum length of 30 amino acids. The mean MSA depth 7797 with a standard deviation of 6282. The maximum depth is 23525 and the minimum depth is 5. The mutations contain either one or two substitutions on small protein. We follow the train-val split introduced in [18] and additionally filter out proteins in our training set that are similar to those in our evaluation datasets. Specifically, we train on 116 proteins with 213,000 total mutations, of which 97, 7,000 are double mutants and 117,000 are single mutants. We hold out cDNA2, validation set of 18 mini-proteins with 22,000 total double mutations. The proteins in our training set have at most 36% sequence similarity to those in cDNA2. A comprehensive analysis of the dataset, experimental assay, and filtering criteria are found in their paper [1].
ProThhrmn [48] contains experimentally measured thermodynamic data AAG for protein mutants. We use PTMul, split containing only higher-order mutations curated by DDGun [44]. There are 88 proteins with 846 total mutations in PTMul. The proteins in training set have most 35% sequence similarity to those in PTMul.
S669 [54] a manually curated dataset of mutant proteins and their AAG measurements. S669 contains mutations with one substitution. There are 94 proteins with 669 total single mutations in S669. The proteins in our training set have at most 30% sequence similarity to those in S669.
ProteinGym [49] a curated set of proteins in which mutations are assessed for some fitness phenotype. These fitness phenotypes vary drastically, from cell growth to expression to fluorescence. The fitness values are always normalized so higher value indicates stronger fitness. For classification metrics, the fitness values are binarized by some threshold either defined manually or arbitrarily at the median value. We also consider a subset of the proteins in ProteinGym in which fitness correlates highly with thermodynamic stability. These phenotypes include thermostability, expression, etc. There are 87 proteins totaling over 1.6 million (potentially higher-order) mutations. The proteins in training set have at most 45% sequence similarity those in ProteinGym.
T2837 [18] is aggregate dataset with mutant proteins and their AAG measurements. T2837 contains mutations with one substitution. There are 129 proteins with 2837 total single mutations in T2837. The proteins in our training set have at most 34% sequence similarity to those in T2837.
A.2 Metrics
We evaluate our model's ability to assess mutations in ProTherm, S669, and ProteinGym. Additionally, we evaluate our model's ability to detect stabilizing mutations in cDNA2.
A.2.1
Mutation Assessment
The model predicts the AAG measurement under all mutations with experimental data. The models evaluated by the following regression and classification metrics. The classification metries are
16

---

useful in applications where it is unnecessary to provide exact AAG and a binary stabilizing or de-stabilizing label (or third label, neutral) is sufficient.
Following the literature, we additionally evaluate reverse mutations (in contrast to the standard "forward" mutations) where the "to" and "from" amino acids are flipped (e.g. PIK K1P). The model takes the mutant sequence as input and is expected to output the negation of the original AAG value when mutating the original native amino acid.
Spearman correlation coefficient (r,) evaluates the model's ability rank the mutations by their AAG values. This metric considers only the relative ordering between mutations and disregards the predicted values. Spearman has its limitations because our datasets contain predominantly destabilizing mutations (AAG > 0). Spearman applied on these datasets overwhelmingly measures the model's ability to rank destabilizing mutations which does not directly translate to identifying the most stabilizing mutations [8, 58].
Root Mean Square Error (RMSE) provides measure of how closely the predicted measurements align with the true measurements.
Area Under the receiver operating characteristic Curve (AUC) is popular classification metric that focuses on positives. The receiver operating characteristic curve is created by plotting precision against recall at various classification thresholds. Precision calculated as the ratio of true positives to the sum of true positives and false positives. Recall is calculated as the ratio of true positives to the sum of true positives and false negatives.
Matthew correlation coefficient (MCC) evaluates model's classification performance under imbalanced datasets. takes into account true positive, true negative, false positive, and false negative predictions to measure the overall performance of a model.
A.2.2 Mutation Detection
The regression metrics for mutation assessment are not directly suitable for protein engineering where the goal is find stabilizing mutations. In addition the classification metries, we evaluate our models using per-protein metrics which evaluate the model's ability to detect a few stabilizing mutations from large number of possible mutations. These metrics are applied to heavily annotated proteins (~ 1000 AAG measurements) with randomly sampled mutations.
Stabilizing Spearman correlation coefficient r, is the Spearman coefficient calculated on the subset of experimentally stabilizing mutations. This mimics the case in protein engineering where there too many stabilizing mutations to test in vitro and the goal is to prioritize the most stabilizing ones. This metric is used in [18, 65].
Normalized Discounted Cumulative Gain (nDCG) [30] quantifies the quality of ranking among mutations by considering both the experimental AAG and predicted rank of each mutation. Specif- ically, nDCG computes the sum of experimentally stabilizing AAGs weighted inversely by their ranks. This metric rewards a model that ranks the best mutations first. We truncate the ranking at K 30 mutations to approximate the number of in-vitro experiments that can typically be run in parallel. This metric is found in [59].
Detection Precision (DetPr) measures the fraction of experimentally stabilizing mutations among the top K 30 predicted mutations by the model. This simulates the success rate a protein engineer might expect when deploying this model. K is set to approximate the number of in-vitro experiments that can typically be run parallel.
A.3 Baselines
Several baselines calibrate performance across datasets. A subset of the critical baselines is introduced briefly (see [54] for more).
Mean is simple baseline in which the mean measurement is predicted for all mutations. This baseline is not conditioned on the protein or mutation and instead uses only the dataset statistic.
Multiple Sequence Alignment (MSA) leverages the statisties from a set of similar protein sequences to rank mutations. Similar protein sequences are fetched, aligned, and indexed at mutation position.
17

---

This yields an amino acid frequency count normalized into a probability distribution. The score is defined as the likelihood ratio between the mutated and native amino acid types. While calibrated for given position and protein, it may not be well-calibrated across proteins.
ESM [38] also computes the likelihood ratio between a mutated and native amino acid type [42]. Unlike MSA, the probability distribution decoded from ESM2 [38], a model pre-trained on masked amino acid prediction. Like MSA. ESM may not be well-calibrated across proteins.
DDGun [44] uses a linear combination of scores derived from evolutionary and structural features to compute AAG.
FoldX [67] leverages empirical force field estimate the change in stability. models the Van Der Waals interactions, solvation and hydrogen bonding interactions, changes in entropy, and other physical values.
PROSTATA [74] sequence model for stability prediction built on the ESM2 model. It takes native and mutant protein sequences as input and produces AAG value.
Stability Oracle [18] is a structure-based graph-transformer model for stability prediction.
Tranception [49] is a specialized autoregressive transformer that uses grouped attention to encourage specialization across heads.
B Additional Results
B.1 Higher-order Mutation Results on PTMul Breakdown
ProTherm Multiple (PTMul) contains 846 mutations, of which 536 are double mutations and 310 mutations with more than two mutations. We breakdown the performance on PTMul into double mutations and higher-order mutations (greater than in Table 6. While our model performs similarly against the additive baseline on double mutations, Mutate Everything outperforms the additive baseline on higher-order mutations.
Double Mutations
More Mutations
Method
r,
AUC
MCC
RMSE
AUC
MCC
RMSE
DDGun
0.28
0.63
0.22
2.23
0.15
0.57
0.12
2.19
DDGun3D
0.29
0.61
0.17
2.25
0.19
0.60
0.17
2.20
Additive Baseline
0.52
0.76
0.30
1.95
0.49
0.79
0.45
2.14
Mutate Everything (ours)
0.50
0.75
0.34
2.08
0.60
0.86
0.58
1.99
Table 6: Results on PTMul splits containing double mutations and more than 2 mutations. Mulate Everything outperforms existing methods, especially when there are multiple mutations.
18

---

Forward
Reverse
Method
AUC
MCC
RMSE
AUC
MCC
RMSE
PROSTATA-IFME [18, 74]
0.53
0.75
0.31
1.75
0.53
0.75
0.32
1.75
Stability Oracle [18]
0.62
0.81
0.39
1.65
0.62
0.81
0.39
1.65
Mutate Everything (Ours)
0.65
0.79
0.39
1.59
0.52
0.72
0.28
1.75
Table 7: Comparisons on Single Mutation Tharmodynamic Stability prediction in T2837. We report regression and classification metrics the forward and reverse datasets. refers to Spearman correlation coefficient. Our training proteins have at most 34% sequence similarity with those in T2837. Mulale Everything state-of-the-art T2837 forward.
Test Set
#Mutations
RMSE
MCC
AUROC
Precision
Recall
Accuracy
T2837 Orig
2835
0.6274
0.6472
1.5877
0.3858
0.7956
0.5615
0.4791
0.7894
T2837 TR
2835
0.5301
0.5208
1.7563
0.2802
0.7229
0.8297
0.8572
0.7517
T2837 Orig TR
5670
0.7229
0.7282
1.6741
0,5482
0,8433
0.7765
0.7701
0.7705
Table 8: Mutate Everything regression and classification metrics on T2837 and its reverse datasets. TR refers to Thermodynamic Reversibility. Ip is Pearson correlation coefficient. is Spearman correlation coefficient. Mutations that failed our data pipeline are excluded.
B.2 Single Mutation Results on T2837
We include our results on the newly proposed T2837 dataset in Table [18]. Mutate Everything
achieves Spearman correlation of 0.65 compared to 0.62 of Stability Oracle, outperforming existing methods on the forward dataset. PROSTATA-IFML [18] fine-tunes PROSTATA [74] on the cDNA proteolysis dataset for more direct comparison. Mutate Everything performs competitively with other sequence models on the reverse dataset even without training reverse mutations. We show our
performance on additional metrics for T2837
Table
19


Reviews:
[{"review_id": "deONIFyPTT", "reviewer": "Reviewer_KhLn", "summary": "This work predicts changes in thermodynamic stability for single or higher-order mutations on top of AlphaFold2 modules. The proposed model leverages linear aggregation of mutational scores on all possible sites in the latent space to decode $\\Delta\\Delta G$ value for deep mutations in parallel. \n", "questions": "- How is the ratio '3%' observed in line 176?\n- Similarly, why $\\Delta\\Delta G<-0.5$ are considered stabilizing? Any reference here?\n- Since the authors use ProteinGym, is the proposed method applicable to amino acid addition and deletion?", "limitations": "Although the authors discussed limitations of the research in Section 6, they focused mainly on the technical limitations. No negative society impact was mentioned. ", "rating": 3, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 2, "strengths": "- The workflow is easy to understand, which proposes a lightweight solution for an important direction of research. \n- Compared to most existing work, the proposed method 'runs a strong backbone only once and a lightweight decoder N times in parallel'.", "weaknesses": "- The innovation is incremental. The main algorithm relies heavily on the existing AlphaFold2 model. \n- While the authors claim the main contribution is that they designed 'a simple, parallel decoding algorithm', this is not the first method that tries to decode all possible mutational scores from a large latent representation (a $L\\times 20$ matrix in this paper). See, for instance, https://arxiv.org/pdf/2304.08299.pdf. \n- Notations are poorly explained. For instance, the meaning of p, t, and superscription A, Y are undefined on page 4. \n- References are missed in many places (for instance, ESM2 in Table 3). Also please define abbreviations properly when they first appear (even when they are used widely). For example, AUC and MCC are on page 6.\n- A conclusion section is missed at the end.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "> “The innovation is incremental. The main algorithm relies heavily on the existing AlphaFold2 model.”\n\nWe experimented with multiple backbones, including AlphaFold2, ESM2, MSA-Transformer (see Table 6a). Among these backbones, AlphaFold2 happens to be the strongest performing backbone. We present a new model paradigm that predicts the effects of all mutations, rather than a single mutation or position as in prior works [1,2,3,4]. We see our contribution as orthogonal to the choice of backbone. \n\nEarlier works found that AlphaFold2 is not sensitive to single mutations [5], concluding “AlphaFold may not be immediately applied to other problems or applications in protein folding.” To the best of our knowledge, we are the first to apply AlphaFold2 for predicting ΔΔG values of mutations.\n\n> This is not the first method that tries to decode all possible mutational scores from a large latent representation (a  L × 20  matrix in this paper). See, for instance, https://arxiv.org/pdf/2304.08299.pdf.\n\nThanks for sharing “Accurate and Definite Mutational Effect Prediction with Lightweight Equivariant Graph Neural Networks'' by Bingxin Zhou, Outongyi Lv, Kai Yi, Xinye Xiong, Pan Tan, Liang Hong and Yu Guang Wang. Per NeurIPS guidelines, this paper is considered contemporaneous work as it was arxiv’ed 34 days before the NeurIPS paper deadline. We are happy to discuss it as such in the final version.\n\n> Notations are poorly explained. For instance, the meaning of p, t, and superscription A, Y are undefined on page 4.\n\nWe decided to explain all the notations in the preliminary (p,t,A,Y are all defined in lines 113-115). We can see that this might be confusing and will work on clearing this up.\n\n> References are missed in many places (for instance, ESM2 in Table 3). Also please define abbreviations properly when they first appear (even when they are used widely). For example, AUC and MCC are on page 6.\n\nThank you, we will add them.\n\n> A conclusion section is missed at the end.\n\nWe omitted a conclusion due to space constraints. We will include one in the additional page allotted for the final version.\n\n> How is the ratio '3%' observed in line 176?\n\nLine 176 states \"Only 3% of the mutation sets in our training set are stabilizing.” This is a dataset statistic. Our training set consists of mutations and their corresponding ΔΔG value. The percentage of mutations with ΔΔG < -0.5 kcal/mol is 3%.\n\n> Why ΔΔG < − 0.5 are considered stabilizing? Any reference here?\n\nWe follow Benevenuta et al. [6] in categorizing mutations with a 0.5 kcal/mol decrease in free energy as stabilizing. We will add the reference to the final version.\n\n> Since the authors use ProteinGym, is the proposed method applicable to amino acid addition and deletion?\n\nWe have not considered insertions and deletions yet. Our primary application of protein engineering modifies an existing protein only slightly to increase thermodynamic stability. Insertions and deletions cause a shift in the entire sequence, leading to global changes in the protein structure.\n\n[1] Benevenuta, S., Pancotti, C., Fariselli, P., Birolo, G., Sanavia, T.: An antisymmetric neural network to predict free energy changes in protein variants.\n\n[2] Li, B., Yang, Y.T., Capra, J.A., Gerstein, M.B.: Predicting changes in protein thermodynamic stability upon point mutation with deep 3d convolutional neural networks.\n\n[3] Umerenkov, D., Shashkova, T.I., Strashnov, P.V., Nikolaev, F., Sindeeva, M., Ivanisenko, N.V., Kardymon, O.L.: Prostata: Protein stability assessment using transformers.\n\n[4] Diaz, D.J., Gong, C., Ouyang-Zhang, J., Loy, J.M., Wells, J.T., Yang, D., Ellington, A.J., Dimakis, A., Klivans, A.R.: Stability oracle: A structure-based graph-transformer for identifying stabilizing mutations.\n\n[5] Pak, M.A., Markhieva, K.A., Novikova, M.S., Petrov, D.S., Vorobyev, I.S., Maksimova, E.S., Kondrashov, F.A., Ivankov, D.N.: Using alphafold to predict the impact of single mutations on protein stability and function.\n\n[6] Benevenuta S., Birolo G., Sanavia T., Capriotti E., Fariselli P. Challenges in predicting stabilizing variations: An exploration. \n"}, {"review_id": "6qTDmecghK", "reviewer": "Reviewer_AB6d", "summary": "The authors are concerned with the task of predicting the effects of single- and double-residue mutations on the thermodynamic stability of a protein. They propose a simple method that involves passing combinations of embeddings from a pretrained model (AlphaFold2 or ESM) to MLPs. Compared to existing approaches, theirs is computationally efficient and conceptually simple. It also achieves good results on various benchmarks.", "questions": "-     \"We compute z(μ) = ft(xp) + ht...\" \n\nIf h^t depends solely on t, why separate f^t and h^t?\n\n-     “Sequence features from the Evoformer and Structure Module are aggregated as input to the decoder.”\n\nWhich sequence features exactly? How are they aggregated?\n\n-     “To evaluate generalization on unseen proteins, we train on cDNA proteins with low similarity to those in all evaluation benchmarks.”\n\nHow exactly?\n\n-     “We train on the cDNA display proteolysis dataset [61], which leverages a high throughput selection assay to extract noisy ∆∆G values for 100 mini-proteins totaling over 100,000 single and double mutations.”\n\nWould be good to see some stats in the main paper. Protein length, MSA depth, etc.\n\n-     “Our model demonstrates exceptional performance in prioritizing stabilizing double mutations over destabilizing ones, achieving a significantly higher normalized discounted cumulative gain of 0.43 compared to 0.25, as well as a superior detection precision of 0.16 compared to 0.10. Our model additionally improves classification metrics MCC and AUC by 0.02 and 0.03 respectively.”\n\nIf I understand this correctly, the baseline referenced here is just a version of your model where the heads at the end are replaced with simple addition? Why is that a meaningful comparison?\n\n-     “While other methods also handle multiple mutations, they adopt a train and test split where unique mutations in the test set are not included in the training process. This inadvertently leads to training and testing on the same set of proteins [12, 28, 55]. PTMul proteins have at most 35% homology to the proteins used in these methods’ training. To fairly evaluate generalization to new proteins, we exclude these inflated comparisons from our study.”\n\nIs this widely known? If not, this’ll require more specific justification, either in the main paper or the supplement.\n\n- You finetune AlphaFold2 for the task. I'm kind of curious how well the method works if you simply train the MLPs with the embedding model frozen.\n\n- Could you elaborate more on how AlphaFold2 was finetuned? There are some nontrivial implementation details here; e.g. how cropping is handled.", "limitations": "Authors include a brief discussion of limitations. They do not address potential negative societal impacts of their work (not that I think that would be warranted in this case).", "rating": 7, "confidence": 4, "soundness": 3, "presentation": 4, "contribution": 3, "strengths": "The method is original, efficient, simple to understand, and works well. It addresses an important question in structural biology that is currently far from a solution. It has several favorable properties; for example, one can hot-swap in better embeddings as new models become available.", "weaknesses": "Some of the evaluations lack some important details, and a few claims seem questionable to me. See the \"questions\" section below. Overall, I think the evaluation for high-order mutations is a little light. In particular, while the model is evaluated on a mixture of double- and triple-residue mutations, individual results for each group are not provided. Given that the model wasn't trained on triple-residue mutations, and a purported strength of this method its ability to cleanly generalize to high-order mutations, this seems like an important omission.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "> While the model is evaluated on a mixture of double- and triple-residue mutations, individual results for each group are not provided… this seems like an important omission.\n\nAgreed. While our model performs similarly against the additive baseline on double mutations, it excels at higher-order mutations. We will include this analysis in the final version.\n\nBelow are results for the double mutations in PTMul\n|                   | Spearman    | AUROC       | MCC         | RMSE        |\n|-------------------|-------------|-------------|-------------|-------------|\n| DDGun             | 0.28        | 0.63        | 0.22        | 2.23        |\n| DDGun3D           | 0.29        | 0.61        | 0.17        | 2.25        |\n| Additive Baseline | 0.52 (0.01) | 0.76 (0.01) | 0.30 (0.01) | 1.95 (0.01) |\n| MEM (ours)        | 0.50 (0.01) | 0.75 (0.01) | 0.34 (0.01) | 2.08 (0.01) |\n\nand triple and more mutations in PTMul\n\n|                   | Spearman    | AUROC       | MCC         | RMSE        |\n|-------------------|-------------|-------------|-------------|-------------|\n| DDGun             | 0.15        | 0.57        | 0.12        | 2.19        |\n| DDGun3D           | 0.19        | 0.60        | 0.17        | 2.20        |\n| Additive Baseline | 0.49 (0.04) | 0.79 (0.04) | 0.45 (0.06) | 2.14 (0.08) |\n| MEM (ours)        | 0.60 (0.01) | 0.86 (0.01) | 0.58 (0.02) | 1.99 (0.02) |\n\n> If h^t depends solely on t, why separate f^t and h^t?\n\nWe tried both and found that adding it improves performance. We believe that separating them benefits training, as f^t is a deep network and h^t is an embedding that encodes each amino acid type separately.\n\n> “Sequence features from the Evoformer and Structure Module are aggregated as input to the decoder.” Which sequence features exactly? How are they aggregated?\n\nThe Evoformer and Structure Module both output LxD representations for the input sequence where L is the sequence length and D is the hidden dimension. This is “Single repr. (r,c)” for Evoformer in AlphaFold2 Figure 1. More precisely, in the OpenFold code, it is `s` at the following lines: https://github.com/aqlaboratory/openfold/blob/1d878a1203e6d662a209a95f71b90083d5fc079c/openfold/model/evoformer.py#L823 and https://github.com/aqlaboratory/openfold/blob/1d878a1203e6d662a209a95f71b90083d5fc079c/openfold/model/structure_module.py#L753C1-L754C1 \n\nThe representations are normalized (LayerNorm) and added together before fed into the decoder. We experimented with the MSA and pair representations but decided not to use them to standardize the inputs from all backbones.\n\n> “To evaluate generalization on unseen proteins, we train on cDNA proteins with low similarity to those in all evaluation benchmarks.” How exactly?\n\nWe compute the sequence similarity between the cDNA proteins and proteins in the validation set. We filter out any cDNA protein from our training set with higher than 30% sequence similarity with any test protein. This ensures that the proteins in the test set are unseen.\n\n> cDNA protein statistics.\n\nGreat idea. The cDNA proteins average 56.1 amino acids in length with a maximum length of 72 and minimum length of 30 amino acids. The mean MSA depth is 7797 with a standard deviation of 6282. The maximum depth is 23525 and the minimum depth is 5. A comprehensive analysis of the dataset, experimental assay, filtering criteria are found in their paper [1].\n\n> The baseline referenced here is just a version of your model where the heads at the end are replaced with simple addition. Why is that a meaningful comparison?\n\nThe impact of two mutations together on ΔΔG can differ from the impact of the two mutations performed separately (epistasis). Our finding that the model outperforms this additive baseline suggests that our model learns the interactions between single mutations.\n\n> Is it widely known that training and testing should have low homology?\n\nYes, proteins with 35% or more sequence similarity are typically highly similar (e.g. same function in related organisms). Prior works have shown that training on proteins with substantial overlap with the test set leads to heavily overestimated performance [2].\nMany prior works build training (Q1744 [3]) and validation splits (s669 [4], t2837 [5]) with low similarity overlap. In this work, we filter our training set to keep low sequence similarity with all our validation sets.\n\n> Frozen embedding model\n\nWe find that fine-tuning the AlphaFold2 backbone improves performance on S669.\n|           | Spearman | AUC  | MCC  | RMSE |\n|-----------|----------|------|------|------|\n| Freeze    | 0.48     | 0.72 | 0.20 | 1.46 |\n| Fine-tune | 0.56     | 0.76 | 0.37 | 1.36 |\n\n> Could you elaborate more on how AlphaFold2 was finetuned?\n\nWe do not subsample the protein sequence during fine-tuning, as the proteins are shorter than the crop size. The MSA sampling is performed randomly.\n\nMany implementation details are hard to explain in plain English. We are happy to publish the code upon acceptance.\n\n[1] Tsuboyama K., Dauparas, J., Chen, J., Mangan N., Ovchinnikov S., Rocklin, G.. Mega-scale experimental analysis of protein folding stability in biology and protein design. \n\n[2] Montanucci L, Savojardo C, Martelli PL, Casadio R, Fariselli P. On the biases in predictions of protein stability changes upon variations: the INPS test case.\n\n[3] Li, B., Yang, Y.T., Capra, J.A., Gerstein, M.B.: Predicting changes in protein thermodynamic stability upon point mutation with deep 3d convolutional neural networks. \n\n[4] Pancotti, C., Benevenuta, S., Birolo, G., Alberini, V., Repetto, V., Sanavia, T., Capriotti, T. and Fariselli, P. Predicting protein stability changes upon single-point mutation: a thorough comparison of the available tools on a new dataset.\n\n[5] Diaz, D.J., Gong, C., Ouyang-Zhang, J., Loy, J.M., Wells, J.T., Yang, D., Ellington, A.J., Dimakis, A., Klivans, A.R.: Stability oracle: A structure-based graph-transformer for identifying stabilizing mutations."}, {"review_id": "ENH3BvQJ6T", "reviewer": "Reviewer_xRxg", "summary": "The authors present the \"Mutate Everything Method\", a simple method that builds on top of protein representations obtained through existing models to predict the effect of single and higher-order mutations. They apply the method to representations from ESM2 and AlphaFold, and show the performance on several downstream tasks.", "questions": "1. The way the method is presented, it seems like the authors claim that it will be possible to reliably predict the effect of single mutations as well as high-order effects. Even though some evidence is presented, it is not conclusive, especially not for the higher-order case. I understand that this depends on e.g. the availability of data sets, but perhaps the presentation can be a bit more nuanced.\n2. Even though the method is presented as being simple, there seem to be a lot of MLPs to train! Can you give more details about these MLPs (i.e. how many layers, size of the layers) and the overall memory requirements?\n3. The paper goes straight from results to \"Limitations\", without any conclusions or discussion. As a reader, this section is clearly missing.\n4. Table 4, 5, and 6 show no standard deviations. Is there a good reason for that? Otherwise, it would be great if those could be added.\n5. The part where you predict the residual between single mutation $\\Delta\\Delta G$s and teh combined effect is a bit vague to me. Why was this beneficial? \n6. Related to the previous point, when you say in line 174-175 that you learn the residual to the sum of experimental $\\Delta\\Delta G$s for the constituent single mutations, what does this mean? Are you giving the model the experimental values of single mutations? Does this mean you always need to have experimental data for single effects if you want to predict combined effects? Or am I misunderstanding the sentence? The wording is a bit confusing.\n7. For clarity, it would be great if you could incorporate $x$, $\\mu$, and $z$ (and perhaps $f$ and $h$) in Figure 2.\n8. What is the dimensionality of $d$?\n9. Some of the baselines only show up in tables without really being explained in the main text or the appendix.\n10. There seem to be some references missing in the first part of the introduction (line 25-30), about the use of machine learning for finding stabilizing mutations. And perhaps more citations can be added to section 2.3 (line 104-111) as well.", "limitations": "The authors describe the limitations of their work briefly but accurately. The only thing missing could be a discussion about the low spearman correlation values (in general).", "rating": 5, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 2, "strengths": "The proposed method is based on a very simple but effective idea, and the paper is written in a clear and mostly comprehensive manner. This seems like a relatively low-resource and fast method to get more use out of existing representations. Additionally, a substantial amount of results is presented over a wide range of benchmarks and metrics.", "weaknesses": "Even though I enjoyed reading this paper and I appreciate the idea that a simple method can make a difference for real downstream tasks, I am not entirely convinced that this method, which esentially consists of training a lot of MLPs on top of existing representations, fits the NeurIPS venue well. It doesn't seem like a substantial contribution to the field of machine learning, and this paper might therefore be better suited for a different venue where it would have more impact. Moreover, the reported spearman correlations to experimental results are very low in some places, and around 0.5 at most, signifying a very weak correlation at best. This begs the question whether this is an unfortunate metric or if there's a deeper issue causing the correlation to experimental results to be weak.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "> Reported spearman correlations to experimental results are very low in some places, and around 0.5 at most, signifying a very weak correlation at best. This begs the question whether this is an unfortunate metric or if there's a deeper issue causing the correlation to experimental results to be weak.\n\nWe agree that the spearman correlation coefficient can be a suboptimal metric. Many datasets contain predominantly destabilizing mutations. Spearman overwhelmingly measures the model’s ability to rank destabilizing mutations [2,3], which does not directly translate to identifying the most stabilizing mutations.\n\nPrior works have proposed using metrics more in line with finding stabilizing mutations. Area under precision-recall curve (AUC), Matthews correlation coefficient (MCC) [1], and normalized Discounted Cumulative Gain (nDCG) [4] assess a model’s ability to identify stabilizing mutations. We report all metrics including Spearman for completeness.\nWe will include a discussion of the metrics in the final version.\n\n> The way the method is presented, it seems like the authors claim that it will be possible to reliably predict the effect of single mutations as well as high-order effects. Even though some evidence is presented, it is not conclusive, especially not for the higher-order case. I understand that this depends on e.g. the availability of data sets, but perhaps the presentation can be a bit more nuanced.\n\nRigorously studying higher-order mutations is hard. Public data is scarce and expensive to collect. The data that does exist does support our conclusions. We express our concerns in the Limitations section in line 285. We are more than happy to add nuance to the presentation where needed. What were the things the reviewer would like to change? \n\n> MLP details and memory requirements\n\nAn adapter maps the backbone-specific hidden dimensionality to D=128 and all subsequent layers operate at D=128.\n\nEach amino acid projection $f^t$ in the amino acid expansion is a linear layer.\n\nThe single mutation decoder $g^1$ is a linear layer.\n\nThe higher-order mutation decoder transforms the previous embedding with a 2-layer MLP. These representations are aggregated and are fed into a 3-layer MLP to predict ddg.\nThe MLPs use LayerNorm and ReLU.\n\nWe train on single and double mutations from small proteins with at most 72 amino acids. This keeps the memory requirement small for the higher-order decoder.\n\nThank you for the feedback, we will report these thoroughly in the paper and release code upon acceptance.\n\n> Missing Conclusion\n\nWe omitted a conclusion due to space constraints. We will include one in the additional page allotted for the final version.\n\n> Table 4, 5, and 6 show no standard deviations. Is there a good reason for that? Otherwise, it would be great if those could be added.\n\nThank you for the suggestion. We will add standard errors for all of our experiments in the final version.\n\n> The part where you predict the residual between single mutation ΔΔGs and the combined effect is a bit vague to me. Why was this beneficial?\n\nWe tried predicting the residual and directly predicting the combined ΔΔG for higher order mutations. We found predicting the residual to be easier to train. We will mention this in the exposition.\n\n> Related to the previous point, when you say in line 174-175 that you learn the residual to the sum of experimental  ΔΔGs for the constituent single mutations, what does this mean? Are you giving the model the experimental values of single mutations? Does this mean you always need to have experimental data for single effects if you want to predict combined effects?\n\nWe do not provide the model with any experimental ΔΔG values. For higher order mutations, our model first predicts the ΔΔG for each constituent single mutation (computationally), then predicts a residual to the sum.\n\n> For clarity, it would be great if you could incorporate x, μ, and  z  (and perhaps f and h) in Figure 2.\n\nGreat idea, thank you. We will add the symbols to the figure.\n\n> What is the dimensionality of d?\n\nWe map a backbone-specific hidden dimension to D=128.\n\n> Some of the baselines only show up in tables without really being explained in the main text or the appendix.\n\nWe explain our baselines in Supplement Section 1.3 and provide citations for existing works. We will try to work this into the final version given the extra content page.\n\n> There seem to be some references missing in the first part of the introduction (line 25-30), about the use of machine learning for finding stabilizing mutations. And perhaps more citations can be added to section 2.3(line 104-111) as well.\n\nThank you. Will do. One citation for the explosion in biological data is “Mega-scale experimental analysis of protein folding stability in biology and protein design” by Kotaro Tsuboyama, Justas Dauparas, Jonathan Chen, Niall M. Mangan, Sergey Ovchinnikov, Gabriel J. Rocklin.\n\n[1] Broom, A., Trainor, K., Jacobi, Z., Meiering, E.M. Computational Modeling of Protein Stability: Quantitative Analysis Reveals Solutions to Pervasive Problems.\n\n[2]  Benevenuta, S., Birolo, G., Sanavia, T., Capriotti, E., Fariselli, P. Challenges in predicting stabilizing variations: An exploration.\n\n[3] Pucci, F., Schwersensky, M., Rooman, M. Artificial intelligence challenges for predicting the impact of mutations on protein stability.\n\n[4] Qiu, Y., Wei, G.W. Persistent spectral theory-guided protein engineering.\n"}, {"review_id": "PFs57VCdfC", "reviewer": "Reviewer_31pU", "summary": "The paper proposes a new method, Mutate Everything Method, to model a protein's thermodynamic stability based on mutations in a protein's sequence. A key distinguishing feature of the method is the ability to perform large amount of parallel evaluations, which significantly speed up computational efficiency. The authors first introduce the general challenge in designing stable protein though mutations, describe thermal stability as a physical metric and introduce the general outline of their method. Next, the authors describe related work in protein engineering, protein structure prediction and protein stability modeling followed by a more precise, mathematical description of the problem setup for the proposed method. Essentially, the Mutate Everything Method relies on taking protein sequence embeddings from pretrained models (AlphaFold, ESM) and then predicting thermal stability of the mutations with lightweight MLP heads with one head corresponding to each mutation. In the case of multiple mutations, the authors aggregate the outputs from each individual mutation head. In their experiments, the authors study their method on a variety of sequence mutation datasets, including both single and higher-order mutations with the general results indicating better modeling and compute performance of the Mutate Everything Method. The order also perform an ablation of representations from different embedding encoders, including AlphaFold, ESM2 and MSA-Transformer.", "questions": "* Can you provide more detail on how you ensure the mapping of the mutation in the embedding space is aligned with the mapping of the mutation in the original sequence space?\n* Can you describe in more detail how your parallel evaluation ensure proper mapping of mutations to the correct protein sequence? I am trying to get a better sense of whether the parallel evaluation are essentially operating by creating a larger \"meta-sequence\" or if there is any type of multi-processing happening.\n* It seems like the method requires a frozen embedding encoder as it is implemented right now. How robust would the method be to finetuning the encoder?", "limitations": "The authors include a section on limitations mainly focused on modeling shortcomings. ", "rating": 8, "confidence": 4, "soundness": 4, "presentation": 3, "contribution": 3, "strengths": "The paper has the following strengths:\n* Originality: The paper proposes a new, pragmatic method for a relevant problem in modeling how mutations affect a protein thermodynamic stability.\n* Quality: The paper performs detailed evaluations of the method on multiple types of tasks, including multiple datasets, and compares to various other methods along both modeling performance and compute efficiency.\n* Clarity: The purpose, goals and details of the method and results are generally well presented.\n* Significance: The paper tackles a relevant problem and show notable performance improvements in multiple settings.", "weaknesses": "The paper could be improved by:\n* Providing more clarity into the parallel evaluation pipeline. Figure 3 is not very clear in how exactly multiple mutation evaluations are processed and aggregated. I suggest adding labels that show how a mutation changes a particular part in the sequence similar to Figure 2. [Clarity]\n* The analysis on homology in Section 5.3 is interesting and adds an important dimension on how to properly design data splits for these kinds of tasks. It would be nice to see a deeper discussion on this; potentially in the appendix due to space constraints. [Quality, Significance]\n\n", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "> Figure 3 is not very clear in how exactly multiple mutation evaluations are processed and aggregated. I suggest adding labels that show how a mutation changes a particular part in the sequence similar to Figure 2. [Clarity]\n\nFigure 3 showcases how we can apply Mutate Everything to decode ΔΔG values in parallel. How the multiple mutation head processes and aggregates is identical to that shown in Figure 2c. We will clarify the head in Figure 3.\n\n> The analysis on homology in Section 5.3 is interesting and adds an important dimension on how to properly design data splits for these kinds of tasks. It would be nice to see a deeper discussion on this; potentially in the appendix due to space constraints. [Quality, Significance]\n\nWe agree that homology plays an essential role in designing data splits. There are many works in the literature that recognize this problem [1] and design more informed dataset splits [2,3,4]. We will elaborate on this discussion in the additional page provided in the final version.\n\n> Can you provide more detail on how you ensure the mapping of the mutation in the embedding space is aligned with the mapping of the mutation in the original sequence space?\n\nOur training objective aligns the embeddings with the input sequence. We have an embedding F[L,t] for each position L and “to” amino acid t. The loss aligns a mutation at position L to amino acid t (and its ΔΔG) with the corresponding embedding F[L,t].\n\n> Can you describe in more detail how your parallel evaluation ensure proper mapping of mutations to the correct protein sequence? I am trying to get a better sense of whether the parallel evaluation are essentially operating by creating a larger \"meta-sequence\" or if there is any type of multi-processing happening.\n\nSingle mutations are modeled exhaustively. The decoder outputs a set of Lx20 single mutation representations from which all single mutation ΔΔGs are decoded. For example, F[10,A] decodes the change in stability when position 10 is mutated to Alanine. During training, the mutation representation is aligned to the experimental ground truth ΔΔG.\n\nAll higher-order mutations are made up of the same set of Lx20 single mutations. In our model, these single mutation representations are indexed and aggregated to represent higher-order mutations. For example, a mutation at position 10 to Alanine and 15 to Cysteine is represented as F[10,A] + F[15,C]. The loss will ensure alignment to experimental ΔΔG.\n\nComputationally, we perform one forward pass through the backbone network. The backbone is the computationally most expensive portion of the network. Computation does not depend on the number of higher order mutations considered. ΔΔGs of higher order mutations are computed as inner sums between embeddings at certain positions, and are thus computationally very fast. They do however scale exponentially with the order of mutation considered (there is a linear number of single point mutations, a quadratic number of two point mutations, a cubic number of three point mutations, …).\n\nThe structure of our model enables us to compute millions of ΔΔG values for one protein with only a single forward of the backbone.\n\n> It seems like the method requires a frozen embedding encoder as it is implemented right now. How robust would the method be to finetuning the encoder?\n\nWe find that fine-tuning the AlphaFold2 backbone improves performance. The current method fine-tunes the backbone.\n|           | Spearman | AUC  | MCC  | RMSE |\n|-----------|----------|------|------|------|\n| Freeze    | 0.48     | 0.72 | 0.20 | 1.46 |\n| Fine-tune | 0.56     | 0.76 | 0.37 | 1.36 |\n\n\n[1] Montanucci L, Savojardo C, Martelli PL, Casadio R, Fariselli P. On the biases in predictions of protein stability changes upon variations: the INPS test case.\n\n[2] Li, B., Yang, Y.T., Capra, J.A., Gerstein, M.B.: Predicting changes in protein thermodynamic stability upon point mutation with deep 3d convolutional neural networks.\n\n[3] Pancotti, C., Benevenuta, S., Birolo, G., Alberini, V., Repetto, V., Sanavia, T., Capriotti, T. and Fariselli, P. Predicting protein stability changes upon single-point mutation: a thorough comparison of the available tools on a new dataset.\n\n[4] Diaz, D.J., Gong, C., Ouyang-Zhang, J., Loy, J.M., Wells, J.T., Yang, D., Ellington, A.J., Dimakis, A., Klivans, A.R.: Stability oracle: A structure-based graph-transformer for identifying stabilizing mutations.\n"}]

Author rebuttal:

Meta-review: This paper proposed an efficient method to perform one and double mutations of proteins in a single pass on top of the AlphaFold module. Experimental results on multiple data sets show the effectiveness of the proposed approach. 

After the discussion, all the reviewers like the simplicity of the proposed approach and agree this is a good contribution in solving a practical and important problem.
Decision: Accept (poster)