Interactions of the Receptor Binding Domain of SARS-CoV-2 Variants with hACE2: Insights from Molecular Docking Analysis and Molecular Dynamic Simulation

Simple Summary Since the onset of the COVID-19 pandemic in late 2019, SARS-CoV-2 has evolved via genetic changes, resulting in numerous variants of concern (VOCs) and interest (VOIs). Using protein-protein docking and dynamics simulation, we examined the interactions of five SARS-CoV-2 variations’ receptor-binding domains with the human angiotensin-converting enzyme 2 (hACE2) receptor in host cells. A comparison of protein-protein docking and dynamics simulations showed that these point mutations significantly altered the structural behavior of the spike (S) protein, affecting RBD binding to hACE2 at the respective sites. Further research is needed to determine whether these changes affect drug–S protein binding and its potential therapeutic impact. Abstract Since the beginning of the coronavirus 19 (COVID-19) pandemic in late 2019, severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) has been evolving through the acquisition of genomic mutations, leading to the emergence of multiple variants of concern (VOCs) and variants of interest (VOIs). Currently, four VOCs (Alpha, Beta, Delta, and Gamma) and seven VOIs (Epsilon, Zeta, Eta, Theta, Iota, Kappa, and Lambda) of SARS-CoV-2 have been identified in worldwide circulation. Here, we investigated the interactions of the receptor-binding domain (RBD) of five SARS-CoV-2 variants with the human angiotensin-converting enzyme 2 (hACE2) receptor in host cells, to determine the extent of molecular divergence and the impact of mutation, using protein-protein docking and dynamics simulation approaches. Along with the wild-type (WT) SARS-CoV-2, this study included the Brazilian (BR/lineage P.1/Gamma), Indian (IN/lineage B.1.617/Delta), South African (SA/lineage B.1.351/Beta), United Kingdom (UK/lineage B.1.1.7/Alpha), and United States (US/lineage B.1.429/Epsilon) variants. The protein-protein docking and dynamics simulation studies revealed that these point mutations considerably affected the structural behavior of the spike (S) protein compared to the WT, which also affected the binding of RBD with hACE2 at the respective sites. Additional experimental studies are required to determine whether these effects have an influence on drug–S protein binding and its potential therapeutic effect.

