key: cord-0278199-0uff0qg2
authors: Aggarwal, Abhishek; Naskar, Supriyo; Maroli, Nikhil; Gorai, Biswajit; Dixit, Narendra M.; Maiti, Prabal K.
title: Mechanistic Insights into the Effects of Key Mutations on SARS-CoV-2 RBD-ACE2 Binding
date: 2021-09-21
journal: bioRxiv
DOI: 10.1101/2021.09.20.461041
sha: 203f2c9e3581af1d47ddf4c4beba3bdb21b4503c
doc_id: 278199
cord_uid: 0uff0qg2

Some recent SARS-CoV-2 variants appear to have increased transmissibility than the original strain. An underlying mechanism could be the improved ability of the variants to bind receptors on target cells and infect them. In this study, we provide atomic-level insight into the binding of the receptor binding domain (RBD) of the wild-type SARS-CoV-2 spike protein and its single (N501Y), double (E484Q, L452R) and triple (N501Y, E484Q, L452R) mutated variants to the human ACE2 receptor. Using extensive all-atom molecular dynamics simulations and advanced free energy calculations, we estimate the associated binding affinities and binding hotspots. We observe significant secondary structural changes in the RBD of the mutants, which lead to different binding affinities. We find higher binding affinities of the double (E484Q, L452R) and triple (N501Y, E484Q, L452R) mutated variants than the wild type and the N501Y variant, which could contribute to the higher transmissibility of recent variants containing these mutations.

The COVID-19 pandemic has caused more than 4.5 million deaths so far. The severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2), the causative agent of COVID-19, attaches to target host cells with its spike protein, called S, by binding to the angiotensinconverting enzyme 2 (ACE2) receptor, mainly expressed in the lungs [1] [2] [3] [4] , and causes the infection. The binding of S to ACE2 triggers conformational changes in S from its metastable prefusion state to a stable post-fusion state 5 . The spike protein is divided into two subunits, S1 and S2, in the prefusion stage. The receptor binding domain (RBD) in S1 is responsible for the binding with ACE2, whereas S2 mediates the subsequent fusion of viral and cell membranes, allowing viral entry 5 . Characterizing the molecular details of the interaction between the RBD and ACE2 plays an important role in understanding the process of SARS-CoV-2 infection of target cells. This understanding may also help explain the improved transmissibility that is seen with SARS-CoV-2 variants that carry mutations potentially affecting RBD-ACE2 binding.

SARS-CoV-2 belongs to the betacoronavirus genus, which includes SARS-CoV and the middle-east respiratory syndrome virus (MERS) 6 . SARS-CoV-2 and SARS-CoV have high sequence similarity, of ~80%. SARS-CoV-2, however, has spread far more than SARS-CoV and that may be due to the higher binding energy of SARS-CoV-2 with ACE2 7, 8 . More recently, several variants of SARS-CoV-2 carrying mutations in RBD, and other regions of S have been reported that have heightened transmissibility and can cause infections with different severity among individuals 9, 10 . The double mutated variant B.1.617 has already caused a second wave of COVID cases in India 11 . There have been several reports from the Indian sub-continent about the potential infection and death rate or severe complications among the B.1.618 triple mutant variants 12 . The underlying molecular mechanism behind the binding of these variants with human ACE2 receptor is still unknown.

In this study, we estimate the binding energy between the ACE2-RBD of the SARS-CoV-2 wild-type, and of its single (N501Y), double (E484Q, L452R) and triple (E484Q, L452R, and N501Y) mutated variants. We have considered the two mutations E484Q and L452R seen in the variant B.1.617 responsible for the second wave in India 13 . We considered these mutations along with B.1.1.7 N501Y variant, one of the early mutations argued to have caused increased transmissibility 14 . We calculate the binding energy differences due to these mutations and elucidated the underlying molecular mechanisms using all-atom molecular dynamics simulations. We examine structural rearrangements of the RBD in the mutated variants leading to these differences. We also employ the ASGB method for a quantitative analysis of the binding mechanisms.

The three-dimensional structures of ACE2 and RBD domains of SARS-CoV-2 were obtained from the protein data bank PDB ID: 6M17 15 . We studied various mutated variants of this structure as listed in Table 1 . The mutated structures were prepared using the CHIMERA software package. Using the xLEaP module of AmberTools18 16 , the protein structures were enclosed in a water box with a buffer of 1.5 nm in all three directions. The systems were then charge neutralized by adding an appropriate number of Na + ions. Additional Na + ions and an equal number of Clions were further added to achieve physiological salt concentration of 150 mM. TIP3P water model 17 was used for solvating the proteins. ff99sb-ildn 18 force field was used to represent the inter-and intra-molecular interactions of the protein atoms. Interactions involving ions were described using the Joung Cheatham parameter set 19 .

