Mutational Scanning and Binding Free Energy Computations of the SARS-CoV-2 Spike Complexes with Distinct Groups of Neutralizing Antibodies: Energetic Drivers of Convergent Evolution of Binding Affinity and Immune Escape Hotspots

1. Introduction

In-depth structural and biochemical studies of the SARS-CoV-2 Spike (S) glycoprotein have provided critical insights into the mechanisms governing viral transmission and immune evasion [1,2,3,4,5,6,7,8,9,10,11,12,13,14,15]. The S glycoprotein is characterized by remarkable conformational flexibility, which plays a pivotal role in its function. This flexibility is particularly pronounced in the S1 subunit, which comprises several essential domains: the N-terminal domain (NTD), the receptor-binding domain (RBD), and two conserved subdomains, SD1 and SD2. The dynamic nature of the S1 subunit enables the S glycoprotein to undergo significant structural rearrangements, facilitating its adaptation to various stages of viral entry, including receptor binding and membrane fusion. The NTD of the S glycoprotein plays a crucial role in the initial attachment to host cells, facilitating the virus’s first contact with the host membrane. In contrast, the RBD is responsible for binding to the angiotensin-converting enzyme 2 (ACE2) receptor, a critical step for viral entry into host cells. Moreover, the RBD transitions between “up” and “down” conformations, modulating its accessibility to both the receptor and neutralizing antibodies. This conformational plasticity not only enhances the virus’s ability to infect host cells but also contributes to its capacity to evade immune detection. The conserved subdomains, SD1 and SD2, contribute to the structural integrity of the S glycoprotein by stabilizing its prefusion conformation. These subdomains also play a pivotal role in orchestrating the transition to the postfusion state, a process essential for membrane fusion and subsequent viral entry [1,2,3,4,5,6,7,8,9,10,11,12,13,14,15]. Together, these structural features underscore the S glycoprotein’s adaptability and highlight its central role in SARS-CoV-2 pathogenesis and immune evasion strategies [10,11,12,13,14,15]. Biophysical studies have further elucidated the thermodynamic and kinetic principles governing the behavior of the S glycoprotein. These investigations have revealed the energy barriers and conformational dynamics that regulate its transitions between functional states. Such insights have deepened our understanding of the S glycoprotein’s role in viral infection and immune evasion, highlighting its adaptability and complexity. For instance, the interplay between the NTD, RBD, and conserved subdomains underscores the finely tuned mechanisms that enable SARS-CoV-2 to balance receptor binding, membrane fusion, and immune escape [16,17,18]. Collectively, these findings provide a comprehensive framework for understanding S glycoprotein’s central role in viral pathogenesis and inform the development of targeted therapeutic and preventive strategies.

The extensive array of cryo-electron microscopy (cryo-EM) and X-ray structures of SARS-CoV-2 Spike (S) protein variants of concern (VOCs) in various functional states, along with their interactions with antibodies has provided significant insights into the virus’s adaptability. These studies have demonstrated that VOCs can induce structural changes in the dynamic equilibrium of the S protein. Such alterations influence the distribution of functional states, impacting the binding affinities of the S proteins with different classes of antibodies and determining the efficacy of these antibodies in neutralizing the virus [19,20,21,22,23,24,25]. These structural changes are of paramount importance because they directly influence how the S protein interacts with the host immune system. For example, certain VOCs can adopt specific conformations that reduce their recognizability by pre-existing antibodies, enabling them to evade immune detection and neutralization. This immune evasion is often achieved through mutations that stabilize the S protein in a closed or shielded state, limiting antibody access to key epitopes, particularly in the RBD. Conversely, other conformational states may expose previously hidden or novel epitopes, creating opportunities for the immune system to recognize and target the virus more effectively. The ability of VOCs to dynamically shift between these conformational states highlights the remarkable adaptability of the S protein and its role in viral fitness. This plasticity not only facilitates immune evasion but also enhances the virus’s capacity to infect host cells by optimizing receptor binding and membrane fusion.

The BA.2.86 variant, a distinct sublineage of the Omicron variant, emerged in mid-2023 and rapidly drew global attention due to its remarkable genetic divergence from previously circulating strains [26,27]. This variant stands out for harboring more than thirty mutations in its S protein compared to its predecessor, BA.2, marking one of the most significant evolutionary leaps since the initial emergence of the Omicron variant from the ancestral Wuhan strain. Such a substantial genetic shift underscores the virus’s continued capacity for rapid adaptation and evolution, posing new challenges for public health and vaccine efficacy. The evolutionary trajectory of BA.2.86 is believed to be driven by a combination of factors, including accelerated mutation rates in chronically infected individuals and recombination events between pre-existing viral lineages [26,27,28]. These processes may have facilitated the accumulation of mutations that confer selective advantages, such as enhanced transmissibility, immune evasion, and receptor binding affinity. This dual ability—strengthening receptor interactions while avoiding immune detection—highlights the variant’s evolutionary sophistication and its potential to undermine existing immunity. Structural studies of BA.2.86 have been instrumental in unraveling the functional implications of its mutations, particularly in the S protein. Key mutations, including R403K, F486P, and R493Q, have been identified as critical determinants of the variant’s altered biological properties [29,30].

The JN.1 variant, which evolved from the BA.2.86 lineage, represented a major step in the ongoing evolution of SARS-CoV-2. Detected in late 2023, JN.1 quickly became the dominant strain due to its increased ability to spread and evade the immune system [31,32,33]. Structural studies showed that key mutations, like L455S in the RBD improved its ability to bind to the ACE2 receptor and avoid neutralizing antibodies [34,35]. KP.2, a descendant of JN.1, emerged with additional mutations, such as R346T and F456L, which further boosted its ability to evade immunity and spread more efficiently [36,37,38]. Similarly, KP.3, another subvariant of JN.1, carried mutations like Q493E and F456L, which worked together to strengthen ACE2 binding and increase resistance to antibodies [39,40,41]. These changes made KP.3 one of the fastest-spreading variants in 2024.. JN.1 subvariants LB.1 (JN.1 + S:S31-, S:Q183H, S:R346T, S:F456L), and KP.2.3 (JN.1+ S:R346T, S:H146Q, S:S31-) which convergently acquired S31 deletion in addition to the above substitutions, have spread as of June 2024 and contribute to immune evasion and the increased relative effective reproduction number [39,40].These changes further enhanced immune evasion and transmissibility [39,40]. Meanwhile, the XEC variant, a recombinant strain combining elements of KS.1.1 and KP.3.3, appeared in mid-2024 with additional mutations like T22N and F59S. These mutations improved its ability to infect cells and evade immune responses, making it a potential candidate to become the next dominant strain [42,43]. The evolution of JN.1, KP.2, KP.3, and XEC variants demonstrated SARS-CoV-2’s ability to adapt through mutations that enhance its spread and immune evasion [36,37,38,39,40,41,42,43]. This ongoing evolution underscores the importance of continuous surveillance and adaptive vaccine strategies to keep pace with the virus.

The rapidly expanding body of structural studies on SARS-CoV-2 antibodies has provided critical insights into the molecular mechanisms underlying their binding competition with the ACE2 host receptor. These studies have revealed that the S protein contains multiple distinct antigenic sites, which can be targeted to achieve efficient cross-neutralization. By synergistically targeting both conserved and variable epitopes within the RBD antibodies can effectively neutralize the virus, even in the face of emerging variants [44,45,46]. Initially, SARS-CoV-2 antibodies were broadly categorized into several main classes based on their binding characteristics. Among these, class 1 and class 2 antibodies are particularly notable because they target epitopes that overlap with the ACE2 binding site. This overlap allows these antibodies to directly block the virus’s ability to attach to and enter host cells, making them potent tools for neutralization [44,45]. To further understand the antigenic landscape of the S protein and the mechanisms of antibody-mediated neutralization, researchers have extensively examined the diversity of binding epitopes in antibody-S protein complexes. These studies have led to a detailed classification of antibodies into distinct categories, each defined by their specific binding sites and neutralization mechanisms [47,48]. The structural principles underlying these classifications, as well as the mechanisms by which antibodies neutralize the virus, have been comprehensively summarized in several reviews [49,50,51,52]. These studies indicated that antibodies targeting conserved epitopes are more likely to retain their efficacy across different variants, while those targeting variable regions may need to be updated to address emerging strains. Moreover, the interplay between antibody binding and the conformational dynamics of the S protein has been a key focus of research. Studies have shown that certain antibodies can lock the S protein in specific conformations, preventing it from transitioning to the open state required for ACE2 binding. Others may destabilize the S protein, disrupting its ability to mediate membrane fusion. These insights have not only advanced our understanding of SARS-CoV-2 neutralization but also provided a framework for designing next-generation antibodies that can effectively combat the virus’s evolving strategies for immune evasion.

By combining biophysical and structural data, researchers have classified RBD-directed antibodies into seven distinct groups, often labeled RBD-1 to RBD-7 [49]. These antibodies can be further divided into four main classes based on how they interact with the RBD. Class I antibodies neutralize the virus by directly competing with ACE2, binding only to the “up” conformation of the RBD. Class II antibodies can bind to both “up” and “down” RBD conformations. Class III antibodies target areas outside the immediate ACE2 binding site, attaching to the RBD in both “up” and “down” states. Class IV antibodies bind to a hidden (cryptic) epitope that is only accessible in the “up” conformation. A similar classification system was developed using 54 human antibodies targeting the RBD, grouping their binding sites into six different classes [53]. Combining antibodies from different classes—especially those targeting both conserved and variable RBD epitopes—can enhance cross-neutralization and reduce the risk of the virus mutating to escape immunity [54,55]. Some highly potent antibodies work in unique ways, such as causing allosteric changes in the S protein. These changes can block receptor binding or disrupt other critical functions, preventing infection without directly interfering with ACE2 recognition. This approach adds another layer of defense against the virus’s ability to adapt and evade the immune system.

