Omicron: A Heavily Mutated SARS-CoV-2 Variant Exhibits Stronger Binding to ACE2 and Potently Escapes Approved COVID-19 Therapeutic Antibodies.

The new SARS-CoV-2 variant of concern "Omicron" was recently spotted in South Africa and spread quickly around the world due to its enhanced transmissibility. The variant became conspicuous as it harbors more than 30 mutations in the Spike protein with 15 mutations in the receptor-binding domain (RBD) alone, potentially dampening the potency of therapeutic antibodies and enhancing the ACE2 binding. More worrying, Omicron infections have been reported in vaccinees in South Africa and Hong Kong, and that post-vaccination sera poorly neutralize the new variant. Here, we investigated the binding strength of Omicron with ACE2 and monoclonal antibodies that are either approved by the FDA for COVID-19 therapy or undergoing phase III clinical trials. Computational mutagenesis and free energy perturbation could confirm that Omicron RBD binds ACE2 ~2.5 times stronger than prototype SARS-CoV-2. Notably, three substitutions, i.e., T478K, Q493K, and Q498R, significantly contribute to the binding energies and almost doubled the electrostatic potential (ELE) of the RBDOmic-ACE2 complex. Omicron also harbors E484A substitution instead of the E484K that helped neutralization escape of Beta, Gamma, and Mu variants. Together, T478K, Q493K, Q498R, and E484A substitutions contribute to a significant drop in the ELE between RBDOmic-mAbs, particularly in etesevimab, bamlanivimab, and CT-p59. AZD1061 showed a slight drop in ELE and sotrovimab that binds a conserved epitope on the RBD; therefore, it could be used as a cocktail therapy in Omicron-driven COVID-19. In conclusion, we suggest that the Spike mutations prudently devised by the virus facilitate the receptor binding, weakening the mAbs binding to escape the immune response.

Introduction

Since its emergence, severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) has been found to continuously evolve and raise new variants of concerns (VOCs) to avoid host hostilities, i.e., evade the host immune response, increase transmission, and aggress the pathogenesis of coronavirus disease 2019 (COVID-19). This host adaptation by the virus has been demonstrated by the rise of VOCs, including Alpha, Beta, Gamma, and Delta variants that weaken the neutralizing efficacy of antibodies (1–4). Most recently, a new strain of the SARS-CoV-2 named Omicron by the World Health Organization has emerged in South Africa (November 24, 2021) and spread worldwide within a short period. Researchers around the globe are racing to determine whether Omicron poses a threat to the immunity induced by the COVID-19 vaccine (5).

Omicron harbors many novel mutations in structural and non-structural proteins, leading to serious concerns over vaccine failure, immune escape (5), and increased transmissibility. More than 32 mutations were found in the Spike protein alone, where 15 of these mutations reside in the receptor-binding domain (RBD), which are vital to both receptor and viral neutralizing antibodies. The non-structural proteins encoded by the ORF1ab contain mutations in the nsp3 (K38R, V1069I, Δ1265, L1266I, A1892T), nsp4 (T492I), nsp5 (P132H), nsp6 (Δ105-107, A189V), nsp12 (P323L), and nsp14 (I42V). Nsp3 (Plpro) and nsp5 (3Clpro, main protease) are proteases that cleave the polypeptide encoded by ORF1a and ORF1ab. 3Clpro and nsp12 [RNA-dependent RNA-polymerase (RdRp)] are primary targets for drugs that block the polypeptide cleaving and viral protein synthesis (6). Using structural models and as confirmed by the preliminary data in a preprint study (7), we found that mutations in nsp5 and nsp12 are not close to the active site and may not hinder the effect of antiviral drugs; nonetheless, these proteins play a vital role in innate immune response (interferon induction), requiring further experimental investigation (6). Omicron also had mutations in the other structural proteins, including Envelope (E) (T9I), Membrane (M) (D3G, Q19E, and A63T), and Nucleocapsid (N) (P13L, Δ31-33, R203K, G204R), further enhancing their infectivity. Since N protein is highly immunogenic (8, 9), these mutations could help escape the host immune response.

