Immune Escape Mechanisms of SARS-CoV-2 Delta and Omicron Variants against Two Monoclonal Antibodies That Received Emergency Use Authorization.

Multiple-site mutated SARS-CoV-2 Delta and Omicron variants may trigger immune escape against existing monoclonal antibodies. Here, molecular dynamics simulations combined with the interaction entropy method reveal the escape mechanism of Delta/Omicron variants to Bamlanivimab/Etesevimab. The result shows the significantly reduced binding affinity of the Omicron variant for both antibodies, due to the introduction of positively charged residues that greatly weaken their electrostatic interactions. Meanwhile, significant structural deflection induces fewer atomic contacts and an unstable binding mode. As for the Delta variant, the reduced binding affinity for Bamlanivimab is owing to the alienation of the receptor-binding domain to the main part of this antibody, and the binding mode of the Delta variant to Etesevimab is similar to that of the wild type, suggesting that Etesevimab could still be effective against the Delta variant. We hope this work will provide timely theoretical insights into developing antibodies to prevalent and possible future variants of SARS-CoV-2.

Coronavirus disease 2019 (COVID-19), caused by the severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2), is spreading globally, with more than 520 million confirmed cases and 6 million deaths as of June 2022.[1] As a member of the coronavirus family, SARS-CoV-2 belongs to the β genus single-stranded RNA viruses, and its surface transmembrane spike glycoprotein (S-protein) can specifically recognize the host cell receptor angiotensin-converting enzyme 2 (ACE2) through the receptor-binding domain (RBD) to mediate virus invasion.[2−4] This process has been identified as a critical step of viral infection in the body, and the RBD has also become a significant target for drug and antibody research.[5,6]

Currently, monoclonal antibodies (mAbs) are promising for treating and preventing COVID-19, and the effectiveness of various mAbs has been validated.[7−11] Among them, Bamlanivimab (LY-CoV555) and Etesevimab (LY-CoV016) have received emergency use authorization (EUA) from the U.S. Food and Drug Administration (FDA) for the treatment of patients with mild to moderate symptoms due to their good response in clinical trials.[12,13] They are both RBD-targeting antibodies and can specifically bind to the RBD and interfere with the regular recognition between the S-protein and ACE2, preventing the virus from invading host cells (Figure A and B).

Figure 1

(A) The action mode of monoclonal antibodies (mAbs) hinders SARS-CoV-2 invasion into host cells. (B) Binding sites of the two mAbs (Bamlanivimab and Etesevimab) to RBD. (C) The mutation residues of the SARS-CoV-2 Delta and Omicron variants in the RBD are green and orange, respectively, and the binding sites of mAbs to the RBD are marked with dashed ellipses.

However, most antibodies are designed against early strains, and the persistent mutation of SARS-CoV-2 poses challenges for antibody efficacy and novel antibody development. As a key site for receptor binding and antibody targeting, mutations in the RBD can directly affect the ability of the virus to invade the human body and change epitopes. Many variants exhibit stronger virulence and infectivity[6,14−16] and even produce immune escape.[17−22] It is urgent to investigate whether the current epidemic lineages can escape existing antibodies.

As the two most popular mutation strains recently, Delta (lineage B.1.617.2) and Omicron (lineage B.1.1.529) have become the focus of many researchers.[17,23−28] Since its appearance in India in October 2020, the Delta variant has spread to more than 185 countries and was listed as a variant of concern (VOC) by the World Health Organization (WHO) on May 11, 2021.[29] There are two mutations, L452R and T478 K, on its RBD (Figure C), which can enhance the affinity between the virus and ACE2 and trigger immune escape.[30−32] Compared with the Delta variant, the Omicron variant appeared much later, but this did not affect its rapid spread. Since the Omicron variant was first detected on November 9, 2021, it was classified as a VOC by the WHO in just half a month, and it has become prevalent in most regions of the world.[29] More importantly, it has an astonishing 15 mutations on its RBD (G339D, S371L, S373P, S375F, K417N, N440K, G446S, S477N, T478K, E484A, Q493R, G496S, Q498R, N501Y, and Y505H) (Figure C). As for other VOCs, only one (N501Y) in the Alpha variant, three (K417N, E484K, and N501Y) in the Beta variant, three (K417T, E484K, and N501Y) in the Gamma variant, and two (L452R and T478K) in the Delta variant, the Omicron variants essentially contain their mutation sites except for L452. This suggests that Omicron variants may have unprecedented infection and immune evasion capabilities, further accelerating the spread of the epidemic.