High-throughput yeast display screening has been employed to map the profiles of RBD escape mutations for 247 human anti-RBD neutralizing antibodies [56,57,58]. Through unsupervised clustering, these antibodies were classified into six distinct epitope groups (A–F) without relying on structural data, yet this classification aligns well with earlier structural studies [44,45,47,53]. Groups A–D in this classification correspond to the RBS A–D neutralizing antibodies described in [53] and overlap with the class 1–2 neutralizing antibodies outlined in [45]. The binding epitopes for antibodies in groups A–D overlap with RBD residues that are critical for ACE2 binding, making them particularly effective at blocking viral entry. Group A and B neutralizing antibodies, such as LY-CoV016 and AZD8895, specifically bind to the RBD in its “up” conformation, which is the state required for ACE2 interaction. In contrast, group C and D antibodies, including LY-CoV555 and REGN-10987, exhibit broader binding capabilities, attaching to the RBD in both “up” and “down” conformations [59]. This flexibility allows them to neutralize the virus more effectively across different structural states. Groups E and F antibodies, which are analogous to class 3 and 4 antibodies in other classifications, target epitopes that do not overlap with the ACE2 binding site [45]. These antibodies neutralize the virus through mechanisms that do not involve direct competition with ACE2, such as stabilizing the RBD in a conformation that prevents its interaction with the receptor or inducing allosteric changes that disrupt viral function. This comprehensive classification system, supported by both functional and structural data, provides a detailed understanding of how different antibodies target the SARS-CoV-2 RBD. It highlights the diversity of neutralizing strategies and underscores the importance of targeting both conserved and variable epitopes to achieve broad and resilient protection against evolving viral variants.

A recent groundbreaking study used advanced methods to analyze the binding sites (epitopes) of antibodies produced after a breakthrough infection with the BA.1 variant in vaccinated individuals. This study identified the escape mutations for 1,640 antibodies targeting the RBD of the SARS-CoV-2 spike protein and grouped them into 12 distinct epitope categories [60]. Groups A–C include antibodies that target the ACE2-binding site, directly blocking the virus from attaching to human cells. Group D antibodies, such as REGN-10987, LY-CoV1404, and COV2-2130, bind to a specific region (residues 440–449) on the RBD and are further divided into subgroups D1 and D2. Groups E and F are split into smaller subgroups (E1–E3 and F1–F3) that cover the front and back of the RBD. Antibodies in groups A, B, C, D1, D2, F2, and F3 can neutralize the virus by competing with ACE2, while those in groups E1, E2.1, E2.2, E3, and F1 do not compete with ACE2 but may still neutralize the virus through other mechanisms [60]. In a follow-up study, this data were further expanded with antibody escape profiles from individuals who recovered from BA.2 and BA.5 infections. Using advanced computational techniques like multidimensional scaling and t-SNE (t-distributed stochastic neighbor embedding), they mapped the epitope groups for an even larger set of antibodies [61]. This study revealed that the BA.5 variant, particularly with the R346T mutation, can escape most antibodies in groups D1, E1, and E2.1/E2.2. Additionally, the K444N mutation helps the virus evade antibodies in group D2. These findings highlight how specific mutations allow the virus to dodge immune responses more effectively [61]. Another study used deep mutational scanning (DMS) to analyze 781 monoclonal antibodies targeting the RBD, isolated from individuals with repeated Omicron infections. By using graph-based clustering to analyze escape mutations, researchers identified 12 major epitope groups on the BA.5 RBD. This study also showed that repeated Omicron infections trigger the production of more Omicron-specific antibodies, which target different RBD epitopes and escape mutations compared to antibodies generated after infection with the original (wild-type) virus [62].

A recent study used advanced high-throughput yeast-display-based DMS assays to explore the molecular mechanisms behind broadly neutralizing antibodies (Abs) produced after XBB and JN.1 variant infections. This work analyzed the escape mutation profiles of 2,688 antibodies, including 1,874 isolated from individuals infected with XBB or JN.1, and grouped them into 22 distinct clusters based on their binding characteristics [63]. Among these clusters, certain groups of antibodies stood out for their ability to neutralize the virus effectively. Groups A1/A2 (Class 1) and B (Class 1/2) compete directly with ACE2, blocking the virus from binding to human cells. Groups D2/D3/D4 and F3 can also compete with ACE2 and show strong neutralizing potential. Some antibody groups do not compete with ACE2 and are less likely to neutralize the virus effectively, including groups E1/E2 (Class 3), E3 (Class 5), and F1 (Class 4). These antibodies target different regions of the RBD and S protein and are less potent in preventing viral entry [63]. The study also revealed that reinfection with the JN.1 variant triggers the production of Omicron-specific neutralizing antibodies with significantly broader neutralizing capabilities compared to infections caused by the XBB or HK.3 variants. Specifically, antibodies in the F3 group were found to play a key role in enabling this broad neutralization, highlighting their importance in combating JN.1 and related lineages [63]. These findings provide valuable insights into how the immune system responds to emerging SARS-CoV-2 variants and how specific antibody groups contribute to protection.

Recent high-throughput studies have proposed a strategy for developing broadly neutralizing antibodies (Abs) against the latest JN.1 variants by targeting non-immunodominant epitopes on the RBD and conserved regions critical for viral function. This approach minimizes the virus’s ability to evolve escape mutations, making it harder for new variants to evade immunity [64]. Several promising antibodies from the E1 group (BD55-3546, BD55-3152, BD55-5585, BD55-5549, and BD55-5840, also known as SA58) and the F3 group (BD55-4637, BD55-3372, BD55-5483, and BD55-5514, also known as SA55) were discovered using this approach [64]. Notably, SA55 demonstrated exceptional neutralizing ability, effectively countering escape mutations found in variants like BQ.1, BQ.1.1, and XBB [64]. These findings, along with other recent studies, highlight how the antibody response to SARS-CoV-2 has evolved in response to the Omicron variant’s antigenic shifts, particularly from XBB to JN.1. Antibodies in groups F3, A1, B, and D3 remain highly effective against JN.1 subvariants, retaining their neutralizing potency. In contrast, antibodies in groups A2, D2, D4, and E1/E2.1 are more likely to be evaded by these variants [62,63,64].

Two recently discovered antibodies, CYFN1006-1 and CYFN1006-2, have shown remarkable ability to neutralize all tested SARS-CoV-2 variants, even outperforming the previously identified broadly neutralizing antibody SA55 [65]. These antibodies bind to epitopes on the outer surface of the RBD that overlap with those targeted by other well-known antibodies like LY-CoV1404, REGN10987 and S309. Their consistent performance across variants highlights their potential as powerful tools for combating the virus’s evolution. In another breakthrough, a yeast-display system combined with a machine learning (ML)-guided approach were used to design and screen a large library of antibody variants leading to the discovery of VIR-7229, a class 1 human antibody that targets the receptor-binding motif (RBM) of the spike protein. VIR-7229 has demonstrated potent neutralization against a wide range of SARS-CoV-2 variants, including EG.5, BA.2.86, and JN.1 [66]. Further advancements in high-throughput DMS assays have identified BD55-1205, a class 1/group A1 antibody with exceptional neutralizing capabilities. BD55-1205 exhibits broad and ultra-high affinity binding against all major SARS-CoV-2 variants evaluated, including XBB, BA.2.86, and JN.1-derived subvariants. Importantly, this antibody also presents a high barrier to viral escape, making it a promising candidate for therapeutic development [67]. The growing body of structural, functional, and biophysical studies has significantly expanded our understanding of how antibodies interact with the SARS-CoV-2 spike protein. These studies have cataloged a wide range of RBD escape mutations, revealing distinct patterns of mutational hotspots that allow the virus to evade antibody binding.

Computer simulations have become an indispensable tool in advancing our understanding of the S protein, its interactions with the ACE2 receptor, and its evasion of neutralizing antibodies at the atomic level [68,69,70,71,72]. These simulations offer unprecedented insights into the structural and energetic factors that govern viral-host interactions and the strategies the virus employs to escape immune detection. By leveraging molecular dynamics (MD) simulations and Markov state models (MSM), we have systematically mapped the conformational landscapes of Omicron variants such as XBB.1 and XBB.1.5, as well as their complexes with ACE2 and neutralizing antibodies [73]. These studies have provided a detailed understanding of how the virus’s structure changes during infection and immune evasion. Mutational scanning and quantitative binding analysis of the XBB spike variants with ACE2 and a panel of class 1 antibodies have offered mechanistic insights that align with experimental observations [74,75]. Notably, epistatic interactions - where mutations at one site influence the effects of mutations at another—have been identified among key binding hotspots like Y501, R498, Q493, L455F, and F456L. These interactions enhance the virus’s ability to bind ACE2 more tightly. At the same time, convergent mutations at F456L and F486P have been shown to play a critical role in enabling the virus to resist a wide range of neutralizing antibodies [74,75]. To explore these mechanisms further, we have integrated AlphaFold2-based atomistic structural predictions with ensemble analyses of S protein-ACE2 complexes for dominant Omicron variants, including JN.1, KP.1, KP.2, and KP.3 [76]. This approach has allowed scientists to investigate how convergent evolutionary hotspots balance the need for strong ACE2 binding with the ability to evade immune responses. Key findings include the identification of binding energy hotspots and epistatic interactions involving mutations at L455, F456, and Q493. These interactions help the virus maintain or restore ACE2 binding affinity while simultaneously facilitating immune escape. Previous studies have also revealed that the S protein can function as an allosteric regulatory machine, using its intrinsic flexibility to regulate binding and immune evasion. By combining molecular dynamics (MD) simulations, ensemble-based mutational scanning of protein stability and binding, and perturbation-based network profiling of allosteric interactions, our studies have examined the mechanisms of the S binding with antibodies [80,81,82]. For example, atomistic modeling of the S protein trimer in complex with antibodies like H014, S309, S2M11, and S2E12 has demonstrated that antibodies can induce specific and functionally relevant changes by modulating the allosteric propensities and collective dynamics of the S protein [80]. These studies have shown that antibody-escaping mutations often target structurally adaptable energy hotspots and allosteric effector centers, which control the protein’s functional movements and communication pathways. Importantly, the escape patterns observed for ultrapotent antibodies are not solely determined by localized changes in binding interactions. Instead, they arise from a complex interplay between structural stability, binding strength, and long-range allosteric communication within the S protein [82]. This suggests that the virus’s ability to evade antibodies is influenced by a delicate balance of multiple energetic factors and interaction networks. Experimental and computational studies have further suggested that the cross-neutralization activity of antibodies against Omicron variants is driven by a trade-off between multiple energetic factors leading to the evolving escape hotspots involved in antigenic drift and convergent evolution [83,84]. These insights highlight the sophisticated mechanisms by which SARS-CoV-2 adapts to immune pressure and underscore the importance of targeting conserved regions or disrupting critical interactions in the design of next-generation vaccines and therapeutics.

