All-Atom Simulations of Human ACE2-Spike Protein RBD Complexes for SARS-CoV-2 and Some of its Variants: Nature of Interactions and Free Energy Diagrams for Dissociation of the Protein Complexes.

The spike protein of SARS-CoV-2 is known to interact with the human ACE2 protein via its receptor binding domain (RBD). We have investigated the molecular nature of this interprotein interaction and the associated free energy diagrams for the unbinding of the two proteins for SARS-CoV-2 and some of its known variants through all-atom simulations. The present work involves generation and analysis of 2.5 μs of unbiased and 4.2 μs of biased molecular dynamics trajectories in total for five explicitly solvated RBD-ACE2 systems at full atomic level. First, we have made a comparative analysis of the details of residue-wise specific interactions of the spike protein with ACE2 for SARS-CoV-1 and SARS-CoV-2. It is found that the average numbers of both direct interprotein and water-bridged hydrogen bonds between the RBD and ACE2 are higher for SARS-CoV-2 than SARS-CoV-1. These higher hydrogen bonded interactions are further aided by enhanced nonspecific electrostatic attractions between the two protein surfaces for SARS-CoV-2. The free energy calculations reveal that there is an increase in the free energy barrier by ∼1.5 kcal/mol for the unbinding of RBD from ACE2 for SARS-CoV-2 compared to that for SARS-CoV-1. Subsequently, we considered the RBDs of three variants of SARS-CoV-2, namely N501Y, E484Q/L452R, and N440K. The free energy barrier of protein unbinding for the N501Y variant is found to be ∼4 kcal/mol higher than the wild type SARS-CoV-2 which can be attributed to additional specific interactions involving Tyr501 of RBD and Lys353 and Tyr42 of ACE2 and also enhanced nonspecific electrostatic interaction between the protein surfaces. For the other two mutant variants of E484Q/L452R and N440K, the free energy barrier for protein unbinding increases by ∼2 and ∼1 kcal/mol, respectively, compared with the wild type SARS-CoV-2, which can be attributed to an increase in the number of interprotein hydrogen bonds for the former and also to enhanced positive electrostatic potential on the RBD surfaces for both of the variants. The successive breaking of interprotein hydrogen bonds along the free energy pathway of the unbinding process is also found out for all five systems studied here.

Introduction

The Severe Acute Respiratory Syndrome Coronavirus 2 (SARS-CoV-2) has been one of the most infectious and deadliest viruses to hit the mankind all over the world.[1−6] SARS-CoV-2, first reported in December 2019,[1−4] is the latest of the three human coronaviruses (CoVs) that have been reported so far. The other two are SARS-CoV (severe acute respiratory syndrome coronavirus, now referred to as SRS-CoV-1), which was first reported in 2003,[7−10] and MERS-CoV (Middle East respiratory syndrome coronavirus), first reported in 2012.[11−13] The SARS-CoV-2 viruses possess spherical structures with diameters ranging from 60 to 140 nm.[4,14] A long RNA polymer is tightly packed at the center of the sphere which is surrounded by a protective nucleocapsid. This is further protected by a membrane envelope made of lipids with inserted spike (S), membrane (M), and envelope (E) proteins. The spike protein consists of an extracellular N-terminus, a transmembrane (TM) domain anchored in the viral membrane, and a short intracellular C-terminal segment.[15,16]

The viral entry of SARS-CoV-2 takes place through its binding to the human angiotensin-converting enzyme 2 (ACE2) receptor via its spike protein. Once the virus interacts with the host cell, extensive structural rearrangement of the spike protein occurs, which allows the virus to fuse with the host cell membrane.[17−19] The spike protein exists in trimeric form with 1273 amino acids per monomer. The monomeric form of the spike protein of SARS-CoV-2 consists of an N-terminus signal peptide (residues 1–13), the S1 subunit (residues 14–685), and the S2 subunit (residues 686–1273)[20,21] (Figure a). The S1 and S2 subunits are responsible for receptor binding and membrane fusion, respectively. In the S1 subunit, there is an N-terminal domain (NTD; residues 14–305), a receptor-binding domain (RBD; residues 319–541) and two subdomains (SD1 and SD2).[1,21−23] The S2 subunit consists of fusion peptide (FP; residues 788–806), heptapeptide repeat sequence 1 (HR1; residues 912–984) and sequence 2 (HR2; residues 1163–1213), transmembrane region (TM; residues 1214–1237), and intracellular domain (IC; residues 1238–1273).[1,20−23]

Figure 1

(a) Residue sequence of different domains of the spike protein of SARS-CoV-2.[20] Here, the domains such as NTD (N-terminal domain), RBD (receptor binding domain), SD1 (subdomain 1), SD2 (subdomain 2), FP (fusion peptide), HR1 (heptad repeat 1), HR2 (heptad repeat 2), TM (transmembrane region), and IC (intracellular domain) are shown. (b) Structure of the RBD domain of spike protein from SARS-CoV-2 showing the secondary structural elements.[27] The RBM part is shown in red color. (c) Structure of RBD bound to ACE2 for SARS-CoV-2.[27] Two subdomains of ACE2 are colored in gray and tan, respectively.

The receptor binding domain (RBD) of the spike protein has been reported to be involved in binding with the host ACE2.[24−30] The RBD is composed of five-stranded antiparallel β-sheets (β1, β2, β3, β4, and β7) connected by α-helices and loops (Figure b),[27] and it runs from residues 319 to 541 of the spike protein for SARS-CoV-2. An extended insertion from residues 438 to 506 between two β-sheets (β4 and β7) contains most of the residues that bind to ACE2 (Figure c)[27] and is referred to as the receptor-binding motif or RBM.

