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

Severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) is the causative agent of the current coronavirus disease 2019 (COVID-19) pandemic. It is known that the receptor-binding domain (RBD) of the spike protein of SARS-CoV-2 interacts with the human angiotensin-converting enzyme 2 (ACE2) receptor, initiating the entry of SARS-CoV-2. Since its emergence, a number of SARS-CoV-2 variants have been reported, and the variants that show high infectivity are classified as variants of concern according to the United States Centers for Disease Control and Prevention. In this study, we performed both all-atom steered molecular dynamics (SMD) simulations and microscale thermophoresis (MST) experiments to characterize the binding interactions between ACE2 and RBD of all current variants of concern (Alpha, Beta, Gamma, and Delta) and two variants of interest (Epsilon and Kappa). We report that RBD of the Alpha (N501Y) variant requires the highest amount of force initially to be detached from ACE2 due to the N501Y mutation in addition to the role of N90-glycan, followed by Beta/Gamma (K417N/T, E484 K, and N501Y) or Delta (L452R and T478 K) variants. Among all variants investigated in this work, RBD of the Epsilon (L452R) variant is relatively easily detached from ACE2. Our results from both SMD simulations and MST experiments indicate what makes each variant more contagious in terms of RBD and ACE2 interactions. This study could shed light on developing new drugs to inhibit SARS-CoV-2 entry effectively.

Introduction

Reported in late 2019, severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) emerged and has rapidly infected people worldwide. As of mid-September 2021, 230 million cases and 4.71 million deaths have been reported globally.[1] Despite worldwide efforts to overcome the current coronavirus disease 2019 (COVID-19) pandemic, the rise of various SARS-CoV-2 variants may deteriorate the efficacy of vaccination and other countermeasures.

The SARS-CoV-2 virus utilizes receptor-binding domain (RBD) of the S1 protein, a part of trimeric spike (S) glycoprotein,[2,3] for viral entry through the RBD interaction with the human receptor angiotensin-converting enzyme 2 (ACE2). Since ACE2 can interact with RBD of both SARS-CoV-2 and SARS-CoV (or SARS-CoV-1, the virus that caused the 2002–2004 SARS outbreak), there have been many studies not only to understand binding interactions between RBD and ACE2 but also to characterize the difference between SARS-CoV-1 and SARS-CoV-2.[4−6]

In September 2020, the Alpha variant, lineage B.1.1.7, was first detected in southeast England and quickly became a populated lineage in the United Kingdom. The variant was subsequently detected in the United States in December 2020.[7,8] The Beta variant, lineage B.1.351, was first detected in South Africa in May 2020 and found in the United States at the end of January 2021.[9] At that time, there was another identified Gamma variant, which is known for lineage P.1,[10,11] in the United States that was initially found in Japan from a traveler from Brazil. In November 2020, the Epsilon variant, lineage B.1.427, was detected in California in the United States.[12] Recently, two additional variants, Kappa (lineage B.1.617.1) and Delta (lineage B.1.617.2), first identified in India at the end of 2020, were detected in the United States.[13] Since the emergence of diverse SARS-CoV-2 variants, Alpha, Beta, Gamma, and Delta variants have been classified as variants of concern by the United States Centers for Disease Control and Prevention (CDC) due to their high infectivity.

Several studies have been performed experimentally and computationally to better understand the highly contagious characteristics of these variants.[14−16] For example, Tian et al. conducted an experimental and computational study to capture the role of the N501Y mutation in Alpha, Beta, and Gamma variants.[14] They suggested that the π–π interactions and π–cation interactions are responsible for the enhanced interactions between RBD and ACE2. However, only the N501Y mutation was examined in their study, although other potentially important mutations have emerged. More recently, Socher et al. performed energy decomposition analysis from molecular dynamics simulations to compare the interaction energies between ACE2 and RBDs of Alpha, Beta, and Gamma variants.[15] They investigated each specific mutation, N501Y, K417N/T, and E484 K, and reported that F486, Q498, T500, and Y505 in RBDs are important residues across viral variants in the RBD–ACE2 interface.

In this study, using all-atom steered molecular dynamics (SMD) simulations and microscale thermophoresis (MST) experiments (see Methods, Supporting Information), we report the differential interactions between human ACE2 and RBD of SARS-CoV-2 of all variants of concern (Alpha, Beta, Gamma, and Delta) as well as two variants of interest (Epsilon and Kappa). The study also provides a better understanding of such differences at the molecular level.

Methods

Computational Methods

