Possible link between higher transmissibility of B.1.617 and B.1.1.7 variants of SARS-CoV-2 and increased structural stability of its spike protein and hACE2 affinity

The Severe Acute syndrome corona Virus 2 (SARS-CoV-2) outbreak in December 2019 has caused a global pandemic. The rapid mutation rate in the virus has caused alarming situations worldwide and is being attributed to the false negativity in RT-PCR tests, which also might lead to inefficacy of the available drugs. It has also increased the chances of reinfection and immune escape. We have performed Molecular Dynamic simulations of three different Spike-ACE2 complexes, namely Wildtype (WT), B.1.1.7 variant (N501Y Spike mutant) and B.1.617 variant (L452R, E484Q Spike mutant) and compared their dynamics, binding energy and molecular interactions. Our result shows that mutation has caused the increase in the binding energy between the Spike and hACE2. In the case of B.1.617 variant, the mutations at L452R and E484Q increased the stability and intra-chain interactions in the Spike protein, which may change the interaction ability of human antibodies to this Spike variant. Further, we found that the B.1.1.7 variant had increased hydrogen interaction with LYS353 of hACE2 and more binding affinity in comparison to WT. The current study provides the biophysical basis for understanding the molecular mechanism and rationale behind the increase in the transmissivity and infectivity of the mutants compared to wild-type SARS-CoV-2.