The solvated systems were initially energy minimized for 5000 steps using the steepest decent algorithm. The minimized structures were then heated from 0 to 300 K within a period of 50 ps with the solute atoms held fixed by a harmonic potential of strength 20 kcal/mol/Å 2 . After heating, the systems were subjected to a 2 ns long unrestrained equilibration. Finally, the equilibrated structures were subjected to 200 ns long production runs in the NVT ensemble.

All the simulations were performed using the AMBER18 simulation package 16 . Periodic boundary condition was employed in all the three directions. SHAKE constraints were used for all bonds containing hydrogen allowing the use of a time step of 2 fs. Langevin thermostat with a collision frequency of 2 ps -1 was used for temperature regulation while isotropic Berendsen barostat with a coupling constant of 0.5 ps was used to regulate pressure in the NPT simulations. A cut-off of 10 Å was used to compute the short-range Lennard-Jones (LJ) interaction, while PME (Particle Mesh Ewald 20 ) was employed for the calculation of longrange electrostatic interaction. The snapshots were visualized using VMD 21 for analysis.

The MMPBSA method 22,23 employed in the MMPBSA.py module of AMBER18 was used to calculate the binding energies of the different ACE2-RBD complexes. The binding free energy difference ( ) of the SARS-CoV-2 RBD and ACE2 complex formation is calculated as 24

where , , and 2 represent the free energies of the RBD-ACE2 complex, individual RBD and individual ACE2 receptor, respectively. Eq. (1) can be decomposed into different interactions and can be written as

Here, is the change in gas-phase molecular mechanics energy, is the entropy change (which has been neglected in this work), and represents the solvation free energy change upon ligand binding. can be further written as

Here , , and represent the change in bonded energies (bond, angle, and dihedral), electrostatic energies, and van der Waals energies upon ligand binding, respectively.

is the sum of the nonpolar solvation energy ( ) and electrostatic solvation energy ( ), computed using the Poisson-Boltzmann method. can be written as,

Here, (=0.00542 kcal/Å 2 ) is the surface tension, while =0.92 kcal/mol, and SASA represents the solvent-accessible surface area of the molecule.

The Alanine Scanning with MM/GBSA (ASGB) method is used on the last 10 ns of the 200 ns long MD simulation of each of the four RBD-ACE2 complex structures. In this method, all the residues lying within the 5 Å range of the protein-protein interaction interface are mutated. We use ASGB method with dielectric constants of 1, 3, and 5 for nonpolar, polar, and charged residues 24-26 , respectively. The binding free energy difference due to a point mutation was calculated by using the following relations,

Here, → ( ) is the change in the binding free energy difference upon ligand binding due to a single point-mutation in RBD, while → is the change in the solvation free energy difference upon ligand binding due to single point-mutation in RBD. represent the difference in electrostatic and van der Waals energy between the mutated residue ( ) in RBD with ACE2 to that of pre-mutated residue ( ) in RBD with ACE2.

The structural stability of wild type, the single (N501Y), double (E484Q, L452R), and triple (N501Y, E484Q, L452R) mutants has been evaluated by analysing the RMSD fluctuations of backbone atoms and Cα atoms (Fig. 1) . Fig. 4 We have calculated the number of hydrogen bonds between RBD-ACE2 for wild type and Using MM-PBSA method, we compute the total binding energy of RBD with ACE2 (Fig. 6) . To gain further insight into the mechanism causing these changes in the binding affinities, we have used the alanine scanning method (ASGB) to predict the hotspot residues in RBD involved in ACE2-RBD binding. The residues Asn709, Tyr710, Ile733, Val744, Glu745, Phe751 and Asn762 show more than 2 kcal/mol energy differences. The reason for high binding energy from these residues is their positioning at the ACE2 binding site. Among these residues, we found that Glu745 contributed more due to its hydrogen bond formation with the Phe10 of ACE2. Further, the contribution from residues 669R, 678K, 710Y, 747F, 754Q, 759Q, 761T and 766Y are higher for the wild type ACE2-RBD (Fig. 7) . Large conformational changes are caused due to the mutations in RBD that increase the number of residues that interact with ACE2 receptor. We find that 747F is a key contributor to the binding energy of all the mutated variants of RBD we considered. The total contribution from the hotspot residues is -20.84 kcal/mol and single and double mutant were identified as -29.94 kcal/mol and -36.1 kcal/mol. The triple (N501Y, E484Q, L452R) mutant variant shows a higher binding affinity than wild type (-33.06 kcal/mol). As shown in Table 2 , the contribution from each residue is found to be higher for the mutated variants than the wild type. Moreover, the double mutated (E484Q, L452R) variant shows the highest binding affinity. We note that ASGB does not consider the repulsive interactions which may lead to a difference between the computed binding energies using MMPBSA and ASGB. Nevertheless, the qualitative results obtained using ASGB agree well with MMPBSA and indicate that the double mutant (E484Q, L452R) RBD spike has the highest binding affinity for ACE2 among the variants considered.