Through the tireless efforts of scientists, some recent studies have evaluated the effects of the Delta and Omicron variants on Bamlanivimab and Etesevimab. The results show that the Delta variant can escape Bamlanivimab but is still inhibited by Etesevimab, while the Omicron variant has obvious immune escape to both antibodies.[33−36] These studies provide experimental guidance for the subsequent use and research of mAbs. However, the molecular mechanism of their interaction and escape remains to be improved, especially the identification of the key sites where the antibody-binding epitope changes due to mutation, which is helpful for the improvement of the design of mAbs.

In our work, molecular dynamics (MD) simulations were used to study the binding and escape mechanisms of the Delta and Omicron variants to Bamlanivimab and Etesevimab. Detailed energetic and conformational analyses were used to determine the origin of the binding differences caused by the mutations. For the binding free energy calculation, the enthalpy term was calculated by the molecular mechanics/generalized Born surface area (MM/GBSA) method,[37,38] and the entropy was obtained by the interaction entropy (IE) method.[39] In addition, alanine scanning combined with the IE (ASIE) method was used to calculate the contribution of individual residues to the binding affinity for key site prediction.[40−42] It is hoped that our work will provide timely theoretical insights into the development of vaccines and antibodies against Delta or Omicron and future novel variants.

To complement the existing experimental results from the molecular level, a total of six systems (Bamlanivimab–WT/Delta/Omicron and Etesevimab–WT/Delta/Omicron) were used for the MD simulations. All systems were in a stable binding state during the simulation with the root mean square deviation (RMSD) fluctuating almost below the 4 Å range (Figure S1A). When binding to Bamlanivimab, both variants exhibited higher flexibility compared to the wild type (WT), especially the Bamlanivimab–Omicron system, which was significantly more flexible at the binding interface than the other two systems (Figure S1B). When binding to Etesevimab, all systems exhibited similar flexibility in the binding interface (Figure S1C). The slightly higher flexibility is reflected in the binding interface of the antibody in the Etesevimab–Omicron complex.

Then, the binding free energy of each system was calculated by MM/GBSA and IE methods (Table S1), and the results were consistent with recently reported experimental evidence. Taking the randomness into account, two repeated simulations were performed additionally for each system, and the binding affinity was calculated by the same method. showing that three independent simulations obtained similar results. That is, the Omicron variant has significantly weaker binding free energy with both mAbs compared to the WT, while the Delta variant only exhibits lower binding affinity to Bamlanivimab, and its binding affinity with Etesevimab has no significant change compared to WT. In addition, the standard error of the mean for the binding free energy of each system was calculated (Table S1), and the results showed that the error was too small compared to the energy change caused by the mutations to affect the conclusion. From each energy term, we found that both the Delta and Omicron variants can significantly weaken the electrostatic interaction between the RBD and mAbs. This may be due to the mutations introducing plenty of positively charged amino acids (L452R and T478K for Delta; N440K, T478K, Q493R, Q498R, and Y505H for Omicron) while reducing the number of negatively charged amino acids (E484A for Omicron) so that more positive charges accumulated in the binding interface (Figure S2), which enhanced the electrostatic repulsion between the RBD and mAbs. In addition, the Omicron variant can weaken the van der Waals (vdW) interaction energy between the RBD and mAbs. This may be due to the obvious changes in the conformation and binding mode caused by the multisite mutation of the Omicron variant, leading to fewer atomic contacts at the binding interface. The heavy-atom contacts between the RBD and mAbs for different distance cutoffs in each system are calculated and support this speculation (Figure S3).

Quantitative calculation initially reveals the energy origin of the affinity change of each system, and the subsequent analysis of essential dynamics was used for the in-depth exploration of mutation-induced changes in the binding mode. The dynamic conformational sampling of each system is projected onto the first two eigenvectors (Figure A), and it is clear that the Bamlanivimab–Delta, Bamlanivimab–Omicron, and Etesevimab–Omicron complexes have a more dispersed distribution than the WT system, suggesting that they have lower structural stability affected by the mutation. Then, the characteristic dynamic fluctuations of each complex are depicted by the maximum and minimum eigenvalues of the lowest frequency principal component (Figure B). The result shows that all complexes are deflected to varying degrees around the interaction interface; the deflection angle is determined by the center of mass of the mAb of the structure corresponding to the maximum/minimum eigenvalue and the center of mass of the RBD. The Omicron variant has the highest rotational degrees of freedom of all systems, while the Delta variant exhibits greater torsion than the WT only when bound to Bamlanivimab. These changes further indicate that the structural stabilities of the Bamlanivimab–Delta, Bamlanivimab–Omicron, and Etesevimab–Omicron systems were reduced, supporting the results of the energy calculations.

