key: cord-1028969-r425ms9q authors: Bhadane, Rajendra; Salo-Ahen, Outi M. H. title: High-throughput molecular dynamics-based alchemical free energy calculations for predicting the binding free energy change associated with the common mutations in the spike receptor-binding domain of SARS-CoV-2 date: 2022-03-08 journal: bioRxiv DOI: 10.1101/2022.03.07.483402 sha: c69efe45d80a22f7643ff5f5414cf2eff0ab4ac5 doc_id: 1028969 cord_uid: r425ms9q The ongoing pandemic caused by SARS-CoV-2 has gone through various phases. From the initial outbreak the virus has mutated several times, with some lineages showing even stronger infectivity and faster spread than the original virus. Among all the variants, beta, gamma, delta and the latest (omicron) are currently classified as variants of concern (VOC) while the remaining are labelled either as variants of interest (VOI) or variants under monitoring (VUM). In this work, we have focused on the mutations observed in important variants, particularly at the receptor-binding domain (RBD) of the spike protein that is responsible for the interactions with the host ACE2 receptor and binding of antibodies. Studying these mutations is particularly important for understanding the viral infectivity, spread of the disease and for tracking the escape routes of this virus from antibodies. Molecular dynamics (MD) based alchemical free energy calculations have been shown to be very accurate in predicting the free energy change due to a mutation that could have a deleterious or a stabilising effect on the protein itself or its binding affinity to another protein. Here, we investigated the significance of six commonly observed spike RBD mutations on the stability of the spike protein binding to ACE2 by free energy calculations using high throughput MD simulations. For comparison, we also used other (rigorous and non-rigorous) binding free energy prediction methods and compared our results with the experimental data if available. The alchemical free energy-based method consistently predicted the free-energy changes with an accuracy close to ±1.0 kcal/mol when compared with the available experimental values. As per our simulation data the most significant mutations responsible for stabilising the spike RBD interactions with human ACE2 are N501Y and L452R. In the current pandemic caused by SARS-CoV-2(1) , the study of viral mutations can give us understanding on the infectivity, pathogenesis, and drug resistance of the virus as well as insight into vaccination efficiency. Further, it may help tracking immune escape routes and consequent emergence of new diseases. During the second wave of COVID-19, SARS-CoV-2 mutations drew much more attention than those of any other virus studied till date. However, it is important to point out that viral mutation is not an unusual phenomenon. For any organism, the mutation rate is the change in genetic information that passes from one generation to the next. The viral mutation rate involves a process that starts from the viral particle's attachment to the host, its entry, release of genetic material, replication of viral particles, assembly and release or escape from the host cell (2) . Viruses show a higher mutation rate than other organisms due to their short generation time and large population size. The mutation rate of viruses also depends on the type of virus. RNA viruses show higher mutation rates than DNA viruses, e.g., HIV virus shows extremely high mutation rate due to its higher replication rate (3) . When HIV and SARS-CoV-2 viruses are compared, despite both being RNA viruses, mutations of SARS-CoV-2 are slower than in HIV. There are also many other factors that could influence the mutation rate of a particular virus. It is also important to understand that not every mutation is beneficial to a virus as most mutations are deleterious, leading to the instability of the viral particles (3) . The highly deleterious mutation gets rapidly purged while neutral or mildly deleterious stays. These mild mutations impact the virus phenotype and could be beneficial for pathogenicity, infectivity, transmissibility and/or antigenicity of virus particles (4, 5) . Since the first outbreak of SARS-CoV-2 in Wuhan, China, the virus has mutated to various lineages. Among these so-called Pango (Phylogenetic Assignment of Named Global Outbreak) (6) lineages some are classified as variants of concern (VOC), variants of interest (https://www.ecdc.europa.eu/en/covid-19/variants-concern). Of special interest in the present work are the mutations occurring in the receptor-binding domain (RBD) of the SARS-CoV-2 spike protein (8) , (9) since that domain plays a major role in the entry of the virus into the host cell by binding to the human angiotensin-converting enzyme 2 (hACE2). For example, a sub lineage of the alpha variant that was first detected in the U.K in December 2020 carries E484K and N501Y mutations in the spike RBD and is 30-50% more infectious than the original virus (10) . The exponential growth in computational power has facilitated the utilisation of in silico methods. Computationally expensive free energy calculations that were previously thought to be impossible can now be performed in a few hours. For an ideal free-energy prediction method, the predictive accuracy should be of the same range as that obtained by experimental methods. In a so-called "alchemical" transformation, an amino acid is transformed from one state to another via a non-physical pathway. This transformation can be either reversible or irreversible and can be performed by an equilibrium or non-equilibrium method, respectively(11). The major advantage of this approach is the prediction accuracy that matches that reached in experiments (12) . In this work, we performed alchemical free energy calculations, utilising high throughput molecular dynamics (MD) simulations to predict the binding free energy change (with respect to binding partner ACE2) upon six spike RBD mutations commonly observed in current VOC and early VOI lineages. We performed altogether 24 equilibrium MD simulations and 1200 non-equilibrium MD simulations to ensure exhaustive conformational sampling of each of the mutants. We also carried out comparative studies with less rigorous computational methods and discuss the results in the light of available experimental data and/or other computational studies. The results of this comparative study may shed further understanding on the impact of these mutations on infectivity and spread of SARS-CoV-2. Most importantly, we show that the relatively fast free energy calculations may provide a reasonable alternative for expensive experimental mutation studies. Selection of mutations for the study: As per the Pango lineage definition, till date, there are more than 1200 lineages found for SARS-CoV-2(13). The B.1 lineages with the D614G mutation in the spike protein observed in early-2020 gave rise to more than 20% increase in transmissibility compared to the original strain of SARS-CoV-2 Wuhan, China. Other variants also appeared, including many lineages with E484K as the common mutation that was found to reduce the binding affinity of antibodies (14) (15) (16) (17) (18) . Figure 1 shows the rough timeline of VOC, VOI and VUM lineages with prominent mutations in the RBD region of the spike protein. The envelope of SARS-CoV-2 has homotrimeric spike glycoproteins that consist of a S1, S2 and a short cytoplasmic domain. To interact with the human angiotensin-converting enzyme 2 (hACE2) the virus utilises the amino acids 331-524 of the S1 domain, that is, the receptor- (formed by antiparallel β1-4 and β7 strands) with short connecting helices (α1-α3) and loops. Between β4 and β7 strands there is an external subdomain (residues 438-506) that is formed by α4 and α5 helices and a flexible loop that connects two short β strands (β5 and β6). This subdomain is also called as the receptor-binding motif (RBM) region since it makes direct contacts with the hACE2 receptor. Except for the R346K mutation that is in the core, right after the α1 helix, all other selected mutations are in the RBM region. L452R is at the β5 strand, T478K and E484K/E484Q mutations are located between the β5 and β6 strands while N501Y precedes the α5 helix (50) . Of note, only residues at positions 501 and 484 have direct interactions with hACE2. Figure 2 shows the location of the mutations selected to be studied here. Currently, there are various tools available for predicting the functional and/or structural impact of mutations on a protein or protein-protein complexes. In this work, we have utilised Schrödinger's Prime/MM-GBSA-based residue scanning (Prime, Schrödinger, LLC, New York, NY, 2021)(51,52) as well as CUPSAT (53) (54) (55) and FoldX(56-59) tools for predicting stabilising or destabilising effect of the selected mutations on spike RBD-hACE2 complex. In addition, we calculated the effect of the mutations on the stability of the spike RBD alone. The major advantage of these calculations is that they are fast and computationally non-costly. CUPSAT and FoldX tools are also available for free. The downside of these predictions is that they are less accurate and, thus more qualitative than quantitative. The predicted stabilising and destabilising effects of the selected spike RBD mutations on the RBD binding interaction with hACE2 and the stability of RBD are summarised in Table 2 . A ΔΔG value ≥ |1.0| kcal/mol was counted as a significant change ("hotspot") for CUPSAT/FoldX predictions and ≥ |3.0| kcal/mol for Prime predictions (60) . Since all the mutations are found in the viable virus variants, it can be assumed that they should not destabilize the actual spike RDB structure drastically. The prediction tools consistently suggest that mutations L452R and T478K are stabilising the binding interaction. Predictions for the other mutations are inconsistent with these tools. According to the CUPSAT prediction for the spike RBD alone, none of the mutations affects the overall stability of the spike RBD protein significantly (only R346K and T478K are predicted to have a slight destabilising effect while the others are stabilising the protein). The Prime/MM-GBSA-based residue scanning results for some mutants (especially N501Y) were observed to significantly depend on the used refinement cut-off distance for the sidechain predictions. The side-chain prediction with backbone sampling or simple backbone minimization (see Methods) was performed using different cut-off values for including the nearby residues in the refinement. Table 2 shows the affinity change (ΔAffinity) predicted by Prime using the backbone sampling with cut-off of 6 Å. This cut-off gave the most reasonable ΔAffinity value for the N501Y mutant as that mutation is known to increase the spike binding affinity to hACE2 (see the results with other cut-off values in Supporting Information, Table S1 ). On the other hand, the change in spike RBD stability was predicted using minimization of the mutated residue as the refinement method (cut-off 0 Å) to avoid possible errors caused by allowing the protein structure to move more than necessary. According to these results, all but the E484 mutants are predicted to increase the spike RBD binding affinity for hACE2 (T478K and N501Y showing the largest effect on the affinity), while E346K and E484K are predicted to have the greatest effect on the spike RBD stability among the mutants. Of note, the E484K mutant has been experimentally shown to exhibit a slightly lesser thermal stability than the WT whereas N501Y mutant is as stable as the WT (43) . Almost all FoldX predictions for the mutants were less than 1.0 kcal/mol in quantity, suggesting an insignificant effect on the hACE2 affinity or spike RBD stability. However, inconsistent with experimental evidence, N501Y was predicted to reduce the spike RBD binding affinity to hACE2 by both FoldX and CUPSAT, and even by Prime with insufficient backbone sampling (Supporting Information, Table S1 ). The FoldX optimised structures of the wildtype (WT) and mutant spike RBD bound to hACE2 were used as input for conventional MD simulations. Root mean square deviation (RMSD) and Information, Figures S3 and S4 ). In addition, the hydroxyl oxygen of Tyr501 may form a hydrogen bond/a dipole-ion interaction with Lys353 or its phenyl ring engages Lys353 in pication interactions (Figures 4a, b) . Further, the trajectory analysis shows that L452R mutation increases hydrogen bonding as well as Coulombic and vdW interactions. This is likely due to the introduction of the charged residue (Arg452), which causes a change in the local environment (Figures 4c-4e) , thus favouring increased electrostatic attraction and/or new polar interactions between the proteins. Interestingly, Figure 4c shows how Arg452 of the spike RBD mutant interacts with the backbone carbonyl of another RBD residue, Tyr449, thus indirectly providing additional stability to the interaction between Tyr449 and Asp38 of hACE2. E484K mutation shows a slight increase in H-bonding and vdW interaction energy but the total interaction energy decreased due to the decrease in the Coulombic interaction energy as compared with the WT. The opposite charge seems to be non-optimal for the interaction as the salt bridge between the WT Glu484 and hACE2 Lys31 is lost upon the mutation. In sum, according to the conventional MD simulations, L452R stabilises the spike RBD-hACE2 interactions while all the other studied mutations seem to destabilise the interactions somewhat. We did not observe a consensus between the conventional MD simulation results and those from the non-rigorous methods except for L452R. In addition, limited sampling during the MD simulations might affect the ability to predict the stability of the interactions upon mutation. Hence to increase the conformational sampling of the mutants, we carried out high-throughput MD simulations, utilizing non-equilibrium alchemical free energy calculations to investigate the effects of these mutations on hACE2 binding. To perform these calculations, hybrid topology files containing information on both WT and the mutated residues were created. Among the selected mutations, only R346K and N501Y are non-charge-changing mutations i.e., before and after mutation the change in the charge for these mutants remains equal to zero. This is in contrast to charge-changing mutations (L452R, T478K, E484Q and E484K) where after mutation the net charge change does not equal to zero. Estimating the free energy difference for such charge-changing mutations required an additional step of performing corrections in the charge and was achieved by a so-called analytical post-simulation charge correction method (61) . Table 3 contains the binding free energy calculated before and after the charge corrections (see also According to the free energy calculations, N501Y is the most significant stabilising mutation affecting the binding of spike RBD to hACE2. This is followed by L452R that shows a mild to moderate stabilising effect. The R346K mutation also shows a weak stabilising effect. The remaining three T478K, E484Q and E484K mutations are predicted to be weakly or negligibly destabilising. We then compared the predicted binding affinity of three mutations with experimentally reported binding affinities ( Figure 5 , Table 1, Supporting Information, Table S3 ). Table 2 ) have been converted to correspond to the convention that negative values are stabilizing and positive values destabilising. The results show that the MD-based alchemical free energy calculations accurately and consistently predicted the affinity when compared with the other prediction methods. The Prime-predicted values also look relatively accurate in Figure 5 , but we need to remember that its reliability for some residues may drastically depend on the choice of a suitable distance cutoff during backbone sampling, which here was guided by the available information from experimental studies. All methods were successful in predicting the effect of the L452R mutation; experimental range of G from -0.50 to -0.35 kcal/mol ( Table S1 ). The viral mutation is an ongoing process and can evolve in response to selective pressures. The SARS-CoV-2 mutation pattern on spike protein observed till date suggests that during the first phase of COVID-19 infections, mutations with increased infectivity overcame the early variant. This can be seen for example from the disappearance of B.1 and B. There are still no reports on the effect of the R346K mutation on the binding interaction with hACE2, but it has been shown that it helps in escaping class 3 antibodies (65, 66 More extensive conformational sampling (i.e., longer MD simulations with several parallel runs) might have given better interaction energies for that mutant, which would also be consistent with the Prime residue scanning results that emphasize the need for sufficient backbone sampling in the environment of the mutation. Notably, the rigorous alchemical free energy calculation method consistently predicted the free energy changes of all mutations with an accuracy close to ±1.0 kcal/mol when compared with the available experimental values. In sum, the results of this study are consistent with the experimental data on that both N501Y and L452R mutations stabilise the binding interactions between spike RBD and hACE2. On the other hand, neither E346K in the spike RBD core nor the mutations E484Q, E484K and T478K in the RBM offer significant advantage on binding to hACE2. Our study shows also that the MD-based alchemical free-energy calculations work accurately and consistently for different types of mutations (including the more challenging charge-changing mutations) in different structural environments. Based on the reports on increased infectivity of alpha and beta variants and the sudden spike in infections in India in April-May 2021, the key mutations in the spike RBD of these variants were selected for the study. The 3D structure of the spike receptor binding domain (RBD) that is complexed with the angiotensin-converting enzyme-2 (ACE2) (PDB ID: 6LZG(49)) was retrieved from the Protein Data Bank (PDB) (70) . The selected structures were processed using the Protein Preparation Wizard of Maestro (Schrödinger Release 2021-2 New York, NY, USA). Briefly, the missing hydrogen atoms were added to the structure and the hydrogen bond network was optimized using PROPKA at pH 7.0. All water molecules beyond 3 Å were removed and a restrained minimization was conducted using the OPLS4(71) force field with a convergence criterion of 0.3-Å root-mean-square deviation (RMSD) for all the heavy atoms. FoldX. The 3D structure of the spike RBD in complex with hACE2 (PDB ID: 6LZG) was stored in the FoldX directory. Using the RepairPDB command, residues having bad torsion angles, steric clashes, or high total energy were identified and repaired. The BuildModel command was then used to create the mutant PDB files and the corresponding WT structures. Finally, using the AnalyzeComplex command the interaction energy between the spike RBD and hACE2 was computed for each mutant and the corresponding WT structure, after which the free energy difference between the WT and the mutants was calculated according to ΔΔG = ΔG mut -ΔG WT (72) . Here, we calculated the stability of the spike RBD separately using just the RBD structure (not the complex) as the input. The output structures from the residue minimization protocol with cut-off 0 Å were used as input for performing the conventional MD simulations. The OPLS4 force field was applied for these calculations. The WT and mutated structures of spike RBD in complex with hACE2 were submitted to a The simulation systems were prepared using the System Builder tool of the Desmond module. The single point charge (SPC) water (74) was chosen as the explicit solvation model. Each system was neutralized using an appropriate number of Na + or Clcounter ions. An orthorhombic simulation box, with Periodic Boundary Conditions (PBC) and a 10-Å buffer space between the solute and the box edge, was used for each system. The simulation systems were minimized and equilibrated before the production simulation in a stepwise manner. After the system relaxation, the production simulation was performed in the NPT ensemble for 100 ns, using a reversible reference system propagation algorithm (RESPA) integrator (73) . The temperature (300 K) was set using the Nosé-Hoover chain thermostat (75) (76) (77) , with a relaxation time of 1.0 ps. The pressure was set at 1.01325 bar with the Martyna-Tobias-Klein barostat(78), using isotropic coupling and a relaxation time of 2.0 ps. Long-range interactions were handled using the U-series method (79) and for short range interactions, a cut-off radius of 9.0 Å was used. The MD trajectories were analyzed using the Maestro built-in Simulation Interactions Diagram and simulation event analysis tool and Microsoft Excel360. Alchemical relative free energy calculations exploit and avoid the need to simulate binding and unbinding events by making use of the fact that the free energy is a state function and can be explained as in Equation 1 (80, 81) ). According to Equation 1 and Figure 3 , the difference in the free energy for example for the L452R mutation can be calculated by performing two independent MD simulations and estimating the free energy change from Leu-452 to Arg-452 in the presence or absence of hACE2, i.e., (∆ ) and (∆ ), respectively. In an alchemical simulation this change from state A to state B is performed by introducing an alchemical progress parameter, lambda (  ) that controls the potential energy function ( ;  ) (q is the x,y,z coordinates of the simulation system). This is achieved by generating and simulating a hybrid topology of amino acids composed of the atoms that represent both A and B states. A subset of the energetic terms in ( ;  ) are modulated such that at  the atoms representing state A are activated and those representing state B are non-interacting "dummy atoms" and vice versa at  . In this way the alchemical transformation is performed and the free energy of mutating A to B in any environment (e.g., bound or unbound state) is computed as per the following Equation 2 (80, 82) . To achieve a good phase-space overlap between the two thermodynamic states, multiple Preparation of Hybrid Topology. The PMX webserver (85, 86) was used to generate the hybrid topology files for the selected amino acids. The pdb2gmx command was executed to reassign correct hydrogen atoms in the input coordinate files. The Amber99SB force field was used (12, 87) and the output files containing the topology information were generated in a Gromacs readable form. Equilibrium and non-equilibrium MD simulations. All MD simulations (equilibrium and non-equilibrium) were carried out on Puhti HPC cluster(88) with GROMACS-2020.5 (89) (90) (91) (92) (93) (94) . The hybrid topologies generated by the PMX webserver were used as an input. For each mutation, we prepared an individual simulation box and performed forward and backward transition with two simulation systems representing the protein in state A and state B (λ= 0 and λ=1, Figure 3 ). Amber99SB force field and the Joung and Cheatham ion parameters were used (95) . Each state was solvated using a dodecahedral water box with simple point charge 216 (SPC216)(74) 3-point water model and neutralized with Na + and Cl − ions at a 0.15 M concentration. Each simulation system was then energy-minimized for 100000 steps using the steepest descent algorithm (96) . This was followed by 1-ns NPT ensemble simulation with positional restraints, using stochastic leap-frog integrator and isotropic Berendsen pressure coupling (97) . The final production MD simulations were performed using the stochastic leapfrog integrator for 50 ns (98) . All bond lengths were constrained using the LINCS algorithm (99) . The constant pressure and temperature condition was maintained using the Parrinello−Rahman pressure coupling(100) at 1 atm and temperature coupling at 298.15 K. A 2-fs time step was used, and the trajectory was recorded at every 10 ps. To treat the long-range electrostatics, Particle Mesh Ewald (PME) method with a direct space cut-off of 1.1 nm and Fourier grid spacing of 0.12 nm were used (101, 102) . The relative strength of the Ewald-shifted direct potential was set to 10 -5 . Van der Waals interactions were smoothly switched off between 1.0 and 1.1 nm. To perform the alchemical transformations after the equilibrium MD simulations, fastgrowth non-equilibrium simulations were performed to estimate the binding free energy (ΔΔG) value (12, 82, 83) . For this, each equilibrium simulation was used as input for a corresponding nonequilibrium run; the first 25 ns of the trajectory were omitted, and the last 25 ns of the simulation were used to generate 50 snapshots from every 100 ps. The non-equilibrium simulation was performed for a total of 5 ns with the lambda value changing continuously in forward and reverse direction from state 0 to 1 and vice-versa at each time step with a frequency of transition of 2x10 -7 . The standard errors of the ΔG estimates were obtained by bootstrapping and reflect the variability in the datasets (trajectories) analyzed. Among the six mutations studied here, four mutations are charge transformative, which means that when alchemical transformation takes place from state A to state B, the total charge present on the amino acid changes. As we use the hybrid topology approach for the alchemical transformation, during the neutralisation stage only one state (state A or B) is sufficiently neutralised, leaving the system with a net charge of the other state (82) . This gives rise to a technical error during the simulation when using the state-of-the-art PME method for handling the Coulombic interactions. Hence, for a non-neutral system, any extra charge is neutralized by the implicit introduction of a uniform background charge (101, 102) . Due to the contribution of this background charge in the free energy difference, an additional correction scheme is needed to correct the finite charge effect (61, 103) . Here, we used an analytical post-simulation correction A Guide to COVID-19: a global pandemic caused by the novel coronavirus SARS-CoV-2. The FEBS journal Cellular and Molecular Life Sciences Pathogenesis of Virus Infections Neutral Theory and Rapidly Evolving Viral Pathogens. Molecular biology and evolution SARS-CoV-2 variants, spike mutations and immune escape A dynamic nomenclature proposal for SARS-CoV-2 lineages to assist genomic epidemiology Tracking SARS-CoV-2 variants Estimated transmissibility and impact of SARS-CoV-2 lineage B.1.1.7 in England SARS-CoV-2 B.1.1.7 sensitivity to mRNA vaccine-elicited, convalescent and monoclonal antibodies. medRxiv : the preprint server for health sciences Estimated transmissibility and impact of SARS-CoV-2 lineage B.1.1.7 in England Alchemical binding free energy calculations in AMBER20: Advances and best practices for drug discovery Accurate and Rigorous Prediction of the Changes in Protein Free Energies in a Large-Scale Mutation Scan Assignment of epidemiological lineages in an emerging pandemic using the pangolin tool. Virus Evolution Comprehensive mapping of mutations in the SARS-CoV-2 receptor-binding domain that affect recognition by polyclonal human plasma antibodies. Cell Host and Microbe Complete Mapping of Mutations to the SARS-CoV-2 Spike Receptor-Binding Domain that Escape Antibody Recognition. Cell Host and Microbe Sensitivity of infectious SARS-CoV-2 B.1.1.7 and B.1.351 variants to neutralizing antibodies Deep Mutational Scanning of SARS-CoV-2 Receptor Binding Domain Reveals Constraints on Folding and ACE2 Binding Evidence of escape of SARS-CoV-2 variant B.1.351 from natural and vaccine-induced sera Transmission, infectivity, and neutralization of a spike L452R SARS-CoV-2 variant Fact Sheet for Health Care Providers Emergency Use Authorization of Bamlanivimab and Etesevimab SARS-CoV-2 Variant Classifications and Definitions Genomic Variations in the Structural Proteins of SARS-CoV-2 and Their Deleterious Impact on Pathogenesis: A Comparative Genomics Approach. Frontiers in Cellular and Infection Microbiology SARS-CoV-2 spike E484K mutation reduces antibody neutralisation. The Lancet Microbe Multiple SARS-CoV-2 variants escape neutralization by vaccine-induced humoral immunity How Dangerous Is the Delta Variant (B.1.617.2)? SARS-CoV-2 and ORF3a: Nonsynonymous Mutations, Functional Domains, and Viral Pathogenesis. mSystems Pathogenicity and Transmissibility of Delta and Lambda Variants of SARS-CoV-2, Toxicity of Spike Protein and Possibilities for Future Prevention of COVID-19 Structural and biochemical rationale for enhanced spike protein fitness in delta and kappa SARS-CoV-2 variants Neutralization of the SARS-CoV-2 Mu Variant by Convalescent and Vaccine Serum R346K Mutation in the Mu Variant of SARS-CoV-2 Alters the Interactions with Monoclonal Antibodies from Class 2: A Free Energy Perturbation Study Structural and functional insights into the major mutations of SARS-CoV-2 Spike RBD and its interaction with human ACE2 receptor Transmission, infectivity, and neutralization of a spike L452R SARS-CoV-2 variant Epitope Classification and RBD Binding Properties of Neutralizing Antibodies Against SARS-CoV-2 Variants of Concern. Frontiers in Immunology SARS-CoV-2 B Mutations L452R and E484Q Are Not Synergistic for Antibody Evasion. The Journal of Infectious Diseases Systemic effects of missense mutations on SARS-CoV-2 spike glycoprotein stability and receptor-binding affinity Epitope Classification and RBD Binding Properties of Neutralizing Antibodies Against SARS-CoV-2 Variants of Concern. Frontiers in Immunology SARS-CoV-2 spike L452R variant evades cellular immunity and increases infectivity Membrane fusion and immune evasion by the spike protein of SARS-CoV-2 Delta variant Crucial Mutations of Spike Protein on SARS-CoV-2 Evolved to Variant Strains Escaping Neutralization of Convalescent Plasmas and RBD-Specific Monoclonal Antibodies Experimental Evidence for Enhanced Receptor Binding by Rapidly Spreading SARS-CoV-2 Variants Predicting Mutational Effects on Receptor Binding of the Spike Protein of SARS-CoV-2 Variants Receptor binding, immune escape, and protein stability direct the natural selection of SARS-CoV-2 variants Anton Van Der Merwe P. Effects of common mutations in the sars-cov-2 spike rbd and its ligand the human ace2 receptor on binding affinity and kinetics. eLife Enhanced binding of the N501Y-mutated SARS-CoV-2 spike protein to the human ACE2 receptor: insights from molecular dynamics simulations SARS-CoV-2 spike protein N501Y mutation causes differential species transmissibility and antibody sensitivity: a molecular dynamics and alchemical free energy study. Molecular Systems Design & Engineering Reduced neutralization of SARS-CoV-2 B.1.1.7 variant by convalescent and vaccine sera Fast prediction of binding affinities of the sars-cov-2 spike protein mutant n501y (UK variant) with ace2 and miniprotein drug candidates Structural and Functional Basis of SARS-CoV-2 Entry by Using Human ACE2 Structure of the SARS-CoV-2 spike receptorbinding domain bound to the ACE2 receptor A hierarchical approach to all-atom protein loop prediction On the Role of the Crystal Environment in Determining Protein Side-chain Conformations CUPSAT: prediction of protein stability upon point mutations Structural analysis and prediction of protein mutant stability using distance and torsion potentials: role of secondary structure and solvent accessibility. Proteins Computational modeling of protein mutant stability: analysis and optimization of statistical potentials and structural features reveal insights into prediction model development The FoldX web server: an online force field Local water bridges and protein conformational stability. Protein Science Electrostatic enhancement of diffusion-controlled protein-protein association: comparison of theory and experiment on barnase and barstar Prediction of water and metal binding sites and their affinities by using the Fold-X force field Applying Physics-Based Scoring to Calculate Free Energies of Binding for Single Amino Acid Mutations in Protein-Protein Complexes Calculating the binding free energies of charged species based on explicit-solvent simulations employing lattice-sum methods: An accurate correction scheme for electrostatic finite-size effects Estimated transmissibility and impact of SARS-CoV-2 lineage B.1.1.7 in England Tracking Changes in SARS-CoV-2 Spike: Evidence that D614G Increases Infectivity of the COVID-19 Virus Evaluating the Effects of SARS-CoV-2 Spike Mutation D614G on Transmissibility and Pathogenicity Mapping mutations to the SARS-CoV-2 RBD that escape binding by different classes of antibodies Identification of SARS-CoV-2 spike mutations that attenuate monoclonal and serum antibody neutralization In vitro data suggest that Indian delta variant B.1.617 of SARS-CoV-2 escapes neutralization by both receptor affinity and immune evasion SARS-CoV-2 neutralizing antibody structures inform therapeutic strategies Escape from neutralizing antibodies by SARS-CoV-2 spike protein variants. eLife The Protein Data Bank OPLS4: Improving Force Field Accuracy on Challenging Regimes of Chemical Space In silico mutagenesis of human ACE2 with S protein and translational efficiency explain SARS-CoV-2 infectivity in different species Scalable algorithms for molecular dynamics simulations on commodity clusters ACM/IEEE Conference on Supercomputing, SC'06 Interaction Models for Water in Relation to Protein Hydration A unified formulation of the constant temperature molecular dynamics methods A molecular dynamics method for simulations in the canonical ensemble Canonical dynamics: Equilibrium phase-space distributions Constant pressure molecular dynamics algorithms The u -series: A separable decomposition for electrostatics computation with improved accuracy Best Practices for Alchemical Free Energy Calculations Insights from the First Principles Based Large Scale Protein Thermostability Calculations Accurate Calculation of Free Energy Changes upon Amino Acid Mutation Protein thermostability calculations using alchemical free energy simulations Efficient estimation of free energy differences from Monte Carlo data A User Friendly Interface for Alchemistry Automated protein structure and topology generation for alchemical perturbations Comparison of multiple amber force fields and development of improved protein backbone parameters 0: A package for molecular simulation and trajectory analysis GROMACS: Fast, flexible, and free Algorithms for highly efficient, load-balanced, and scalable molecular simulation GROMACS 4.5: A highthroughput and highly parallel open source molecular simulation toolkit Tackling exascale software challenges in molecular dynamics simulations with GROMACS High performance molecular simulations through multi-level parallelism from laptops to supercomputers Determination of alkali and halide monovalent ion parameters for use in explicitly solvated biomolecular simulations The Origin of the Method of Steepest Descent Molecular dynamics with coupling to an external bath A Leap-Frog Algorithm for Stochastic Dynamics Polymorphic transitions in single crystals: A new molecular dynamics method Particle mesh Ewald: An N·log(N) method for Ewald sums in large systems A smooth particle mesh Ewald method An overview of electrostatic free energy computations for solutions and proteins Accurate Calculation of Relative Binding Free Energies between Ligands with Different Net Charges The authors declare no competing interests.