Next, we have compared our estimated binding energy data with available experimental values. 

There are no conflicts of interest to declare. Hotspot residues of ACE2 binding to wild type, N501Y single mutated, E484Q, L452R double mutated, and N501Y, E484Q, L452R triple mutated RBD obtained using the ASGB calculation.

Full-Genome Evolutionary Analysis of the Novel Corona Virus (2019-NCoV) Rejects the Hypothesis of Emergence as a Result of a Recent Recombination Event. Infection

Characterization of the Receptor-Binding Domain (RBD) of 2019 Novel Coronavirus: Implication for Development of RBD Protein as a Viral Attachment Inhibitor and Vaccine

Clinical Features of Patients Infected with 2019 Novel Coronavirus in Wuhan

A Novel Coronavirus from Patients with Pneumonia in China

The Potential Role of Procyanidin as a Therapeutic Agent against SARS-CoV-2: A Text Mining, Molecular Docking and Molecular Dynamics Simulation Approach

Silico Approaches to Detect Inhibitors of the Human Severe Acute Respiratory Syndrome Coronavirus Envelope Protein Ion Channel

Structure of the SARS-CoV-2 Spike Receptor-Binding Domain Bound to the ACE2 Receptor

Genomic Characterisation and Epidemiology of 2019 Novel Coronavirus: Implications for Virus Origins and Receptor Binding

Early Empirical Assessment of the N501Y Mutant Strains of SARS-CoV-2 in the United Kingdom

Molecular Mechanism of the N501Y Mutation for Enhanced Binding between SARS-CoV-R s Spike Protein and Human ACE2 Receptor

Convergent Evolution of SARS-CoV-2 Spike Mutations, L452R, E484Q and P681R, in the Second Wave of COVID-19 in Maharashtra

Triple Mutant Bengal Strain (B.1.618) of Coronavirus and the Worst COVID Outbreak in India

SARS-CoV-2 Variant B.1.617 Is Resistant to Bamlanivimab and Evades Antibodies Induced by Infection and Vaccination

The N501Y Spike Substitution Enhances SARS-CoV-2 Transmission

Structural Basis for the Recognition of SARS-CoV-2 by Full-Length Human ACE2

The Amber Biomolecular Simulation Programs

Comparison of Simple Potential Functions for Simulating Liquid Water

Improved Side-Chain Torsion Potentials for the Amber Ff99SB Protein Force Field. Proteins: Structure, Function, and Bioinformatics

Determination of Alkali and Halide Monovalent Ion Parameters for Use in Explicitly Solvated Biomolecular Simulations

Particle Mesh Ewald: An N·log(N) Method for Ewald Sums in Large Systems

VMD: Visual Molecular Dynamics

MMGBSA as a Tool to Understand the Binding Affinities of Filamin-Peptide Interactions

The MM/PBSA and MM/GBSA Methods to Estimate Ligand-Binding Affinities

Calculation of Hot Spots for Protein-Protein Interaction in P53/PMI-MDM2/MDMX Complexes

Computational Alanine Scanning with Interaction Entropy for Protein-Ligand Binding Free Energies

Quantitative Analysis of ACE2 Binding to Coronavirus Spike Proteins: SARS-CoV-2 vs. SARS-CoV and RaTG13

Experimental Evidence for Enhanced Receptor Binding by Rapidly Spreading SARS-CoV-2 Variants

Differential Interactions Between Human ACE2 and Spike RBD of SARS-CoV-2 Variants of Concern

We thank DST, India for computational support through TUE-CMS computational facility.