Figure 2

(A) Two-dimensional eigenvector projections and (B) characteristic dynamic fluctuations of the Bamlanivimab–RBD and Etesevimab–RBD complexes. The WT, Delta, and Omicron systems are shown in blue, green, and orange, respectively.

Following the assessment of the global motion patterns, the local conformational transitions at the binding interface are discussed via changes in the distances of each residue between the RBD and the antibody (Figure S4). We define ΔD as the difference between the average distance for the Cα atoms of each residue in the equilibrium phase of the MD simulation and the original distance in the initial structure. Negative and positive values indicate closer or farther motion between residues, colored blue and red, respectively.

In the Bamlanivimab–RBD system, the distance change of the Delta variant is significantly different from that of the WT (Figure S4A and B). In the regions corresponding to 444–455@RBD (the region where the L452R mutation is located) and 488–501@RBD, the Delta variant has closer contact with the light chain and residues 100–110 of the heavy chain of Bamlanivimab than the WT, but the distance from residues 28–65 of the heavy chain of Bamlanivimab is farther than that of the WT. Through the structural observation of the Bamlanivimab–RBD complex (Figure S5A), we found that 100–110@Bamlanivimab is an independent region extending from the heavy chain, located between the light chain and the RBD. The Delta variant is close to this region while being away from the main part of the heavy chain (residues 28–65), which may reduce the contact area of the heavy chain with the RBD, resulting in a more unstable binding state.

In addition, residues 444–455 and 488–501 of the Delta variant were close to the light chain of Bamlanivimab, but the minor changes did not play a significant role compared to the greater distance between them. For the region near the T478K mutation, the area closest to the light chain of Bamlanivimab on the RBD, the Delta variant showed a tendency to move away from the light chain in this region. This would impair the binding of the two, further reducing the potency of the antibody. As for the Bamlanivimab–Omicron complex, it showed more red patches in the image than the WT (Figure S4C), especially in the region corresponding to the RBD and the antibody heavy chain, indicating that more residues separated from each other at the binding interface. This may result in fewer atomic contacts between the RBD and the antibody (Figure S3A), thereby exacerbating the escape of the Omicron variant to Bamlanivimab.

The Etesevimab–Delta and Etesevimab–WT systems showed similar distance changes (Figure S4D and E), which was consistent with the results of energy calculations that their binding free energies are similar. This implies that the Etesevimab work with the virus is less affected by the Delta variant and can still exert a neutralizing effect. As for the Etesevimab–Omicron system (Figure S4F), more large-area red patches appear compared to the WT. Although the part of region corresponding to 453–493@RBD appears blue, most of the other residues show a tendency to move away from the binding interface, resulting in fewer atomic contacts between Etesevimab and the Omicron RBD (Figure S3B). Through structural analysis of the complex (Figure S5B), we found that 453–493@RBD is a long loop region close to Etesevimab, while other regions of the RBD show a trend away from Etesevimab. This may cause these regions away from the antibody to drag 453–493@RBD and, thus, prevent it from approaching Etesevimab in subsequent movements, further promoting the escape of the Omicron variant to Etesevimab.

Next, hydrogen bond network analysis was used to further evaluate the differences in binding modes resulting from the mutation-induced conformational changes (Figure A and B). Bamlanivimab–Delta and Bamlanivimab–Omicron have more low-occupancy hydrogen bonds (occupancy <70%) but fewer high-occupancy hydrogen bonds (occupancy ≥70%) compared to the WT, especially for the Omicron variant, which has only one stable hydrogen bond. This may be since the conformational transition induced by the mutation interferes with the normal contact between Bamlanivimab and the RBD, thereby reducing the stability of the hydrogen bond network. Interestingly, more stable hydrogen bonds were formed between residues S494@RBD and E102@Bamlanivimab in the Bamlanivimab–Delta complex compared to the WT. According to the above analysis, this is due to the fact that residue E102 is located in the region where the heavy chain of Bamlanivimab extends (100–110@Bamlanivimab), and the Delta variant has closer contacts to this region. But with the separation between the RBD and the Bamlanivimab heavy chain main part, the occupancy of hydrogen bonds between residues E484@RBD and R50, R96@Bamlanivimab decreases. The Delta variant still has fewer high-occupancy hydrogen bonds compared to the WT; this results in the reduction of the binding stability between the RBD and Bamlanivimab.

