Binding Interactions between Receptor-Binding Domain of Spike Protein and Human Angiotensin Converting Enzyme-2 in Omicron Variant.

The emergence of new SARS-CoV-2 Omicron variant of concern (OV) has exacerbated the COVID-19 pandemic because of a large number of mutations in the spike protein, particularly in the receptor-binding domain (RBD), resulting in highly contagious and/or vaccine-resistant strains. Herein, we present a systematic analysis based on detailed molecular dynamics (MD) simulations in order to understand how the OV RBD mutations affect the ACE2 binding. We show that the OV RBD binds to ACE2 more efficiently and tightly predominantly because of strong electrostatic interactions, thereby promoting increased infectivity and transmissibility compared to other strains. Some of the OV RBD mutations are predicted to affect the antibody neutralization either through their role in the S-protein conformational changes, such as S371L, S373P, and S375F, or through changing its surface charge distribution, such as G339D, N440K, T478K, and E484A. Other mutations, such as K417N, G446S, and Y505H, decrease the ACE2 binding, whereas S447N, Q493R, G496S, Q498R, and N501Y tend to increase it.

Since the outbreak of the COVID-19 pandemic, its pathogen SARS-CoV-2 has continuously mutated and evolved, resulting in the emergence of major variants of concern (VOC). These VOC have been observed to alter the virus characteristics, such as infectivity, transmissibility, antigenicity, and pathogenicity.[1] The most recently identified SARS-CoV-2 VOC is the Omicron variant (OV) (B.1.1.529), which has quickly become the dominant strain.[2,3] It has the highest number of amino acid (AA) mutations of any known SARS-CoV-2 VOC, with over 30 mutations in the spike (S) protein, of which 15 are in the receptor-binding domain (RBD),[4] the main target for vaccine and treatment developments.[5−8] The presence of many mutations in OV S-protein has raised concerns about elevated transmissibility, immunological escape, and vaccine and treatment failures.[9−15] OV has been identified to contain several key mutations observed also in other SARS-CoV-2 VOC, such as K417N, E484A, and N501Y, that change the virus sensitivity to neutralization or increase the infectivity.[16] Moreover, it contains many novel mutations that have not been observed previously, and their biological effects are largely unknown. Because the binding between RBD and the angiotensin-converting enzyme-2 (ACE2) facilitates viral entry, initiating the infection process, the fundamental understanding of how the OV RBD interacts with ACE2 is pivotal for understanding the viral infection mechanism and its evolution, as well as for therapeutic development of effective means to reduce its spread.

The present study aims to investigate how OV mutations affect the binding strength between RBD and ACE2 by highlighting the role of each mutation, its underlying mechanism, and the pertinent binding driving forces. All-atom molecular dynamics (MD) simulations in explicit solvent have been implemented in order to study the dynamics and binding mechanism of the RBD-ACE2 system, followed by the molecular mechanics (MM) generalized Born surface area (MM-GBSA) method to predict the binding affinity and the binding profile. The results are compared with the previously reported analysis of the unmutated, wild type (WT) RBD-ACE2 system[17] in order to assess the effect of each mutation and the nature of its interaction based on the per-residue and pairwise decomposition schemes. MD simulations have proven valuable to investigate the dynamic and binding processes of the RBD-ACE2 complex in WT and many previous VOC.[18−25]

We implement a procedure, previously developed for both Alpha and Beta VOC, to build a computational model of the Omicron RBD-ACE2 system.[17] We briefly summarize it as follows. First, the interface structure of RBD-ACE2 complex (PBD ID:6M0J) is used as a template to create the OV RBD-ACE2 system, including all 15 RBD mutations:[26] G339D, S371L, S373P, S375F, K417N, N440K, G446S, S477N, T478K, E484A, Q493R, G496S, Q498R, N501Y, and Y505H.[4] The last 10 mutations are subsequently exposed to the ACE2 receptor as shown in Figure .

Figure 1

Omicron RBD-ACE2 interface system. (a) Front view of unbound Omicron variant (OV) RBD with all 15 mutations shown in blue and (b) top view. (c) Bound OV RBD-ACE2 model.