The human angiotensin converting enzyme 2 (ACE2) is ellipsoid shaped and possesses largely α-helical structure.[31−34] The structure is divided into approximately equal-sized subdomains 1 and 2. The two subdomains are not simply the N- and C-terminal halve; rather, the polypeptide chain crosses between the two subdomains five times (Figure c). Subdomain 1 consists of the residues 19–102, 286–431, and 521–580, while subdomain 2 contains the residues 103–285, 432–520, and 581–612. The subdomain 1 of ACE2 is mainly responsible for binding with the RBD of spike protein from SARS-CoV-2. The RBD interacts with the first long helix (H1: residues 21–53) and loop residues 82–83 at the end of the second long helix of ACE2, and also with a β-hairpin loop (residues 352–354) and an α-helix loop (residues 325–330) of subdomain 1 of ACE2.[27] Both the spike protein and ACE2 also contain several glycans.[4,35−38] These glycans act as a shield to the host immune response.[39−41] Also, several studies have revealed that these glycans play important roles in developing favorable conformation of the receptor binding domain (RBD) for binding with ACE2.[35,39]

The spike proteins of SARS-CoV-2 and SARS-CoV-1 possess about 80% sequence identity.[17,18,20] Both the spike proteins function in a similar manner in terms of their human-infection and pathogenic mechanisms.[24−29] The structures of the RBD-ACE2 complexes for spike proteins from SARS-CoV-1 and SARS-CoV-2 can be greatly superimposed on each other. The superimposed structures of the RBD-ACE2 complexes are shown in Figure S1 of the Supporting Information. The primary sequence alignment of the RBD domain of spike proteins from SARS-CoV-1 and SARS-CoV-2 shows that the RBD is almost conserved (Figure S2). Therefore, the favorable binding interactions between the RBD and ACE2 is a prerequisite for viral infection caused by both SARS-CoV-1 and SARS-CoV-2.

Several mutated variants of SARS-CoV-2 have also been reported over past two years.[42−50] These mutations increase the transmissibility of the infection and have also been reported to show increased resistance toward vaccines.[51] In particular, the N501Y mutation is present in the RBD domain where the spike makes initial contact with the host ACE2 protein. The N501Y, or the Alpha variant, was first sequenced in April 2020 and is linked with a SARS-CoV-2 variant that is an independent lineage from B.1.1.7.[44,45,48] Another variant of the lineage B.1.617, the so-called Kappa variant, has double mutations (L452R and E484Q) in the RBD domain.[49] The N501Y and L452R mutations are also present, respectively, in the Beta and Delta variants of the virus. There is another variant (N440K) of the lineage B.1.36 where Asn440 of RBD domain is mutated with Lysine.[50] We note that the N440K mutation is also there as one of the multiple mutations present in the Omicron variant of the SARS-CoV-2 virus.[52]

The binding of RBD with human ACE2 has been studied in many recent studies through computational means.[53−77] Although these studies have greatly contributed to our understanding of the RBD-ACE2 binding process, many details of the molecular mechanism of the binding, the process of unbinding and associated free energy diagrams, and their variations with mutation still remain open. For example, while the experimental work of ref (20) reported that the receptor binding domains of SARS-CoV-2 and SARS-CoV-1 bind with ACE2 with similar affinities, the biomechanical force measurements and computational studies of refs (30 and 53−55) reported that SARS-CoV-2 binds with ACE2 with higher affinity than SARS-CoV-1. Some of the existing studies[56,57,59,60,69] have also looked at effects of mutations on RBD-ACE2 binding through computational means but many issues still need to be resolved. The work of refs (56 and 57) concluded a greater flexibility of the RBD for variants of SARS-CoV-2 which possibly lead to more stable RBD-ACE2 complex. However, the work of ref (59) found similar binding affinities for all the variants studied in that work. The experimental work of refs (78 and 79) reported a significantly greater binding affinity for some of the variants of SARS-CoV-2 which is also supported by calculations of electrostatic potentials.[80] An accurate estimate of the binding affinities at room temperature could be obtained from calculations of free energy diagrams of the binding or unbinding of the RBD-ACE2 complexes for different RBDs through full atomistic simulations of the solvated proteins. The present study makes a contribution toward this end.

In the current work, we have investigated the complexes of human ACE2 and RBD of spike proteins of SARS-CoV-2, some of its variants, and also SARS-CoV-1 by using all-atom molecular dynamics simulations. Specifically, we have calculated the free energy changes during the unbinding of RBD from ACE2 and also the nature of interprotein interactions, both direct and water mediated, of the RBD-ACE2 interfaces for the different systems considered here. Since the spike protein is a potential drug target and also acts as antigen of the virus, the present study will be useful for design of new inhibitors and vaccines for prevention of the disease caused by SARS-CoV-2.[81,82] An important objective of the current study has been to study how the mutations enhance the infectivity of SARS-CoV-2. The present work involves generation and analysis of 2.5 μs of unbiased and 4.2 μs of biased molecular dynamics trajectories for five solvated RBD-ACE2 systems at full atomic level with an average system size of about 230 000 atoms.

We have organized the rest of the paper as follows. In Section , the details of model preparation and the protocols of molecular dynamics simulations are presented. The computational details of biased simulations involving umbrella sampling method for calculations of the free energy diagrams of the unbinding process are also included here. In Section , the current simulation results of the interprotein structure and interactions of the RBD-ACE2 complexes, the free energy diagrams of the unbinding of RBD from ACE2, and also a comparison of the free energy diagrams of protein unbinding for SARS-CoV-1, SARS-CoV-2, and also three variants of SARS-CoV-2 are presented. The results of sequential breaking of interprotein hydrogen bonds along the free energy pathway of the unbinding process are also discussed in this section for all five systems studied here. Finally, a brief summary and conclusions of the present work are included in Section .