At the genomic level, the region downstream of the 5 -untranslated region (UTR), which encompasses ORF1a and ORF1b, accounts for 67% of the total genome, and encodes viral replicase and protease. The remaining genomic region preceding the 3 -UTR possesses four ORFs that encode S, E, M, and N structural proteins, as well as 9 or 10 interspersed ORFs that correspond to accessory proteins [8][9][10].
The spike (S) glycoprotein is a critical component of viral infection. It adheres to the host cell's surface receptor, human angiotensin-converting enzyme 2 (hACE2), allowing viral cellular entry via endosome formation and/or plasma-membrane fusion [11,12]. The S protein (1273 amino-acid residues) exists in a trimeric prefusion form, which comprises an amino (N)-terminal signal peptide (SP) (residues 1-13), the S1 subunit (residues 14-685), and the S2 subunit (residues 686-1273) ( Figure 1). It is thought that the host furin protease is responsible for the cleavage of the S protein into its S1 and S2 subunits [13]. The S1 subunit contains an N-terminal domain (residues 14-305) and a receptor-binding domain (RBD; residues 319-541). The S2 subunit contains a fusion peptide (FP) (residues 788-806), heptapeptide repeat sequence 1 (HR1) (residues 912-84), HR2 (residues 1163-1213), a transmembrane (TM) domain (1213-1237 residues), and a C-terminal cytoplasmic domain (residues 1237-1273). S1 and S2 are responsible for binding with the host-cell receptor and membrane fusion, respectively [14,15]. After entering the cell, the virus releases its genomic RNA into the cytoplasm. Both the 5 ORF1a and ORF1b are immediately translated by host-cell ribosomes, forming precursor polypeptides that are referred to as pp1a and pp1ab. These then undergo autoproteolysis, forming 16 enzymatic nsps, which are assembled into a three-dimensional (3D) supramolecular enzymatic complex known as RNA-dependent RNA polymerase (RdRp). RdRp binds to the +ssRNA genome, forming a replicationtranscription complex (RTC), which mediates these processes. The RTC activity results in the synthesis of sub-genomic mRNAs, whose translation produces a multitude of structural and accessory proteins [16,17].
Since the beginning of the COVID-19 pandemic in late 2019, SARS-CoV-2 has been evolving through the acquisition of genomic mutations, leading to the emergence of multiple specific variants of concern (VOCs) and variants of interest (VOIs) (Figure 2). More recently, several VOIs and VOCs with the potential for increased transmissibility and virulence have been identified, which may enhance disease severity, as well as show resistance to the prevailing vaccination program worldwide [18][19][20][21]. The US government's SARS-CoV-2 Interagency Group (SIG) and the European Centre for Disease Prevention and Control (ECDC) (https://www.ecdc.europa.eu/en/covid-19/variants-concern (accessed on 28 August 2021)) regularly evaluate new evidence on variants discovered through epidemic intelligence, rules-based genomic variant screening, or other scientific sources on a regular basis. They classified SARS-CoV-2 variants into three categories: variants of interest (VOI), variants of concern (VOC), and variants of high consequence (VOHC). They defined VOI as a genetic variant associated with altered receptor binding, reduced neutralization by previous infection or vaccination antibodies, reduced treatment efficacy, potential diagnostic impact, or predicted increase in transmissibility or disease severity. This includes Eta (B. The global scientific community is rigorously following the three most rampantly spreading SARS-CoV-2 VOCs identified in the United Kingdom, South Africa, and Brazil [22][23][24]. The first SARS-CoV-2 VOC, B.1.1.7 (also known as 20I/501Y.V1, VOC202012/01, and the UK variant) was identified in the United Kingdom in December 2020 (Kirby, 2021). Within a few months, this variant became prevalent in the United Kingdom and has since expanded to 114 other nations worldwide [25,26]. This variant has around 17 mutations in the S protein alone, including a crucial nonsynonymous mutation (N501Y) at position 501 in the RBD in which asparagine (N) has been replaced with tyrosine [27].
Mutations in the RBD can help to enable strong affinity and binding capacity to hACE2, leading to higher transmissibility [28,31,32]. Furthermore, these mutations may result in a reduction in antibody neutralization [33,34]. This could have a major impact on the efficacy of existing vaccines. This has the potential to have a significant impact on the efficacy of the existing vaccines [17,35]. As a result, global monitoring of the continuing genomic changes in SARS-CoV-2 is critical for identifying areas linked with drug resistance and vaccine evasion in order to create successful antiviral medicines. Several drug-repositioning studies and compound evaluations have been conducted to discover novel antiviral medicines against SARS-CoV-2 utilizing experimental and theoretical/computational methodologies [36][37][38].
Furthermore, genetic differences in the active/binding site region of molecular targets can have a significant impact on the binding mechanism and affinity of the lead compounds. Consequently, this will be of relevance to the efficacy of previously investigated promising candidates. The current study modelled the structure of the S protein of newly found variants and examined novel antiviral medications, utilizing a variety of computational drug design methodologies, with the aim of identifying a viable treatment for COVID-19.

Results and Discussion
The SARS-CoV-2 that causes COVID-19 continues to mutate. The majority of variants created through alterations in amino acids in the RBD have been found to be less infectious [39], but certain variations investigated have been resistant to some neutralizing antibodies [40]. A mutation that occurred outside the RBD region (D614G) was reported to be more infectious, although there was no evidence that this variant was resistant to neutralizing antibodies [39]. Mutations that occur in the RBD region are likely to have an impact on the attachment of the virus to hACE2 [41]. Therefore, the present study focused on analyzing the interaction between several RBD variants with hACE2 through molecular docking and MD simulation studies. Along with the RBD Gamma, Delta, Beta, and Alpha variants, the United States (US/lineage B.1.429/Epsilon) variant and the WT were examined in this study.