Figure 3

(A) Number and occupancy of the hydrogen bonds formed at the interaction interface. (B) Stable hydrogen bonds (occupancy ≥70%) of each complex, marked with green dashed lines.

Both the Etesevimab–Delta and Etesevimab–Omicron systems have more low-occupancy hydrogen bonds than the Etesevimab–WT system. For high-occupancy hydrogen bonds, the stable hydrogen bond originally formed between residues R457@RBD and S53@Etesevimab in the WT system has decreased occupancy in the Delta system, but this phenomenon is compensated by the hydrogen bond formed between residues R403@RBD and Y92@Etesevimab so that Etesevimab–Delta and Etesevimab–WT have the same number of high-occupancy hydrogen bonds. As for the Omicron variant, the two hydrogen bonds involving residues R403@RBD and Y473@RBD become more unstable compared to the WT. This implies that the Delta variant may have less effect on the hydrogen bond network between Etesevimab and the RBD, while the more unstable hydrogen bond network exhibited by the Omicron variant may facilitate its escape from Etesevimab.

Then, the contributions of binding free energy for these residues near the binding interface were calculated using the ASIE method (Figure and Table S2), and they are generally considered important sites for mutation-induced changes in the binding affinity of the system. We define these residues with an absolute value of the binding free energy difference from the WT system ≥1 kcal/mol as hot-spot residues. In the Bamlanivimab–Delta complex, mutated residue 478 has a relatively obvious energy change, which is due to its mutation from neutral threonine to positively charged lysine, enhancing the electrostatic repulsion between the antibody and the RBD. Residues Q493 and S494 in Bamlanivimab–Delta are also more sensitive to mutation, which may be related to their proximity to the 100–110@Bamlanivimab region.

Figure 4

Binding free energy contribution of the important sites of the (A) Bamlanivimab–RBD and (B) Etesevimab–RBD complexes obtained by the ASIE method. The energy difference between the mutant system and the WT system (ΔΔG) is projected on the protein structure and residues with |ΔΔG| ≥ 1 are labeled in the images.

For the Bamlanivimab–Omicron system, it has more hot-spot residues than Bamlanivimab–Delta, probably due to the larger number of mutation sites of the Omicron variant. According to the results, these hot-spot residues mainly come from those that involve the mutation of charged amino acids, such as K417N, N440K, T478K, E484A, Q493R, and Q498R. The residues mutated to positively charged amino acids can induce attenuated electrostatic interaction energy, whereas mutations to negative do the opposite. In addition, residues V483 and F490 exhibit significant energy reduction; the former is mainly due to weakened electrostatic energy and enhanced polar solvation energy, and the latter comes from weakened vdW energy.

Etesevimab–Delta and Etesevimab–WT have similar binding modes, differing only slightly in residues 403 and 478. This may be since the charged residues introduced by the mutation are not close to the binding interface and do not cause a significant effect on the binding affinity. For the Etesevimab–Omicron complex, residues 417 and 478 have obvious energy changes among the mutated residues, both of which have weakened electrostatic interaction energy after mutation, which is not conducive to their interaction with the mAbs. Other residues involved in the mutation also have changes in the electrostatic interaction energy, but they are offset by the polar solvation energy, resulting in no obvious change in the affinity of these residues with the antibody. Furthermore, residues E406, R408, Y421, and N460 show obvious sensitivity to mutation. They clustered in a similar region, suggesting that large conformational changes may have occurred in this region, leading to significant energy changes in nearby residues.

As the epidemic progresses, the Omicron variant gradually replaces other mutant strains as the overwhelmingly dominant variant. To verify the reliability of the above findings and to further elucidate the escape mechanism of Omicron variants to mAbs, we extend the MD simulation of the Bamlanivimab/Etesevimab–Omicron systems to 1.5 μs. For comparison, the WT systems are also extended to the same time scale.

It is found that the RMSD of each system consistently fluctuated within a reasonable range and eventually stabilized around 4 Å (Figure S6A), indicating that the binding mode of mAbs–RBD stabilized and the systems had converged. The Bamlanivimab/Etesevimab–Omicron complexes have higher flexibility of residues near the binding interface than the WT systems (Figure S6B and C), suggesting that the Bamlanivimab/Etesevimab–Omicron systems had undergone a pronounced conformational shift that may trigger an unstable binding pattern between the RBD and the mAbs; this point is also confirmed by the following computation.