Computational Details

The structure of the spike receptor binding domain (RBD) bound with human angiotensin-converting enzyme 2 (ACE2) for SARS-CoV-2 is taken from crystal structure of PDB code 6M0J.[27] In this crystal structure, the RBD domain has 229 residues from Thr333 to Gly526 and the ACE2 enzyme has 603 residues from Ser13 to Asp615. The experimental structure also has some N-acetyl-β-glucosaminide (NAG) glycans linked with both ACE2 and RBD.[27] Among them, four glycans are linked to residues Asn53, Asn90, Asn322 and Asn546 of ACE2. One NAG is linked to the Asn343 residue of RBD. The structure of ACE2-RBD complex with the glycans is shown in Figure S3 for SARS-COV-2. In the current simulations, we have considered only those glycans that are present in the experimental crystal structure.[27] The full glycosylation of the ACE2-spike protein complex is missing in the current simulations. It has been shown both computationally and experimentally that the glycosylation plays an important role in the binding of SARS-CoV RBD with ACE2.[28,35,83,84] In the context of the current work, we note that the interactions of the RBD with the ACE2 is probably not specific; hence, it may not significantly affect the results for the differences between the different variants of SARS-CoV-2 or between SARS-CoV-1 and SARS-CoV-2 spike proteins. The system also contains one Zn2+ ion bound in the ACE2 protein.[27] The parameter and topology files are generated using the CHARMM-GUI web server.[85] The AMBER-FF14SB force field is used to prepare the topology files[86] for this system and also of other systems described later. The current system is neutralized by adding 24 Na+ ions and solvated using TIP3P[87] water. We added 63 712 water molecules to solvate the protein complex so that the full simulation system contains 203 806 atoms in total. The same water model was also used to solvate protein complexes of other systems of the current study which are described below.

The structure of the spike RBD domain of SARS-CoV-1 bound with human ACE2 is taken from the crystal structure of PDB code 3SCI.[88] In this structure, the RBD domain has 228 residues from Cys323 to Glu502 and the ACE2 enzyme has 603 residues from Ser13 to Asp615. No glycan is linked with the structure. The system also contains one Zn2+ ion bound in the ACE2 protein.[88] Again, the parameter and topology files are generated using CHARMM-GUI web server[85] and the AMBER-FF14SB force field is used to prepare the topology files.[86] The system is neutralized by adding 23 Na+ ions and the protein complex is solvated using 71 118 water molecules. The full simulation system contains 225 629 atoms.

The structure of RBD of the N501Y variant of SARS-CoV-2 bound to human ACE2 is prepared from the crystal structure of 6M0J.pdb[27] by mutating Asn501 with tyrosine in the RBD domain and generating the topology files by using the CHARMM-GUI web server.[85] We added 24 Na+ ions to neutralize the system and 70,981 water molecules to solvate the protein complex. The full simulation system contains 225 620 atoms. The structures of RBDs of the double mutant (E484Q and L452R) and N440K variants of SARS-CoV-2 bound to human ACE2 are also prepared in a similar manner from the crystal structure of PDB code 6M0J.[27] The double mutant system is neutralized by adding 22 Na+ ions and the protein complex is solvated using 77 834 water molecules so that full simulation system contains 246 117 atoms. The N440K mutant system is neutralized by adding 23 Na+ ions and the protein complex is solvated using 77 834 water molecules. In total, there are 246 179 atoms in this simulation system.

We performed fully atomistic unbiased molecular dynamics simulations of all five systems by using the AMBER18 suite of programs.[89] The simulations were carried out using the pmemd.cuda code in AMBER18, which greatly accelerates the explicit solvent atomistic simulations.[90,91] The following protocols were followed for equilibration of all the systems. First, the energy of all water molecules was minimized for 2000 steps (1000 steps of steepest descent and 1000 steps of conjugate gradient minimization) while the rest of the system was restrained. This was followed by 2000 steps (1000 steps of steepest descent and 1000 steps of conjugate gradient minimization) of energy minimization of the entire system without any constraint. The system is then gradually heated to 300 K through simulation for 500 ps. Subsequently, isothermal–isobaric (NPT) simulation was run for 1 ns with gradual relaxation of positional restraints of water molecules. This step was followed by simulation for a 1 ns run, where the entire system was relaxed under NPT condition without any restraint. The temperature was maintained at 300 K. The Langevin thermostat with a collision frequency of 1.0 ps–1 was used. The pressure was maintained at 1 bar with a relaxation time of 2 ps. Subsequently, each system was further equilibrated for 5 ns under NVT condition. All simulations were carried out with a time step of 1 fs. The shake algorithm was used to restrain the bonds linking the heavy atoms and hydrogen atoms. The electrostatic interactions were calculated using the particle-mesh Ewald sum method. A nonbonded cutoff of 10 Å was used in calculations of the Lennard-Jones interactions and real-space part of electrostatic interactions of the Ewald summation method. In the production phase of each system, the simulation was run for 500 ns under NVT condition at 300 K for calculation and analysis of various properties. Since the interactions between RBD and ACE2 are dominated by specific hydrogen bonds and nonspecific electrostatic interactions, we have made detailed calculations of these interactions from the simulation trajectories. The electrostatic potential surfaces are calculated for an equilibrated configuration for each system using the web server APBS-PDB2PQR software[92,93] based on the Poisson–Boltzmann equation using the Adaptive Poisson–Boltzmann Solver (APBS)[92] and visualized using the PyMOL software tool.[94] Both direct interprotein and water-bridged hydrogen bonds between the protein residues are investigated to understand the roles of hydrogen bonds in the binding and unbinding of the RBD-ACE2 complexes for all the systems studied here.