Analysis of the Modeled RBD Structures
The 3D structures of all studied RBDs were modeled and minimized using the SWISS-MODEL web server. The results were validated using the PROCHECK service, and are presented in Ramachandran plots, as shown in Figure 2 and Table 1. The WT RBD (PDB ID: 6M0J) crystal structure served as a control. The Ramachandran plot was divided into four types of areas by ProCheck: most favored, additional allowed, generously allowed, and disallowed [42]. If the percentage of non-glycine residues in the disallowed region was <15%, the protein structure was considered to be of high quality: the lower the percentage, the higher the quality of the protein structure [43,44]. Based on the results of the analysis, five models of the structure of RBD variants could be accepted based on the non-glycineresidue data in the disallowed region. The percentage of amino-acid residues in the most favored region and the disallowed region indicated the quality of the structure based on its geometry: the greater the percentage of amino-acid residues in the most favored region and the lower the percentage of residues in the disallowed region, the better the quality of the structure.

Molecular Docking Analysis
Protein-protein docking uses the biomolecular interaction between two molecules to determine the strength of a complex. Protein-protein docking simulation of hACE2 and RBD variants was performed on the HDOCK webserver. The results are visualized as illustrated in Figures 3-8. The Biovia Discovery Studio was used to generate 3D visualizations, while LigPlot+ and PDBsum were used to generate 2D visualizations.
As shown in Table 2, the present analysis of the interaction revealed the formation of H-bonds between hACE2 and residues in the active site of RBD Gamma: specifically Tyr449 Tyr449  Tyr449  Tyr449  Tyr449  Tyr449  Tyr41  Thr500  Thr500  Thr500  Thr500  -Thr500  Gln42  Gly446  -Gly446  Gly446  Gly446  Gly446  Tyr449  Tyr449  Tyr449  Tyr449  Tyr449  Tyr449  ----Gln498  -Tyr83 Asn487 Considering the three residues required for the formation of the hACE2-RBD interaction, the Gamma and Beta variants had hydrogen bonds on the Gln493 residue, while Phe486 had hydrophobic contact. The H-bond was located on the Gln493 residue in the Delta and Epsilon variants, while the hydrophobic contacts were located on Phe486 and Asn501. In the Alpha and WT variants, the three residues showed hydrophobic contact alone. Of the six variants analyzed, four showed H-bonding at the Gln493 residue while two variants showed hydrophobic contact. Six variants made only hydrophobic contact with the Phe486 residue, whereas four variants made hydrophobic contact with Asn501. Salt bridges were formed during the interaction of hACE2 with RBD WT Asp30-Lys417

Molecular Dynamics Simulation of Protein-Protein Complexes
The stability and pattern of binding of protein-protein interactions (PPIs) change according to physiological circumstances and time [49]. To provide insight into the dynamic status and estimation of different bond forms between two proteins, long-term molecular dynamics (MD) simulation of the selected protein-protein binding complexes (hACE2 and RBD) were performed. While examining the structural stability of the protein complexes during the 100-ns MD simulation, it was discovered that the stable structure was preserved, as illustrated in Figure 9. The Alpha variant was the least stable complex, with a fluctuation of up to 4 Å. By contrast, the results for the Rg in all structures indicated that structural stability did not deteriorate during the MD simulation. MD simulations revealed that SARS-CoV-2 VOCs interacted more strongly with the host receptor hACE2. This finding was consistent with prior research indicating that the N501Y mutation, which affected the primary contact residue in RBD, increased the viral attachment proclivity for hACE2 [50,51]. The combination of mutations S477N/E484K, E484K/N501Y, and K417T/E484K/N501Y increased viral affinity for hACE2 and also antibody resistance [52].
As illustrated in Figure 10A, the most flexible regions during the 100-ns simulation were the terminal parts of S proteins. While no significant change was observed in the ACE2 protein compared to the WT during the simulation, we observed significant conformational changes in the S protein, especially in the Delta variant. As seen in Figure 10C,D, we observed residual fluctuations in areas 1, 2, 3, 4, and 6. This might have signified that the immune response was delayed, as the antibody recognition sites may also have changed. Recent studies have shown that vaccine efficacy was lower in the Delta variant, which was consistent with our findings [53]. While no significant fluctuation was observed in other variants, it was clearly seen in residues in the range of 477 to 480 in the Gamma variant.
As the hydrogen-bond stability of the complexes was examined during the 100-ns MD simulation, the Tyr83-side/Asn487-side, Gln493-side/Glu35-side, Thr500-side/Asp355side, and Gly502-main/Lys353-main hydrogen bonds, which played a role in the binding of the S to hACE2 proteins, were found to be critical in all variants (Table 3). This was consistent with the discovery of Khan et al. [31] that hACE2-RBD formed hydrogen bonds between Glu35-Gln493 and Lys353-Gly502. However, it was observed that hydrogen-bond formation between Tyr505-side/Glu37-side occurred at a high rate in WT and the Epsilon variant. The Ser19-side and Ala475-main H-bonds, which played a role in the formation of the complex, were effective in the WT, Gamma, and Delta variants, but not in the others. Ashwaq et al. [54] reported similar findings.
Notably, the hydrogen bond between Tyr449-side and Asp38-side was ineffective in the WT but highly effective in the other variants. While a hydrogen bond between Tyr501side and Asp38-side was observed in the Gamma variant, this interaction did not occur in the others. A hydrogen bond was established between Gln498-side and Asp355-side in the Delta variant; However, this formation was not observed in the other variants.