In addition, Omicron had multiple mutations in the Spike protein, which are associated with increased infectivity and antibody evasion. Out of 32 mutations, half of them hold the potential to dampen the potency of therapeutic antibodies and enhance the ACE2 binding. Omicron has also been shown to infect triple-vaccinated individuals who have received BNT162b2 jabs (10). Here, we conducted molecular modeling and mutational analyses to delineate how the new variant enhances its transmissibility and escapes against the FDA-approved Spike-neutralizing COVID-19 therapeutic antibodies. Our results may provide new insights into therapeutic management against the infection caused by Omicron.

Results

Mutations in the Omicron RBD Strengthen the Spike–ACE2 Interaction

Omicron is unique among the previously reported SARS-CoV-2 VOCs, showing multiple mutations in Spike and other genes. According to the unrooted phylogenic analysis using the global ~4,000 full-genome SARS-CoV-2 sequences from the Global Initiative on Sharing All Influenza Data (GISAID), Omicron stands distant from other VOCs ( ). A full-length trimeric 3D model was constructed by substituting the respective amino acids of previously reported reference (Wuhan strain, PDB ID: 7VNE) structure into Omicron. There are three deletion sites in the N-terminus domains (NTD) and at least 15 substitutions in the RBD region. Omicron Spike also harbors mutations that were reported in the previous VOCs such as K417, T478, E484, and N501 ( ). Of these, at least 11 mutations are involved in ACE2 binding ( ), substantially affecting their binding affinity ( ). In addition, Omicron Spike, compared with the prototype SARS-CoV-2, has three deletions, i.e., Δ69-70, Δ143-145, and Δ211, and one highly charged insertion, i.e., ins214EPE at 214 positions.

Figure 1

(A) Phylogeny of the Omicron and annotation of the mutation in Spike protein. The unrooted phylogenetic tree was constructed from the Nextstrain servers. Wuhan-Hu-1/2019 strains were taken as a reference sequence. (B) The full-length Delta and Omicron Spikes were built to annotate the relative (not exact) positions of the mutations on the surface map of Spike. (C) The amino acids mutated in the RBD of Omicron are shown concerning the ACE2 interface. Residues are colored according to the electrostatic map of the WT strain. The respective Omicron mutations are depicted in the panel below the RBD surface map.

Table 1

RBD binding interface residues of ACE2 and therapeutic antibodies.

ACE2	CT-p59	Sotrovimab	Etesevimab	Bamlanivimab	AZD1061	AZD8895	Casirivimab	Imdevimab
Lys417	Arg403	Asn334	Arg403	Tyr449	Arg346	Lys417	Lys417	Arg346
Gly446	Tyr449	Leu335	Lys417	Gly482	Lys444	Ala475	Tyr453	Asn440
Glu484	Asn450	Glu340	Asp420	Glu484	Tyr449	Gly476	Ser477	Lys444
Asn487	Tyr453	Asn343	Tyr421	Gly485	Asn450	Ser477	Glu484	Tyr449
Gln493	Glu484	Thr345	Leu455	Phe486	Glu484	Thr478	Phe486	Gln498
Gly496	Phe486	Arg346	Asn460	Gln493		Asn487	Tyr489
Gln498	Gln493	Lys356	Tyr473	Ser494		Tyr489	Leu492
Thr500	Ser494		Ala475	Tyr495		Gln493
Gly502	Tyr505		Asn487
Tyr505			Gln493

Bold residues are shared by ACE2 and mAbs on the RBD interface.