Convergent evolution has driven the emergence of SARS-CoV-2 lineages that independently acquire mutations at key amino acid residues, including R346, K444, N450, N460, F486, F490, Q493, and S494. This phenomenon reflects the strong selective pressure exerted by immunity from prior infections or vaccinations, pushing the virus to adapt in ways that enhance its ability to infect host cells while evading immune detection. Notably, Omicron sublineages have evolved to exploit charged interactions at hotspots like R346, N440K, and K444, which help balance the virus’s need for strong ACE2 binding with its ability to escape neutralizing antibodies [85]. Studies have highlighted the importance of electrostatic interactions as a critical thermodynamic force in the binding of the S protein to the ACE2 receptor [86,87,88]. Our computational studies of the Omicron variants RBD binding with ACE2 showed the favorable binding contributions provided by RBD residues K378, R403, K424, K440, K444, K460, N477, K478 that are determined by strong electrostatic interactions mediated by lysine residues, which is the result of significant accumulation of positively charged substitutions interacting with the negatively charged ACE2 binding interface [89,90]. Strategies that target conserved electrostatic hotspots or disrupt critical charged interactions may be more resilient to viral evolution. For example, antibodies engineered to bind to regions with minimal electrostatic variability, such as the stem helix of the S2 subdomain of the S protein, have shown broad neutralizing activity against multiple variants [91,92].

Together, experimental and computational studies have provided compelling evidence that the cross-neutralization activity of antibodies against Omicron variants is governed by a complex and delicate balance and trade-off of multiple energetic factors and interaction contributions. These factors are linked to the evolving escape hotspots associated with antigenic drift and convergent evolution, which enable the virus to adapt to immune pressure while retaining its capacity to infect host cells. Moreover, the evolutionary trade-offs that shape the virus’s ability to balance immune evasion and ACE2 binding are complex and often antibody-dependent. While some mutations may enhance immune escape, they could also impose fitness costs that limit the virus’s transmissibility or replication efficiency. Understanding these trade-offs is critical for predicting the evolutionary trajectory of SARS-CoV-2 and designing interventions that can effectively target emerging variants. One of the key challenges lies in disentangling the contributions of individual interactions to the overall binding energy and stability of the S protein-antibody complex. For instance, how do changes in electrostatic interactions at a single residue, such as K444, influence the binding of different antibody classes? How do hydrophobic effects and hydrogen bonding contribute to the stability of the S protein in its open and closed conformations? Addressing these questions requires integrating high-resolution structural data with detailed dynamic and energetic analysis of the S protein binding with diverse groups and classes of antibodies. Although significant progress has been made in understanding the principles of cross-neutralization and immune evasion, the molecular and energetic details that govern these processes often lack quantifiable analysis of interactions and relative contributions of RBD residues.

In this study, we conducted a comprehensive structure-based mutational scanning of the receptor-binding domain (RBD) residues and performed binding free energy computations for the S-RBD complexes with a panel of neutralizing antibodies. These antibodies span four major classification groups—A, B, C, and D—each targeting distinct binding epitopes on the RBD. The simulated RBD-antibody complexes included: Group A: CB6/LY-CoV016 (etesevimab) [93], Group B: AZD8895 [94] and REGN1033 [95], Group C: LY-CoV555 [96], and Group D: CoV2-2130/AZD1061 [94], REGN10987 [95], and LY-CoV1404 [97]. Coarse-grained CABS model [98,99,100] was used for simulations of the cryo-EM structures of the SARS-CoV-2 S complexes with the panel of antibodies. Additionally, for comparison we performed atomistic MD simulations for S-RBD complex with group A LY-CoV016 antibody. Coarse-grained simulations are combined with mutational profiling of the RBD residues in complexes with the antibodies to enable quantitative comparison with the experimental data on antibody escape. To rigorously analyze the binding affinities of these S-Ab complexes, we employed the Molecular Mechanics/Generalized Born Surface Area (MM-GBSA) approach. This method allowed us to compute binding free energies and perform residue-based energy decomposition, providing detailed insights into the contributions of individual RBD residues to antibody binding. Using a combination of rapid mutational scanning with the simplified energy model and accurate binding free energy analysis using MM-GBSA, we examined the contributions of key RBD residues, with a special emphasis on convergent mutation hotspots. These hotspots, which have emerged repeatedly across Omicron sublineages, reflect the selective pressure exerted by immunity from prior infections or vaccinations.

Our findings reveal that while electrostatic interactions often play a central role in mediating antibody resistance, certain convergent mutational hotspots—such as R346 and K444—are also associated with critical binding affinity hotspots for antibodies targeting conserved RBD epitopes. We show that mutations at convergent mutational hotspots R346 and K444 can alter the electrostatic landscape of the RBD, reducing the binding affinity of some antibodies while maintaining or even enhancing the binding of others. This dual role highlights the complex interplay between immune evasion and the preservation of viral fitness. Our findings offer useful insights into the delicate balance between hydrophobic and electrostatic interactions that govern the binding of neutralizing antibodies to the S-RBD. By elucidating the role and contributions of electrostatic interactions for the RBD residues this study provides better understanding of the mechanisms governing antibody resistance. These findings underscore the importance of targeting conserved regions of the RBD to achieve broad and durable protection against evolving SARS-CoV-2 variants. Our findings underscored on atomistic level the role of targeting conserved epitopes—regions of the RBD that remain relatively unchanged across variants—to enhance the durability and effectiveness of antibodies. By avoiding interfaces with positively charged RBD residues, which are often subject to mutational changes, broadly neutralizing antibodies can maintain their neutralizing capacity even as the virus evolves. This strategy is particularly important in the context of SARS-CoV-2 rapid antigenic drift, which has led to the emergence of numerous variants with distinct immune evasion properties.

2. Results

2.1. Structural Analysis of S-RBD Binding with Four Classes of Antibodies A-D

Our study commenced with a detailed structural analysis of the S-RBD in complex with four distinct groups (A–D) of neutralizing antibodies. These antibodies exhibit diverse binding mechanisms and epitope specificities, which directly influence their neutralization potency and susceptibility to viral escape mutations. Group A antibodies, exemplified by CB6/LY-CoV016 (etesevimab), are characterized by their ability to bind exclusively to the RBD in the ‘up’ conformation, which is essential for their neutralization activity [93]. The binding epitope of LY-CoV016 spans a broad region of the RBD, encompassing residues 403, 405, 408, 409, 415, 417, 420, 421, 453, 455, 456, 457, 458, 459, 460, 473, 474, 475, 476, 477, 487, 489, 490, 493, 494, 498, 500, 502, 503, 504,505 (Figure 1A,D, Table S1). This extensive epitope significantly overlaps with the ACE2 binding interface, enabling the antibody to effectively block viral entry into host cells. The binding of LY-CoV016 to the RBD is mediated by a combination of charged residues, hydrophobic contacts, and salt-bridge interactions. Key residues involved in these interactions include R403, K417, L455, F456, Y473, and N487 (Figure 1A,D, Table S1). Notably, the salt-bridge interaction between K417 of the RBD and a negatively charged residue in the antibody is critical for binding affinity. Interestingly, this interaction is disrupted by the K417N mutation, which is a major escape mechanism employed by the virus to evade neutralization [59,93]. LY-CoV016 utilizes both its heavy and light chains to engage a wide epitope on the RBD (Table S1). However, escape mutations are predominantly driven by changes in RBD residues that interact with the heavy-chain complementarity-determining region (CDR). Key RBD positions where mutations confer escape from LY-CoV016 binding include K417, D420, L455, F456, Y473, A475, N487, and G504 [59]. Our findings underscore that the RBD residues involved in the interfacial contacts with LY-CoV016—but not critical for RBD stability or folding—are likely candidates for escape hotspots. These residues, including K417, D420, L455, F456, Y473, A475, N487, and G504, are frequently mutated in emerging variants to evade antibody-mediated neutralization. This observation emphasizes the importance of targeting conserved regions of the RBD that are less prone to mutational changes.

The neutralizing antibodies in group B incude AZD8895 [94] and REGN1033 [95], that bind to the left shoulder of the RBD (Figure 1B,E). Our structural analysis showd that the binding epitope resides for AZZD8895 include K417, L455, F456, K458, Y473, A475, G476, S477, T478, P479, C480, V483, E484, G485, F486, N487, C488, Y489, and Q493 positions on the RBD (Figure 1B,E, Table S2). Despite a significant number of the epitope positions, G485 and especially F486 and N487 residues are involved in the largest number of interaction contacts with AZD8895 (Figure 1B,E, Table S2). This is consistent with the experimental mutational escape data showeing that group B neutralzing antibodies are very sensitive to the changes at the F486, N487 and G476 sites [59]. However, these major targeting sites of AZD8895 are critically involved in ACE2 binding, and can be less proned to escape. Instead, even though fewer contacts are formed by AZD8895 with S477, T478 and E484 (Table S2).

REGN10933 interacts with RBD residues 403, 406, 417, 421, 453, 455, 456, 473, 475, 476, 477, 478, 484, 485, 486 487, 488, 489, 490, 493, 494, and 498 (Table S3). The largest numbr of contacts are formed with K417, F456, E484, F486, N487 and Y489 (Table S3). The pattern and compositions of interfacial contacts with the binding epitope is similar for group B antibodies AZD8895 and REGN10933. However, we noticed that for REGN10933 there are more contacts with Y453, L455 and Q493 (Table S3). Interestingly, these small differences can have an impact on the escape mutational profile as for REGN10933, these positions emerged among moderate escape hotspots while no escape is provided in these three positions for AZD8895 [59]. Overall, however, the evolution of escape positions tend to minimize mutations in conserved positions such as Y453 but may leverage mutations in L455 and Q493. The latter RBD sites emerged as prominent sites of mutations in the new wave of Omicron variants including JN.1. KP.2 and KP.3 where L455S and Q493E enabled significant changes in the immune escape profile. Interestngly, two Ab cocktails – AZD8895/AZD1061 and REGN10933/REGN10987 consisting of antibodies from group B (AZD8895, REGN10933) and group D (AZD1061, REGN10987) provided complementary binding to RBD and demonstrated significant potential in preventing mutational escape. REGN10933 and REGN10987 have escaped by different mutations as mutations at F486 escaped neutralization only by REGN10933, whereas mutations at K444 escaped neutralization only by REGN10987 [101].

Group C neutralizing antibodies LY-CoV555 (bamlanivimab), exhibit a unique ability to bind to the RBD in both the ‘up’ and ‘down’ conformations, which contributes to their higher neutralization potency compared to antibodies in other groups [59,96]. This dual-conformation binding capability allows Group C antibodies to effectively target the RBD regardless of its structural state, enhancing their ability to block viral entry into host cells. The binding epitope of LY-CoV555 encompasses a broad range of RBD residues, including Y351, Y449, L452, L455, F456, T470, I472, N481, G482, V483, E484, G485, F486, N487, Y489, F490, L492, Q493, and S494 (Figure 1C,F, Table S4). Many of these residues are conserved hydrophobic sites that play critical roles in RBD function and are less prone to mutational changes due to their structural and functional importance. Among these, V483, E484, G485, and F486 form a significant number of contacts with the antibody, making them key contributors to the binding interaction. Experimental studies have shown that these RBD positions, particularly E484, exhibit the highest mutation escape scores, indicating their importance in antibody evasion [59]. For example, the E484K mutation, in Beta variant, and the E484A mutation in many Omicron sublineages, are well-documented escape mechanisms that significantly reduce the binding affinity and neutralization efficacy of Group C antibodies. These mutations alter the electrostatic and structural properties of the RBD, enabling the virus to evade immune detection. The susceptibility to escape mutations, particularly at key residues like E484, emphasizes the need for targeting conserved epitopes and designing antibodies that can withstand viral evolution.