The binding affinity of the Omicron variant to Bamlanivimab/Etesevimab remains significantly weaker than that of the WT (Figure A). The difference mainly arises from the dramatic diminution in the electrostatic and vdW energies (Table S3); the former is caused by the large number of positively charged residues introduced through the mutation, and the latter comes from the unstable binding mode between the mutated RBD and the mAbs with fewer atomic contacts (Figure B). These are consistent with the 100 ns conclusions.

Figure 5

(A) Binding free energy calculated by the MM/GBSA and IE methods. (B) Average number of heavy-atom contacts between the RBD and mAb when the distance cutoff is set to 4–10 Å. (C and D) Superposition of the average structure of the complexes with the initial structure for Bamlanivimab/Etesevimab–WT/Omicron systems in the equilibrium phase of the 1.5 μs MD simulation.

The lowest frequency principal component is then used to describe the characteristic dynamic fluctuations of each system (Figure S7). Their dynamic characteristics are found to be similar to the results above; that is, the mAbs deflect on the axis of the binding interface, and the Omicron variant exhibits greater rotational freedom than the WT, suggesting unstable binding of the Omicron variant to the mAbs compared to the WT.

Subsequently, the distance changes between residues near the binding interface are calculated for further elaboration of the escape trend of the Omicron variant toward Bamlanivimab/Etesevimab (Figure S8). Apparently, the distribution of the distance changes of WT systems is similar to that of the 100 ns simulation, indicating that the WT is still stably binding to the mAbs. However, a large dark red patch appears in the Omicron variant systems, suggesting that it is away from the mAbs, resulting in less atomic contact and hydrogen bonding between the two (Figure S9).

Interestingly, although the Omicron variant exhibits a more unfavorable binding mode to the mAbs compared to the WT, there are still a few blue patches (Figure S8). To visualize the escape pattern of the Omicron variant toward mAbs, the average conformation of the equilibrium phase is used for comparison with the initial structure (Figure C and D). We find that the average structures of the Bamlanivimab/Etesevimab–WT systems are similar to the initial structure with minor differences, which are not sufficient to affect the sustained stabilization of the mAbs against the WT. In contrast, for the Bamlanivimab/Etesevimab–Omicron complexes, there is a significant deflection from the initial structure. The conformational change leads to an open hinge-like angle between the binding interface of the mAbs and the RBD, with one side of the mAbs gradually moving away from the RBD while the other side is temporarily immobile or slightly closer to the RBD. This explains why the Omicron variant shows a large number of red patches but still retains a few blue patches (Figure S8). However, with a gradual increase of the opening angle, the contact between mAbs and the RBD becomes less and less; the RBD will get rid of the mAbs’ grasp on it, thus producing the immune escape.

In summary, we dissect the impact of recently circulating SARS-CoV-2 Delta and Omicron mutant strains on Bamlanivimab and Etesevimab, two RBD-targeting mAbs already approved for therapy. The results show that the Omicron variant appears to be more prone to immune escape than the Delta variant, from both the energetic and conformational perspectives. The Omicron variant carries multiple mutations that, on the one hand, weaken the electrostatic interaction between the RBD and the mAbs by introducing a large number of positively charged amino acids. On the other hand, it triggers a significant conformational deflection resulting in less contact between the RBD and the antibody, thus reducing the binding stability of both. As for the Delta variant, the introduction of two positively charged amino acids in the RBD also results in a weakening of the electrostatic interaction energy between the RBD and the mAbs, but this change is balanced by the polar solvation energy. The real reason for the reduced affinity of the Delta variant for Bamlanivimab is that the mutation causes the RBD to be close to the extended region of the antibody (100–110@banranivirumab) but away from the main part of the heavy chain (28–55@banranivirumab). This change increases the rotational freedom of Bamlanivimab binding to the RBD and reduces the number of stable hydrogen bonds between the two, leaving the complex in an unstable binding state. In addition, the binding mode of the Delta variant to Etesevimab is similar to that of the WT, which may be because the mutation site is far away from the binding interface and the presence of the mutation does not cause significant conformational changes. Our work reveals the escape mechanism of mainstream SARS-CoV-2 variants against two mAbs that received EUA at the molecular level, which provides timely theoretical insights for the improvement of existing antibodies and the development of novel antibodies.