Subsequently, we used the umbrella sampling method[95] to calculate the free energy diagrams for the unbinding of RBD from human ACE2. The separation between the centers of mass (COMs) of RBD and ACE2 along the normal (x) direction is used as the reaction coordinate. The COM separation is gradually varied from 40 to 70 Å with an increment of 0.25 Å. A harmonic force constant of 10 kcal mol−1 Å−2 is used for the restraining potential in each window. Each umbrella sampling window was first equilibrated for 2 ns and then run for another 5 ns for the production purpose. The reweighted equilibrium probability distributions for all the umbrella sampling windows and the free energy profiles were obtained through the Weighted Histogram Analysis Method (WHAM)[96] using the code of ref (97). The error estimates of the free energy profiles were performed over 100 rounds of bootstrapping analysis.[98] We used the “Monte Carlo Bootstrap Error Analysis” technique in the code[97] for error estimation.

Results and Discussions

Contact Maps between RBD and ACE2

We have analyzed the contacts between heavy atoms of RBD and ACE2 with a cutoff distance of 3.5 Å for SARS-CoV-1, SARS-CoV-2 and the three mutated systems considered in this work. The distributions of the number of contacts between the heavy atoms of ACE2 and RBD with cutoff distance of 3.5 Å for SARS-CoV-1 and SARS-CoV2 (wild type) and its variants are presented in parts a and b of Figure . The fluctuations of the number of contacts between the heavy atoms of ACE2 and RBD along the simulation trajectories are shown in Figure S4 of the Supporting Information for both SARS-CoV-1 and SARS-CoV-2 and its variants. The average number of contacts is ∼42 for SARS-CoV-2 and ∼37 for SARS-CoV-1 which means the RBD and ACE2 are more closely bound for SARS-CoV-2. The contact numbers for the three variants of SARS-CoV-2 (Figure b) are also found to be higher than that of SARS-CoV-1. The average contact number (44) is found to be the highest for the N501Y mutant, which shows a tighter binding of the RBD-ACE2 complex for this variant. The distributions of the interface distances between ACE2 and RDB for the current systems are presented in parts c and d of Figure , and their fluctuating values along the simulation trajectories are shown in Figure S5. In these calculations, the interface heavy atoms of ACE2 that lie within 3.5 Å of an RBD heavy atom were selected. The interface heavy atoms of RBD were also identified in a similar manner and, subsequently, the center-of-mass distances between the interface heavy atoms of ACE2 and RBD were calculated. The average interface distance appears to be ∼4.2 Å for SARS-CoV-2 and ∼4.6 Å for SARS-CoV-1. This again shows that the RBD-ACE2 complex for SARS-CoV-2 is more closely bound than that for SARS-CoV-1. The average interface distance is also found to be ∼4 Å for all three variants of SARS-CoV-2, with the minimum distance being somewhat lower for the N501Y variant which means a tighter binding of RBD of this variant to ACE2.

Figure 2

Distributions of the number of contacts between the heavy atoms of RBD and ACE2 with a cutoff distance 3.5 Å for (a) SARS-CoV-1 and SARS-CoV-2 and (b) three variants of SARS-CoV-2. Distributions of interface distances between ACE2 and RBD are shown for (c) SARS-CoV-1 and SARS-CoV-2 and (d) three variants of SARS-CoV-2.

We have also calculated the contact maps between the residues of RBD and ACE2 for all the five systems and the results are presented in Figure . These maps show the binding regions of RBD and ACE2 that interact with each other. The results of these contact maps show that there are mainly two regions of human ACE2 protein which effectively interact with the RBD of spike protein. These are residues Ser19 to Glu87 of α-helix 1 and α-helix 2 (highlighted using boxes named as A, B, and C) and residues Val318 to Phe400 of β-sheet, β-turn, and α-helix (highlighted using boxes named as D and E). The contact maps also show that there are three regions of RBD which take part in the interactions with ACE2. These are residues 470 to 510, 435 to 460, and 395 to 415 for SARS-CoV-2 and its variants. The results of these contact maps show that mainly the RBM part (residues 438 to 506) of RBD of the spike protein interacts with ACE2. In parts a–e of Figures , we have highlighted the interactions with boxes (named A, B, C, D, E, and F) for RBDs of all the wild type and mutant systems studied here. It is found that there are mainly five interacting regions except SARS-CoV-1 where four interacting regions are found in its RBD. No interaction is found in the region named C for SARS-CoV-1 (Figure a). Also, the number of contacts is higher for SARS-CoV-2. We also observed an additional interacting region, referred to as F for the N501Y variant (Figure c). The interacting regions of RBDs of the other two variants are found to be the same as found for the wild type SARS-CoV-2.

Figure 3

Contact maps between the residues of RBD and ACE2 for (a) SARS-CoV-1, (b) SARS-CoV-2, (c) SARS-CoV-2 (N501Y), (d) SARS-CoV-2 (E484Q/L452R), and (e) SARS-CoV-2 (N440K). The cutoff distance is 3.5 Å for the heavy atoms of RBD and ACE2.

Coulomb Potentials on the Surfaces of RBD and ACE2: Nonspecific Electrostatic Interactions