We monitored the relative binding strength of RBD–ACE2 complexes of both prototype and Omicron strains using a protein design strategy and calculating binding affinity and stability changes in terms of relative change in energies. We observed that the individually substituted residues had a slight effect on the local stability of the RBD–ACE2 complexes ( ). However, by performing endpoint molecular mechanics generalized born surface area (MMGBSA) binding free energy calculation, we could demonstrate a substantial increase in the binding affinity by T478K, Q493K, and Q498R, leading to an overall increase in the binding affinity of the RBDOmic with ACE2 (ΔGWT = −64.65 kcal/mol < ΔGOmic = −83.79 kcal/mol; ). We also investigated the change in electrostatic potential of the RBDOmic relative to that of RBDWT because the five residues in the RBM region of RBD are mutated from the polar to the positively charged residues (i.e., N440K, T478K, Q493K, Q498R, and Y505H). Surprisingly, we could observe that the electrostatic energy of ACE2–RBDOmic was double as that of ACE2–RBDWT, which in turn doubled the polar solvation free energies of the ACE2–RBDOmic ( ). Per residues, energy distribution suggests that mutations in RBDOmic directly participate in the binding, enhancing the binding strength of amino acids in the same network ( ). To validate our finding, we performed molecular dynamics simulations (MDS) using GROMACS (11) and calculated their binding free energies using the widely acceptable and more authentic molecular mechanics Poisson–Boltzmann surface area (MMPBSA) approach (12). As expected, we observed that RBDOmic, compared with RBDWT, had over 2.5 (BFE = −2,642.5 kJ/mol) times stronger binding affinity toward ACE2 (BFE = −951.9 kJ/mol). In addition, the electrostatic potential of RBDOmic was increased by ~1.5 times due to the polar-to-positive amino acids substitution ( ). In addition, energy perturbation per amino acid could confirm that the four amino acids, i.e., N440K, T478K, Q493K, and Q498R, directly contribute to the change of the total energy and the electrostatic potential, whereas K417N and E484A compensate the energy change ( ). Among 15 substituted amino acids, K417N and Y505H exhibited a slight reduction in binding energy due to the breakage of salt bridges between K417 of the RBD and D30 of ACE2; nonetheless, this breakage was compensated by the salt bridge between E35 of ACE2 and Q493K substitution in RBDOmic ( , right panel). Although the simulation time was short (20 ns), the root mean square deviation (RMSD) of the RBD proteins was not much different ( ). However, the number of hydrogen bonds between RBDWT–ACE2 showed a transient shift from high (N = ~7.5) to low (N = ~5) to high (N = ~7.5). This effect was not observed in RBDOmic–ACE2, and the hydrogen bond number remained consistent (N = ~7.5) ( ). Taken together, we suggest that Omicron binds ACE2 with greater affinity, partly explaining its increased transmissibility.

Figure 2

Relative effect of mutations in Omicron RBD on the ACE2 binding. (A) Effect of 15 individual mutations on the binding and stability of RBDOmic–ACE2 was monitored relative to that of RBDWT–ACE2. (B) The binding free energies (measured through MMPBSA) as consequences of all 15 mutations at once were monitored for both RBDOmic–ACE2 and RBDWT–ACE2. (C) Per-residue energy contribution was monitored, and the hotspots of RBD were labeled. The change in the hydrogen bond network of the selected hotspots is shown at the right. (D, E) Root mean square deviation and hydrogen bonds at the RBD–ACE2 interface as a function of time are displayed for both RBDOmic–ACE2 and RBDWT–ACE2 complexes.

Mutations in the RBDOmic Deteriorate the Binding of Therapeutic Antibodies and Garble Their Epitopes on RBD

To evaluate whether the mutations that strengthen the RBDOmic–ACE2 interaction affect the RBD-targeting COVID-19 therapeutic antibodies, we constructed structural models of eight monoclonal antibodies (mAbs) bound to RBDOmic. The antibodies like CT-p59, developed by Celltrion from the peripheral blood mononuclear cells (PBMCs) derived from a convalescent plasma of COVID-19 patients, and sotrovimab are used as solo COVID-19 therapeutics undergoing phase III clinical trials (13, 14). The other six mAbs were approved for COVID-19 therapies on an emergency basis (15), which are used as cocktail therapy to tackle the immune escape by the newly acquired mutants ( ).

