key: cord-1034562-3ld9ux5y
authors: Mandal, Nabanita; Padhi, Aditya K.; Rath, Soumya Lipsa
title: Molecular Insights into the Differential Dynamics of SARS-CoV-2 Variants of Concern (VOC)
date: 2021-10-22
journal: bioRxiv
DOI: 10.1101/2021.10.22.465272
sha: bc7aaf86a89a85fb61f09ee247005c46a6c90c7c
doc_id: 1034562
cord_uid: 3ld9ux5y

Severe Acute Respiratory Syndrome Coronavirus 2 (SARS-CoV-2) has affected the lives and livelihood of millions of individuals around the world. It has mutated several times after its first inception, with an estimated two mutations occurring every month. Although we have been successful in developing vaccines against the virus, emergence of variants has enabled it to escape therapy. Few of the generated variants are also reported to be more infectious than the Wild-type. The World Health Organization (WHO) has prioritized the variants into Variants of Concern (VOC) and Variants of Interest (VOI) to focus on the variants that pose a threat to public health. In this study, we compare the structural and dynamic attributes of all the RBD/ACE2 complexes for the reported VOCs, namely, Alpha, Beta, Gamma, and Delta through atomistic simulations. Results indicated that the orientation and binding energies of the VOCs were different from the Wild-type. Specifically, we observed that the Receptor Binding Domain (RBD) in B.1.351 (Beta) and B.1.617.2 (Delta) underwent a relative rotation to make the complex more compact. Protein dynamics, however, show that their fluctuations were similar to the Wild-type. It was also reflected in the calculation of the total interaction energies. Overall, it was observed that electrostatic interactions play a major role in the formation of the complexes. Detailed residue level energetics revealed that the most prominent changes in interaction energies were seen particularly at the mutated residues which were present at RBD/ACE2 interface. We found that B.1.167.2 (Delta) is one of the most tightly bound variants of SARS-CoV-2 with dynamics similar to Wild-type. High binding affinity of RBD towards ACE2 is indicative of an increase in the viral infectivity. Therefore, the intrinsic details presented in this study would be useful for the design and development of effective therapeutic strategies for the emerging variants of the virus.

Severe Acute Respiratory Syndrome Coronavirus 2 (SARS-CoV-2), is one of the largest known RNA viruses with a single-stranded RNA ranging between 26,000 to 32,000 bases [1] . Most of these RNA viruses are prone to mutations [2, 3] . Moreover, RNA viruses have a higher mutation rate compared to the DNA viruses, thereby, reflecting a higher replication fidelity of the DNAdependent DNA polymerases over that of the RNA-dependent RNA polymerases [4] .

Additionally, positive-sense single stranded RNA viruses, such as SARS-CoV-2, have a much higher mutation rate than the negative-sense single stranded RNA viruses [5] . The estimated rate of mutation reported is two per month for SARS-CoV-2 [6] . Mutation is also one of the primary generators of diversity among the genomes, including the viral genomes [6, 7] . Thus, we are observing an emergence of variants of SARS-CoV-2 since the first incidence of the Coronavirus disease 2019 (COVID-19) [8] .

Although SARS-CoV-2 originated in 2019, it has undergone roughly about 82062 mutations according to the GISAID database [9] . Particularly, mutation in the gene which encodes the Spike glycoprotein has been reported to be very high [10] . It is also well documented that those mutations which are found to affect the glycosylation of viral proteins affect the viral life-cycle [11] .

Mutations at those residues which are used for recognition by antibodies make the virus resistant to antibody mediated neutralization [12] . In Spike protein, mutations can affect the transmission rate of the virus and also the disease outcome [13] . The first reported mutation of the Spike protein was D614G mutation, which was found to enhance the SARS-CoV-2 transmission [13] . The N501Y mutation increased the affinity of the Spike protein for its receptor, Angiotensin 5 Converting Enzyme 2 (ACE2), thereby increasing the chances of viral transmission [14] . Mutation E484K is known to contribute to the evasion of antibody neutralization [15] . D796H and H655Y mutations that are present in the Spike protein, are associated with reduced affinity towards the neutralizing antibodies [16] . The World Health Organization (WHO) has recently assigned different labels for the generated variants of SARS-CoV-2. They can be broadly separated into two categories, namely, the Variants of Concern (VOC) and the Variants of Interest (VOI) [17] .

VOC has increased transmissibility and severity of (COVID- 19) The SARS-CoV-2 has several structural and non-structural proteins which assist the virus from initial attachment, entry to replication and other vital functions [18, 19] . One of the largest and prominent structural proteins of SARS-CoV-2 is the Spike protein. This protein lies on the external viral membrane and helps in initial attachment of the virus with the host receptors [19] .