A fully-glycosylated SARS-CoV-2 RBD and ACE2 complex was obtained from the COVID-19 Protein Library in the CHARMM-GUI Archive (6vsb_1_1_1_6vw1.pdb).[17] The complex includes 6 N-linked glycans: five glycans in ACE2 (Asn53, Asn90, Asn103, Asn322, and Asn546) and one glycan in RBD (Asn343). For system generation, parameter setup, and corresponding mutations, we utilized CHARMM-GUI Solution Builder.[18,19] From the WT RBD structure, each variant was modeled with the following mutations: Alpha (N501Y), Beta (K417N, E484K, N501Y), Gamma (K417T, E484K, N501Y), Epsilon (L452R), Kappa (L452R, E484Q), and Delta (L452R, T478K). The CHARMM36(m) force field[20,21] for protein and carbohydrates with TIP3P water model[22] was used with 0.15 M of K+ and Cl– ions for mimicking physiological conditions. The system size was determined to be large enough (about 190 Å × 190 Å × 190 Å) to have the proteins solvated enough when they are fully detached. The total number of atoms is approximately 550,000.

The overall simulation details are nearly identical to our previous work.[4] NAMD simulation software[23] was used for the pulling simulations with the COLVARS method. As an initial condition, the SARS-CoV-2 RBD and ACE2 complex structures were aligned along the X-axis, and the center of mass (COM) of each protein was calculated to apply the external force on the proteins. The effective force acting on the COMs of both proteins can be calculated through the following equationwhere k is the spring constant, v the moving speed of the spring potentials (also called dummy atoms), R(t) the current position of the selected protein COM, and the COM-COM unit vector. This force enables the spring-connected proteins to move in the opposite directions to pull away two proteins. The moving speed of proteins was set to 0.5 Å/ns along the X-axis, and a spring constant of 5 kcal/mol/Å2 was applied to the COM of each protein to have both proteins move along the X direction and restrict moving along the Y and Z directions. For better statistics, 20 independent simulations for each system were performed (140 systems total, 20 replicas of 7 variants) with at least 40 ns of each simulation run. The pulling simulations stopped when RBD and ACE2 were completely detached from each other.

The van der Waals interactions were switched off smoothly over 10–12 Å using a force-based switching function.[24] The electrostatic interactions were calculated by the particle-mesh Ewald method with a mesh size of 1 Å.[25] To constrain bond lengths involving hydrogen atoms, the SHAKE algorithm was used.[26] The simulation time step was set to 4 fs with the hydrogen mass repartitioning method.[27,28] Equilibration simulations were performed with the NVT (constant particle number, volume, and temperature) ensemble where positional and dihedral restraints were employed. The restraint was gradually decreased during the equilibration simulations. The NPT (constant particle number, pressure, and temperature) ensemble was then applied for the production runs, where the Langevin piston method[29] was used for the pressure control. The simulation temperature was set to 303.15 K with the Langevin damping control method.

Experimental Methods

The recombinant human ACE2 protein (GenBank accession: AF291820.1, Sino Biological 10108-H08H; Wayne, PA) was labeled with RED-NHS (second Generation) dye using the Monolith Protein Labeling Kit (NanoTemper Technologies, MO-L011, München, Germany). Labeled ACE2 (5 nM, final concentration) was mixed with the RBD proteins (WT or variants, 2-fold diluted in a 15-step starting from 1.5 to 4 μM) in a PBS buffer supplanted with 0.1% Pluronic F-127. The RBD proteins include WT (ACRObiosystems, SPD-C52H3, Newark, DE, GenBank accession: QHD43416.1), Alpha (ACRObiosystems, SPD-C52Hn), Beta (ACRObiosystems, SPD-C52Hp), Epsilon (Sino Biological, 40592-V08H28), Kappa (Sino Biological, 40592-V08H88), and Delta (Sino Biological, 40592-V08H90). All the recombinant proteins used in this study were produced in HEK293 cells and presumably fully glycosylated. The mixed RBD + ACE2 samples were separately loaded into 16 premium glass capillaries (NanoTemper Technologies, MO-K025). The 16 capillaries were then placed in the reaction chamber in the order of concentration. MST measurements were conducted on a Monolith NT.115 instrument (NanoTemper Technologies) at 20% excitation power at 24 °C. The measurement was repeated at least three times. Kd calculations were performed using the MO Affinity Analysis software (NanoTemper Technologies).

Results and Discussion

Almost All Variants Show Increased Interactions with ACE2