Group D antibodies included CoV2-2130/AZD1061 (cilgavimab) [94], REGN10987 (imdevimab) [95], and LY-CoV1404 [97] exhibit unique binding mechanisms and epitope specificities that contribute to their neutralization potency and therapeutic potential. The distinct orientation allows group D antibodies target a specific loop formed by RBD residues 440–449, which is critical for their binding and neutralization activity (Figure 2) The binding epitope of AZD1061 includes residues T345, R346, N439, N440, L441, D442, S443, K444, V445, G446, G447, Y449, L450, L452, E484, F490, L492, Q493, S494, and P499 (Figure 2A,C, Table S5). Among these, G446S has been shown to substantially impair the neutralization capability of AZD1061 against Omicron variants [59]. REGN10987 forms contacts with RBD positions 346, 439, 440, 441, 443, 444, 445, 446, 447, 448, 449, 450, 498, 499, 500, and 501 (Table S6). The largest number of contacts are formed by K444, V445, Y449 and P499 (Table S5). LY-CoV1404 (bebtelovimab) targets the same binding epitope but forms contacts with an extended set of residues, including T345, R346, N439, N440, L441, D442, S443, K444, V445, G446, G447, N448, Y449, N450, Q498, P499, T500, N501, G502, V503, and Q506 (Figure 2B,D, Table S7). This antibody interacts with the entire stretch of residues 439–450 and also engages key ACE2 binding interface sites such as T500 and N501, as well as adjacent residues 503–506. LY-CoV1404 binds to a region that overlaps with the ACE2-interacting site of the S protein, making it accessible in both the ‘up’ and ‘down’ conformations of the RBD. Although LY-CoV1404 is classified as a group D or class 2 binder, the structural location of its epitope is closer to that of the canonical class 3 binder, S309. The binding epitope of LY-CoV1404 is similar to that of imdevimab (REGN10987), but it has a significantly larger contact surface area on the RBD (584 Ų compared to 343 Ų for REGN10987) [95,97].

The electrostatic Interactions are formed by R346 and K444 with complementary residues in the antibody, contributing to binding stability, while the hydrophobic contacts are provided by L441, V445 and P499. The salt-bridge interaction between K444 of the RBD and a negatively charged residue in the antibody is particularly important for binding affinity. However, this interaction can be disrupted by mutations such as K444N, which is a known escape mechanism employed by the virus. REGN10987 is susceptible to escape mutations at K444, V445, and G446 that are critical hotspots for escape mutations [59]. Structural binding mechanism for this group of antibodies involves a network of electrostatic, hydrophobic, and hydrogen-bonding interactions that contribute to its high binding affinity and neutralization potency.

To summarize, a comparative structural analysis of binding for the four groups of antibodies offers several important observations. First, group A antibodies bind exclusively to the ‘up’ conformation, while Groups B, C, and D can bind to both ‘up’ and ‘down’ conformations, providing broader neutralization potential. Second concerns epitope specificity. While group A and B antibodies target overlapping epitopes with the ACE2-binding site, group C antibodies focus on conserved hydrophobic regions, and group D antibodies target the 440–449 loop. All groups rely on a combination of electrostatic, hydrophobic, and hydrogen-bonding interactions, but the key specific residues vary. For example, group A antibodies depend heavily on K417, while group D antibodies rely on K444 and G446. Although all four groups are susceptible to escape mutations, group A is particularly sensitive to K417N, while group C to E484K, and group D to K444N. Groups A, B, C, and D antibodies exhibit distinct binding mechanisms and epitope specificities, their susceptibility to escape mutations highlights the challenges of key factors and contributions of individual residues to the evolving landscape of immune escape.

2.2. CG-CABS Simulations and Collective Dynamics Reveal Role of Hinge Sites as Positions of Antibody Escape

We performed CG-CABS and all-atom MD simulations for CB6/LY-CoV016 (Supporting Information, Figure S1). CABS trajectories were subjected to atomistic reconstruction and refinement thereby allowing for a direct comparative analysis with MD simulations performed in an explicit solvent environment. Using a comparison of CG-CABS and MD simulations, we verified the reliability of the proposed simplified and less time-consuming simulation model. A comparison of the RMSF profiles obtained from CABS-based reconstructed trajectories and all-atom MD simulations showed a highly similar conformational dynamics profile for the RBD residues LY-CoV016 (Supporting Information, Figure S1). CG-CABS trajectories featured more stable profiles of the S-RBD complexes, while MD simulations highlighted a more dynamic nature of the RBD regions. Both simulation models reproduced stability of the conserved core of the RBD antiparallel β strands (β1 to β4 and β7) (residues 354-358, 376-380, 394-403, 431-438, 507-516) and a particularly significant stabilization of β-sheets (β5 and β6) (residues 451-454 and 492-495) that anchor the RBM region to the central core (Supporting Information, Figure S1). The observed similarities of the dynamic profiles suggested that CG-CABS simulations accompanied by atomistic reconstruction can be used to provide accurate and affordable simulation approach for quantifying flexibility of the S-RBD complexes with the panel of antibodies.

We also characterized collective motions and determined the hinge regions in the SARS-CoV-2 S-RBD complexes using principal component analysis (PCA) of trajectories derived from CG-CABS simulations using the CARMA package [102]. The local minima along these profiles are typically aligned with the immobilized in global motions hinge centers, while the maxima correspond to the moving regions undergoing concerted movements leading to global changes in structure [77,80,103,104]. The low-frequency ‘soft modes’ are characterized by their cooperativity and there is a strong relationship between conformational changes and the ‘soft’ modes of motions intrinsically accessible to protein structures [77,80,103,104]. Hinge sites are regions within a protein that enable the protein to undergo conformational changes necessary for its biological functions. When hinge positions in spike-antibody complexes correspond to stable sites, several important effects can occur. Stable hinge positions may limit the flexibility of the spike protein and hinder the protein’s ability to undergo necessary conformational changes, potentially affecting its function in viral entry and immune evasion. Stability at hinge positions can influence how antibodies bind to the spike protein. If these positions are too rigid, it might reduce the effectiveness of antibody binding, as the protein may not adopt the optimal conformation for interaction with the antibody. These sites are characterized by a balance between flexibility and rigidity, which is essential for maintaining structural integrity while allowing functional dynamics. In the context of the SARS-CoV-2 spike protein and its interactions with neutralizing antibodies, we will examine a pivotal role of hinge positions in coordinating the dynamics of antibody binding and viral escape.

For LY-CoV016 Ab, the hinge positions are distributed along the RBD interface and among these functionally important positions are D405, E406, T415, K417, D420, Y421, L455, F456, A475, Q493 and N501 residues (Figure 3A). Only a small fraction of these sites corresponded to the binding energy hotspots (Y421, Y453, L455, and F456). It should be pointed out that mutational sensitivity maps for CB6 binding yielded only minor binding energy changes upon substitutions in D405, E406, K417, D420, A475, and N501 residues. The major experimentally determined escape sites are K417, D420, N460, Y473, and A475 [59] and these sites correspond to flexible hinge sites that control dynamics of antibody binding. The results highlight the importance of hinge positions in coordinating the dynamics of antibody binding and viral escape. We found that key antibody-escaping mutations can target these regulatory hinge sites, which control collective motions and allosteric interactions within the RBD. While only a small fraction of these sites correspond to binding energy hotspots, mutations at flexible hinge positions can significantly alter the effectiveness of antibody binding by disrupting conformational dynamics and reducing the antibody’s ability to bind effectively. Although mutations in these sites are not lethal for binding affinity, they can alter effectiveness of antibody as the protein may not adopt the optimal conformation for interaction with the antibody. Some mutations at hinge positions can lead to the occlusion of epitopes, which are the specific parts of the antigen that antibodies recognize and bind to. This occlusion can prevent antibodies from binding effectively, allowing the virus to escape immune detection. The results showed that key antibody-escaping mutations may often target regulatory hinge positions that coordinate collective motions and allosteric interactions.

The slow mode profile for group B AZD8895 complexed with RBD revealed a broad and deep minimum associated with RBD positions 473-490 (Figure 3B) as this stretch of residues is involved in the S-RBD binding. By projecting the sites of immune escape S477, T478, E484, G485, F486 and N487 on the slow mode profile, it appears that these positions are aligned wih the wider hinge region. The binding epitope of group C LY-CoV555 encompasses a broad range of RBD residues, and several conserved hydrophobic sites including V483, G485, and F486 along with E484 position are associated with the escape mutations. Experimental studies have shown that these RBD positions, particularly E484, exhibit the highest mutation escape scores, indicating their importance in antibody evasion [59]. The slow mode profile for LY-CoV555-RBD complex featured local minima and hinge site near RBD residues 483-486 (Figure 3C). Here again we observed that the major antibody escape center is functionally important not only for modulating binding affinity but also global movements of the antibody around the RBD. Structural maps of the slow mode profiles illustrated the overall pattern of formation hinge sites of the RBD-antibody complexes for groups A,B and C (Supporting Information, Figure S1).