Figure 3

Mutations in the Omicron RBD distort the epitopes of therapeutic mAbs. (A–D) Crude epitopes of seven selected mAbs are shown on the RBD. Antibodies used as cocktails are labeled with their sponsors. All variable light chains are colored yellow or orange and variable heavy chains are colored red. (E) Changes in the binding affinity of the RBDOmic–mAbs relative to RBDWT–mAbs are shown. The binding energies were calculated through endpoint MMGBSA. (F) The binding free energies (measured through MMPBSA) as consequences of all 15 mutations at once were monitored for RBDOmic–etesevimab and RBDWT–etesevimab. (G) The binding free energies (measured through MMPBSA) as consequences of all 15 mutations at once were monitored for RBDOmic–CT-p59 and RBDWT–CT-p59.

Since mAbs in their respective cocktail therapy regimen do not share overlapping epitopes on the RBD, except etesevimab and bamlanivimab (sponsored by AbCellera) where the light chain variable domains show a slight clash ( and ) and are capable of neutralizing the virus independently, we investigated the change in their interface and binding strength of the mAbs with RBDOmic individually. Remarkably, we found a substantial drop in the total binding energies of bamlanivimab and CT-p59 when bound to RBDOmic ( ). The total binding energy is the sum of the four energies (listed in ). We could see that vdW (Van der Waals potentials) and SA (solvation free energy) energies did not affect the binding strength; however, electrostatic potentials (ELE) had a significant shift in the calculated energies. The magnitude of the change in ELE energies was similar to that of RBDOmic–ACE2; however, the effect was opposite, i.e., in RBDOmic–ACE2 ELE energies favor the binding, whereas RBDOmic–mAbs ELE opposed the binding strength ( ). All the mAbs showed a significant drop in ELE energies, but AZD1061 (AstraZeneca) showed a slight drop (RBD-AZD1061 = −204.24 kcal/mol > RBDOmic-AZD1061 = −112.35 kcal/mol). To validate the relative change of total binding energies and the shift of the electrostatic potential of RBDOmic–mAbs, we analyzed the MDS of the two complexes (i.e., RBDOmic–etesevimab and RBDOmic–CT-p59) and calculated their BFE using MMPBSA. Interestingly, we found that the change of the electrostatic potential of both complexes exhibited a similar trend for both MMGBSA and MMPBSA approaches. The MMPBSA values were calculated as the statistical outcome of 100 frames extracted from the 20-ns MDS trajectory ( ). The total binding energies for RBDOmic–mAbs complexes were substantially higher compared with the RBDWT–mAbs, suggesting that Omicron can escape both etesevimab and CT-p59 (regdanvimab) ( ).

Table 2

The binding energies of RBDOmic and RBDWT with seven therapeutic mAbs are listed.

Sponsor	mAbs	VDW	ELE	GB	SA	Total
AbCellera and Eli Lilly	Etesevimab	−131.41	−179.48	243.28	−17.15	−84.76
	Omic-Ete	−147.42	65.04	19.31	−18.16	−81.24
	Bamlanivimab	−95.81	−13.97	33.7	−13.01	−89.1
	Omic-Bam	−108.74	602.59	−532.27	−13.12	−51.54
AstraZeneca	AZD8895	−84.37	−30.35	71.27	−10.46	−53.92
	Omic-Az95	−96.77	76.73	−24.25	−11.62	−55.91
	AZD1061	−94.96	−204.24	238.04	−11.49	−72.64
	Omic-Az61	−98.83	−112.35	149.92	−12.29	−73.55
Regeneron	Imdevimab	−76.1	−27.36	74.73	−9.67	−38.4
	Omic-Imd	−81.32	−5.84	64.23	−10.15	−33.08
	Casirivimab	−103.52	−174.81	228.04	−14.66	−64.94
	Omic-Cas	−113.25	−55.14	119.9	−14.42	−62.91