Surprisingly, all the mutations that were observed in the VOC were primarily located in the RBD 6 ranging from residue 333-527 of the Spike protein [20] . The RBD of the SARS-CoV-2 Spike protein interacts with a human ACE2 receptor. This receptor is found on the lung alveolar epithelial cells and plays a primary role in protection against lung injury in humans [21] . Several studies have shown that the difference in the sequence of RBD between SARS-CoV-1 and SARS-CoV-2 have increased the binding affinity of the SARS-CoV-2 RBD towards ACE2 [22] . It is therefore essential to learn how these mutations impact the association of the RBD with the ACE2.

Very recently, the cryo-EM structure of the Delta Variant has been solved [23] , however, the differences in structures and binding are yet to be unraveled. In the present study, we used molecular modelling tools to model the RBD domains of all the reported VOCs. Subsequently, we compared the generated models of the VOCs in the RBD/ACE2 complex with the Wild-type by using extensive molecular dynamics (MD) simulations and identified several key features which result in differential activity of the VOCs. This study provides mechanistic and molecular insights of the VOCs and would prove crucial for understanding the structure-function relationship as well as in the development of effective therapeutic strategies for the emerging variants of SARS-CoV-2.

The crystal structure of the Spike protein RBD associated with ACE2 was taken from Protein Data Bank (PDB ID: 6LZG) as the starting structure [24] . The RBD has residues ranging from Tyr333 to Pro527. Although 13 residues were missing from the N-terminal of ACE2, the N-terminal 7 residues do not directly interact with the RBD [24, 25] , hence they were not modelled. This structure was considered as the Wild-type system. The variants namely, P.1(Alpha), B.1.1.7(Gamma), B.1.351(Beta) and B.1.167.2 (Delta) were generated by mutating specific residues in Wild-type Spike protein after aligning the RBD sequences shown in Figure S1 .

Modeller 10.1 molecular modelling suite was used to generate the new models based on the Wildtype RBD/ACE2 template [26] . To understand the effect of mutations on the structure and dynamics of the Spike-ACE2 complex, we performed all-atom MD simulations of Wild-type (WT) and four variants of Spike-RBD/ACE2. We also checked the initial structures of the generated variants by observing the distribution of their phi-psi angles and other stereochemical properties by using the PROCHECK [27] suite of programs.

Atomistic MD Simulation was carried out using Gromacs MD Simulation package [28] using CHARMM36 force field parameters [29] . Each system was subjected to energy minimization in steepest descent and then in conjugate gradient for 2000 steps. After initial relaxation, a cubic simulation box consisting of three-site TIP3P water molecules and neutralizing ions was created for the systems [30] . The box dimensions were 10 x 10 x 10 Å. For charge neutralization 24 ions were randomly placed by replacing the corresponding solvent molecules. Subsequently, energy minimization and thermalization was performed to avoid any bad contacts which might have been created due to the mutations and addition of water and ions. Periodic boundary condition was implemented during simulation. The systems were gradually heated from 0 to 310 K for 200 ps.

Then the systems were equilibrated at 310 K in NVT ensemble using modified Berendsen 8 thermostat [31] for about 500ps and then equilibrated in NPT ensemble using 1 atmospheric pressure using Parrinello-Rahman barostat [31] for 1 ns. A time step of 2 fs was used for all the equilibration and subsequent production runs. After the convergence of potential energy and density, production simulation was carried out for the Wild-type and VOCs for 100 ns in NPT where the coordinates were saved at the interval of every 1000 ps. Particle-mesh Ewald method was used to treat the long-range electrostatic interactions [32] . VMD and Pymol were used for visualization of the trajectories [33, 34] . All the analyses were carried out using Gromacs tools [35] .

The binding energy between RBD and ACE2 for WT and VOCs [36] . In this methodology, the binding energy of the target-ligand or protein-protein is typically defined as

where ΔGprotein, ΔGcomplex, and ΔGligand represent the total free energies of the complex, the ligand, and the protein also separately in the solvent, respectively. Further, the free energy of the separate entity is represented as

where EMM stands for the average molecular mechanic's potential energy in the vacuum, Gsolvation denotes the free energy of solvation. TS stands for the entropic augmentation to the free energy in a vacuum, here S and T denote the entropy and temperature, respectively. The EMM consists of nonbonded and bonded terms, including torsion, bond angle and electrostatic (Eelec) and the Van der waal (Evdw) interactions, respectively & the solvation free energy, Gsolvation takes both electrostatic and non-electrostatic (Gpolar and Gnonpolar) components. The binding free energy for the complexes was calculated from 50 snapshots over the last 10 ns of the simulation trajectories.