Binding Free Energy of Protein-Protein Complex Simulation Trajectory
As the binding free energies of the S protein variants with hACE2 were examined, we observed the strongest result in the Delta variant (Table 4, and Figure 11). Recent investigations have demonstrated that the Delta variant is more effective than others at binding to hACE2 [55]. Additionally, this variation appears to spread at a considerably faster rate than the other three VOCs [25]. Residue 478, which is located in the flexible loop, is most likely to come into contact with ACE2. This appears to strengthen the interaction between the Delta variant and ACE2, which accounts for the dramatic increase in the proportion of infectivity of this variant [56]. The interaction between RBD WT and hACE2 exhibited the lowest binding affinity. This implied that the existence of mutations in RBD resulted in a greater binding affinity for hACE2. Laffeber et al. [57] demonstrated experimentally that RBD carrying the N501Y mutation had a sevenfold greater affinity for the hACE2 receptor than WT RBD. However, it was reported that mutations at the position K417N/T decreased the binding affinity [58]. Furthermore, a K-to-N mutation significantly reduced the binding affinity between N417/Y501-RBD and ACE2 when compared to the Y501-RBD to ACE interaction [59].  Figure 11. MM-PBSA binding free energy estimation of Spike-ACE2 complexes.
As the residual contribution to the binding free energies was examined, no significant changes were observed in the S protein variants, but significant differences were observed in the hACE2 protein ( Figure 12). We observed that many residues of the S protein had lower binding free energies than the hACE2 residues, with the most significant change being observed in the interaction of the Delta variant. In particular, we calculated that Asp30, Glu35, and Glu37 in hACE2 were significantly separated from the others as the residues with the lowest binding free energy in the Delta variant. We observed that the His34 residue had high binding free energy in the Alpha, Epsilon, and Beta variants, while the value for Glu75 was significantly lower in the Gamma, Delta, and Beta variants. Chakraborty et al. [60] discovered that His34 of ACE2 had the second highest energetic contribution (4.68 kcal/mol) when compared to Arg403 of RBD via a direct hydrogen bond. While it was observed that Lys417, one of the S protein residues, was ineffective in the Gamma and Beta variants, we found that Leu452 in the Delta and Epsilon variants had a relatively low binding free energy of <200 kJ/mol and was largely dissociated from the other variants. In the WT RBD, Lys417 formed a salt bridge with Glu30 of RBD [31]. Salt bridges, like disulfide bonds, can act as keystone interactions [61].