Celltrion	CT-p59	−105.66	−7.38	49.61	−13.75	−77.18
	Omic-CT-p59	−114.85	352.19	−278.48	−14.27	−55.4

mAbs, monoclonal antibodies; VDW, Van der Waals potentials; ELE, electrostatic potentials; GB, polar solvation potentials (generalized born model); SA, non-polar contribution to the solvation free energy calculated by an empirical model.

Next, to evaluate which mutations are mainly involved in weakening the RBDOmic–mAbs interactions, we calculated per amino acid energy perturbation for CT-p59 and bamlanivimab when bound to RBDWT and RBDOmic. We observed that two hotspots, i.e., R96 in CDRL3 and R50 in CDRH2 of bamlanivimab, established highly stable salt bridges with the E484 of RBDWT, losing their binding entirely upon E484A mutation in RBDOmic ( ). In addition, E102 and R104 in CDRH3 showed a 50% reduction in binding energies. Similarly, the hotspots in CDRL1 and CDRH3 lost their bindings due to the mutations of E484A, Q493K, and Y505H in RBDOmic. By contrast, N501Y slightly strengthens the binding, establishing a hydrogen bond with D57 in CDRH2 ( ). Next, we extended our search and used two more SARS-CoV-2 neutralizing antibodies, C102 and C105, isolated from the convalescent plasma of a single donor for their RBDOmic affinity (16). Both mAbs bind to an overlapping epitope on the RBM and predominantly utilize the heavy chain variable domain (VH) CDRs as main paratopes ( ). Like other RBD-binding mAbs, both C102 and C105 exhibited a drastic reduction in binding affinity driven by a significant change in the electrostatic potential ( ). To investigate how sotrovimab retains its neutralization efficacy, we constructed the RBDOmic–sotrovimab model and found that sotrovimab binds to a highly conserved epitope on the RBD and, among 15 mutations in the RBDOmic, faces only G339D mutation. The retained Omicron neutralization in pseudovirus assay (17) may indicate that the salt bridges between the CDRH3 and RBD may override the clash between RBDOmic D339 and Y100 in CDRH3 that could potently destabilize the RBDOmic–sotrovimab interaction ( ). These results suggest that mutations in the Omicron Spike are precisely designed by the virus, facilitating receptor binding but hindering antibody binding simultaneously and that antibodies recognizing conserved epitope on the Spike of SARS-CoV-2 variants could be used as pan-variant therapeutics. Overall, the escape of Omicron from a large pool of antibodies, especially those approved by the FDA after undergoing extensive clinical trials and safety measures, raises serious concerns about the efficacy of therapeutic mAbs in Omicron-infected patients.

Figure 4

Per-residue changes in the binding affinity of RBD–mAbs were monitored and the hotspots on CDRs of (A) bamlanivimab and (B) CT-p59 are labeled. The change in the hydrogen bond network of the selected hotspots is shown at the right.

Discussion

To this end, it is well known that SARS-CoV-2 is rapidly evolving and makes at least two mutations per month in its genome (18, 19). The virus is capable of adapting to the host environment by increasing transmissibility and evading immune response, as exemplified by the continuous rise of VOCs (20, 21). Although tremendous efforts have been made in vaccine development and COVID-19 therapeutics, including mAbs and COVID-19 pills by Merck, the emergence of VOCs has raised concerns over the efficacy of neutralizing antibodies (21, 22). Even though these variants had a limited number of mutations, they successfully escaped the immune response, at least partly if not entirely. Omicron harbors four or five times more mutations in the Spike protein than other SARS-CoV-2 VOCs and raise more serious concerns (see ).