All the systems were stable during this period.

To understand the intermolecular interactions formed between RBD and ACE2 for WT and all the VOCs, contacts (hydrogen bonds and salt bridges) were computed and analyzed from the last 10 ns MD simulated trajectories using GetContacts [37] . The hydrogen bonds were shown in a clustergram to make the interpretation clear for visualization.

The initial coordinates of the RBD/ACE2 complex were taken from the crystal structure (PDB ID:

6LZG) [24] . As shown in Figure 1 and Figure (data not shown). Very recently, the structure of the B.1.167.2 has been deposited in the protein data bank (PDB ID: 7V8B) [23] . We have made a comparison between the structures generated after simulation with that of the Cryo-EM structure in the later parts of the article. The

Ramachandran plot of the generated models did not show drastic changes as expected from a static model ( Figure S3 ). Notably, these mutations were found to occur close to the RBD/ACE2

interface. It is therefore fascinating to investigate how Spike-ACE2 interaction and dynamics might be affected due to mutations.

To begin with, we ran atomistic MD simulations for WT, B. however its presence is also seen in other organs [38] . Further reviews on its role in various metabolic pathways are described elsewhere [39] . Structurally, ACE2 is an alpha helical protein [40] , where the crystal structure and our models that used the crystal structure as template, contain residues ranging from 19-615 of the ACE2 protein. The N-terminal helical part of the protein is the primary site of interaction with the RBD of the Spike protein of SARS-CoV-2. The overall binding mode of ACE2 with SARS-CoV-1 and SARS-CoV-2 is known to be largely similar [40, 41] . The WT RBD shows that the N-terminal domain along with residues Glu329, Asn330 and For analysis, the PCA was performed on the backbone Cɑ atoms of RBD and ACE2 separately [43] . Figure S4 shows the comparison of motion of the proteins along the first principal components where we could observe potential protein dynamics. From Figure S4 it is evident that landscape was then constructed as a function of the PC1 and PC2 coordinates. The highly stable protein conformation is shown in red, other low energy states are colored either in blue, green, or cyan. In the WT complex, the ACE2 was confined to a single cluster whereas, RBD explored two separate clusters (Figure 4a ). This indicates that RBD can exist in two different conformations after being bound to ACE2. Although not much changes could be seen in the binding of RBD, the difference in clusters can be mainly attributed to the loop dynamics (Figure 4f ). When we compared the variants, we found that a greater number of conformational states of the ACE2 in The trace values are correlated with the total variance in the values of eigenvectors where higher 14 values indicate more variation [44] . For ACE2 in WT, P. To explore the rationale behind the differential dynamics, we used the MM/PBSA to calculate the total binding energies of the protein complexes (Table 1) . We used the last 10 ns of simulated trajectories of all the systems for MM/PBSA based binding energy calculations, where the systems were mostly stable. Table 1 

Interfacial residues of proteins play a significant role in the association as well as in governing stability of the protein complexes. In the earlier studies by Spinello et al., [45] they compared the interface of SARS-CoV-1 with SARS-CoV-2 and found substantial differences in the interaction energies between residues of ACE2 and RBD. Here, we further calculated residue-wise contribution towards the total binding energy of the complex for both ACE2 and RBD by using MM/PBSA. Several of the residues were found to show drastic differences in their binding energies. We compared the energy of those residues that contributed >10 kJ/mol towards the binding energy of the complex. The interaction energies of RBD show radical changes in values for the P.1 complex. Here, almost all of the residues, except Glu484, contribute negligibly towards protein binding (Figure 5a ). We further noticed that E484 mutates into K484 in both P. Subsequently, we checked if complementary changes take place in the ACE2 receptor of the complexes. Upon comparing the interaction energy of the ACE2 protein residues with RBD, calculated from the last 10ns of the trajectory (Figure 5b ), we found that in the WT complex except The orientation of the RBD with respect to ACE2 was checked to observe any conformational changes that might take place due to mutation. It was seen that at the interface region Gamma moved away and both Beta and Delta moved closed towards the ACE2 receptor. This also shows that the interacting residues might be more closely associated in Beta and Delta systems. The RBD region is shown as cyan bar in (b), where the RBM is highlighted in orange. 

Severe Acute Respiratory Syndrome Coronavirus-2 (SARS-CoV-2): An Update