The interactions between the RBD domain and human ACE2 enzyme are dominated by hydrogen bonded interactions between specific residues and also nonspecific electrostatic interactions between the protein surfaces.[27,30,53] In this subsection, we investigate the nonspecific electrostatic interactions by calculating the electrostatic potentials on solvent accessible surfaces of ACE2 and RBD for all the five systems (Figure ). The electrostatic potential on the surface of ACE2 protein is found to be highly negative (Figure a). The electrostatic potential on the H1 helix, which mainly interacts with the RBD of both SARS-CoV-1 and SARS-CoV-2, is also negative (Figure b). The electrostatic potentials on various parts of the RBD surfaces of the five systems are found to be mainly positive with variable charge distributions (Figure c–g). The positively charged surface regions of RBD favor its binding with the negatively charged ACE2. The electrostatic potentials on protein surfaces show that there are more patches of negative regions on the RBD surface for SARS-CoV-1 than SARS-CoV-2 (Figure c,d) which means a relatively weaker binding of RBD with ACE2 for SARS-CoV-1. For the N501Y mutant of SARS-CoV-2, an additional positive patch is observed on its RBD surface (Figure e) which favors its binding with ACE2. For the double mutant E484Q/L452R variant of SARS-CoV-2, the mutation increases the positive electrostatic potential patches on RBD as shown in Figure f. An increase in the regions of positive electrostatic potential is also observed for the N440K variant of SARS-CoV-2 (Figure g). Thus, the mutations in RBD can increase the regions of positive electrostatic potential on its surface and enhance its binding with ACE2 which primarily has negative potential on the receptor parts of its surface. We note that the current results of opposite electrostatic potentials of RBD and ACE2 calculated for their molecular dynamics equilibrated configurations are in overall agreement with similar observations made earlier for the protein structures taken from protein data bank,[80] Monte Carlo simulation[53] and homology models.[54,99]

Figure 4

Electrostatic potentials on the solvent accessible surfaces of ACE2 and RBD proteins obtained by using the APBS method.[92] The color code is varied from −1 (red) to +1 (blue).

Specific Residue-wise Interactions between RBD and ACE2

The residues of RBD at the ACE2-RBD interface are mostly conserved for SARS-CoV-1 and SARS-CoV-2 (Figure a). The RBD residues of SARS-CoV-2 that can interact with ACE2 are Arg403, Asn439, Val445, Gly446, Tyr453, Leu455, Phe456, Tyr473, Ala475, Gly476, Phe486, Asn487, Tyr489, Gln493, Tyr495, Gly496, Gln498, Thr500, Asn501, Gly502, Val503, and Tyr505 (Figure a). Similarly, the interacting residues of RBD of SARS-CoV-1 at the interface are Lys390, Arg426, Ser432, Thr433, Tyr440, Phe442, Leu443, Phe460, Pro462, Asp463, Phe472, Asn473, Tyr475, Asn479, Tyr481, Gly482, Tyr484, Thr486, Thr487, Gly488, Ile489, and Tyr491 (Figure a). It is found that both direct residue–residue and water bridged residue–water–residue hydrogen bonds are present between RBD and ACE2 (Figure b,c). We have calculated the numbers of both types of hydrogen bonds between the residues of RBD and ACE2 for SARS-CoV-1, SARS-CoV-2 and its variants along the simulation trajectories (Figures S6 and S7), and the corresponding distributions are shown in Figure a–d. The average number of direct hydrogen bonds is found to be 7 for SARS-CoV-1 and 8 for SARS-CoV-2. A number of water-bridged hydrogen bonds are also found to be present between the residues of RBD and ACE2 (Figures , parts b and d). The average number of water-bridged hydrogen bonds for SARS-CoV-1 is found to be 4, and that for SARS-CoV-2 is 5.

Figure 5

(a) Superimposition of RBD structures showing the interface residues for SARS-CoV-1 and SARS-CoV-2. (b) Snapshot at the end of the simulation trajectory to show direct hydrogen bonds between the residues of RBD and ACE2 for SARS-CoV-2. (c) Snapshot at the end of the simulation trajectory to show water-bridged hydrogen bonds between the residues of RBD and ACE2 for SARS-CoV-2.

Figure 6

Distributions of (a) direct interprotein and (b) water-bridged hydrogen bonds for SARS-CoV-1 and SARS-CoV-2. Distributions of (c) direct interprotein and (d) water-bridged hydrogen bonds for the three variants of SARS-CoV-2.

We also calculated the percentage occupancies of various residue-specific hydrogen bonds (Tables S1–S5) for both direct interprotein and water-bridged hydrogen bonds. The percentage occupancies of some of the direct interprotein hydrogen bonds are also shown in parts a and b of Figure . The results show that the percentage occupancies of residue-wise hydrogen bonds are higher for SARS-CoV-2 than SARS-CoV-1 for many of the hydrogen bonds. The percentage occupancies of water-bridged hydrogen bonds are also found to be higher for many of such hydrogen bonds for SARS-CoV-2 than SARS-CoV-1. It is found that not only the average number of hydrogen bonds between RBD and ACE2 is higher, the average percentage occupancies of the interprotein hydrogen bonds are also greater for SARS-CoV-2 than SARS-CoV-1. The results of Figure b show that the occupancies of RBD-ACE2 hydrogen bonds are generally comparable or higher for the mutant types than the wild type SARS-CoV-2. The interactions between various residues of RBD and ACE2 for the mutant variants of SARS-Cov-2 are discussed below in more details.

Figure 7