Of particular interest is a comparison of slow mode profiles for group D antibodies that displayed a common deep minimum around positions K444, V445 and G446 as well as second minimum around T500 (Figure 3D-F). Although these antibodies commonly target RBD residues 440–449, which is critical for their binding and neutralization activity, group D antibodies form favorable interactions with other RBD sites. Despite important role of 440-449 for immune escape, several sites are particularly important with many escape mutations for K444, V445, and G446 [59]. Strikingly, K444 and V445 correspond to the exact hinge positions in the complexes (Figure 3D-F), thus pointing to the important role of these residues as modulators of the RBD-antibody movements and dynamics. Structural maps of the slow mode profiles for group D antibody complexes with RBD showed the increased rigidity around residues 444-446 along essential mobility modes (Supporting Information, Figure 2). The binding mode for group D antibodies appeared to enable significant mobility of the RBM residues in functional movements that would allow for modulation of the RBD flexibility and exposure of the RBD to ACE2 interactions (Supporting Information, Figure S2). While Group D antibodies do not directly target all ACE interacting residues, their binding epitopes partially overlap with the ACE2-binding site, particularly at residues K444, V445, and Y449. We could see that K444 and V445 are critical for both antibody binding and ACE2 interactions. Mutations at these sites (e.g., K444Q, V445A) can disrupt antibody binding while also affecting ACE2 affinity, though to a lesser extent. Y449 residue is involved in hydrophobic interactions with both ACE2 and Group D antibodies, making it a key site for competitive binding. AZD1061 and REGN10987 can partially block ACE2 binding and their competitive potential is moderate compared to antibodies that directly target the ACE2-binding interface. LY-CoV1404 has a stronger competitive potential due to its engagement with key ACE2-binding residues T500 and N501. By binding to these critical sites and forming two primary hinge positions ( 444-446 and 500-501) LY-CoV1404 can effectively block ACE2 binding and also reduce global functional movements of the RBD thus limiting the dynamic potential of the RBD to optimize interactions with the ACE2 receptor.

Hinge regions are characterized by their flexibility, which allows the RBD to transition between different conformational states. This flexibility is essential for the spike protein’s function, such as binding to the ACE2 receptor or adopting the optimal conformation for antibody interaction. However, this same flexibility makes hinge regions vulnerable to mutations that can alter the protein’s dynamics and disrupt antibody binding. The relatively flexible nature of K444 and V445 in mediating both hydrophobic and electrostatic contacts make them ideal for immune escape as mutations at hinge positions can cause the RBD to adopt alternative conformations that occlude epitopes, preventing antibodies from binding effectively. Mutations in these hinge sites can disrupt the collective motions and allosteric interactions of the RBD, reducing the antibody’s ability to neutralize the virus. Moreover, mutations at hinge positions can perturb antibody binding without compromising the virus’s ability to infect host cells, providing an evolutionary advantage. Since the virus is under constant evolutionary pressure to evade immune detection while maintaining its ability to infect host cells, hinge region can provide an ideal target for escape mutations because they allow the virus to balance these competing demands. The flexibility of hinge regions enables the virus to rapidly adapt to immune pressure by evolving mutations that disrupt antibody binding.

2.3. Mutational Profiling of Protein-Antibody Binding Interfaces

To provide a systematic comparison, we constructed mutational heatmaps for the RBD interface residues of the S complexes with the antibodies. We began by analyzing the mutational heatmap for LY-CoV016, a representative group A antibody. The heatmap revealed several strong binding hotspots where mutations cause significant destabilization of the antibody-RBD complex (Figure 4A). Y421 is surrounded by Y33, G54, G55, and S53 of the antibody’s heavy chain, forming a network of interactions critical for binding. Y453 makes a moderate number of contacts with the antibody and plays a key role in stabilizing the complex. L455 forms hydrophobic interactions with the antibody, and F456 makes a large number of contacts with M101, V98, N32, S31, S53, P100, L99, and Y33 of the heavy chain, making it a dominant binding hotspot (Figure 4A). Y489 also contributes to binding, though to a lesser extent than the other hotspots. The mutational heatmap analysis highlights the trade-off between binding affinity and immune evasion. While residues like Y421, Y453, L455, and F456 are critical for antibody binding, they are less likely to mutate due to their role in RBD stability. Instead, the virus targets more flexible residues that are less critical for ACE2 binding but sensitive to electrostatic interactions. In particular residues like K417, D420, and N460 are more flexible and sensitive to electrostatic interactions, making them prime targets for escape mutations [59].

Mutational heatmap reproduced a general trend in escaping hotspots and mutations, particularly highlighting K417N mutation that disrupts a critical salt-bridge interaction, significantly reducing binding affinity (Figure 4A). Our analysis indicated that L455S, which alters hydrophobic contacts can weaken antibody binding as well as F486V that induced considerable binding loss leading to escape. In this context, N460K combined with F486V diminishes neutralization efficacy for the BA.5 variant [59,60,61,62]. A similar pattern was seen in the mutational heatmap for group B AZD8895 where the main binding hotspots are F486, Y489, F456 and to a lesser extent G476 and N487 where for the latter large destabilization can occur upon N487E and N487K mutations (Figure 4B). For LY-CoV555, the largest destabilization changes are associated with mutations of F486, including F486D ( ΔΔG = 5.01 kcal/mol), F4856E ( ΔΔG = 4.61 kcal/mol), F486K ( ΔΔG = 4.4 kcal/mol) and F486N ( ΔΔG = 4.1 kcal/mol) (Figure 4B). Additionally, there are highly destabilizing mutations F456P (ΔΔG = 2.42 kcal/mol), N487E (ΔΔG = 2.28 kcal/mol). The experimental mutational escape data showed that group B neutralzing antibodies are very sensitive to the changes at the F486, N487 and G476 sites [59]. However, these major targeting sites of AZD8895 are critically involved in ACE2 binding, and are less proned to escape. S477, T478 and E484 are found to tolerant to mutational changes leading to moderate changes of ΔΔG < 1.0 kcal/mol but some mutations S477N ( ΔΔG = 1.19 kcal/mol), T478K ( ΔΔG =1.3 kcal/mol) and E484A (ΔΔG = 1.14 kcal/mol) are more deatbilizing and are less critical for ACE2 binding. Indeed, mutations in these posiitions S477N/T478K/E484A can result in a significant escape from goup B antibodies with the Omicron variant owing likely to the altetered pattern of electrostatic interactions [59].

Similar mutational heatmap was obtained for another group B antibody REGN1093 (Supporting Information, Figure S3A). The key binding hotspots are K417, Y421, Y453, L455, F456, F486, Y489 and Q493 where mutations in F456 and F486 are especially destabilizing. These results showed that REGN10933 has an expanded range of strong binding hotspots as compared to AZD8895 which could arguably produce a greater repertoire of escape hotspots (Figure 3B, Supporting Information, Figure S3A). This is exactly what was found in the experimental studies [59] in which the escape profile for REGN10933 included K417, Y453, L55, Y472, T476, F486, N487, Y489 and Q493 positions. We specifically highlighted the results of mutational scanning for Omicron and particularly JN.1/KP.3 mutational sites (D339H, K356T, R403K, K417N, V445H, G446, N450D, L452W, L455S, F456L, N460K, N481K, E484K, F486P, Q493E) (Supporting Information, Figure S3B). It can be seen that the largest destabilization changes occur upon mutations L455S, F456L, F486P, Q493E. These results are consistent with experimental data showing that the newly emerging Omicron variants become highly resistant to group antibodies including REGN10933 [59]. While REGN10933 has potent neutralization against earlier variants, emerging variants like JN.1 and KP.3 may exhibit reduced susceptibility to REGN10933 due to specific mutations in the RBD. Indeed, JN.1 a descendant of the BA.2.86 lineage and carries additional mutations, including L455S, F456L, and Q493E. KP.3 is a sublineage of JN.1 and includes additional mutations such as R346T, F456L, and Q493E. Our results indicated that mutations L455S, F456L and Q493E can disrupt antibody binding, suggesting that JN.1 and KP.3 may exhibit reduced susceptibility to REGN10933.

LY-CoV555 (bamlanivimab) is a group C neutralizing antibody and mutational scanning of the RBD revealed that certain residues are particularly sensitive to mutations, leading to significant destabilization of the antibody-RBD complex (Figure 4C). The largest destabilizing changes were observed for F490D (ΔΔG = 3.4 kcal/mol), F490E (ΔΔG = 3.12 kcal/mol) and F490N (ΔΔG = 2.93 kcal/mol) where mutations introduce polar or charged residues, disrupting the hydrophobic interactions critical for binding. Other large destabilizing mutations included V483S (ΔΔG = 1.94 kcal/mol) that replaces a hydrophobic valine with a polar serine, weakening the binding interface as well as Q493S (ΔΔG = 1.89 kcal/mol) that disrupts hydrogen-bonding interactions, reducing binding affinity. Mutations at E484 such as E484D (ΔΔG = 1.81 kcal/mol), E484A (ΔΔG = 2.02 kcal/mol) and E484N (ΔΔG = 1.73 kcal/mol) alter the electrostatic landscape of the RBD, reducing the antibody’s binding affinity (Figure 4C). The mutational scanning results align with experimental data showing that residues V483, E484, F486, and F490 are highly susceptible to escape mutations [59]. These positions exhibited the highest mutation escape scores, making them prime targets for viral adaptation. Although V483, F486, and F490 are conserved hydrophobic sites critical for RBD stability and ACE2 binding, the virus can still employ some mutations at these positions to evade antibodies. However, such mutations often come at a cost to viral fitness, as they can destabilize the RBD or reduce ACE2 binding affinity. While residues like V483, F486, and F490 are critical for antibody binding, their role in RBD stability and ACE2 binding limits the virus’s ability to mutate these sites without compromising its infectivity. As a result, the virus tends to target second-tier binding hotspots, such as E484, which are less detrimental to ACE2 binding but still induce significant loss in antibody binding affinity. Understanding the binding escape mechanisms of LY-CoV555 provides valuable insights for therapeutic design, suggesting that antibodies engineered to target conserved hydrophobic residues like V483, F486, and F490 may be more resilient to viral evolution, as these sites are less prone to mutational changes.

For AZD1061 (cilgavimab), another group D antibody, the primary binding affinity hotspots are located at residues V445, Y449, and F490 (Figure 4D). These residues form critical interactions with the antibody, contributing to its high binding affinity and neutralization potency. Specifically, V445 and F490 are involved in hydrophobic interactions, while Y449 participates in both hydrophobic contacts and hydrogen bonding, stabilizing the antibody-RBD complex. In contrast, residues such as K444, G446, G447, and P499 play a somewhat less significant role in binding, though they still contribute to the overall interaction network. Despite their secondary role, K444 and G447 exhibit the highest experimental escape scores [59], making them dominant sites for viral escape mutations. This pattern likely arises because V445, Y449, and F490 are not only critical for antibody binding but also essential for RBD stability and ACE2 binding. As a result, mutations at these sites could compromise viral fitness, making them less prone to evolutionary changes. In contrast, residues like K444 and G447 are more flexible and less critical for RBD function, allowing the virus to exploit these positions for immune evasion without significantly affecting its ability to infect host cells.