To substitute the AAs of WT with mutated AAs of OV, we used the Dunbrack backbone-dependent rotamer library[27] implemented by UCSF Chimera,[28] which also allows for a careful control of backbone dihedral angles. For example, the torsion angles of K417N and N501Y mutations from the available RBD Beta VOC have been used in our OV model,[29] as well as the angles of T478K mutation in the Delta variant.[30] Second, the resulting OV RBD-ACE2 complex is solvated with 27000 explicit water molecules together with the appropriate number of ions (1 Zn2+, 1 Cl–, and 22 Na+) using the TIP3P explicit water model, implemented in AMBER,[31,32] with the most recent ff14SB force field used for parametrizations of intermolecular and intramolecular interactions in the OV complex.[33] Third, the same MD steps of minimization, heating, equilibration, and two MD production runs were reproduced, but this time the MD runs were extended over 500 ns (1 μs in total).[17] For each run, 2500 frames are extracted from the entire simulation and used for the binding free energy (BFE) postprocess analysis. Finally, the MM-GBSA method[17,34−37] was applied to compute the BFE and to quantify the complete binding profile,[17] allowing for per-residue and pairwise BFE decompositions to identify the role and the nature of interaction for each mutation in the OV model.

Figure a shows the root-mean-square deviation (RMSD) and their frequency distributions for two replicate MD simulations of the OV. Overall, the complex achieves stable interfacial interactions as shown by small RMSD fluctuations. We find that the OV RBD-ACE2 complex has a relatively RMSD smaller than that previously reported of the WT from 100 ns MD simulations (average from both runs is 2.32 vs 2.53 Å of WT).[17] The 15 OV RBD mutations thus lead to a more stable interfacial complex with ACE2. Surprisingly, the inset figure depicts a higher overlap between two MD runs. Even though the simulation time is only 500 ns, this is a good sign that our simulations are reproducible and convergent. This is even more pronounced by the value of the root mean squared inner product (RMSIP) of 0.82 for the first two principal components (PC1 and PC2) (Figure S1 in the Supporting Information). This behavior is also evident from time evolutions and distribution of BFE (Figure S2).

Figure 2

Stability of the OV RBD-ACE2 complex and slight conformation changes of mutated residues at RBD core. (a) Time course and distribution of the root mean square deviation (RMSD) of the heavy atoms of OV RBD-ACE2 complex in both MD runs. (b) RMSF for 15 residues that are mutated in OV and compared with WT.

In comparison with WT, the OV mutations are seen to only slightly change the RBD structure. This is even more pronounced when their root-mean-square fluctuations (RMSFs) are compared, as shown in Figure S3. The average RMSF of mutated RBD in OV is 1.49 Å vs 1.29 Å of WT RBD, indicating that the OV RBD is relatively less rigid than WT RBD, thus allowing for larger fluctuations. Figure b compares the RMSF of only the mutated AAs in OV vs WT. In particular, D339, L371, P373, and F375 are relatively more flexible than their WT counterparts. This is consistent with a recent study based on the cryo-EM S-protein structure, showing that the L371, P373, and F375 significantly alter the conformation and mobility of RBD,[12] thus demonstrating that the computational MD simulations can confirm and reproduce what is observed experimentally.

The MM-GBSA method has been applied to calculate the BFE of the OV RBD-ACE2 system at 310 K (37 °C), with neutral pH (7.4) and 0.15 M uniunivalent NaCl salt concentration. Table presents the BFE and its decomposition for OV and in comparison with WT.[17] It shows that the OV RBD binds ACE2 more strongly than WT, with relative binding energy of −1.67 kcal/mol, consistent with recent experimental and computational studies.[11−13,38−40] Interestingly, our predicted BFE value of −14.53 kcal/mol for OV is close to the Alpha BFE value of −14.7 kcal/mol.[17] The complete thermodynamic decomposition listed in Table shows the Coulomb electrostatic interaction (ΔEele) of OV to be more than twice that of WT. The increase in ΔEele of OV is mainly the result of five AAs at the RBM changing from polar to positively charged residues (N440K, T478K, Q493R, Q498R, and Y505H). The electrostatic component is thus seen as the main reason behind the more effective binding of ACE2 and OV RBD, which may also elucidate the reason behind the highly contagious nature of OV (Figure S4). The electrostatic interaction has been shown to be the primary source of increasing the binding of SARS-CoV-2 to ACE2 compared to SARS-CoV as well as enhancing the binding affinity of VOC to ACE2.[23−25] ΔEele of OV, however, creates a higher desolvation energy (ΔGGB) that is indispensable in the formation process and cannot be avoided (Table ). On the other hand, the van der Waals interaction (ΔEvdW) plays a key role in stabilizing and governing the RBD–ACE2 interaction as well as enhancing their binding by gaining −2.71 kcal/mol as compared to WT.