Average occupancies of direct interprotein hydrogen bonds for (a) SARS–CoV-I and SARS-CoV-2 and (b) three variants of SARS-CoV-2, namely N501Y, E484Q/L452R, and N440Y.

SARS-CoV-2 (N501Y)

In this variant of SARS-CoV-2, the RBD of the spike protein is modified by mutating Asn501 to Tyr501. In the wild type SARS-CoV-2, the nearest residues of Asn501 of RBD are Tyr42 and Lys353 of ACE2. For the wild-type SARS-Cov-2, the Asn501 residue does not make good interactions with residues Tyr42 and Lys353 of ACE2 (Figure a). In the mutated system, the average number of hydrogen bonds between RBD and ACE2 is found to be 8 which is the same as in wild type SARS-CoV-2 (Figure a). However, the average number of water-bridged hydrogen bonds is now 7 which is higher than that for the wild type SARS-CoV-2 (Figure b). For the wild-type protein, the hydrogen bond occupancy of the Asn501 and Lys353 residues is found to be only 0.85%. The Asn501 also does not make any effective interaction with Tyr42 of ACE2. In the mutated system where Asn501 is replaced by Tyr501, the hydrogen bond occupancy between Tyr501 and Lys353 increases to 4.52%. Besides, the Tyr501 residue of RBD also participates in a T-shaped π–π interaction with residue Tyr42 of ACE2[100] (Figure b). An analysis of the minimum distances between Asn501-Tyr42 for the wild type and Tyr501-Tyr42 for the mutated type of SARS-CoV-2 reveal that the average minimum distance is ∼6 Å for Asn501-Tyr42 for the wild type, and it is reduced to ∼3 Å for Tyr501-Tyr42 for the mutated system (Figure c).

Figure 8

(a) Interactions between Asn501 of RBD with Tyr42 and Lys353 of ACE2 for SARS-CoV-2. (b) Interactions between Tyr501 of RBD with Tyr42 and Lys353 of ACE2 for the N501Y variant of SARS-CoV-2. (c) Variation of minimum distances between Asn501 and Tyr42 for SARS-CoV-2, and between Tyr501 and Tyr42 for the N501Y variant of SARS-CoV-2 (N501Y).

SARS-CoV-2 (E484Q/L452R)

The average number of hydrogen bonds between RBD and ACE2 for this double mutant variant is 9 which is higher than that for the wild type SARS-CoV-2. The average number of water-bridged hydrogen bonds between RBD and ACE2 is found to be 8 which is also higher than that found for the wild type SARS-CoV-2. The Glu484 residue at RBD of the wild type weakly interacts with Lys32 of ACE2 with a hydrogen bond occupancy of 2.23%. When this residue is mutated with Gln484, the hydrogen bond occupancy between Gln484 and Lys32 becomes 0.2%; i.e., the hydrogen bonding interaction between these two sites is further weakened. The Leu452 residue of the wild type protein is not present at the interface and does not participate in direct interaction with any residue of ACE2. However, these double mutations increase the overall positive electrostatic potential on the RBD and allows it to interact more effectively with ACE2 which possesses overall negative electrostatic potential on its surface. This enhanced electrostatic interaction, in turn, increases the occupancies of some of the already existing hydrogen bonds leading to an increase in their average numbers. This can also be seen from Figure . Thus, for the double mutant variant of SARS-CoV-2, enhanced electrostatic interactions primarily contribute to the stronger binding of RBD with ACE2.

SARS-CoV-2 (N440K)

The Asn440 residue of RBD is generally located away from any residue of ACE2. The average number of hydrogen bonds between RBD and ACE2 for this variant is 7 which is lower than the wild type SARS-CoV-2. However, the average number of water-bridged hydrogen bonds between RBD and ACE2 is 6 which is slightly higher than the wild type SARS-CoV-2. When the Asn440 residue is mutated with Lysine, it also increases the positive electrostatic potential on the RBD due to the positively charged nature of the lysine residue. This increased positive electrostatic potential on the RBD surface of the N440K variant can also be seen from Figure g, and this enhanced positive potential gives rise to a stronger binding of the RBD-ACE2 complex.

Unbinding of RBD from ACE2: Free Energy Diagrams along the Dissociation Coordinate

The free energy diagrams of unbinding of RBD from ACE2 were calculated by using the umbrella sampling method for SARS-CoV-1, SARS-CoV-2, and the three variants of SARS-CoV-2 (N501Y), SARS-CoV-2 (E484Q/L452R) and SARS-CoV-2 (N440K). These calculations show how the free energy changes when RBD is slowly pulled apart from ACE2 so as to move the RBD-ACE2 complex from its initial bound to the final unbound or dissociated state in aqueous medium. We have considered the center-of-mass separation between the two proteins as the relevant collective variable or the reaction coordinate of the unbinding process. The free energy diagrams show up to which separation any effect of binding between the two proteins still exists and also provide a quantitative estimate of the free energy barrier for unbinding of RBD from ACE2 under their room temperature solvated conditions and how it changes with mutation.

