The effect of the multiple mutations in Omicron RBD on its binding to human ACE2 receptor and immune evasion: an investigation of molecular dynamics simulations

SARS-coronavirus-2 (SARS-CoV2) Omicron variant (B.1.1.529) is of great concern to the world due to multiple mutations that may have an impact on transmissibility and immune evasion. Compared to the wild type (WT), there are 15 mutations in the Omicron receptor-binding domain (RBD), 10 of which are in the receptor-binding motif (RBM), where the host angiotensin-converting enzyme 2 (ACE2) interacts directly with. As a comparison, the currently dominant variant Delta (B.1.617.2) only has 2 mutations (L452R and T478K) or an additional E484K mutation in the RBM. As many as 15 mutations in Omicron RBD make it very hard to predict whether the mutations would increase the binding affinity to ACE2, particularly considering that 10 mutations crowded in the RBM. To understand the combinatorial mutation effect on Omicron RBD binding to ACE2 and potential immune evasion, we calculated the binding affinities of the WT/Delta/Omicron RBDs to ACE2 and antibodies with 600 ns molecular dynamics simulations for each system. We found that Omicron RBD has slightly weaker ACE2 affinities than WT RBD (-29.39 ± 2.96 Kcal/mol vs. -33.13 ± 3.26 Kcal/mol), and much lower affinities than Delta RBD (42.76 ± 2.38 Kcal/mol). Further analysis revealed that Omicron N501Y increase ACE2 binding but Q493K and Q498R decrease ACE2 binding. In addition, Omicron RBD might escape the launched monoclonal antibodies (mAbs) Etesevimab and clinical BD368-2 but may still sensitive to the launched mAbs Bebtelovimab.

The current global epidemiology of SARS-CoV-2 is characterized by a predominance of the Delta variant. Of 839 119 sequences uploaded to GISAID database with specimens collected in the last two months, more than 99.8% are Delta variant [2]. Liu et al. found that the binding affinity between R452/K478-RBD of Delta (B.1.617) variant and ACE2 is ~2 times higher than that of wild type (WT) RBD and ACE2 (8.3 nM) [3]. Therefore, the viral mutation may have high risk of increasing transmissibility and virulence or decreasing of effectiveness of therapeutics, making it more difficult to control the epidemic situation effectively.
SARS-CoV-2 variant Omicron (B.1.1.529) was newly announced by the World Health Organization (WHO) on November 26, 2021. Astonishingly, this variant carries as many as 30 single point mutations, 3 deletion mutation and one insertion mutation on spike protein, the main antigenic target of many monoclonal antibodies [4]. Among the mutations, 15 of them are on receptor-binding domain (RBD), especially 10 mutations are in the receptor-binding motif (RBM), which can interact directly with the angiotensin-converting enzyme 2 (ACE2) entry receptor on host membranes when the virus infects host cells [5]. In comparison, other detrimental variants have less than 5 RBD mutations [6]. For example, Delta (B.1.617.2) only has 2 mutations (L452R and T478K) or an additional E484K mutation in the RBM. Therefore, it is speculated that Omicron variant may significantly impact the effectiveness of current prophylactic and/or therapeutic drugs and binding affinity to ACE2. Consequently, Omicron mutant has aroused wide concern and many countries have taken measures on entry restrictions to prevent its rapid spread.
Bloom et al. systematically changed every single amino acid in the RBD of the spike protein and determine the effects of the substitutions on ACE2 binding [7]. By checking the 15 Omicron RBD mutations, it was found that 9 single mutations (S371L, S373P, S375F, K417N, G446S, E484A, G496S, Q498R, Y505H) should decrease the binding affinity to ACE2 while 6 mutations (G339D, N440K, S477N, T478K, Q493K, N501Y) might increase the binding affinity to ACE2. However, different from Alpha, Beta, Gamma and Delta variants, Omicron have much more mutations, making the predictions for ACE2 binding and immune evasion very difficult due to the combinatorial mutations, particularly considering that 10 mutations crowded in the RBM interacting with each other.
To understand the combinatorial mutation effect on Omicron RBD binding to ACE2 and potential immune evasion, we calculated the binding affinities of the WT/Delta/Omicron RBDs to ACE2 and antibodies by 600 ns molecular dynamics (MD) simulations for each system to shed light on the detailed implications of SARS-CoV2 Delta and Omicron variant on COVID-19 epidemiology, severity, effectiveness of public health and social measures, or other relevant characteristics. We found that Omicron RBD has slightly weaker ACE2 affinities than WT RBD (-29.39 ± 2.96 Kcal/mol vs. -33.13 ± 3.26 Kcal/mol), and much lower affinities than Delta RBD (-42.76 ± 2.38 Kcal/mol). Further analysis revealed that Omicron N501Y increase ACE2 binding but Q493K and Q498R decrease ACE2 binding. In addition, Omicron RBD might escape the launched monoclonal antibodies (mAbs) Etesevimab and clinical BD-368-2 but may still sensitive to the launched mAbs Bebtelovimab.