To gain molecular insight into the difference of all variants that are classified as variants of concern (Alpha, Beta, Gamma, and Delta) and two additional variants of interest (Epsilon and Kappa), pulling force analysis was performed on each RBD–ACE2 complex (Figure A) as a function of distance (D) between the COMs of RBD and ACE2 proteins. Our fully-glycosylated S RBD–ACE2 complex model (Figure B, C) was employed for the pulling simulation.[17] As shown in Figure A, most variants have increased force profiles than WT except for the Epsilon variant, indicating that the variants have strengthened interactions with ACE2. It should be noted that the amount of average force at D = 53 Å shows a good match with our previous WT study,[4] where we utilized only the N-linked glycan (N-glycan) core structure for all N-glycans. In this study, we used the most probable N-glycan structures (Figure C) that are larger than the core structure.

Figure 1

(A) Average force profiles of WT (red), Alpha (blue), Beta (orange), Gamma (sky blue), Epsilon (green), Kappa (pink), and Delta (gray) variants as a function of the distance between the centers of mass of RBD and ACE2. (B) Initial snapshot of WT. Residues subjected to each mutation are shown as solid sticks (N501, K417, E484, L452, and T478). RBD and ACE2 are, respectively, colored in light gray and yellow. All N-glycans, water, and ions are hidden for clarity. (C) Initial snapshot of WT with clockwise 90° rotation along the normal from (B). All N-glycans are depicted in different colors. Any other residues, water, and ions are not shown for clarity.

The Alpha Variant Could Have the Highest Chance of RBD–ACE2 Interaction

Figure A shows that at D = 53 Å, the Alpha variant clearly requires the highest initial force to pull the RBD–ACE2 complex in the opposite direction. The difference can be explained in Figure B, a two-dimensional contact map between RBDAlpha and ACE2 at D = 53 Å, where RBD Y501 presents increased interactions with ACE2 Q42, Y41, and D38. Such contacts are decreased or even lost in the case of RBDWT or RBDEpsilon lacking the N501Y mutation (Figure A, D). To quantify the contact frequency between RBD residue 501 (N501 for WT, Epsilon, Kappa, and Delta; Y501 for Alpha, Beta, and Gamma) and ACE2, the number of heavy atom contacts was calculated (Figure A). The contact was counted if RBD residue 501 positioned within 4.5 Å of heavy atoms of key interacting residues of ACE2 protein. Notably, Y501 of Alpha, Beta, and Gamma variants retain more contacts (about 40%) than N501 of WT, Epsilon, Kappa, and Delta variants. As shown in Figure B and C, Alpha Y501 is located closer to ACE2 Y41 and K353 than WT N501 at D = 53 Å, and thus, it has the π–π and π–cation interactions with neighboring Y41 and K353, which is in accordance with the recent cyro-EM study.[30] On top of the Y501–ACE2 interactions, RBDAlpha also contains the highest amounts of contacts with ACE2 N90-glycan (Figure S3), which could be the reason why it has been reported as the most common lineages by June 19, 2021, among the estimated proportions of SARS-CoV-2 lineages according to the CDC,[31] although this study considers only single RBD out of the trimeric SARS-CoV-2 S protein.

Figure 2

Two-dimensional contact maps at D = 53 Å. (A) Interacting residue pairs between RBDWT and ACE2. RBD residues subjected to mutation are shown in colored boxes at the bottom: (B) blue for Alpha, (C) orange for Beta, and (D) green for Epsilon. The contact frequency is numbered with colors from light blue to dark blue. Dark red and yellow colors on the map, respectively, represent increased and decreased interactions between RBD and ACE2 upon mutations.

Figure 3

(A) The average number of contacts between RBD residue 501 and ACE2. (B, C) Representative snapshots at D = 53 Å of (B) Alpha variant and (C) WT. (D) Average number of contacts between RBD residue 417 and ACE2 and (E, F) their interacting residue pairs at D = 53 Å of (E) Beta and (F) Alpha variants. (G) Average number of contacts between RBD residue 478 and ACE2 and (H, I) key interaction pairs at D = 78 Å of (H) Delta and (I) Epsilon variants. The overall color scheme is the same as in Figure , and each mutated residue in each variant is shown using the same colors (i.e., red for WT, blue for Alpha, orange for Beta, green for Epsilon, and gray for Delta). Interacting residues are depicted as solid sticks, and residues losing their interactions are shown as transparent sticks. RBD and ACE2 are presented in light gray and yellow, respectively.

Beta and Gamma Variants Display Weaker RBD–ACE2 Interactions than the Alpha Variant, Yet They Are Still Stronger than WT