For another group D antibody REGN10987 (imdevimab), the most significant destabilization of the antibody-RBD complex is caused by mutations at V445, particularly V445A (ΔΔG = 2.5 kcal/mol), V445D (ΔΔG = 2.65 kcal/mol), and V445E (ΔΔG = 2.26 kcal/mol) (Figure 4E, Supporting Information Figure S3C,D). Additionally, large destabilization changes are observed for mutations at K444, such as K444Q and K444E, as well as for mutations at P499 (Figure 4E). These residues are part of the major escape sites for REGN10987, which are dominated by K444, V445, G446, and G447 [59]. We also analyzed binding free energy changes associated with BA.2.86, JN.1, KP.2 and KP.3 mutations for REGN10987 binding (Supporting Information Figure S3D). The binding free energy changes showed significant loss of binding upon N440K, V445H, G446S mutations along with somewhat smaller but still appreciable change upon L455S mutation. These key mutational changes present in JN.1, KP.2 and KP.3 variants, particularly N440K and V445H may induce enhanced immune escape from REGN10987. Interestingly, REGN10933 and REGN10987 exhibit distinct escape pathways. Mutations at F486 specifically escape neutralization by REGN10933, while mutations at K444 escape neutralization by REGN10987. Hence, while L455S, F456L, F486P, Q493E mutations are detrimental for REGN10933, mutations N440K, K444Q, V445H can cause escape from REGN10987. Mutations escaping both antibodies include combination of F486V and K444Q, E484K + K444Q and L455S + F456L: These mutations, observed in the JN.1 and KP.3 variants, can disrupt binding to both antibodies by altering the structural conformation of the RBD. The experimental evidence showed that XBB variant that carries mutations V445P and F486S may enhance immune escape from REGN10933 and REGN10987. Similarly, our energetic analysis is consistent with the evidence that JN.1 and KP.3 variants that carry L455S, F456L and Q493E mutations can collectively reduce the binding affinity of both antibodies.