Preparation of mAb-S protein and ACE2-S protein complexes
The Omicron spike trimer was modelled by SWISS-MODEL Server in Alignment mode [8]. The Omicron homology model with the RBD up was chosen for further analysis. The Omicron spike were superimposed to a WT spike/ACE2 complex (PDB ID: 6vw1 [9]) to create an Omicron-ACE2 complex. We retrieved 3 marketed or clinical RBD-specific antibodies bound to S protein from the Protein Data Bank. The Omicron RBD domain containing residue 334-526 were truncated from the full-length S protein.
In order to get the intact structures for WT/Delta RBD and antibodies, missing residues in flexible loops were modeled using SWISS-MODEL. Delta RBD model was created by PyMOL2.5 [10] to yield K417N and E484K on the basis of 7vvs [11]. The Delta model in our simulations have 4 mutations: K417N, L452R, T478K, and E484K.

System preparation
Protonation states were assessed using H++ 3.2 [12,13] at pH 7.4. A cubic explicit water box described using the TIP3P model was used to solvated the complex system, which was extended by 10 Å from the solute. An atmosphere of 150 mM NaCl was also included in all simulations. The generated models were parametrized using amber ff1+4SB force fields [14] for protein. Subsequently, the parameter files created by tleap in Amber18 [15] were converted to gromacs format. About, 5000 steps of minimization including 2500 steps of steepest descent minimization and 2500 steps of conjugate gradient minimization were performed to remove bad contacts during the energy minimization phase. Equilibration in NPT ensemble was run at 1.0 bar and 300 K for 500 000 steps at 2 fs/step. Gromacs2020.2 [16,17] software package was used to run the minimization, equilibration simulations with position constraints (1 kcal/mol/Å 2 ) on protein.

Molecular dynamics (MD) simulations
Mdrun module in Gromacs2020.2 was used to perform 200ns MD production simulations at 300 K, 1 bar for all complexes. Temperature and pressure were controlled by Langevin thermostat [18] and a Nosé-Hoover Langevin barostat [19,20]. Bonds involving hydrogen atoms were fixed by the SHAKE algorithm [21]. The cutoff distance applied for van der Waals interactions was 10 Å. All simulations were performed using particle-mesh Ewald (PME) for long-range electrostatic interactions [22].
Mdconvert [23] was used to convert the trajectories to amber format. Cpptraj module in Amber18 was used for trajectory processing and analysis.

Binding free energy calculation
Binding free energy (ΔG) of RBD-antibody or ACE2-RBD complexes was calculated by MM/GBSA [24] method using GB OBC model(igb = 5)with a salt concentration of 150 mM. 750 snapshots extracted evenly from 50-200ns trajectories were used for binding free energy calculation. In this study, the internal and external dielectric constants were set to be 1.0 and 78.5 separately. The free energy decomposition analysis was carried out using an internal program with idecomp = 1.

Results and Discussion
Hard to instinctively predict the effect of the crowed multiple mutations in Omicron RBD on its binding to human ACE2 receptor and immune evasion There are 1273 amino acids in the WT spike (UniProt ID: P0DTC2), which is composed of S1 (residues 13 -685) and S2 (residues 686 -1273)/S2' (residues 816 -1273). The