Table 1

Decomposition of BFE (kcal·mol–1) of RBD-ACE2 Complex in OV and WT

energy	OV (SEM)	WT (SEM)a	ΔΔGb
ΔEvdW	–92.81 (0.1)	–90.10 (0.1)	–2.71
ΔEele	–1486.5 (0.65)	–700.92 (0.6)	–785.58
ΔEMM	–1579.31 (0.66)	–791.03 (0.7)	–788.28
ΔGGB	1534.19 (0.64)	748.49 (0.6)	785.7
ΔGSA	–13.27 (0.01)	–13.21 (0)	–0.06
ΔGsol	1520.92 (0.64)	735.28 (0.6)	785.64
–TΔS	43.86	42.89	0.97
ΔGbind	–14.53 (0.1)	–12.86 (0.1)	–1.67

All values are taken from ref (17), and SEM is the standard error of the mean.

ΔΔG = ΔGOV – ΔGWT.

To gain a better understanding of the nature and impact of each mutation on the BFE in RBD-ACE2, in terms of per-residue fractions, decomposition schemes have been implemented and are shown in Figure S5 for all OV and WT single AAs. Figure a shows the per-residue BFE decomposition change for the 15 mutations between OV and WT (ΔΔG = ΔGOV – ΔGWT), while their AA–AA interaction pair maps are displayed in Figures b. These mutations are divided into three groups based on their influence on binding: neutral, decreased, and increased binding.

Figure 3

BFE decompositions of RBD-ACE2 complex in OV and WT. (a) Change in per-residue BFE decomposition (ΔΔG = ΔGOV – ΔGWT). The 15 OV mutations are classified based on their ACE2 binding properties as neutral, increased, or decreased binding. (b) Pairwise BFE decomposition for only the mutated OV residues (blue characters on y-axis) that form pairs with ACE2, compared with their WT (black). (c) Details of the interactions for the five important mutations in OV (Q493R, G496S, Q493R, N501Y, and Y505H) and compared with their WT. The black dashed lines represent possible hydrogen bonds or salt-bridges, while the red dashed lines represent hydrophobic interactions.

Although mutations outside of RBM and away from the interface, such as G339D, S371L, S373P, and S375F, do not affect the RBD–ACE2 binding, they may still give rise to other biological consequences. For instance, substituting neutral G339 with highly negatively charged D339 changes the surface charge distribution of RBD (Figure ), which may impair the binding between RBD and antibody. Indeed, a recent study found that this mutation has a slightly higher escape fraction for Sotrovimab, a human neutralizing monoclonal antibody (mAb).[41] Aside from that, changing G339 to D339 results in a longer side chain, which probably varies the local intramolecular interactions, particularly with the N343 glycosylation site.[12] We did not include glycans in our simulation because earlier work suggested that the RBD of the S-protein had far less glycans than the S-protein itself and did not interact directly with ACE2.[42] Mutating the polar residue (S) to the hydrophobic residues, L at 371, P at 373, and F at 375, forms a unique cluster that changes the biochemical properties of this RBD region in ways not previously observed in any other strains. This again allows OV to escape from the class 4 antibodies and some other antibodies from class 1, 2, and 3.[11,12] Surprisingly, our MD simulations revealed that S371L, S373P, and S375F mutations are more flexible and induce a conformational change in RBD, suggesting that they have a higher chance to evade antibody recognition.

Figure 4

Comparison of electrostatic potential surface of RBD in mutated Omicron (top) and WT (bottom). From left to right, the first four shapes are of the RBD front view, rotated sequentially by 90°, while the fifth one represents the RBD top view. Important residues that changed their electrostatic potential surfaces are labeled. Negative-potential residues are shown in red, near-neutral residues in white, and positive-potential residues in blue.