The free energy diagrams are shown in parts a and b of Figure . The free energy difference between the bound and unbound states of the RBD-ACE2 complex is found to be 10.20 kcal/mol for SARS-CoV-1. The experimental free energy difference from the equilibrium constant K for SARS-CoV-1 is estimated to be 10.17 kcal/mol.[101] It is clear from the free energy diagrams that this free energy difference also acts as the free energy barrier for the unbinding process from initial bound to the final unbound state. The calculated free energy barrier for RBD-ACE2 unbinding for SARS-CoV-2 is found to be 11.75 kcal/mol which can be compared with the experimental free energy difference of 11.35 and 12.24 kcal/mol obtained from K values reported in ref (27) and ref (20), respectively, for SARS-COV-2. The above-mentioned calculated values are without standard state corrections. In the later part of this section and in the Supporting Information, we have discussed standard state corrections to the unbinding free energies. As can be seen, even after standard state corrections, the agreements with experiments are reasonably good for both SARS-CoV-1 and SARS-CoV-2. We also note in this context that the experimental results of free energy differences are deduced from measured equilibrium dissociation constants (K) of ACE2 and the full trimeric spike protein[20,27,101] The contributions of the non-RBD parts of the spike protein including that of the glycosylation[35] on parts other than the RBD are missing in the current free energy calculations since only the glycosylated RBD part, rather than the full trimeric spike protein, is considered in the current calculations to study its binding with ACE2. Also, as discussed in Section , full glycosylation of ACE2 is missing in the current calculations. We also note that the calculated increase in the free energy barrier by 1.55 kcal/mol for SARS-CoV-2 is in very good agreement with the corresponding increase of about 1.5 kcal/mol found through a combination of measurements of biomechanical forces and transition state theory for dissociation of RBD-ACE2 complexes for SARS-CoV-1 and SARS-CoV-2.[30]

Figure 9

Free energy diagrams for the unbinding of RBD from ACE2 for (a) SARS-CoV-1 and SARS-CoV-2 and (b) N501Y, E484Q/L452R, and N440K variants of SARS-CoV-2. The free energies were obtained through biased simulations using the umbrella sampling method. R is the separation between the COMs (center of mass) of the two proteins which has been used as the dissociation coordinate in the umbrella sampling calculations. The error estimates shown on the free energy profiles were calculated over 100 rounds of bootstrapping analysis.[98]

For the N501Y variant of SARS-CoV-2, it is found that this mutation significantly affects the binding of RBD with ACE2 (Figure b). The mutation of Asn501 with Tyrosine gives rise to two additional interactions between RBD and ACE2. The first one is a weak hydrogen bonding interaction of Tyr501 with Lys353 of ACE2 and the second one is its T-shaped π–π interaction with Tyr42 of ACE2. The estimated barrier of unbinding is found to be 15.65 kcal/mol which shows an increase of free energy barrier by around 4 kcal/mol for this mutated system. The higher barrier for the dissociation of RBD-ACE2 complex for this so-called Alpha variant is consistent with the results of a recent steered molecular dynamics study[69] which showed that this variant (N501Y) requires higher initial force to pull RBD from ACE2 than the other variants. The mutations of E484Q and L452R also lead to an increase in the free energy barrier for unbinding of the RBD-ACE2 complex. As discussed before, the mutation of Glu484 to Gln484 increases the overall positive potential on RBD, hence enhances its binding with ACE2. Although Arg452 in the mutated system does not directly interact with any residue of ACE2, the positively charged arginine increases the electrostatic interaction with ACE2. The estimated free energy difference between the bound and unbound states is 13.75 kcal/mol which is ∼2 kcal/mol higher than that for the wild type SARS-CoV-2. For the N440K variant of SARS-CoV-2, it is found that the mutated residue is not involved in any direct interaction with ACE2. However, as discussed before, the positive charge of the lysine residue enhances the electrostatic interaction of RBD with ACE2, which leads to an increase in the free energy barrier for unbinding by ∼1 kcal/mol compared to that for the wild type SARS-CoV-2. The standard state corrections to the unbinding free energies are also calculated for all the systems by following the method described in refs.[102,103] and the results are included in the Supporting Information (Table S6).

Successive Breaking of Hydrogen Bonds during Dissociation of RBD-ACE2 Complex

We have investigated the successive breaking of interprotein hydrogen bonds during the dissociation of RBD-ACE2 complex for all the systems studied here. In particular, we have looked at which interprotein hydrogen bonds are broken at which separation between the COMs (center of mass) of RBD and ACE2 and where such breakings appear on the free energy diagrams along the dissociation coordinate. The results of these calculations are shown in Figure . Further details of these results are included in Table S7 in the Supporting Information. The results show that, for SARS-CoV-1 (Figure a), most of the interprotein hydrogen bonds are broken by the COM separation of ∼50 Å, only two hydrogen bonds still survive after 55 Å, and the last hydrogen bond breaks when the reaction coordinate reaches the value of 60 Å. For SARS-CoV-2 (Figure b), however, only one hydrogen bond is broken by COM separation of 50 Å, five hydrogen bonds still survive after 55 Å, and as many as three hydrogen bonds finally break at COM separation of 62 Å. Thus, although the COM separation is around 47 Å for the stable equilibrium configuration of the RBD-ACE2 complex where maximum number of interprotein hydrogen bonds are present, some of these hydrogen bonds still exist even when the COM of RBD is pulled by 15 Å from ACE2 so as to reach a value of 62 Å for the reaction coordinate of the unbinding process. This happens because as the RBD-ACE2 complex dissociates along the reaction coordinate, some conformational changes of the two proteins also take place so as to retain their hydrogen bonds and maximize their binding to the extent possible during the unbinding process. Similar results are also found for the three variants of SARS-CoV-2 (Figures c–e). In Figure , we have shown configurations of the RBD-ACE2 complex during the dissociation process for SARS-CoV-2. The snapshots are taken from the umbrella sampling simulations for different values of the dissociation coordinate. Similar snapshots of the RBD-ACE2 complexes during their dissociation from umbrella sampling simulations are shown in Figures S8–S11 for SARS-CoV-1 and the three variants of SARS-CoV-2.

Figure 10