Omicron RBD doesn't enhance binding affinity to ACE2
We performed 200ns all-atom molecular dynamics simulation for WT, Delta and Omicron variant RBD-ACE2 complex in triplicate and the mean value of three repetitions for each measurement was used for analysis. According to the predicted binding free energy (Table 1), the ΔGs of ACE2-RBDDelta in three repeating experiments are all about 10 kcal/mol higher than that of ACE2-RBDWT. The average number of ΔG of RBDDelta to ACE2 (-42.76 ± 2.38 kcal/mol) demonstrated stronger binding than that of that of RBDWT (-33.13 ± 3.26 kcal/mol). The result is in good agreement with existing experimental findings, proving that our method is reliable to a certain extent [3]. By contrast, the difference of ΔG between ACE2-RBDOmicron and ACE2-RBDWT is only 3.74 kcal/mol ( Table 1), suggesting that although there are a high number of mutations in S protein RBD, they don't increase affinity to ACE2 and even slightly weaken the binding.

Binding free energy decomposition
In order to understand the detailed implications of every mutation on the binding to ACE2, we conducted binding free energy decomposition based on all nine 50-150ns trajectories. As indicated in Figure 2A

The role of N501Y in Omicron variants
The binding free energy decomposition of residues ( Figure 2) shows that energy contribution of N501 in RBDWT is -2.12 ± 1.13 kcal/mol, while the energy contribution of Y501 in RBDOmicron is -6.97 ± 0.30 kcal/mol (a stronger attraction). By analyzing molecular dynamic trajectories, we found that Y501RBD and Y41ACE2 could form a strong π-π stacking interaction ( Figure 3A). As shown in Figure 3B, the center of mass distance between Y501RBD and Y41ACE2 is around 5.5 Å throughout the 200ns trajectory.

The role of Q493K and Q498R in Omicron variants
The common feature of Q493K and Q498R mutation is that both of them mutates from electrically neutral Q to positively charged amino acids (K and R). By binding free energy decomposition, the energy contributions of both K493 and R498 are decreased.
Specifically, the energy contribution of K493 in RBDOmicron decreased about 1.90 kcal/mol to Q493 in RBDWT, and the energy contribution of R498 in RBDOmicron decreased about 4.15 kcal/mol to Q498 in RBDWT. As shown in Figure 4, Q493 in RBDWT could form tight interactions with K31 and E35 in ACE2. By calculating distances of K31ACE2-Q493RBD and K31ACE2-K493RBD ( Figure 4A), we found that K31ACE2 moves away from K493RBD, which may be due to the repulsion of lysine with the same electrical properties. As shown in Figure 5C, K353ACE2 also moves away from R498RBD. Hence, we speculate that the repulsion between positively charged residues is the main cause of the decreased binding free energy in ACE2-RBDOmicron. Differently, K493RBD forms a tighter interaction with E35ACE2 ( Figure 4B

Omicron variant may cause immune evasion from some antibodies
To evaluate the binding affinity between three launched or clinical mAbs and two variant RBD, MM/GBSA calculations were carried out with RBDWT as control based on each MD runs lasting 200ns in triplicate. As shown in Table 2, the relative binding free energy to Etesevimab and BD-368-2 between WT and Delta RBD is -25.44, -17.21 kcal/mol, which means this variant may impair the effectiveness of these two antibodies.
In comparison, the difference in binding affinity to Bebtelovimab is only -4.81 kcal/mol, suggesting that the potency of this antibody is less likely to be affected by the mutant.    transmissibility. In addition, spike trimer may function differently from a single RBD.
Lastly, we calculate the binding energies between Omicron/WT/Delta RBD and mAbs and found that Omicron RBD might escape the launched Etesevimab and clinical BD-368-2 but may still sensitive to the launched Bebtelovimab.

Data Availability
The data that support the findings of this study are available from the corresponding authors upon reasonable request.