Introduction
The SARS-CoV-2 (Severe Acute Respiratory Syndrome-Coronavirus detected in December 2019 in the Wuhan province of China has caused the COVID-19 pandemic. As of April 22, 2021, there are more than 143,184,614 confirmed cases, and 3,047,322 people have lost their lives (https://covid19.who.int/). The SARS-CoV-2 belongs to the family of beta corona virus, the same class of viruses that created the pandemic in the past such as SARS-CoV and MERS 1,2 . SARS-CoV-2 possesses a large single-stranded RNA as genetic material and has four main structural components, namely, Envelope protein, Spike protein, Membrane protein and Nucleocapsid 3,4 . The main structural element that enables this virus to attach to the host receptor is the spike glycoprotein, and it also gives the crown-like appearance to the virus, hence named as Coronavirus [5][6][7] . The Spike glycoprotein of SARS-CoV-2 attaches to the human Angiotensin Converting Enzyme (hACE2) receptor and then activated by another human enzyme, Transmembrane Protease Serine (TMPRSS2), to entry inside the host cells 7,8 . Since Spike is the primary target receptor for the entry and the main virulence factor of the virus, various therapeutic drugs and vaccines are being made and tested against it 9 . Although multiple medications such as remdesivir or hydroxychloroquine lopinavir and ritonavir have been recommended by the World Health Organization (WHO) against COVID-19, their efficacy is still the topic of debate [10][11][12]  However, the SARS-CoV-2 cases are still increasing at an alarming rate all over the globe, and the primary rationale behind it is the rapid accumulation of mutations in the SARS-CoV-2.
In the last few months, variants of the SARS-CoV-2 have been reported. Some of them are the Variant of Concern (VOC), which have increased the infectivity or have the potential of immune scape. Almost all the VOCs reported till now have the mutations in the Spike glycoprotein of the virus, which has increased the binding affinity of the virus to hACE2 or has acquired immune scape potential 13,14 . The Lineage B.1.1.7 or 20I/501Y.V1 was detected in the United Kingdom in October 2020. This variant has increased the transmissibility by 40-80%, which has been partially correlated with N501Y mutation in Receptor Binding Domain (RBD) of Spike protein 15 (Figure 1A). In December 2020, B.1.351 was detected in the South African population, which could infect more younger people and had three primary 4 mutations in the RBD of Spike protein, namely, N501Y, K417N and E484K 16 . Similarly, the lineage P.1 detected in January 2021 in the Brazilian population had three mutations of concern in Spike RBD, namely, N501Y, K417T and E484K 14,17 . In our previous study, we had reported that N501Y mutation could enhance the ACE2 affinity and possibly confer resistance towards the antibodies 14   In the recent study it was shown that E484 involved in the hydrogen bonding with S2M11 neutralizing antibody and mutation at E484Q/K abrogated the S2M11 mediated neutralization 19 . The B.1.618 (triple mutant), recently detected in the four Indian states (Maharashtra, Delhi, West Bengal and Chhattisgarh), has been characterized by the deletion of TYR145 and HIS146 as well as E484K and D614G mutation in the Spike glycoprotein 20 ( https://covlineages.org/ ). Although the sudden spike in COVID-19 cases in India is believed to be because of these two new variants and that might be due to the higher binding affinity towards hACE2 and immune escape ability, there are no biochemical and biophysical studies reported on them yet.
In this study, our aim was to study the thermodynamic effects of the mutations in the RBD region of the Spike glycoprotein interacting with hCAE2 and compare that with the wildtype.
We accordingly studied two crucial variants B.1.1.7 and B.1.617, which caused the spike in 5 the COVID19 cases in various countries including India. These two variants possess significant mutations in the RBD domain of the Spike glycoprotein and have a higher infectivity rate. Therefore, to study and compare the dynamics, interactions and binding free energy of Spike protein variants with hACE-2 at the molecular level, we performed the classical Molecular Dynamic (MD) simulations.

MD Simulations
The X-ray crystal structure of SARS-CoV-2, Spike RBD bound with hACE2 was retrieved from Protein Data Bank (PDB) having PDB ID 6M0J. Along with wildtype, two mutants of Spike protein were created, namely, B.1.1.7 (N501Y) and B.1.617 (L452R and E484Q) using the Maestro suite of Schrodinger software 21 . All the three structures were then pre-processed for missing side chains, deleting waters, the addition of hydrogens, hydrogen bond optimization and restrained minimization using the Protein preparation wizard of Schrodinger software 21 . The prepared mutated structures were then simulated for 50ns and the last frame structure was taken for further studies. The following protocol was adopted for the MD simulations of all three prepared structures -each system was solvated with the TIP3P water model in an orthorhombic periodic boundary box. To prevent interaction of the protein complex with its own periodic image, the distance between the complex and the box wall was kept 10 Å. The system was then neutralized by the addition of appropriate number of Na+/Cl− ions depending on the protein-ligand complex using OPLS3e forcefield. Then the energy of the prepared systems was minimized by running 100 ps low-temperature (10K) Brownian motion MD simulation (NVT ensemble) to remove steric clashes and move the system away from an unfavourable high-energy conformation. Further, the minimized systems were equilibrated in seven steps in NVT and NPT ensembles using the "relax model system before simulation" option in the Desmond Schrodinger suite 21  The energy contribution of the mutated residues was then compared with wildtype residues.
The Prime module of Schrodinger software was used for the calculation of the energy contribution of the individual residues. The solvent model used here was Surface Generalized Born (SGB), with variable dielectric enabled, the internal dielectric was 1.00 and solvent dielectric was 80.00 21 .
The following equation was used for the calculation of prime energy of individual residues: Total energy = Covalent total + Non-Bonded total + Other -SGB14|SGB-Torsinol Here, other energy = SGB self, Nonpolar, Hydrogen bond, Packing, self-contact.

Results and Discussion
The mutant spike proteins have a better binding affinity with hACE2 in comparison to wildtype.
The wildtype (WT) Spike-hACE2 complex, along with the prepared and equilibrated B. and WT_Spike-ACE2 (2.82 ± 0.84 Å) as shown in Figure 2A. When RMSF of the simulated complexes were analysed, it was found that the residues number ARG355 to PHE400 of spike protein was more flexible, however, no significant difference was found when WT was compared with both the variants. Also, the fluctuation in the mutant residues was not high, although, in the case of B.1.617, it was found that residue VAL445 had more fluctuations than WT and N501Y mutant ( Figure 2B)   It was noted that the number of hydrogen bonds was similar in case of WT and B.1.617. We further analysed the significant residues, to find out which of them have greater than 30% of the occupancy of hydrogen bond throughout the simulation. It was found that B.1.1.7 and B.1.617 had more residues interaction than WT. In the case of WT, only three residues (TYR453, THR500 and GLY502) of Spike protein were making hydrogen bond with hACE2 for more than 30% of simulation time. In comparison, in the B.1.1.7 Spike mutant, five residues (LYS417, AlA475, ASN487, THR500 and GLY502) were involved, and in B.1.617, there were six residues (LYS417, TYR449, ASN487, TYR489, THR500 and GLY502) of spike protein that had significant hydrogen bond interactions with hACE2 ( Figure S1). It was observed that THR500 and GLY502 were the critical residues interacting significantly with hACE2 in all three complexes. None of the mutated residues in B.1.617 (L452R, E484Q) or N501Y were found to be interacting significantly with hACE2. Hence, it was essential to investigate if these mutated residues of Spike protein had interaction with any other Spike residues or any interaction with hACE2 for any fraction of time. To analyse the changes in the interaction due to mutation, we extracted the three structures at the 50ns interval from all the three simulated complexes. It was found that in the case of B.1.617, in the 50 th ns frame, neither ARG452 nor GLN484 was involved in any polar contact with other residues, while in 100 th ns and 150 th ns frame, it was found that GLN484 was making hydrogen bond contact with SER349 and ASN450 of Spike protein itself, while in WT Spike protein, GLU484 was making hydrogen bond only with SER439. So, it was observed that there was an increase in intra-chain interaction of Spike protein due to mutation of E484Q. Similarly, when B.1.1.7 variant was compared with WT, it was found that due to ASN to TYR mutation at 501 th residue, there was increase in the hydrogen bonding with LYS353 of hACE2 (Figure 3). The increase in the intra-chain interaction in case of B.1.617 indicated that it may interfere in the human antibodies interaction with Spike protein, which may lead to higher infectivity rate.
While, in case of B.1.1.7 the increase in hydrogen bond contact with hACE2 indicated the higher binding affinity of this mutant with hACE2 in comparison to WT.
The MM/GBSA binding free energy has been earlier reported in literature to correlate with the binding affinity between the complexes 21,24 . Therefore, to assess the affinity of Spike  an increase of energy due to mutation as well as change in intra-chain interaction, which may lead to stabilization and change in the orientation of Spike protein. Recently it was also reported about the abrogation of interactions with neutralizing antibody S2M11 due to E484K mutation 19

Conclusion
In this study, we performed the MD simulations and have compared the binding energy, interactions, and change in dynamics of three protein complexes, namely WT, B.1.617 (L452R, E484Q) and N501Y mutant. We showed that mutants have a higher number of significant hydrogen bond interactions with hACE2, and the binding free energy of the mutants is also higher in comparison to WT. In B.1.617, the mutations were energetically favourable as was also observed in terms of intra-chain residue interactions. In B.1.1.7, the mutation leads to an extra interaction with hACE2 (LYS353) in comparison to WT. The detailed molecular level interaction dynamics of Spike-hACE2 interactions and the predicted increased structural stability of its spike protein and hACE2 affinity can be possibly linked to higher transmissibility of B.1.617 and B.1.1.7 variants of SARS-CoV-2