Are RNA Viruses Candidate Agents for the Next Global Pandemic? A Review

Why are RNA virus mutation rates so damn high?

Emerging SARS-CoV-2 mutation hot spots include a novel RNA-dependent-RNA polymerase variant

Coronavirus envelope protein: current knowledge

SARS-CoV-2 variants, Spike mutations and immune escape

COVID-19: emergence and mutational diversification of SARS-CoV-2

Characteristics of SARS-CoV-2 and COVID-19

Tracking Changes in SARS-CoV-2 Spike: Evidence that D614G Increases Infectivity of the COVID-19 Virus

Ojha Amrita, Peng Haiyong. SARS-CoV

Spike-protein D614G mutation increases virion Spike density and infectivity

The Impact of Mutations in SARS-CoV-2 Spike on Viral Infectivity and Antigenicity

Antibody resistance of SARS-CoV-2 variants B.1.351 and B.1.1.7. Nature

D614G mutation in the SARS-CoV-2 spike protein enhances viral fitness by desensitizing it to temperature-dependent denaturation

Signatures in SARS-CoV-2 Spike protein conferring escape to neutralizing antibodies

SARS-CoV-2 Spike E484K mutation reduces antibody neutralization. The Lancet Microbe 2021

Structural basis of fitness of emerging SARS-COV-2 variants and considerations for screening, testing and surveillance strategy to contain their threat. medRxiv preprint 2021

The origins and potential future of SARS-CoV-2 variants of concern in the evolving COVID-19 pandemic

Neutralizing monoclonal antibodies for treatment of COVID-19

Structural and functional properties of SARS-CoV-2 Spikeprotein: potential antivirus drug development for COVID

Large-scale analysis of SARS-CoV-2 spike-glycoprotein mutants demonstrates the need for continuous screening of virus isolates

Uhal Bruce D. ACE2, Much More Than Just a Receptor for SARS-COV-2

Structure of the SARS-CoV-2 Spike receptor-binding domain bound to the ACE2 receptor

Local refinement of SARS-CoV-2 S-Delta variant (B.1.617.2) RBD and Angiotensin-converting enzyme 2 (ACE2) ectodomain

Structural and Functional Basis of SARS-CoV-2 Entry by Using Human ACE2

Structural basis for the recognition of SARS-CoV-2 by full-length human ACE2

Comparative Protein Structure Modeling Using MODELLER

Stereochemistry and Validation of Macromolecular Structures

High performance molecular simulations through multi-level parallelism from laptops to supercomputers

CHARMM: The Biomolecular Simulation Program

All-Atom Empirical Potential for Molecular Modeling and Dynamics Studies of Proteins

All-Atom Empirical Potential for Molecular Modeling and Dynamics Studies of Proteins

Particle mesh Ewald: An N·log(N) method for Ewald sums in large systems

The PyMOL Molecular Graphics System

VMD: visual molecular dynamics

All-Atom Empirical Potential for Molecular Modeling and Dynamics Studies of Proteins

Open Source Drug Discovery Consortium, Lynn A. g_mmpbsa--a GROMACS tool for high-throughput MM-PBSA calculations

Uncovering patterns of atomic interactions in static and dynamic structures of proteins

Body Localization of ACE-2: On the Trail of the Keyhole of SARS-CoV-2

Stochastic fluctuations in metabolic pathways

Role of angiotensin-converting enzyme 2 (ACE2) in COVID-19

The SARS-CoV-2 Spike Glycoprotein Biosynthesis, Structure, Function, and Antigenicity: Implications for the Design of Spike-Based Vaccine Immunogens

Emergence and spread of SARS-CoV-2 lineage B.1.620 with variant of concern-like mutations and deletions

Principal component analysis: a review and recent developments

The application of principal component analysis to drug discovery and biomedical data

Is the Rigidity of SARS-CoV-2 Spike Receptor-Binding Motif the Hallmark for Its Enhanced Infectivity? Insights from All-Atom Simulations

Structure-Function Analysis of Resistance to Bamlanivimab by SARS-CoV-2 Variants Kappa, Delta, and Lambda

Ten emerging SARS-CoV-2 Spike variants exhibit variable infectivity, animal tropism, and antibody neutralization

We are grateful to the COVID-19 HPC Consortium for providing resources and helping researchers work for a noble cause.

NM and SLR designed research; NM performed research; NM and AKP analyzed data; NM, SLR and AKP wrote the paper.

Supporting Information Available. Includes Figures S1-S6 and Table 1 . This information is available free of charge via the internet

The authors declare no competing financial interests.