We used the previously available structural data of Spike RBD-binding antibodies, Spike itself, and Spike–ACE2 complexes and constructed the mutant Omicron Spike. Omicron Spike contains some of the mutations reported in the previous VOCs. In particular, D614 enhances the receptor binding by increasing its “up” conformation and the overall density of Spike protein at the surface of the virus (23, 24). In addition, five of the amino acids within the RBD region are mutated from polar to positively charged residues (K, R, or H) that paradoxically enhance the receptor binding and weaken the Spike neutralizing interactions ( and ). The RBD mutation has been mapped to predict the neutralization escape from REGN-COV2 (a cocktail of REGN10933 and REGN10987) and LY-CoV016. Such mutations have already been found in patients with persistently infected COVID-19 since late 2020 or early 2021 (25). Single E406W mutation can lead to the viral escape from both antibodies in REGN-COV2, whereas F486K has been reported to escape REGN10933. N440K and K444Q can escape against REGN10987, while K417N, N460T, and A475V can successfully escape against LY-CoV016 antibody (AbCellera), currently approved by the FDA for COVID-19 therapy (25). Unfortunately, Omicron has mutations or amino acids adjacent to those predicted to escape the neutralization of antibodies.

Among the investigated antibodies here, we suggest that AZD1061 may be able to retain Omicron neutralization ( ). During our study, two research groups investigated the neutralization escape of Omicron from the same set of antibodies. They reported consistent results for the mAbs sponsored by Regeneron and AbCellera. Nonetheless, they demonstrated that the Omicron pseudovirus neutralization findings differed in AZD1061 (cilgavimab) and AZD8895 (tixagevimab). In support of our results, Planas et al. have shown that AZD1061 (cilgavimab) and the cocktail (AZD1061+AZD8895), but not AZD8895 alone, could retain the Omicron pseudovirus neutralization (17).

On the other hand, Cao et al. have demonstrated that AZD8895, but not AZD1061, binds Omicron with lower affinity and slightly neutralizes the Omicron pseudovirus at a very higher concentration (26). However, they did not examine the effect of the AZD1061+AZD8895 cocktail on the Omicron. The above studies suggested that sotrovimab (VIR-7831/GSK4182136/S309) holds the promising neutralization efficacy against the Omicron pseudovirus. Thus, we constructed an RBDOmic–sotrovimab model to investigate how sotrovimab retains its neutralization efficacy. Sotrovimab could bind to the highly conserved epitope on the RBD, and among the 15 mutations in the RBDOmic, it faces only G339D mutation.

Although the pathological manifestations are thus far reported to be mild, the threat of Omicron is global, and it is also quite clear that the new variant is more transmissible than Delta. Vaccination has been reported to significantly drop the COVID-19 infection of emerging VOCs, including Delta (27). However, the sera from convalescent subjects infected with different variants of SARS-CoV-2 including Alpha, Beta, Gamma, and Delta and vaccines were found to be ineffective against Omicron. Nonetheless, immunity boosted by the third dose of vaccines or vaccinees infected by the Delta strain has shown effect against Omicron (28). Another study has also reported some preliminary data about the ineffectiveness of vaccines against Omicron. The sera from individuals with the 5-month post-vaccination with Pfizer (29) or AstraZeneca vaccine have failed to inhibit Omicron (30, 31). Similarly, the sera from 6- to 12-month post-infection individuals could not neutralize the new variant. Nevertheless, 5- to 31-fold lower neutralization of Omicron compared with Delta has been reported by boosters and the previously infected vaccines (17). Thus, Omicron can escape against the therapeutic and vaccine-elicited antibodies. Our study and the preliminary data from other studies consistently suggest that a cocktail of Evusheld (AstraZeneca mAbs) and sotrovimab (GSK, S203 mAb) could effectively neutralize the Omicron.

Methods

Model Construction and Optimization

