In this study, we evaluated the ability of the MD+FoldX method to predict SARS-CoV-2 RBD escape mutations using a comprehensive deep mutational scanning dataset. We focused on 19 Ab/Ag systems with structural information of the complexes and escape fractions from the yeast-based DMS technique. Our analysis revealed a positive correlation between predicted affinity and experimental escape fraction. We found that fine-tuning predicted affinity cutoffs using empirical data allowed for more accurate escape mutation identification, underscoring the importance of tailoring cutoffs to specific Ab/Ag interactions rather than adopting a one-size-fits-all approach. Our results demonstrate the potential of MD+FoldX to identify significant affinity-driven escape mutations, some of them already present in variants of concern and interest, and also highlight the need for the development of more accurate methods.This study includes 19 Ab/Ag systems across four distinct Ab classes, all directed against the SARS-CoV-2 S-RBD. This dataset covers Abs targeting similar epitopes, i.e. Abs within the same class, as well as those targeting different structural epitopes, i.e. Abs from different classes. Analysis of the escape fraction data for all the systems showed that only a small proportion of mutations exhibit escape (high escape fractions), consistent with the prevailing understanding that the majority of mutations have negligible effects on antigenicity58. In addition, our analysis of binding affinity predictions for all systems showed that the majority of the values were destabilizing (\(\Delta \Delta G\) > 0 kcal/mol); this is consistent with previous studies in a series of globular proteins using the FoldX algorithm where approximately 70% of mutations were destabilizing59 and with experiments where random mutations tend to be destabilizing60. Notably, our investigation encompasses a substantial dataset comprising 4443 single amino acid mutations across the 19 antibody-antigen systems, each accompanied by experimental escape fraction data. In comparison, related studies have examined significantly fewer mutations: Tandiana et al. validated 27 mutations in different antibody/hen egg white Lysozyme complexes with experimental data61, Gonzalez et al. validated 253 mutations in Ab/Ag complexes27, Miller et al. validated 114 mutations in 27 Ab/Ag systems29, Beach et al. validated 8 mutation in Ab/Ag systems30, and Sharma et al. validated 22 mutations in Ab/Ag systems54. Validating in silico predictions of affinity-driven antibody escape mutations remains challenging due to the scarcity of large-scale experimental data given that obtaining binding affinity data for a large number of mutations is costly and time-consuming. However, our research shows that the extensive escape fraction data from Bloom’s large-scale studies34,35,36 provides a promising alternative for validating this type of prediction.In general, our method was able to identify significant variants of concern and interest such as R346K, L452R, E484K/A/Q, F486F/V, K417N/T, and N501Y. Specifically, the escape mutations E484K/A/Q, identified using the MD+FoldX methodology for LY-CoV555 system, were documented across multiple SARS-CoV-2 variants of concern, including the Alpha (B.1.1.7), Beta (B.1.351), Eta (B.1.525), Mu (B.1.621), and Omicron BA.1 and BA.2 (B.1.1.529) lineages, as well as the Kappa (B.1.627.1) lineages16. Similarly, the mutation F486P was observed in the Omicron (XBB.1.5) subvariant62. E484K that substitutes a positively charged residue with a negatively charged one was demonstrated to escape LY-CoV555 Ab reducing neutralization35. Interestingly, S2X259, a broad sarbecovirus neutralizer, exhibits an escape profile confined to a single mutation, G504D, as revealed by DMS and in vitro escape selection experiments37. The fact that we were able to predict these mutations suggests that our approach has the potential to uncover key escape mechanisms that could undermine current antibody-based treatments and vaccines. This insight provides valuable data for ongoing SARS-CoV-2 research, particularly in developing effective therapeutic antibodies and improving future vaccine formulations. Moreover, understanding these mutation patterns helps researchers monitor viral evolution, ultimately contributing to global efforts in preventing and managing COVID-19 and similar infectious diseases.One primary obstacle we faced was devising a suitable approach to accurately link binding affinity with escape fractions. Here, we assume that the experiment was performed under equilibrium conditions and hence the binding affinity is related to the experimental escape fraction through the Hill relation48,49. We found that the Pearson and Spearman correlation coefficients were consistently positive via the Hill relation. Our mean Pearson correlation coefficient of 0.52 is comparable to the accuracy observed in previous MD+FoldX studies, where a correlation coefficient of 0.39 was reported for predicted versus experimental \(\Delta \Delta G_{\text {bind}}\) in Ab/Ag systems27. In contrast, our mean Spearman correlation coefficient was 0.25, also consistent with the limited performance of other similar methods studying Ab/Ag systems, where Spearman correlation coefficients rarely exceed 0.2827,63. We also confirmed that using a multi-frame approach (MD+FoldX) improves correlation coefficients compared to relying on a single structure (FoldX), as observed in our previous studies27 but with a larger and systematic set of data from DMS experiments. We attribute this improvement to the relaxation of the initial structure, especially if it lacks high resolution, and to the incorporation of variability in the initial structure that is used for mutant construction by FoldX (Figs. S9, S10 of the Supplemental Material).While \(\Delta G_{\text {bind}}\) serves as a valuable metric for assessing binding strength, its correlation with the escape fraction may be influenced by a multitude of complex factors inherent to the binding process: (1) \(\Delta G_{\text {bind}}\) provides a thermodynamic view of affinity, but escape fractions also include kinetic factors such as binding and unbinding rates. Thus, while high \(\Delta G_{\text {bind}}\) coupled with low kinetic barriers could result in high escape fractions (indicating a strong correlation), the presence of high kinetic barriers alongside high \(\Delta G_{\text {bind}}\) can obscure this relationship due to impeded escape despite weak binding53,64. (2) Protein folding variations can significantly affect epitope accessibility and stability. Altered epitopes due to folding changes can reduce antibody affinity, thereby affecting the likelihood of escape. In particular, the yeast-DMS approach excludes mutants with significant misfolding33. Conversely, predictive models often do not calculate antigen folding stability; when they do, their accuracy can be variable, posing additional challenges to accurately predicting escape mutations19,20,21,22,23. (3) The presence of multimeric structures, such as trimers in spike proteins, could introduce cooperative binding effects that could affect binding strength and escape. In yeast-DMS experiments, antibodies bind to monomeric yeast-expressed RBDs, and therefore cannot fully capture mutational effects on spike-trimer conformation or effects on antibodies that bind quaternary epitopes33. (4) Differences in the glycosylation patterns of yeasts compared to human cells may influence antibody recognition of glycan-rich epitopes. In particular, N-linked glycans on yeast-expressed proteins are more mannose-rich than those on mammalian-expressed proteins33. In this study, our simulations do not include glycans and hence our prediction accuracy would suffer for such an epitope. (5) Finally, variations in the structural dynamics of the Ab/Ag complexes, the influence of solvent effects like bridging water molecules, and the presence of allosteric effects or cooperative binding events, are potentially significant factors that are very challenging to account for in simulations. Indeed, all of the factors listed above underscore the challenges of directly correlating binding energies with escape fractions and highlight the intricate interplay between structural, kinetic, and biological phenomena in antibody-antigen interactions.Given the factors listed in the previous paragraph that could impact the correlation between predictions and escape fraction data, we will consider three distinct cases. Case 1: Antibodies bind to monomeric non-glycan epitopes. Here, \(\Delta G_{\text {bind}}\) is expected to correlate with f via the Hill relation (Eq. (1)), meaning that significantly larger \(\Delta G_{\text {bind}}\) in mutants compared to wildtype lead to higher likelihood of escape. Case 2: Antibodies bind glycan-free epitopes that span several monomeric units. In this case, the correlation between \(\Delta G_{\text {bind}}\) and escape fractions f may not be strong since the experimental setup includes only monomeric spike protein leading to larger escape fractions (\(f>0.5\)). Case 3: Antibodies bind glycan-rich epitopes. In this case, \(\Delta G_{\text {bind}}\) may not correlate well with escape fractions f since yeast cells produce different glycans. In our study, systems are predominantly consistent with Case 1, including AZD8895, S2E12, S2H14, REGN10933, C105, LY-CoV555, S2D106, S2H13, REGN10987, LY-CoV1404, CoV2-2130, S2H97, S2X259, C002, and S2X35. Here, \(\Delta G_{\text {bind}}\) could be expected to correlate with f via the Hill relation. Positive high Pearson correlation coefficients associated with these systems (Pearson = \(0.28{-}0.68\)) indicate that larger \(\Delta G_{\text {bind}}\) are associated with high escape fractions. By contrast, C121 is consistent with Case 2, Abs targeting glycan-free quaternary epitopes11. In this system the spike trimer structure has 2 Abs bound, the first one (C121up) interacts primarily with 1 RBD in the “down” conformation but also has minor contact with 1 RBD in the “up” conformation. A second Ab (C121down) interacts primarily with 1 RBD in the “down” conformation and has minor contact with a second RBD in the “down” conformation. Although the correlation between \(\Delta G_{\text {bind}}\) and escape fractions f are not expected to be consistent, \(\Delta G_{\text {bind}}\) values could still serve as useful indicators of potential escape events. Systems C110 and C135 are consistent with Case 3, both featuring a glycosylated epitope containing an N343 glycan11. Finally, the C144 system is a hybrid between cases 2 and 3; its epitope spans multiple monomers and includes a glycosylated site (N434). However, its primary interaction is with one of the RBDs that is not glycosylated, and may explain the reasonable correlation found for this system.In this study, we found that adjusting predicted affinity cutoffs based on empirical data was crucial for accurately identifying escape mutations. Overall, there was significant variability in the optimized cutoff values across different systems and classes. This suggests that the optimal cutoff for predicting escape mutations varies depending on the specific Ab/Ag system, indicating that a one-size-fits-all approach may not be suitable for all systems and highlighting the complexity and specificity of Ab/Ag interactions. The challenge of selecting an appropriate cutoff value is further emphasized by the varied approaches used in different studies to pinpoint SARS-CoV-2 immune hot spots. Chang et al.20 used a \(\Delta \Delta G = 0\) kcal/mol cutoff to classify binding stability, with negative values indicating stabilization and positive ones indicating destabilization. Huang et al.19 developed an immune-escaping score by merging binding free energy measurements with variant frequency, using variants of concern and interest to determine mutation frequency cutoffs. Mauria et al.22 refined binding affinity and escape fraction cutoffs using receiver operating characteristic curve analysis, setting \(\Delta \Delta G = – 0.7\) kcal/mol and an escape fraction of 0.001753. In this study, we fixed the escape fraction cutoff and concentrated on fine-tuning the \(\Delta G_{\text {bind}}\) horizontal affinity cutoff. As observed in our results, significant variability among the systems was found, however, the distribution of optimized cutoffs tends to center around our classic cutoff of 2.0 kcal/mol showing that it could serve as a reasonable starting point for initial predictions. It may be particularly applicable in the absence of system-specific data to guide the optimization of the cutoff value.The data reveals Class 3 Abs possess a minimal tolerance for energy changes, corresponding to lower escape cutoffs, whereas Class 4 Abs exhibit a higher tolerance, apparently requiring more substantial energy changes to allow mutations to escape. We can consider several factors that contribute to the difference in the optimized cutoff for these extreme cases. Firstly, the distribution of the experimental escape fractions differs between Class 3 and Class 4 (Figs. S1, S5–S8 in the supplemental material). In class 4, most values are clustered near zero, with a few near 1. In contrast, class 3 values are less clustered at the extremes and have more (uniformly distributed) points in between. Secondly, FP and FN are on average higher for class 3 than for class 4 (33 compared to 8 FP and 28 compared to 13 FN). Thirdly, the distribution of predicted values is flattened towards zero in class 3 compared to class 4 (Figs. S2, S5–S8 in the Supplemental Material). In summary, we can expect that in a scatter plot where FP and FN are lower, extreme clusters are more easily distinguished from each other, and have a higher rate of change of affinity along the escape fraction, resulting in a higher optimized energy cutoff. On the other hand, DMS revealed that higher binding affinity reduces the total number of viral escape mutations65. One hypothesis could be related to the affinity of antibodies in different classes, which could result in different levels of response and consequently different cutoff values. In class 4, where fewer mutations are distributed near 1, this indicates that fewer mutations interfere with binding. For example, S2X35 from class 4 has a binding affinity for SARS-CoV-2 RBD of 0.2 nM34, whereas REGN10987 from class 3 has an affinity of 45 nM66. Experimentally, 42 mutations are considered escape (f > 0.5) for S2X35 compared to 103 for REGN10987. Predicted escape mutations are 16 for S2X35 and 36 for REGN10987. This suggests that fewer mutations lead to escape for class 4, and tighter binding requires higher energy to disrupt binding, thus requiring higher optimized energy cutoffs for escape prediction.It is crucial to weigh the clinical implications of the cutoff selection carefully: an overly stringent cutoff may overlook potential escape mutations, whereas too lenient a cutoff might overpredict escape, potentially triggering unnecessary concern or interventions. Hence, the optimized cutoff must ensure that predictions are both accurate and practically useful for guiding therapeutic strategies and vaccine design. Aiming for the ambitious goal of predicting escape mutations solely through \(\Delta G\) values, adopting a proposed stepwise approach could prove instrumental. We propose the use of system-dependant optimized cutoff whenever possible. If data for optimizations of the particular system is lacking, we can formulate three distinct lists per class based on varying cutoff values to cater to different scenarios. For an optimistic approach, a cutoff based on Q1 in Table 3 is recommended. This threshold may encompass some non-escape mutations but is aimed at generating a concise list suitable for preliminary investigations. The second approach adopts a cutoff based on the median, designed to offer balanced results, making it appropriate for a wide range of cases. Lastly, a more stringent cutoff based on Q3 in Table 3 is suggested for scenarios prioritizing the identification of highly probable escape mutations. While this may not capture all potential escape mutations, it is particularly valuable in therapeutic antibody development, where the focus is on mutations with significant implications for treatment efficacy. In scenarios where systems diverge from those evaluated in this study and lack tailored data for cutoff optimization, the established classic cutoff of 2.0 kcal/mol may provide a practical baseline for initial predictions.In our analysis, we emphasized positive results, since accurate prediction of escape cases is crucial. To evaluate the ability of the method to correctly identify positive cases (mutations that are more likely to be escape mutations due to binding disruption) while considering possible errors in the method (false positives), we compute the precision. This is the proportion of correctly identified escape mutations out of all positive predictions. For Case 1 systems, high precision is possible, provided that the \(\Delta G_{\text {bind}}\) predictions are accurate, as this is the case where Abs binding glycan-free monomeric epitope as in the experiments. In Case 2 systems, \(\Delta G_{\text {bind}}\) values could still serve as useful indicators of potential escape and high precision is possible, assuming accurate \(\Delta G_{\text {bind}}\) predictions. In Case 3 systems, however, \(\Delta G_{\text {bind}}\) is likely to be inaccurate, resulting in low escape prediction precision. Overall, precision varied across the systems: seven exhibited high precision (Precision \(>80\%\)), 14 showed moderate precision (Precision \(>50\%\)), and six had low precision (Precision \(<50\%\)). This highlights the method’s effectiveness in predicting escape mutations but also emphasizes the need for more accurate and efficient binding-free energy estimation methods. Systems classified under Case 1 showed precision ranging from 23 to 100%. Case 2 systems had a precision range of 45% to 90%. As expected, Case 3 systems displayed relatively low precision in predicting escape mutations, ranging from 25 to 43%. However, we also analyzed a more global metric of performance (Fig. S11 of the Supplemental Material). It showed a good classification ability of the model, with the area under the receiver operating characteristic curve (AUROC) values ranging from 0.75 to 0.91 for class 1 systems. This indicates that the model performs well in distinguishing between true positives and true negatives across different thresholds. However, while the ROC curve and AUCROC are useful metrics, they can provide an overly optimistic view of model performance, especially in the context of an unbalanced dataset. The ROC curve considers both true positive rates (TPR) and false positive rates (FPR) at different thresholds. In an unbalanced dataset where the majority of samples are negative, the FPR can be very low due to the high number of true negatives, making the ROC curve appear better than it actually is.Future studies should explore how interactions between mutations, known as epistatic effects, could cause non-additive impacts on binding affinity and immune escape. This research would be crucial for understanding the emergence of antibody escape mutations in clinical variants that contain multiple mutations. While our pipeline could be adapted for multiple mutations, the goal of this study was to evaluate its performance using empirical escape fraction data. The bulk of this data is single mutations, a case where MD+FoldX has been proven to excel at prediction for antibody-antigen systems compared to similar methods but not explored in depth. To our current knowledge, there has not been extensive benchmarking done on fast methods like FoldX to predict the effects of higher-order mutations and their ability to capture epistatic effects. Nor is there currently a large corpus of empirical data on systematic higher-order mutations. Additionally, the FoldX scoring function was primarily trained and tested on single mutations47. In addition, our study suggests that there is a need for methods to more accurately estimate binding free energy changes upon mutation at a reasonable computational speed. Furthermore, the impact of antibody structural variability and flexibility on computational predictions warrants further investigation, especially given the challenges posed by lower-resolution structures. Although FoldX provides more accurate free energies for high-resolution crystal structures (less than 2.6 Å)31, its effectiveness with cryo-EM-derived and low-resolution crystal structures requires thorough evaluation. While affinity plays a key role in antibody escape, it does not encompass all aspects of escape potential. Machine learning methods could significantly enhance the detection of escape mutations. These models trained using escape fractions from DMS experiments, could incorporate additional relevant features beyond mere binding affinity.