LY-CoV1404 (bebtelovimab) is a potent group D neutralizing antibody and mutational scanning showed that the binding energy hotspots for LY-CoV1404 are centered around residues V445, P499, T500, and K444. These residues play a crucial role in stabilizing the antibody-RBD complex through hydrophobic interactions, electrostatic contacts, and hydrogen bonding. The largest destabilizing mutations are observed for V445 including V445E (ΔΔG = 4.24 kcal/mol), V445D (ΔΔG = 4.09 kcal/mol) and V445S (ΔΔG = 3.17 kcal/mol) that disrupt the favorable hydrophobic interactions (Figure 4F). Other hydrophobic hotspots are P499 and T500 (Figure 3D). Of notice are mutations in K444 that engages in electrostatic interactions, with mutations K444D (ΔΔG = 1.93 kcal/mol) and K444E (ΔΔG = 1.64 kcal/mol) significantly reducing binding affinity (Figure 4F). The experiments showed that amino acid substitutions at K444, V445, and G446, as well as some substitutions at P499 and T500, can lead to escape from LY-CoV1404 neutralization [59]. The most important mutations such as V445P found in the XBB variant, causes significant destabilization (ΔΔG = 3.69 kcal/mol) and is a strong escape modification K444T present in the BQ.1 variant, also leads to significant destabilization (ΔΔG = 1.78 Mutations like K444D and K444E disrupt electrostatic interactions, reducing binding affinity without compromising RBD stability or ACE2 binding. Interestingly, not all mutations at these sites lead to escape. For example, V445I (ΔΔG = 1.08 kcal/mol), V445L (ΔΔG = 0.64 kcal/mol) and G446D (ΔΔG = 0.5 kcal/mol) mutations only moderately affect antibody binding and do not result in immune escape, highlighting the nuanced nature of viral adaptation (Figure 4F). The analysis further highlights the trade-off between antibody binding and viral escape. While residues like V445, P499, and T500 are critical for antibody binding, they are also important for RBD stability and ACE2 binding, limiting the virus’s ability to mutate these sites. In contrast, charged residue K444 is less critical for ACE2 binding and prone to escape mutations due to flexibility and appreciable contribution to the antibody binding.

These findings underscore the importance of understanding the specific escape mechanisms for individual antibodies and the potential of combination therapies to enhance neutralization breadth and durability. By targeting multiple epitopes, combination therapies like REGN-COV2 can mitigate the risk of escape mutations and provide more robust protection against emerging variants. All three antibodies target V445, highlighting its importance in binding. LY-CoV1404 and REGN10987 also target K444, while AZD1061 emphasizes Y449 and F490. LY-CoV1404 and REGN10987 share similar escape pathways, with mutations at K444 and V445 being dominant, while AZD1061 is more susceptible to mutations at K444 and G447, which are less critical for RBD stability. The comparative analysis of LY-CoV1404, REGN10987, and AZD1061 reveals both shared and distinct binding and escape mechanisms. While all three antibodies target critical residues like V445 and K444, their specific binding hotspots and escape pathways differ. These differences highlight the importance of combination therapies and engineered antibodies in addressing the challenges posed by viral evolution.

2.4. MM-GBSA Analysis of the Binding Affinities

We utilized conformational ensembles derived from CG-CABS simulations to compute the binding free energies of RBD-antibody complexes using the MM-GBSA (Molecular Mechanics/Generalized Born Surface Area) method. This approach allowed us to perform a detailed binding free energy decomposition, examining the energetic contributions of individual RBD epitope residues. Through this analysis, we identified the binding hotspots that are critical for antibody binding and quantified the roles of van der Waals interactions and electrostatic interactions in the binding mechanism. In the MM-GBSA calculations, we investigated whether the binding affinities and contributions of the major binding hotspots are primarily determined by van der Waals interactions or electrostatic interactions. Additionally, we explored whether positions of immune escape could be associated with binding hotspots where different energetic contributions act synergistically, leading to significant loss of binding upon mutations. Through residue-based binding free energy decomposition, we identified the key residues contributing to the antibody’s binding affinity and explored the mechanisms underlying susceptibility to escape mutations. For group A LY-CoV016 antibody the residue decomposition analysis revealed that the most favorable total binding energies are associated with the following residues: N487 (ΔG = -4.71 kcal/mol), A475 (ΔG = -4.43 kcal/mol), F456 (ΔG = -3.93 kcal/mol), R403 (ΔG = -3.42 kcal/mol), N460 (ΔG = -3.13 kcal/mol) and K417 (ΔG = -2.83 kcal/mol) (Figure 5A). The most favorable van der Waals interactions are formed by F456, L455, Y421, Y489, and K417. These residues participate in hydrophobic contacts, which are essential for stabilizing the antibody-RBD complex (Figure 5B). The largest electrostatic contributions come from K417, followed by D420, E484, and R403. These residues form charged interactions and hydrogen bonds with the antibody, enhancing binding affinity (Figure 5C). Interestingly, K417 exhibits both favorable van der Waals and electrostatic interactions, which can act synergistically to contribute to the antibody’s binding affinity. While the total binding energies are more favorable for residues like N487 and F456, the electrostatic interactions at K417 play a critical role in determining the antibody’s susceptibility to escape mutations.

The experimental escape scores are consistent with our computational findings, identifying K417, N460, A475, D420, F456, and N487 as the dominant escape hotspots [59]. These residues are frequently mutated in emerging variants, enabling the virus to evade neutralization. The binding and escape mechanisms of LY-CoV016 highlight the importance of electrostatic interactions and hydrophobic contacts in determining antibody binding affinity and susceptibility to escape mutations. While residues like N487 and F456 contribute significantly to binding, K417 plays a dual role, with both van der Waals and electrostatic interactions driving its importance as an escape hotspot.

For AZD8895, a representative group B antibody, the largest total binding energies for AZD8895 are associated with F486 (ΔG = -10.0 kcal/mol), N487 (ΔG = -5.07 kcal/mol), T478 (ΔG = -3.99 kcal/mol), G476 (ΔG = -3.93 kcal/mol), S477 (ΔG = -2.7 kcal/mol and F456 (ΔG = -2.04 kcal/mol) (Figure 5D). The most favorable van der Waals interactions are formed by F486, Y489, F456, T478, and N487. These residues participate in hydrophobic contacts, which are essential for stabilizing the antibody-RBD complex (Figure 5E). The largest electrostatic contributions come from E484, D405, E406, D429, and D442 (Figure 5F). However, the favorable electrostatic interactions at these sites are largely offset by desolvation penalties, resulting in moderate total binding energies. The results showed that F486 and N487 have the strongest binding energies but are less prone to mutations due to their critical role in RBD function and ACE2 binding. At the same time, second-tier binding hotspots T478 (ΔG = -3.99 kcal/mol) and S477 (ΔG = -2.7 kcal/mol) have favorable binding energies due to both van der Waals and electrostatic interactions (Figure 5D-F). While their contributions are weaker than those of F486 and N487, they are more flexible and less critical for RBD stability, making them prime targets for escape mutations.

Our results suggest that the virus can exploit multiple second-tier binding hotspots, such as T478 and S477, to evolve mutations that collectively reduce antibody binding. This strategy allows the virus to evade neutralization without compromising its ability to bind ACE2 [59]. The binding and escape mechanisms of AZD8895 highlight the importance of hydrophobic interactions and electrostatic contributions in determining antibody binding affinity and susceptibility to escape mutations. While residues like F486 and N487 are critical for binding, their role in RBD stability and ACE2 binding limits the virus’s ability to mutate these sites. Instead, the virus targets second-tier binding hotspots, such as T478 and S477, to evolve escape mutations that collectively reduce antibody binding.

The binding energy analysis revealed several critical hotspots for group C antibody LY-CoV555, including a cluster of nearby residues: E484, S494, F486, Q493, V483, F490, and G485, as well as a stand-alone peak for Y449 (Figure 6A). E484 position corresponds to the deepest and most significant binding energy peak, contributing ΔG = -12.03 kcal/mol. E484 engages in both electrostatic interactions (ΔG = -3.95 kcal/mol) and van der Waals interactions, making it a dominant binding hotspot. S494 residue contributes significantly to binding, with ΔG = -7.64 kcal/mol, though less than E484. The residue decomposition showed that F486, V483, Y449, G485, and F490 contribute favorably through van der Waals interactions, stabilizing the antibody-RBD complex (Figure 6B), while the strongest electrostatic interactions are observed at E484, followed by D442, E406, D405, and D420 (Figure 6C). The binding energy hotspots identified in our computational analysis are in line with experimental escape profiles. Indeed, E484 is the dominant escape hotspot, with mutations like E484K and E484A significantly reducing binding affinity and neutralization efficacy. V483, F486, F490, Q493, and S494 also serve as escape positions, though to a lesser extent than E484 (Figure 6A-C).

These residues are predicted as secondary binding hotspots, contributing to the antibody’s overall binding affinity. The results indicate that the predicted binding hotspots are strong indicators of potential escape centers for LY-CoV555. This is particularly evident for E484, which not only contributes the most favorable binding energy (Figure 6A-C) but also exhibits the highest escape mutation scores. While residues like E484 dominate the binding landscape, secondary hotspots such as V483, F486, F490, Q493, and S494 also play critical roles in both binding and immune evasion (Figure 6A-C). The virus exploits these hotspots to evolve mutations that reduce antibody binding while maintaining its ability to infect host cells.

MM-GBSA computations for group D antibodies AZD1061 (cilgavimab) (Figure 6D-F), REGN10987 (imdevimab) (Figure 7A-C) and LY-CoV1404 (bebtelovimab) (Figure 7D-F) provided interesting information regarding common and unique binding hotspots and corresponding escape profiles. For AZD1061 (Cilgavimab) the key energy hotspots are R346 (ΔG = -9.23 kcal/mol) which is the most pronounced binding hotspot, contributing significantly to binding affinity (Figure 6D). Other hotspots of binding are K444 (ΔG = -5.94 kcal/mol) with strong van der Waals (ΔG = -4.55 kcal/mol) (Figure 6E) and electrostatic interactions (Figure 6F). Y449 (ΔG = -5.15 kcal/mol): and V445 (ΔG = -3.57 kcal/mol) contribute mostly through hydrophobic interactions (Figure 6E).

For REGN10987 (Imdevimab) the identified hotspots are N440, V445, T500, Y449, N448, P499, L441, and G446 (Figure 7A) with the strongest Van der Waals interactions for V445, N440, T500, P499, L441, and Y449 (Figure 7B). The electrostatic Interactions favor D364, E484, E340, D405, E406, D398, E471, D40 where the residue distribution showed many RBD residues with favorable contributions, while man other RBD sites displayed highly unfavorable electrostatic contacts (Figure 7C). Overall, the electrostatic component is quite noisy for REGN10987 (Figure 7C)and most of the binding hotspots (N440, V445, N448) revealed complex balance and synergistic effects of hydrophobic and electrostatic contributions.

For LY-CoV1404 (Bebtelovimab) K444 (ΔG = -8.2 kcal/mol) and V445 (ΔG = -5.9 kcal/mol) are the strongest binding hotspots (Figure 7D), with significant contributions from both van der Waals and electrostatic interactions (Figure 7E,F). T500, N440, P499, G447, N450, and L441 also contribute to binding, though to a lesser extent. A comparative analysis of energy contributions for group D antibodies showed that AZD1061 and LY-CoV1404 rely heavily on K444 and V445 for binding, with strong contributions from both van der Waals and electrostatic interactions (Figure 6D-F, 7D-F), whereas REGN10987 exhibits a more balanced binding profile, with significant contributions from multiple residues, including N440, K444, V445, and P499 (Figure 7A-C). These results agree with and explain the experimental escape pathways where K444 is a dominant escape hotspot for all three antibodies, with mutations like K444Q and K444T significantly reducing binding affinity. V445 and G446 are also critical escape sites, particularly for LY-CoV1404 and AZD1061. LY-CoV1404 binds to a highly conserved region of the RBD and the antibody forms robust interactions with key residues, K444 and V445, which are critical for binding but less likely to mutate due to their role in RBD stability and ACE2 binding. As a result, LY-CoV1404 is less susceptible to escape mutations that have nullified the activity of other antibodies. For example, mutations like K444Q and V445A only moderately affect binding, whereas they significantly reduce the efficacy of other antibodies.

To summarize, MM-GBSA binding free energy analysis provided detailed insights int the binding mechanisms, energetic contributions, and escape pathways of these antibodies. Some interesting generic observations pointed to the energetic contributions, where van der Waals interactions dominate binding for antibodies driven by F456, F486, and V445 while electrostatic interactions mediated by K417, E484, and K444 are the main determining factors of binding for antibodies. We also noticed that escape hotspots and escape mutations can often but not always emerge in positions associated with the strongest binding affinity centers. While the predicted binding hotspots may be strong indicators of potential escape centers, this mechanism can be adapted to exploit second-tier binding hotspots that are less critical for vital RBD functions such as RBD stability and ACE2 binding. While each group exhibits unique binding hotspots and escape vulnerabilities, common themes emerge, such as the role of key residues and synergistic effects of van der Walls and electrostatic interactions for K417, F486, K444 and Q493 in driving binding mechanisms and evolutionary escape paths. In this context, it is striking to realize that convergent evolution mutations tend to emerge in positions R346, K444, N450, N460, F486, F490, Q493, and S494 where mutations can modulate binding with different groups of antibodies without dramatically compromising ACE2 affinity.

3. Discussion

The results of the MM-GBSA binding free energy analysis, combined with residue-based decomposition, provide a comprehensive understanding of the binding mechanisms and escape pathways of SARS-CoV-2 RBD-antibody complexes for groups A-D. These findings offer critical insights into the molecular determinants of antibody binding affinity, the roles of van der Waals and electrostatic interactions, and the evolutionary strategies employed by the virus to evade neutralization. First, the MM-GBSA analysis revealed that both van der Waals, and electrostatic interactions play crucial roles in stabilizing the RBD-antibody complexes, though their relative contributions vary across antibody groups and specific residues. Hydrophobic contacts, mediated by residues such as F456, F486, and V445, are essential for stabilizing the antibody-RBD complexes. These interactions are particularly dominant in antibodies like LY-CoV016 (Group A) and AZD8895 (Group B), where residues F456 and F486 form strong hydrophobic interactions with the antibody. Importantly, mutational scanning and binding calculations highlighted the role of electrostatic interactions that are mediated by key hotspots K417, E484, and K444 contributing significantly to binding through hydrogen bonds and salt bridges. For example, K417 in LY-CoV016 exhibits both favorable van der Waals and electrostatic interactions, highlighting its dual role in binding and immune evasion. Similarly, E484 in LY-CoV555 (Group C) and K444 in AZD1061 and LY-CoV1404 (Group D) are critical electrostatic hotspots. The synergistic effects of van der Waals and electrostatic interactions at key residues K417, F486, and K444 underscore their dual role in binding and immune evasion. Mutations at these sites disrupt both types of interactions, leading to significant reductions in binding affinity and susceptibility to escape mutations.

These findings align with experimental data, which identify these residues as dominant escape hotspots. The residue-based decomposition analysis identified critical escape hotspots and provided mechanistic insights into how mutations at these sites reduce antibody binding. Residues such as K417, E484, and K444 are frequently mutated in emerging variants (e.g., K417N, E484K, K444Q) due to their dual roles in binding and immune evasion. These mutations disrupt both van der Waals and electrostatic interactions, leading to significant reductions in binding affinity. We found that group B AZD8895, REGN10933 employ hotspots F486, N487, and G476 to form the strongest interactions. While these residues are critical for ACE2 binding and less prone to mutation, second-tier hotspots like S477, T478, and E484 are more flexible and serve as prime targets for escape mutations. This group demonstrates how the virus can exploit less critical residues to evade neutralization without compromising ACE2 binding. At the same time, group D AZD1061, REGN10987, LY-CoV1404 escaped by mutations at K444 and V445 demonstrating the importance of targeting conserved epitopes to develop broadly neutralizing antibodies. Using antibodies that target distinct epitopes ( particularly Group B and Group D antibodies exemplified by REGN10933/REGN10987 pair) in combination can reduce the likelihood of immune evasion, as the virus would need to evolve multiple simultaneous mutations to escape neutralization.

Our results dissected the mechanisms underlying group-specific binding and escape profiles. The analysis revealed distinct binding and escape profiles for each antibody group, reflecting their unique epitope targeting and interaction mechanisms. Group A LY-CoV016 relies heavily on residues like N487, F456, and K417 for binding. K417, in particular, plays a dual role, with both van der Waals and electrostatic interactions driving its importance as an escape hotspot. Group B AZD8895 featured F486 and N487 as critical binding hotspots, but their role in RBD stability limits the virus’s ability to mutate these sites. Instead, the virus targets second-tier hotspots like T478 and S477 to evolve escape mutations. E484 is the dominant binding and escape hotspot for Group C LY-CoV555 with mutations E484K significantly reducing neutralization efficacy. Secondary hotspots like V483, F486, and S494 also contribute to binding and immune evasion. Our results showed that group D (AZD1061, REGN10987, LY-CoV1404) shared K444 and V445 as critical binding hotspots for AZD1061 and LY-CoV1404, with strong contributions from both van der Waals and electrostatic interactions. REGN10987 exhibits a more balanced binding profile, with significant contributions from multiple residues. K444 is a dominant escape hotspot for all three antibodies, though LY-CoV1404 is less susceptible to escape mutations due to its targeting of a highly conserved region. These insights have important implications for the design of next-generation antibodies and therapeutic strategies to combat emerging SARS-CoV-2 variants. By targeting conserved epitopes and leveraging combination therapies, it may be possible to develop more resilient and broadly neutralizing antibodies that can withstand viral evolution.

4. Materials and Methods

4.1. Structure Preparation

The crystal and cryo-EM structures of the RBD-antibody complexes are obtained from the Protein Data Bank [105]. Hydrogen atoms and missing residues were initially added and assigned according to the WHATIF program web interface [106]. The structures were further pre-processed through the Protein Preparation Wizard (Schrödinger, LLC, New York, NY) and included the check of bond order, assignment and adjustment of ionization states, formation of disulphide bonds, removal of crystallographic water molecules and co-factors, capping of the termini, assignment of partial charges, and addition of possible missing atoms and side chains that were not assigned in the initial processing with the WHATIF program [106]. The missing loops in the cryo-EM structures were also reconstructed using template-based loop prediction approach ArchPRED [107]. The side chain rotamers were refined and optimized by SCWRL4 tool [108]. The protonation states for all the titratable residues of the Ab and RBD proteins were predicted at pH 7.0 using Propka 3.1 software and web server [109,110]. The protein structures were then optimized using atomic-level energy minimization with composite physics and knowledge-based force fields implemented in the 3Drefine method [111,112].

4.2. Coarse-Grained Molecular Simulations

Coarse-grained (CG) models are computationally effective approaches for simulations of large systems over long timescales. In this study, CG-CABS model [98,99,100] was used for simulations of the SARS-CoV-2 S complexes with antibodies. In this model, the amino acid residues are represented by Cα, Cβ, the center of mass of side chains and another pseudo atom placed in the center of the Cα-Cα pseudo-bond [98,99,100]. The position of Cα atoms is confined to a cubic lattice of a grid equal to 0.61 Å. The position of the side chain is dependent on the Cα−Cα−Cα angle of the main chain and the amino acid type. We used CABS model as this is a high-resolution knowledge-based coarse-grained force field that is based on potentials of the mean force obtained from statistical analysis of known protein structures and structural correlations of solved protein structures [98,99,100]. The sampling scheme of the CABS model used in our study is based on Monte Carlo replica-exchange dynamics and is modeled as a long random sequence of small local moves of individual amino acids in the protein structure as well as moves of small fragments consisting of two and three residues. The default settings were used for CG-CABS simulations in which soft native-like restraints are imposed only on pairs of residues with the distance between their C_α atoms smaller than 8 Å and both residues being part of the same secondary structure elements. No additional custom-designed distance restraints were applied to the simulation scheme. CABS-flex standalone package dynamics implemented as a Python 2.7 object-oriented package was used for fast simulations of protein structures [100]. A series of independent CG-CABS replica-exchange simulations were performed for each of the systems studied. In each simulation, the total number of cycles was set to 10,000 and the number of cycles between trajectory frames was 100. The conformational ensembles were subjected to MODELLER-based all-atom reconstruction including hydrogen atoms to produce atomistic models of simulation trajectories [113,114].

4.3. All-Atom Molecular Dynamics Simulations

All-atom MD simulations were performed for CB6/LY-CoV016 allowing for a direct comparative analysis of CG-CABS with MD simulations. Using a comparison of CG-CABS and MD simulations, we verified the reliability of the proposed simplified and less time-consuming simulation model. NAMD 2.13-multicore-CUDA package [115] with CHARMM36 force field [116] was employed to perform 1µs all-atom MD simulations for the RBD-Ab complexes. The structures of the complexes were prepared in Visual Molecular Dynamics (VMD 1.9.3) [117] and with the CHARMM-GUI web server [118,119] using the Solutions Builder tool. Hydrogen atoms were modeled onto the structures prior to solvation with TIP3P water molecules [120] in a periodic box that extended 10 Å beyond any protein atom in the system. To neutralize the biological system before the simulation, Na⁺ and Cl⁻ ions were added in physiological concentrations to achieve charge neutrality, and a salt concentration of 150 mM of NaCl was used to mimic physiological concentration. All Na⁺ and Cl⁻ ions were placed at least 8 Å away from any protein atoms and from each other. MD simulations are typically performed in an aqueous environment in which the number of ions remains fixed for the duration of the simulation, with a minimally neutralizing ion environment or salt pairs to match the macroscopic salt concentration [121].

First, minimization was performed for 100,000 steps with all the hydrogen-containing bonds constrained and the protein atoms fixed. In the second stage, minimization was performed for 50,000 steps with all the protein backbone atoms fixed and for an additional 10,000 steps with no fixed atoms. The equilibration is done for 1 ns by gradually increasing the system temperature in steps of 20 K, increasing from 10 K to 310 K, and at each step maintaining a restraint of 10 kcal mol⁻¹ Å⁻² on the protein C_α atoms. After the restraints on the protein atoms were removed, the system was equilibrated for an additional 10 ns. Long-range, non-bonded van der Waals interactions were computed using an atom-based cutoff of 12 Å, with the switching function beginning at 10 Å and reaching zero at 14 Å. The SHAKE method was used to constrain all the bonds associated with hydrogen atoms. The simulations were run using a leap-frog integrator with a 2 fs integration time step. The ShakeH algorithm in NAMD was applied for the water molecule constraints. The long-range electrostatic interactions were calculated using the particle mesh Ewald method [122] with a cut-off of 1.0 nm and a fourth-order (cubic) interpolation. The simulations were performed under an NPT ensemble with a Langevin thermostat and a Nosé–Hoover Langevin piston at 310 K and 1 atm. The damping coefficient (gamma) of the Langevin thermostat was 1/ps. In NAMD, the Nosé–Hoover Langevin piston method is a combination of the Nosé–Hoover constant pressure method [123] and piston fluctuation control implemented using Langevin dynamics [124,125]. An NPT production simulation was run on equilibrated structures for 1µs keeping the temperature at 310 K and a constant pressure (1 atm).

4.2. Binding Free Energy Computations: Mutational Scanning and Sensitivity Analysis

We conducted mutational scanning analysis of the binding epitope residues for the S RBD-Ab complexes. Each binding epitope residue was systematically mutated using all substitutions and corresponding protein stability and binding free energy changes were computed. BeAtMuSiC approach [126,127] was employed that is based on statistical potentials describing the pairwise inter-residue distances, backbone torsion angles and solvent accessibilities, and considers the effect of the mutation on the strength of the interactions at the interface and on the overall stability of the complex. The binding free energy of protein-protein complex can be expressed as the difference in the folding free energy of the complex and folding free energies of the two protein binding partners:

$$\Delta G_{bind}=G^{com}-G^A-G^B$$

(1)

The change of the binding energy due to a mutation was calculated then as the following:

$$\Delta \Delta G_{bind}=\Delta G_{bind}^{mut}-\Delta G_{bind}^{wt}$$

(2)

We leveraged rapid calculations based on statistical potentials to compute the ensemble-averaged binding free energy changes using equilibrium samples from simulation trajectories. The binding free energy changes were obtained by averaging the results over 1,000 and 10, 000 equilibrium samples for each of the systems studied.

2.3. Binding Free Energy Computations

We calculated the ensemble-averaged changes in binding free energy using 1,000 equilibrium samples obtained from simulation trajectories for each system under study. Initially, the binding free energies of the RBD-Ab complexes were assessed using the MM-GBSA approach [128,129]. Additionally, we conducted an energy decomposition analysis to evaluate the contribution of each amino acid during the binding of RBD to antibodies [130,131].

The binding free energy for the RBD-Ab complex was obtained using:

$$\mathsf{\Delta }{\mathrm{G}}_{\mathrm{bind}}={\mathrm{G}}_{\mathrm{RBD}-\mathrm{AB}}-{\mathrm{G}}_{\mathrm{RBD}}-{\mathrm{G}}_{\mathrm{AB}}$$

(3)

$$\mathsf{\Delta }{G}_{bind,MMGBSA}=\mathsf{\Delta }{E}_{MM}+\mathsf{\Delta }{G}_{sol}-T\mathsf{\Delta }S$$

(4)

where ΔE_MM is total gas phase energy (sum of ΔE_internal, ΔE_{electrostatic}, and ΔEvdw); ΔGsol is sum of polar (ΔG_GB) and non-polar (ΔG_SA) contributions to solvation. Here, G_RBD–AB represent the average over the snapshots of a single trajectory of the complex, G_RBD and G_AB corresponds to the free energy of RBD and Ab protein, respectively.

The polar and non-polar contributions to the solvation free energy is calculated using a Generalized Born solvent model and consideration of the solvent accessible surface area [132]. MM-GBSA is employed to predict the binding free energy and decompose the free energy contributions to the binding free energy of a protein–protein complex on per-residue basis. The binding free energy with MM-GBSA was computed by averaging the results of computations over 1,000 samples from the equilibrium ensembles. First, the computational protocol must be selected between the “single-trajectory” (one trajectory of the complex), or “separate-trajectory” (three separate trajectories of the complex, receptor and ligand). To reduce the noise in the calculations, it is common that each term is evaluated on frames from the trajectory of the bound complex. In this study, we choose the “single-trajectory” protocol, because it is less noisy due to the cancellation of intermolecular energy contributions. This protocol applies to cases where significant structural changes upon binding are not expected. Hence, the reorganization energy needed to change the conformational state of the unbound protein and ligand are also not considered. Entropy calculations typically dominate the computational cost of the MM-GBSA estimates. Therefore, it may be calculated only for a subset of the snapshots, or this term can be omitted [133,134]. However, for the absolute affinities, the entropy term is needed, owing to the loss of translational and rotational freedom when the ligand binds. In this study, the entropy contribution was not included in the calculations of binding free energies of the RBD-antibody complexes because the entropic differences in estimates of relative binding affinities are expected to be small [133,134]. MM-GBSA energies were evaluated with the MMPBSA.py script in the AmberTools21 package [135] and gmx_MMPBSA, a tool to perform end-state free energy calculations from CHARMM and GROMACS trajectories [136].

5. Conclusions

This study provides a comprehensive and detailed exploration of the molecular interactions between the SARS-CoV-2 S-RBD and neutralizing antibodies from four distinct groups (A-D). By integrating structural analysis, mutational scanning, and MM-GBSA binding free energy calculations, we have uncovered critical insights into the mechanisms of antibody binding, the roles of key residues in stabilizing these interactions, and the evolutionary strategies employed by the virus to evade immune detection. The residue-based decomposition analysis identified critical escape hotspots and provided mechanistic insights into how mutations at these sites reduce antibody binding. Residues such as K417, E484, and K444 are frequently mutated in emerging variants (e.g., K417N, E484K, K444Q) due to their dual roles in binding and immune evasion. These mutations disrupt both van der Waals and electrostatic interactions, leading to significant reductions in binding affinity. Our results dissected the mechanisms underlying group-specific binding and escape profiles. The analysis revealed distinct binding and escape profiles for each antibody group, reflecting their unique epitope targeting and interaction mechanisms. The study highlights the phenomenon of convergent evolution, where mutations at key residues (e.g., R346, K444, N460, F486, Q493) emerge independently in multiple variants to modulate antibody binding without compromising ACE2 affinity. The study also highlighted the importance of hinge sites in coordinating the dynamics of antibody binding and viral escape. Hinge positions such as K417, E484, and K444 play a critical role in modulating the conformational dynamics of the RBD. Mutations at these sites can disrupt collective motions and allosteric interactions, reducing the antibody’s ability to bind effectively. The flexibility of hinge regions makes them ideal targets for escape mutations, as they allow the virus to rapidly adapt to immune pressure without compromising RBD stability or ACE2 binding.

A comprehensive comparative energetic analysis of group D antibodies suggested that particularly LY-CoV1404, have a strong potential to compete with ACE2 for binding to the RBD due to their engagement with key ACE2-binding residues and their high binding affinity. While AZD1061 and REGN10987 also exhibit competitive binding, their efficacy is more moderate compared to LY-CoV1404. The ability of these antibodies to block ACE2 binding and prevent viral entry makes them valuable tools for therapeutic development. However, the emergence of mutations at critical residues like K444 and V445 highlights the need for continuous adaptation of antibody-based therapies to address evolving viral threats. The findings from this study may have useful implications for the design of next-generation antibodies and therapeutic strategies. Our results suggested that targeting conserved residues like K444 and V445, as seen in Group D antibodies, can lead to the development of broadly neutralizing antibodies that are less susceptible to escape mutations. Identifying critical binding hotspots and escape pathways can also guide the rational design of antibodies with improved binding affinity and resilience to viral evolution.