Preparation of the Macromolecules
The three-dimensional (3D) crystal structure of WT SARS-CoV-2 S RBD bound with hACE2 was retrieved from the Protein Data Bank (https://www.rcsb.org/ (accessed on 26 May 2021)) with PDB ID 6M0J. The 3D structures of the RBD variants were modeled and minimized using the SWISS-MODEL web server (https://swissmodel.expasy.org/ (accessed on 28 May 2021)) [62,63]. The structural accuracy of the RBD protein model was analyzed using the PROCHECK server (https://saves.mbi.ucla.edu/ (accessed on 28 May 2021)) [64]. The analytical results are presented in a Ramachandran plot.

Molecular Dynamics (MD) Simulations
All simulations were run via GROMACS 2021.2 software using the Leap-Frog integration at 2-fs intervals [69]. The MD simulation system was built under periodic boundary conditions (PBCs) with the 'rhombic dodecahedron' scheme. The distance from the protein complex to the corner of the cube was set to 1.2 nm. Amber99SB-ildn was chosen as the force field and 'TIP3P' was preferred as the water model [70]. The system was neutralized by adding 0.15 mM NaCl. Energy minimization was done by the Steepest Descent algorithm in 50,000 steps (minimization stopped when the maximum force was <10.0 kJ/mol). To bring the MD system to the equilibrium phase, 100-ps NVT (constant number of particles, volume, and temperature) and 1-ns NPT (constant number of particles, pressure, and temperature) simulations were performed. In the NVT and NPT stages, all bonds and heavy atoms were restricted by the LINCS (LINear Constraint Solver) algorithm. In the NVT phase, the V-rescale coupling algorithm was the preferred temperature coupling algorithm. In the NPT phase, V-rescale was preferred as the temperature coupling algorithm, and the Parrinello-Rahman with isothermal compressibility was preferred as the pressure coupling algorithm. The temperature and pressure were set to 310 K and 1 atm, respectively. In the MD production phase, unlike the NPT phase, the atomic restrictions were removed and a 100-ns MD simulation was carried out. The Verlet algorithm was used as the cutoff scheme. The Particle-Mesh Ewald (PME) method was preferred for long-range interactions. The cutoff values of the electrostatic and Van der Waals interactions were both set to 1.2 nm.

Analysis of MD Simulations
The root-mean-square deviation (RMSD) and hydrogen bond analyses were performed with the VMD program (ref). The root-mean-square fluctuation (RMSF) and radius of gyration (Rg) analyses were performed with the Gromacs RMSF version 2019.2 and gyrate tools, respectively. All analyses were plotted using GraphPad Prism version 9.1.2 for Windows (GraphPad Software, San Diego, CA, USA; www.graphpad.com (accessed on 10 June 2021)).

Free-Energy Calculations
The binding free energies were calculated using the molecular mechanics Poisson-Boltzmann surface area (MMPBSA) approach, which is an end-point method [71]. The binding free energies, including the entropy contribution, were obtained by taking 100 snapshots at 100-ps intervals in the last 10 ns. The polar component of desolvation was calculated using PB models. In the PB calculations, the partial charges of the proteins were taken from the forcefield parameters. The solvent-accessible surface area (SASA) was preferred as a non-polar contribution. The vacuum electrostatic dielectric constant and the solvent dielectric constant were set to 2 and 80, respectively. The g_mmpbsa tool was used for the MMPBSA calculation [72]. The binding free-energy calculations were made according to the following equation: ∆G Solv = ∆G polar + ∆G nonpolar (2)

Conclusions
The five variants of the RBD of SARS-CoV-2 S-proteins had distinct mutations at their binding sites with hACE2. The present study investigated the PPIs of mutant RBD (K417N, E484K, and N501Y in B.1.351 (Beta); L452R and E484Q in B.1.617 (Delta); K417T, E484K, and N501Y in P.1 (Gamma); L452R in B.1.429; and N501Y in B.1.1.7 (Alpha)) with hACE2. The increased binding affinity of B.1.617 (IN), P.1 (BR), and B.1.351 (SA), (CA) strain with hACE2 when compared to WT indicated the possibility of high transmissibility and rapid spread. Calculation of the binding free energy of the protein-protein complex revealed that the Delta variant had the lowest value, indicating that it bound to hACE2 more strongly than the others. Furthermore, it was predicted that mutations in the B.1.617 (IN)/Delta variant, which caused more residual fluctuations than other variants, may play a role in antibody escape. The MD simulation analysis demonstrated a stable interaction between RBD variants and hACE2. This is expected to have an effect on drug discovery efforts by elucidating residues that could be targeted for disrupting this interface. Overall, this study laid the groundwork for the development of antibodies, vaccines, and drugs against new SARS-CoV-2 variants.

Data Availability Statement:
The authors confirm that the data supporting the findings of this study are available within the article. Raw data were generated on webservers provided in the Section 3. Derived data supporting the findings of this study are available from the corresponding author on request.