Our findings based on BFE decompositions show that the E484A, T478K, and N440K mutations would exhibit the same pattern of neutral binding (Figure ). Interestingly, E484K mutation has also been observed in Beta and Gamma VOC and another variant of interest (VOI), and it has been identified as an immunodominant spike protein residue. Therefore, E484A in OV is expected to greatly reduce the susceptibility of many mAbs, which is fully consistent with the recent studies.[11,12] In our MD simulations, we observed that the E484A mutation eliminates the weak E484:K31 in WT, but it reduces the destabilization that is stemming from possible electrostatic repulsion between E484 of WT and E35 of ACE2 when changing to A484. Therefore, this mutation has no impact on BFE, unlike the E484K mutation in Beta, which increases the interaction marginally.[17] T478K has also been seen in the Delta variant and has no direct influence on the OV RBD-ACE2 interface network.[43] Finally, the N440K has been linked to an increase in antibody neutralization resistance for some antibodies.[11,44] Importantly, the exchange of amino acids at these sites (440, 478, and 484) contributes to the emergence of a different electrostatic surface of RBD, which may play some role in attractive electrostatic interactions with the negatively charged surface of ACE2 (Figure ). In this context, ab initio quantum methodologies offer a more accurate description of the partial charge distributions for these charged residues.[17,45−50]

Similar to the Beta variant, the K417N mutation reduces the OV RBD-ACE2 binding because of the loss of the strong ionic pair with D30 on ACE2.[17] While this result is consistent with previous observations,[29,51] we note that the K417N mutation could affect also the way the RDB clamps the ACE2.[19] Certain mAbs, such as Etesevimab and Casirivimab, have been shown to be affected by K417N.[52] Similarly, the G446S and Y505H reduce the binding, but we suspect that this is because of the impact of other mutations like N501Y rather than the mutations themselves. In our early study, we observed that the N501Y in the Alpha and Beta VOC significantly enhances the interactions with ACE2, but it eliminates the hydrogen bond (HB) between the G446 and Q42 of ACE2 and reduces the strength of the Y505:E37 pair (Figure c).[17] The same is true in OV, even though the residues at 446 and 505 site are now S and H, and as a result, their overall contributions to BFE in OV are smaller than in WT.

On the other hand, S477N, Q493R, G496S, Q498R, and N501Y mutations confer additional strength to the ACE2 binding, while N477 (OV) forms new HBs with S19 on ACE2. Interestingly, our findings indicate that the main source of increased binding is due to the formation of two new strong salt bridges between R493 and R498 of OV and E35 and D38 on ACE2 with pair interaction strengths of −11 and −6.5 kcal/mol, respectively, which was not observed in the WT (Figure b,c). According to our findings, the dominant mutations in the OV that enhance electrostatic interactions to ACE2 are Q493R and Q498R. This finding is in full agreement with experimental results.[53,54] R493 (OV) can also form a salt bridge with D38 with pair interaction strength of −4.5 kcal/mol, but it also makes a weak unfavorable pair with K353 and loses the pair with K31 that exists in the WT. R498 (OV) retains the pairings with Y41 and Q42, but they are stronger than Q498 (WT). However, it loses the pair with K353 because of the N501Y mutation.[17] Unlike G496 of WT, the aliphatic hydroxyl group of S496 (OV) also forms new HBs with D38 on ACE2 with pair interaction strength of −4.3 kcal/mol and maintains the HB with K353 (Figure b,c). Just like in Alpha and Beta VOC, Y501 forms more pairings than N501 of WT (Figure b,c).[17] The number of pairs between Y501 of OV and ACE2 residues is the same as in Alpha and Beta, but their interaction strengths change, particularly the Y501:D38 pair (−0.35 of OV vs −4 kcal/mol of Alpha or Beta[17]).

In conclusion, comprehensive MD simulations have been performed to investigate the effects of OV RBD mutations on ACE2 human cell receptor binding. Our results shed light on the critical roles of the new OV mutations in the conformational changes of RBD (S371L, S373P, and S375F) and changes in its electrostatic potential surface (N440K, T478K, E484A, Q493R, Q498R, and Y505H). As a result, the overall electrostatic attraction between the positively charged OV RBD and the negatively charged ACE2 receptor is twice as strong as in the WT, leading to increased OV contagiousness as well as to the escape from neutralizing antibodies. Our analysis also shows that the loss of some interactions caused by K417N, G446S, and Y505H is completely compensated by the formation of new pairs from S477N, Q493R, G496S, Q498R, and N501Y, resulting in an overall stronger binding of OV’s RBD-ACE2.