The force profiles of Beta and Gamma variants at D = 53 Å present weaker maximum forces than the Alpha variant, albeit they show higher forces than WT at the same distance (Figure A). As shown in Figure B and C, Alpha and Beta variants include the N501Y mutation, while the Beta variant involves two additional mutations, K417N and E484 K. Clearly, compared to WT or Epsilon, Y501 of the Beta variant has increased interactions (colored in dark red box) with ACE2 D38, Y41, and Q42, similar to the Alpha variant. However, it entails decreased contact frequency (shown as the yellow box) between RBDBeta N417 and ACE2 D30/H34, as well as RBDBeta K484 and ACE2 K31, which could explain why Beta has relatively weaker interactions than Alpha. The Gamma variant also shows decreased contact numbers similar to Beta due to its K417T mutation (Figure S2A). The only difference between Gamma and Beta is the K417 mutation, i.e., K417T vs K417N. Figure D compares the number of contacts of residue 417 of all variants that are in contact with heavy atoms of key interacting residues of ACE2. While all other variants containing K417 (i.e., WT, Alpha, Epsilon, Kappa, and Delta) display some RBD–ACE2 contacts from 50 to 60 Å, few interactions were found for the Beta variant. The side chain-shortening mutation from lysine to asparagine could have an impact on the RBD–ACE2 interface, resulting in fewer interactions at the same distance (Figure E, F). Interestingly, T417 of Gamma shows almost no interaction because threonine is even shorter than N417 of Beta. The weakened interactions of RBDBeta N417 and RBDGamma T417 could make them less contagious than the Alpha variant, while the N501Y mutation still allows them to have a strong enough potential to interact with ACE2. Recently, Gobeil et al.[32] and Bhattarai et al.[33] observed similar results from Alpha and Beta variants, and Barton et al.[34] reported corresponding results from Alpha, Beta, and Gamma through different approaches, indicating that our SMD results are reliable enough to investigate the RBD–ACE2 interface of the variants. The weakened interactions of both RBD residues 417 and 484 with ACE2 possibly provide an ability for the virus to escape from neutralizing antibodies targeting the RBD–ACE2 interface, allowing them to transmit more. This could explain why/how the Gamma variant took the second-highest portion by June 5, 2021, among the estimated proportions of SARS-CoV-2 lineages, provided by the CDC.[31]

L452 Mutation of Epsilon Variant Destabilizes RBD Itself, Causing Weakened Interactions with ACE2

Although most variants show similar maximum forces around D = 53 Å, the Epsilon variant shows decreased forces with more fluctuations than other variants (Figure A). The two-dimensional contact map in Figure D confirms its distinct interactions at D = 53 Å, as it shows the least number of contacts between RBDEpsilon and ACE2 (the yellow box represents deceased interactions). For example, K353 residue of all other variants actively interacts with ACE2 Q493, Q496, Q498, T500, N/Y501, G502, and Y505 (Figure A–C and Figure S2A–C). K353 of Epsilon, however, lost its contact with corresponding residues by at least 50%. To investigate the mechanism behind such a big difference, the contact analysis in between RBD residues was performed, where the influence of the L452R mutation was examined by checking its contacts with surrounding residues, L450 and L492 (Figure S4). Interestingly, mutated R452 interacts more with L450 (Figure S4C) and less with L492 (Figure S4A) simultaneously. Note that L450 and L492 are positioned in different β-strands (Figure S4B, D colored in green and orange, respectively), and the L452R mutation makes the RBD–ACE2 interface unstable by shortening each β-strand (i.e., the length of interacting β-strands of Epsilon variant is decreased by almost half). Because of such an unstable RBD structure, the Epsilon variant appears to be detached from ACE2 easier than WT. Indeed, K353 of the Epsilon variant lost contact with ACE2 Q498 and Y505 at D = 55 Å (Figure S4D), but WT holds their interactions at the same distance (Figure S4B). The Epsilon variant has been shown to reduce the neutralization potency of several antibodies in a way that it reduces sensitivity to the antibodies.[35] This result indicates that the Epsilon variant has impacted the world by decreasing the antibody sensitivity, not increasing direct RBD–ACE2 interactions. According to the CDC, as of June 29, 2021, the Epsilon variant deescalated from variants of concern and became a variant of interest since its considerable decrease in lineage proportion in the United States.

Position of T478 Mutation in Delta Variant Could Be Responsible for Its High Infectivity