Successive breaking of interprotein hydrogen bonds during the dissociation of RBD-ACE2 complex for SARS-CoV-1, SARS-CoV-2, and its three variants considered in this study. R is the separation between COMs (cente of mass) of the two proteins.

Figure 11

Snapshots of RBD-ACE2 complex for SARS-CoV-2 taken from umbrella sampling simulations for different values of the reaction coordinate during its dissociation. Here, RBD and ACE2 are colored in green and gray, respectively. R is the separation between COMs (center of mass) of the two proteins.

Summary and Conclusions

We have investigated the interprotein interactions and free energy diagrams for the dissociation of ACE2-RBD complexes for SARS-CoV-2, SARS-CoV-1, and also three different mutant variants of SARS-CoV-2 in aqueous medium through all-atom simulations. Altogether, we have performed all-atom molecular dynamics simulations of five fully solvated RBD-ACE2 systems of average system size of ∼230 000 atoms and total run length of 6.7 μs. The main driving force for the formation of the RBD-ACE2 complex seems to come from nonspecific electrostatic interactions between the two proteins and also from specific hydrogen bonding interactions. The electrostatic potential on the ACE2 surface is negative whereas that on RBD is found to be mostly positive, which favors their binding. The double mutation in the spike protein of SARS-CoV-2 increases the overall charge of RBD by +2 which gives rise to a stronger binding with ACE2. For the N440K variant, although the mutant residue Lys440 does not directly interact with any residue of ACE2, it increases the positive electrostatic potential at RBD surface which, in turn, enhances the stability of RBD-ACE2 complex.

The average numbers of both direct and water-bridged hydrogen bonds are found to be higher for SARS-CoV-2 than SARS-CoV-1. For the N501Y variant of SARS-CoV-2, the average number of direct hydrogen bonds remains the same as for the wild type protein, but the number of water-bridged hydrogen bonds is found to be higher for the mutant variant. The Tyr501 in this variant makes a weak direct hydrogen bond with Lys353 of ACE2 and also participate in a T-shaped π–π interaction. These interactions combined with favorable electrostatic interactions between the two protein surfaces increases the binding strength of the RBD-ACE2 complex for this variant. For the double mutant variant of E484Q/L452R, the numbers of both direct hydrogen bonds and water-bridged hydrogen bonds are found to be higher, and they contribute to a significantly increased binding strength of the RBD-ACE2 complex for this variant of SARS-CoV-2. The average number of direct hydrogen bonds is found to be less for the N440K variant than the wild type, however there is an increase in the number of water-bridged hydrogen bonds for this variant. This increased number of water-bridged hydrogen bonds along with favorable electrostatic interactions between the two protein surfaces increase the binding strength of the RBD-ACE2 complex for this variant.

We have also calculated the free energy diagrams for the process of unbinding of RBD from ACE2 for all the five systems using the umbrella sampling method. It is found that the free energy barrier for unbinding is ∼1.5 kcal/mol higher for SARS-CoV-2 than SARS-CoV-1. The calculated binding free energies of −10.2 and −11.75 kcal/mol for SARS-CoV-1 and SARS-CoV-2, respectively, are found to be in good agreement with those obtained from experimental results.[27,101] The standard state corrections to the unbinding free energies are also calculated for all the systems studied here. We note that the free energy diagrams for unbinding of RBD-ACE2 complexes using umbrella sampling simulations at atomistic level are calculated here for the first time for SARS-Cov-2 and its three variants, namely N501Y, E484Q/L452R, and N440K. In all cases, the binding free energies for the mutated systems are found to be higher than the wild type SARS-CoV-2. The free energy barrier for the dissociation of RBD-ACE2 complex increases to ∼15.65 kcal/mol for the N501Y variant. As discussed before, additional interactions of Tyr501 of the mutant RBD with Lys353 of ACE2, participation in new T-shaped π–π interaction with Tyr42 of ACE2 and also electrostatic interactions between the two oppositely charged protein surfaces increase the free energy barrier for the unbinding of RBD from ACE2. For the double mutant (E484Q/L452R) variant, the main dominating interactions involve electrostatic and hydrogen bonding interactions with ACE2 which give rise to an increased free energy barrier of 13.75 kcal/mol for the dissociation of RBD-ACE2 complex. For the N440K variant also, an increased free energy barrier is found for the unbinding process due to enhanced electrostatic interactions between the RBD and ACE2 proteins. We have also looked at successive breaking of interprotein hydrogen bonds along the free energy pathway of the dissociation of the RBD-ACE2 complex for all five systems studied here. Specifically, we identified the hydrogen bonds, which are the last ones to break during the dissociation process, and also looked at conformational changes of the proteins that take place during the unbinding process. Thus, the present study provides detailed information on the molecular nature of RBD-ACE2 interactions and the molecular processes that are involved in the dissociation of those heterodimer protein complexes for SARS-CoV-1 and SARS-CoV-2 and also for three mutant variants of SARS-CoV-2.

We finally note that the current umbrella sampling calculations were based on a single reaction coordinate without employing any orthogonal restraints. Inclusion of orthogonal restraints help in better convergence of free energies with respect to positional and conformational sampling along the reaction pathway.[103−105] In the current calculations, simulations in each umbrella sampling window were run for 7 ns with the first 2 ns left out of free energy calculations. In Figure S12, we have shown free energy profiles obtained from 3 and 5 ns of production simulations for each umbrella sampling window, and the results are found to be very similar. Nevertheless, the convergence of free energies with respect to conformational sampling of RBD and ACE2 along their unbinding reaction pathway is an important issue which requires further study in future.