For the full-length trimeric Spike, a previously reported PDB ID: 7VNE was used to rebuild the Omicron Spike protein using a Swiss-Model server (32). Other structures used in this study are listed as follows: RBD–ACE2 (PDB ID: 6MOJ), RBD–etesevimab (PDB ID: 7C01), RBD–bamlanivimab (PDB ID: 7KMG), RBD–CT-p59 (PDB ID: 7CM4), RBD–AZD1061 (PDB ID: 7I7E), RBD–AZD8895 (PDB ID: 7I7E), RBD–casirivimab (PDB ID: 6XDG), RBD–imdevimab (PDB ID: 6XDG), and RBD–sotrovimab (PDB ID: 7R6X). For constructing the mutant RBD, free BIOVIA Discovery Studio Visualizer was used (http://www.accelrys.com). All complexes were solvated with TIP3P water cubic box of dimension boundaries extended to 10 Å from protein atoms and neutralized with counter ions, Na+/Cl−, wherever needed. The neutralized systems were energy minimized in GROMACS 2019.6 (33) using CHARMM37 force field (34) and steep descent algorithm. For endpoint binding free energy calculations, the HawkDock server was utilized (35).

Molecular Dynamics Simulation

To calculate the binding free energies of RBDOmic with ACE2 and antibodies, we utilized the GROMACS package for the generation of trajectories and the MMPBSA tool for the free energy perturbation. Each system was solvated in a dodecahedron box filled with TIP3P water model and neutralized by adding counter ions (Na+/Cl−). The neutralized systems were energy minimized as stated above. Next, a two-step equilibration was set up under constant temperature (NVT) and constant pressure (NPT) of 0.2 ns, and the systems were equilibrated. The temperature and pressure were coupled with v-rescale (modified Berendsen thermostat) and Berendsen, respectively (36). The long-range electrostatic interactions were computed by utilizing the particle mesh Ewald algorithm (37). Each system was simulated for 20 ns with all constrains removed. For the calculation of RMSD and hbonds, MD trajectories were converted by removing the jumps and translational and rotational motions using −pbc nojump and −fit rot+trans flags under the trjconv tool in GROMACS. For MMPBSA, every 20th frame was extracted from a 20-ns trajectory in a separate trajectory.

Binding Free Energy Calculation

The MMPBSA (12) approach is best suited for calculating binding free energies of the ligands bound to the same target. Here, RBD is the main target, whereas mAbs and ACE2 are considered as ligands. GROMACS is equipped with g_mmpbsa tool which was used for the calculation of binding free energies. The topology of each system was generated through the older version of GROMACS (v 5.0) as the MMPBSA package has not been updated by the developer till now. The binding energies were calculated according to the equations described in our previous study (38).

Computational Tools Used in This Project

For protein structure visualization, VMD (39), PyMOL (https://pymol.org), and Chimaera Chimera (40) packages were used. For electrostatic surfaces, isolation of the proteins, APBS and APBSrun plugins in PyMOL and VMD were utilized. The interfaces of RBDWT and RBDOmic with mAbs and ACE2 were analyzed by the online server PDBePISA (v1.52) (41), and the binding contribution of individual amino acids was determined. For endpoint binding free energy calculations, the HawkDock server was utilized (35). The hotspot results were validated through the DrugScorePPI web server (42). The unrooted phylogenetic tree was constructed from the Nextstrain (43) servers using ~4,000 full-length SARS-CoV-2 sequences from the GISAID (44) database with reference to Wuhan-Hu-1/2019 as a reference sequence.

Data Availability Statement

The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

Author Contributions

MS and HW contributed toward conceptualization of the project and designed the methodology. MS and HW wrote the original manuscript draft. HW supervised the study and provided funding acquisition. All authors contributed to the article and approved the submitted version.

Funding

This research was supported by grants from the National Research Foundation of Korea (NRF) funded by the Ministry of Science and ICT (MSIT) (NRF-2017M3C9A6047620, NRF-2019R1A5A2026045, and NRF-2017M3A9B6061509) and grant from the Korea Health Industry Development Institute (KHIDI) funded by the Ministry of Health & Welfare, Republic of Korea (HI21C1003). In addition, this study was also supported by KREONET (Korea Research Environment Open NETwork), which is managed and operated by KISTI (Korea Institute of Science and Technology Information).

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.