Newly reported Kappa and Delta variants involve the same L452R mutation as Epsilon, but each variant contains an additional mutation: E484Q (Kappa) or T478 K (Delta). Two-dimensional contact maps (Figure S2B, C) display that Kappa and Delta variants have almost identical interaction patterns to WT between ACE2 K353 and RBD residues (i.e., Q493, Q496, Q498, T500, N/Y501, G502, and Y505). The Delta variant, interestingly, shows distinct features that are not found in other variants. Upon the T478 K mutation, it requires the highest force for the RBD–ACE2 complex to be completely dissociated at D = 78 Å (Figure A). To see what makes the difference, the numbers of contacts between RBD and ACE2 were calculated (i.e., residue 478 and heavy atoms of selected key interacting residues of ACE2). As shown in Figure G, RBDDelta exclusively makes more contacts with ACE2 than other variants. Figure H shows that Delta K478 retains contacts with ACE2 P84 and M82 at D = 78 Å, but Epsilon T478 already lost such interactions. The contacts of residue 478 are observed from the pulling simulations, but in terms of virus entry, it is possible that residue 478 located in the flexible loop first has a chance to contact ACE2. The stronger interactions of Delta K478 with ACE2 could explain why the proportion of Delta variant has dramatically increased with high infectivity. Recently, Baral et al. reported that a subtle reorientation of G496 in Delta induces stronger β-strand interactions and that it could be due to the L452R mutation.[36] It should be noted that both Kappa/Delta and Epsilon share the L452R mutation. Although our results are accordant with recently published studies,[30,32,33,35,36] the reason why Kappa/Delta and Epsilon behaviors are distinctive remained to be further studied, and it might stem from the limitation in our model, as we only employed the L452R mutation in RBD for the Epsilon variant without a D614G mutation. At the time, the Delta variant became the current variant of concern, and it took the highest portion among the estimated variant proportions as of July 3, 2021, according to the CDC.[31]

Microscale Thermophoresis Study Confirms Findings from MD Simulation

To validate the SMD simulation results, we conducted an experimental protein binding assay using MST. MST detects molecular binding kinetics based on the thermophoretic movement of molecules induced by a microscopic temperature gradient inside a glass capillary generated by an infrared laser.[37] MST has been used for detecting viral protein–receptor interactions,[38] including SARS-CoV-2 S proteins.[39] In our assay, human recombinant ACE2 was fluorescently labeled, and various RBD variants were titrated in a 2-fold fashion and mixed with the ACE2. The MST signal was first converted to saturated fraction data and subsequently fitted to a first-order 1:1 binding kinetics model using the manufacturer’s software (Figure S5). The binding affinities of ACE2 and RBDWT were detected to be 27.5 ± 4.8 nM (Figure ). This value is in agreement with a reported Kd range of 5–40 nM measured by surface plasmon resonance.[40] Importantly, our MST data indicate that the Alpha variant binds ACE2 with a 2.3-fold higher affinity (11.8 ± 0.8 nM) than WT. The rest of the variants show slightly different affinities from WT. Beta and Delta variants display approximately 20%–30% higher affinities than WT, and the Epsilon variant shows a 15% lower affinity than WT. In Figure , Kd values from MST experiments were directly compared with the FWT/F ratio from the SMD simulations, where FWT and F are the maximum pulling forces of WT and each variant around D = 53 Å (Figure A). Our MST affinity data are consistent with the SMD simulation data, indicating Alpha and Epsilon variants possess the strongest and weakest binding to ACE2, respectively.

Figure 4

Binding affinities between RBD variants and ACE2 and its comparison with the simulation results. Kd is obtained from microscale thermophoresis experiments. FWT/F is a ratio, where FWT and F are the respective maximum pulling force of WT and of each variant obtained from the SMD simulations.

Conclusions

This study characterizes differential interactions between ACE2 and RBD of all variants that the United States Centers for Disease Control and Prevention classifies as variants of concern and variants of interest. The results indicate that the Alpha variant requires the highest force for initial separation from ACE2, followed by Beta and Gamma variants or the Delta variant. K417N/T mutations of Beta and Gamma appear to make the RBD–ACE2 interactions less strong compared to the Alpha variant. In addition, the Epsilon variant is relatively more easily dissociated from ACE2 than others due to its destabilized RBD structure upon the L452R mutation. The Delta variant specifically shows stronger interactions with ACE2 than other variants at a relatively far distance between RBD and ACE2. The MST experiments show consistent results with the simulation results, where Alpha and Epsilon variants display the strongest and weakest binding to ACE2, respectively.

SARS-CoV-2 variants have been evolving by changing their structures so that they can either strengthen the interactions with the human receptor, i.e., ACE2, or escape from neutralizing antibodies by altering their structures targeting the RBD–ACE2 interface, highlighting their complex behaviors. We hope this study provides valuable information that distinguishes important features of all variants and their interactions with ACE2 and sheds light on developing new drugs to inhibit SARS-CoV-2 entry effectively.