High‐throughput molecular simulations of SARS‐CoV‐2 receptor binding domain mutants quantify correlations between dynamic fluctuations and protein expression

High-throughput molecular simulations of SARS-CoV-2 receptor binding1 domain mutants quantify correlations between dynamic fluctuations and2 protein expression.3 Victor Ovchinnikov1, a) and Martin Karplus1, 2, b) 4 1)Harvard University, Department of Chemistry and Chemical Biology, Cambridge, MA,5 021386 2)Laboratoire de Chimie Biophysique, ISIS, Universit´e de Strasbourg, 67000 Strasbourg,7 France8 (Dated: 26 September 2024)9 Prediction of protein fitness from computational modeling is an area of active research in rational protein design. Here, we investigated whether protein fluctuations computed from molecular dynamics simulations can be used to predict the expression levels of SARS-CoV-2 receptor binding domain (RBD) mutants determined in the deep mutational scanning experiment of Starr et al.1 Specifically, we performed more than 0.7 milliseconds of molecular dynamics (MD) simulations of 557 mutant RBDs in triplicate to achieve statistical significance under various simulation conditions. Our results show modest but significant anticorrelation in the range [-0.4,-0.3] between expression and RBD protein flexibility. A simple linear regression machine learning model achieved correlation coefficients in the range [0.7,0.8], thus outperforming MD-based models, but required about 25 mutations at each residue position for training. Keywords: molecular dynamics, deep mutational scanning, coronavirus, linear regression10 Introduction The COVID pandemic of 2019 resulted in the development of effective vaccines for the11 SARS-CoV-2 coronavirus. However, random mutations in the coronavirus spike protein, high viral trans-12 mission rates, and evolutionary pressure exerted by antibodies continue to cause the emergence of escape13 variants, for which the antibodies induced by the vaccines have reduced affinity.2Because the descendants of14 the original COVID-19 strain are likely to become endemic,3there is a need to develop advanced coronavirus15 vaccines capable of eliciting immunity to different strains, ideally including some that are yet to emerge.16 To address this challenge using rational antigen design, a detailed understanding of the determinants of17 coronavirus infectivity is needed. For rational vaccine design, it is important to be able to predict by com-18 putational methods how amino acid mutations in the receptor binding domain (RBD) of the coronavirus19 spike affect its stability, or binding affinity to the angiotensin converting enzyme receptor 2 (ACE2), and20 a)Electronic mail: ovchinn[email protected] b)Electronic mail: [email protected] 2 Simulation # Solvation #Mutants Duration (ns) cP(p-value) cS(p-value) Correct(%)† RBD simulations to model lMF I 1 cube 55 110 -0.36(7.8e-3) -0.42(1.5e-3) 62.6 2 shell 557 110 -0.22(1.7e-7) -0.29(5.5e-12) 59.5 3 shell 55 1010 -0.30(2.4e-2) -0.31(2.2e-2) 60.5 4 shell‡557 110 -0.34(2.1e-16) -0.35(2.3e-17) 61.8 RBD/ACE simulations to model lKD 5∗shell 557 110 -0.17(4.4e-5) -0.27(8.8e-11) 59.4 5§– – – 0.01(7.7e-1) -0.14(8.3e-4) 54.6 TABLE I. Summary of MD simulations performed in this study. For modeling lMF I, the correlations are between the protein residue RMSFs computed from simulation and the experimental lMF I; For modeling lKD, the correlations are between the protein residue RMSFs computed from simulation and the experimental lKD;†Indicates the proportion of all modeled pairs for which the modeled expression order within the pair matched the experimental order. ‡These simulations were performed with conformational restraints to the initial structure. ∗Both the RBD and ACE2 were simulated in complex, but the average residue RMSF was computed using only the RBD atoms. §The data used to compute the statistics are from simulation 5, but the average residue RMSF was computed over both RBD and ACE protein atoms. to antibodies. For example, stabilization of antigenic spikes in the prefusion conformation4has been used21 to improve vaccine efficacy. While the majority of information needed for rational protein mutagenesis22 derives from experiments, molecular dynamics (MD) simulations are increasingly able to provide dynamic23 information on temporal and/or spatial scales of protein motion that are smaller than the experimental24 resolution.25 Here, we studied whether conventional MD simulations can be used to estimate relative expression fitness26 of RBD mutants. We chose a random subset of 557 RBD sequence mutants from the total set of ≃54,000 in27 the deep mutational scanning (DMS) study of Starr et al.5that were approximately uniformly distributed28 in the level of expression, and performed MD simulations in triplicate for at least 110ns for each mutant.29 The DMS dataset provides the logarithm of the mean fluorescence intensity (lMFI), which was used in the30 assay as a proxy for protein expression.5 31 Results The results of the simulations are described below, and summarized in Table I and Fig. 1.323334 The specific mutants chosen are listed in a supplementary data file.35 First, we chose a small set of 55 mutants, immersed each mutant in a cubic box of explicit water solvent,36 and carried out 110-ns MD simulations in triplicate for each mutant. The initial 10ns were performed37 with harmonic restraints as part of simulation equilibration, and the remaining 100ns were used to compute38 statistics for correlation. As a proxy for mutant fitness, we computed (i) root-mean-square distance (RMSD)39 of the Cαprotein atoms from the simulation structures to the corresponding initial structures, (ii) average40 root-mean-squared fluctuation (RMSF) of protein Cαatoms, and (iii) average RMSF of the centers-of-mass41 (COMs) of all residues; these metrics were averaged over the triplicates. The three metrics were chosen42 based on the assumption that mutants that express poorly do so because of a destabilization of the folded43 3 4 6 8 10 12 log MFI (experiment) 0.8 0.9 1 1.1 1.2 <RMSF>(Ang) cP=-0.35519(7.8e-03) cS=-0.42143(1.5e-03) A Correct(%)=62.6 4 6 8 10 12 log MFI (experiment) 0.8 0.9 1 1.1 1.2 1.3 1.4 1.5 <RMSF>(Ang) cP=-0.21956(1.6e-07) cS=-0.28755(5.5e-12) B Correct(%)=59.5 4 6 8 10 12 log MFI (experiment) 0.5 1 1.5 2 2.5 3 <RMSF>(Ang) cP=-0.3038(2.4e-02) cS=-0.30844(2.2e-02) C Correct(%)=60.5 4 6 8 10 12 log MFI (experiment) 0.52 0.54 0.56 0.58 0.6 0.62 <RMSF>(Ang) cP=-0.33834(2.1e-16) cS=-0.34798(2.3e-17) D Correct(%)=61.8 FIG. 1. Scatter plots of RMSFs from simulations vs. logarithm of experimental MFI for various simulation sets A. Simulations performed in a water cube; B. Simulations in a water shell; C. 1µs-long simulations in a water shell; D. Simulations in a water shell with conformational restraints. conformation of the wild-type (WT/Wuhan) strain, which results in higher overall distance from, and/or44 larger fluctuations around, the WT structure. We note that this hypothesis predicts negative correlations45 between expression level and RMSDs or RMSFs, as is observed in most of our data. The three measures46 gave similar results (Tab. S1), so that in Table I and Fig. 1 we only show results with the RMSF-COM47 metric.48 The Pearson and Spearman (rank) correlation coefficients between the RMSF-COM and the lMFI are49 −0.36 and −0.42, respectively, indicating relatively low, but significant, correlation (respective p-values are50 7.8e-3 and 1.5e-3). Because in protein design it is often of interest to predict mutants with improved binding51 or fitness, we also computed the fraction of predictions that correctly rank any two mutants, relative to52 the lMF I data, at ∼63%. Since a completely uninformed model would predict this value to be 50%, the53 ‘enrichment’ obtained using MD simulations is 63/50 = 1.26, or 26%.54 4 0 20 40 60 80 100 t(ns) -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 cP, cS A 0 20 40 60 80 100 t(ns) -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 cP, cS B 0 200 400 600 800 1000 t(ns) -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 cP, cS C 0 20 40 60 80 100 t(ns) -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 cP, cS D FIG. 2. Time-evolution of the Pearson and Spearman correlation coefficients between RMSDs and lMF I (red and green *, respectively), and the average corelation between the RMSDs of simulation replicates (black ◦and blue ▽, respectively), for various mutant sets. The legends for panels A–D are as in Fig. 1. The above results were obtained from MD simulations in a cubic box of explicit solvent, which con-55 tain ∼70K atoms, and simulation speeds of ∼230 ns/day (on a workstation equipped with an NVIDIA56 RTX2080Ti graphics processor and a Ryzen 9 3900X CPU). Therefore, about 12 hours of simulation were57 needed for each mutant. Although simulations can be run in parallel on a supercomputer, it is desirable58 to explore more approximate MD methods that reduce simulation cost. To do this, we performed the59 next round of simulations using the shell solvation model developed by us.6With this model, each mutant60 protein is placed in a shell of solvent with thickness 8.5˚ A, and solvent evaporation is prevented using a61 half-harmonic potential. Further details can be found in Methods and in Ref. 6. The resulting simulation62 systems consisted of ∼10K atoms, and yielded ∼700 ns/day, i.e.,>3×faster than the simulation in the63 cubic box. The increase in speed was due mainly to the reduced number of atoms, but also to the related64 fact that a shorter nonbonded cutoff (9˚ A) could be used. We note that the solvation shell model has known65 artifacts, such as artificial surface tension, and a small but controllable departure from equilibrium.6 66 With the increased speed, reduced storage requirements, and access to a supercomputer, we were able67 to generate more reliable statistics. We chose 557 mutants, which included the previous set of 55, and68 were uniformly distributed on the mutational landscape. For these 110-ns solvation shell simulations, the69 Pearson and Spearman correlation coefficients between the RMSF-COM and the lMF I were −0.22 and70 −0.29, respectively, corresponding to a decrease in the predictive ability relative to the cubic solvation cell71 (see Fig. 1). The p-values of 1.6e-7 and 5.5e-12, respectively, were much lower than in the previous case,72 due to the larger simulation set size. The fraction of correct predictions was 60%, which corresponds to an73 enrichment factor of 20% over a random model. Thus, the use of a simplified solvation treatment resulted74 in deterioration of the fitness predictions.75 To determine the effect of MD simulation length on the prediction quality, we extended the solvation shell76 simulations of the 55 mutants described before to 1µs, and computed the RMSDs of the RBD structures for77 the three sets of simulations (see Fig. 2) vs. time. Although one might expect that longer simulation lengths78 5 would result in better predictions, e.g., because the additional time would facilitate structural relaxation79 toward low free energy conformations,7the results are generally at odds with this expectation. This may80 be related to the work of Raval et al.8, who found that long MD simulations do not necessarily improve81 the refinement of protein structures built by homology modeling. For all three sets of simulations, the82 highest correlations (ignoring the negative sign) occur for t≤∼ 50ns,i.e. relatively early in the simulation83 trajectories (Fig. 2A-C).84 The best correlation to the lMFI data corresponds to the fully solvated (cubic box) simulations at85 t∼55ns, at which time the Spearman coefficient is ∼ −0.5. For the solvation shell simulations, the best86 correlation is found at t≤25ns, with the Spearman coefficient of ∼ −0.33 and −0.44 for the sets of 110ns87 and 1010ns simulations, respectively. For the other simulations, the correlations tend to worsen gradually88 with time, though not necessarily monotonically.89 While the physical origins for the gradual deterioration in the correlations over the simulation times are90 unclear, the observation that the motions most predictive of fitness occur early in the simulations are of91 practical use. First, it indicates that very long simulations are not necessarily beneficial. Second, the early92 simulation samples are generally in closer proximity to the starting simulation structures (here, the Xray93 crystal structure of the WT RBD9to which point mutations were applied). Sampling in the vicinity of the94 starting structure can be increased with the use of harmonic restraints.95 To test whether such restraints would improve the fitness predictions, we repeated the 557 simulations96 using the shell model with a weak best-fit harmonic restraint on the RMSD between the simulation and the97 starting structure, applied to the heavy atoms for the entire 100ns. The results of this set of simulations98 show an improvement in both Pearson and Spearman correlations (−0.43 and −0.35, respectively), and99 in the enrichment (61.8%) over the unrestrained simulations. The temporal evolution of the correlation100 between the RMSD and lMFI decreases slowly, but monotonically in the restrained simulations, indicating101 that the longer restrained equilibration increases the number of predictive simulation samples.102 Overall, the present MD simulations of RBD mutants show a modest ability to predict expression fitness103 (absolute Spearman correlation of 0.3−0.4 and an enrichment factor of 20%–25% over a random model. It104 is noteworthy that each prediction requires on the order of 100ns of simulation time, which can be generated105 in ∼4 hours on a single GPU workstation. Moreover, different simulations can be run in parallel.106 The DMS dataset of Starr et al.5also contains measurements of the dissociation constant KDfor the RBD107 mutants and the ACE2 receptor. To check how the quality of our fitness predictions from MD simulation108 would compare to predictions of KD, we performed unbiased simulations in triplicate of 559 mutant RBDs109 6 in complex with ACE2 using the solvation shell model, and computed the correlations between the average110 RMSFs and the negative logarithm of KD(lKD). We note that these 559 mutants are different from the111 557 chosen for computing expression correlations, because they were selected to sample uniformly the space112 of experimental lKDrather than lMFI (see Methods). The results are shown in Table I and Fig. S1.113 These simulations were substantially more computationally costly because they involved simulating ACE2114 in addition to RBD, and because our focus in this study is on fitness of the RBD. Therefore, we did not115 perform any simulations in the cubic box.116 The computed correlations between lKDand the RMSFs are significantly lower than those for the117 lMF I simulations (Table I and Fig. S1). When the fluctuations are averaged over both the RBD and118 ACE2 residues, the Pearson and Spearman correlation coefficients are 0.01 and −0.14, respectively. The119 p-value for the latter is 8.3e-4, indicating that the the relative ranking between mutants is significant, de-120 spite providing only a small advantage over a random model (an enrichment of 9.2%). If the RMSFs are121 averaged only over the RBD, the respective correlation coefficients improve to −0.17 and −0.27, and the122 fraction of correctly ordered predictions increases to 59.4% (an enrichment of 18.8%).123 A possible explanation of the correlations between RMSFs and lKD’s is that higher RBD fluctuations124 reflect a larger entropy penalty that has to be overcome to achieve strong binding. The fact that the125 correlation improves when only the RBD residues are considered suggests that such an entropy loss is126 associated specifically with ordering the mutant RBD residues. To investigate this hypothesis, we computed127 an approximation to the conformational entropy of the RBD in the ACE2-bound and unbound states using128 quasiharmonic analysis10,11 in coarse-grained residue coordinates (see Supporting Material for details). It is129 known that quasiharmonic entropy typically overestimates the true conformational entropy (see, e.g., Tyka,130 Clarke, and Sessions 12). However, because the the RBD mutants have nearly identical structures, the131 entropy values corresponding to different mutants are expected to have similar systematic errors, allowing132 a meaningful relative comparison. Scatter plots of the temperature-scaled RBD entropy in the bound133 and unbound states vs. lKDare shown in Fig. S2. The correlations between the negative lKDand the134 computed entropies were found to be in the range [-0.28, -0.24], with p-values of 7e-8 or less, indicating135 modest but statistically significant correlations between the ACE2-RBD binding free energy and entropy136 of the RBD, both in the bound and the unbound states (see Fig. S2). These calculations therefore support137 our interpretation of the anticorrelation between RMSFs and lKD. We note that it is technically possible138 to compute free energies of binding between the RBD and ACE2 exactly within the framework of classical139 statistical mechanics and restrained MD simulations.13 Since such simulations would very likely require long140 7 simulation times (e.g., several µs per complex), they are currently too computationally costly for routine141 high-throughput screening.142 As done in the above case for modeling fitness, we also computed the correlation of lKDand the RMSD143 between the simulation and starting structures of the RBD/ACE2 complexes as a function of time (Fig. S3).144 The best correlation, corresponding to a Spearman coefficient of ∼ −0.14, was at ∼42ns of simulation time.145 Overall, unbiased MD simulations of RBD/ACE2 complexes using the solvation shell model to accelerate146 the calculations, resulted in smaller correlations, relative to the lMFI data, between residue fluctuations of147 the RBD and the experimental lKDdata. Although simulating the complexes in a cubic box could increase148 the correlations, as in the lMFI cases above, we have not done so because of the high computational cost149 involved (∼31K atoms for the shell model vs. ∼255K atoms for the cubic box).150 Fitness prediction methods grounded in physics-based simulations, such as the MD approach evaluated151 here, have the advantage that they generally do not require tuning to fit a specific protein, since they are152 derived from general principles, e.g., an energy function that is applicable to any protein without highly153 specialized chemical modifications.14–17 Their main disadvantages are that (1) some structural information154 is required, e.g., the Xray crystal structure used here, (2) the predictions may not be very accurate, e.g., if155 they are based on statistical samples that require long simulations to converge (though this is not observed156 here, as described above),18 or (3) they fail to incorporate some important factor that affects fitness,157 e.g., the experimental protein expression system has a protease that tends to cleave particular amino-acid158 sequences,19 independently of the thermodynamic stability of the folded protein.159 By contrast, prediction methods based on machine learning (ML) do not require protein structures (al-160 though some algorithms benefit from them20,21), but rely on a preexisting data set for model training.22,23 161 If the training data set is sufficiently representative of the set for which predictions are desired, ML models162 can automatically correct for systematic biases.163 To compare the prediction quality of our MD-based fitness metrics with that from an ML model, we164 parametrized a simple sequence-based linear regression model to compute either RBD mutant fitness (lMFI)165 or binding affinity to ACE2 (lKD), as described in Methods.166 The main results of the ML model are shown in Figs. 3 and S4 for lMF I and lKD, respectively, and167 summarized in Table II. With 50% of the data split between the training and testing sets, the Pearson and168 Spearman correlation coefficients for modeling lMFI are ∼0.75 and ∼0.78, respectively. The error bars169 for the correlation coefficients are generally less than ∼0.005, and represent twice the standard deviation,170 or 95% confidence limits. They are shown in Fig. 3C, in which the various model performance metrics171 8 4 6 8 10 12 experiment (log10[MFI]) 4 6 8 10 12 model Training set (50%) cP=0.758 cS=0.78639 RMSE=0.82135 A Correct(%)=78.9 4 6 8 10 12 experiment (log10[MFI]) 4 6 8 10 12 model Test set (50%) cP=0.7486 cS=0.77672 RMSE=0.8366 B Correct(%)=78.9 2 4 6 10 30 50 70 90 % data used for training 0 0.5 1 1.5 2 2.5 C cP TRAIN cP TEST cS TRAIN cS TEST RMSE TRAIN RMSE TEST qcorrect TRAIN qcorrect TEST FIG. 3. Results of the linear regression model to predict lMF I expression. A & B: scatter plots of model predictions vs. experiment; A: training set; B: testing set; 50% of the data was assigned to each set. In A, ten equispaced contours of the 2D joint probability density function (PDF) of modeled and experimental lMF I’s are drawn, spanning the range of the PDF; the majority of the lMF I values fall near ∼10 which cannot be seen from the scatter plot alone because of the large amount of the data. C. Model performance metrics as a function of the percent data used for training; for each percent value, the remaining data were used for testing; qcorrect represents the fraction of correct predictions of pairwise rank (see text). The training data were chosen randomly, and the error bars in C correspond to twice the standard deviation computed from ten training realizations. Training set Test set Model cPcSRMSE Correct(%) cPcSRMSE Correct(%) lM F I 0.76 0.79 0.82 79.5 0.75 0.78 0.83 78.9 lKD0.73 0.77 1.4 78.8 0.72 0.76 1.4 78.4 TABLE II. Summary of linear regression model results. The training set contained 27192 sequences, corresponding to 50% of the DMS dataset, chosen randomly with uniform probability; the remaining 27192 sequences comprised the test set. p-value for the correlation coefficients are not explicitly shown, because they were zero to within machine precision, which was due to the large number of samples. The statistical uncertainties in the correlation coefficients are indicated in Fig. 3. RMSE denotes the root-mean-square error in the units of experiment. are displayed as a function of the percentage of data used for training vs. testing. The similarity of the172 correlations for the training and testing sets (as long as the training sample size is above 10% of the total173 sample size) indicates that the model shows no siginficant overfitting to the training data, due in part to174 the simplicity of the model.175 Although far more sophisticated and accurate models can be parametrized, such as natural language,22 176 gaussian process,20 or deep learning23 models, our present objective is a qualitative comparison of the177 accuracy of MD-based vs. ML-based models, for which the modeling used here is sufficient. If enough178 training data are available, which in our case corresponds to about 5K mutations, or about 5000/200=25179 mutations for each residue in the spike RBD (ca. residues 331–531), the present ML models provide superior180 predictions of both fitness and binding affinity.181 We note that the model used here is linear in residue space, which implies that the contribution to the182 lMF I or lKDof each residue is additive. However, residues, and therefore mutations, obviously interact183 with each other via the protein structure, which is manifested through epistatic shifts,1,5 and can be modeled184 by parametrizing additional residue interaction weights24 (couplings in the Potts model formalism),25 rather185 9 than only the residue weights, as done here. This obvious drawback of the present model is illustrated in186 Figs. S5 and S6, which shows that the ML model performance is best for single-residue mutants, and187 gradually deteriorates for higher-order mutants. However, the overall correlation between experiment and188 model of >0.7 for both fitness and binding affinity indicates that even simple linear regression models189 outperform molecular dynamics-based criteria, if sufficient model training data are available.190 Conclusion We have investigated the accuracy of two related structural metrics obtained from classical191 MD simulations, RMSD and RMSF, in predicting changes in (1) fitness of residue mutants of the SARS-192 CoV2 RBD and (2) binding of the RBD to the ACE2 receptor. We quantified the predictive ability of MD193 using Pearson and Spearman correlations between RMSFs and RMSDs computed from >0.7ms of total194 simulation time with the deep mutagenesis data of Starr et al.1(Table I).195 The computed correlations with the logarithm of the experimental mean fluorescence intensity are modest,196 in the range 0.3 – 0.4, with the highest correlations observed for the RBDs solvated in a cubic box. The197 less computationally intensive simulations with RBDs solvated in a solvent shell with restraints to the198 starting coordinates of the protein heavy atoms produced a Spearman correlation of ∼0.35. Interpreting199 the simulation data in terms of the percentage of correct ranking of all mutant pairs, we observe at most200 62.6% correctness, which indicates a modest advantage over a random prediction of 50% (a 25% enrichment201 over the random model). We emphasize that the correlations to MD data computed here, although modest,202 are statistically significant, as indicated by the p-values that accompany all correlation values in Table I.203 In particular, the simulations of the larger mutant set (557 mutants) yielded p-values that were all below204 2e-7, with some p-values being zero within the limit of machine precision (∼2e-16). However, although205 the statistical significance of the correlations is robust, demonstrating the robustness required testing the206 large number of mutants. A practical consequence is that obtaining useful predictions with the RMSD and207 RMSF metrics could require a large number of samples. Thus, these metrics are unlikely to be helpful in208 situations where, e.g., a small set of mutations are to be tested.209 The correlations between the fluctuations in the RBD/ACE2 complex and the experimental lKD, are210 lower, being optimal at −0.27 for the Spearman correlation, when only the RBD is considered (i.e. ACE2211 fluctuations are ignored). Because lKDis linearly proportional to the RBD/ACE2 binding free energy212 (bFE), we speculated that the contribution to the bFE reflected in the RMSFs corresponds to the entropic213 cost of ordering the RBD that is necessary for binding. This hypothesis appears to be supported by quasi-214 harmonic analysis of the MD trajectories in coarse-grained residue coordinates, which showed a negative215 correlation between classical qhasiharmonic entropy of the RBD in isolation, as well as in complex with216