Comparative Analysis of Conformational Dynamics and Systematic Characterization of Cryptic Pockets in the SARS-CoV-2 Omicron BA.2, BA.2.75 and XBB.1 Spike Complexes with the ACE2 Host Receptor: Confluence of Binding and Structural Plasticity in Mediating Networks of Conserved Allosteric Sites

Abstract

In the current study, we explore coarse-grained simulations and atomistic molecular dynamics together with binding energetics scanning and cryptic pocket detection in a comparative examination of conformational landscapes and systematic characterization of allosteric binding sites in the SARS-CoV-2 Omicron BA.2, BA.2.75 and XBB.1 spike full-length trimer complexes with the host receptor ACE2. Microsecond simulations, Markov state models and mutational scanning of binding energies of the SARS-CoV-2 BA.2 and BA.2.75 receptor binding domain complexes revealed the increased thermodynamic stabilization of the BA.2.75 variant and significant dynamic differences between these Omicron variants. Molecular simulations of the SARS-CoV-2 Omicron spike full-length trimer complexes with the ACE2 receptor complemented atomistic studies and enabled an in-depth analysis of mutational and binding effects on conformational dynamic and functional adaptability of the Omicron variants. Despite considerable structural similarities, Omicron variants BA.2, BA.2.75 and XBB.1 can induce unique conformational dynamic signatures and specific distributions of the conformational states. Using conformational ensembles of the SARS-CoV-2 Omicron spike trimer complexes with ACE2, we conducted a comprehensive cryptic pocket screening to examine the role of Omicron mutations and ACE2 binding on the distribution and functional mechanisms of the emerging allosteric binding sites. This analysis captured all experimentally known allosteric sites and discovered networks of inter-connected and functionally relevant allosteric sites that are governed by variant-sensitive conformational adaptability of the SARS-CoV-2 spike structures. The results detailed how ACE2 binding and Omicron mutations in the BA.2, BA.2.75 and XBB.1 spike complexes modulate the distribution of conserved and druggable allosteric pockets harboring functionally important regions. The results are significant for understanding the functional roles of druggable cryptic pockets that can be used for allostery-mediated therapeutic intervention targeting conformational states of the Omicron variants.

1. Introduction

The enormous wealth of structural and biochemical studies has revealed mechanisms of the SARS-CoV-2 viral spike (S) glycoprotein consisting of the flexible amino (N)-terminal S1 subunit experiencing functional motions and structurally rigid carboxyl (C)-terminal S2 subunit. Conformational transformations of the SARS-CoV-2 S protein between the distinct closed and open states are driven by global movements of the S1 subunit consisting of an N-terminal domain (NTD), the receptor-binding domain (RBD) and two structurally conserved subdomains SD1 and SD2 which coordinate their dynamic changes to mediate functional responses of the S protein to binding with the host cell receptor ACE2 and antibodies [1,2,3,4,5,6,7,8,9,10]. The cryo-EM and X-ray structures of the SARS-CoV-2 S variants of concern (VOCs) in various functional states and complexes with antibodies detailed molecular mechanisms, diversity of binding epitopes and allosteric communications between functional regions that underlie binding with different classes of neutralizing antibodies. The cryo-EM structures of the S Omicron BA.1 trimer in the open and closed forms revealed that the dominantly populated conformation is the closed state with all the RBDs buried, leading to ‘conformational masking’ that may prevent antibody binding and neutralization at sites of receptor binding [11,12,13]. The thermal stability of the S-BA.1 and S-BA.2 protein ectodomains was examined using differential scanning fluorimetry (DSF) assays, revealing the reduced stability of the BA.1 RBD and increased stability of the BA.2 RBSD more stable than BA.1 but less stable than the Wu-Hu-1 [14,15,16,17]. The cryo-EM structures of the S Omicron protein ectodomain [18,19] showed that the S Omicron BA.1 proteins may preferentially adopt an open 1 RBD-up state predisposed for receptor binding. Moreover, there is greater variability of the S-BA.2 trimers seen in the unbound and ACE2-bound forms [18,19] may represent the molecular basis for higher transmission and infection rates of the BA.2 Omicron sublineage compared to the BA.1 variant. Cryo-EM structures of the S-BA.2 complexes with human ACE2 revealed a reorganization of the highly antigenic NTD regions, resulting in the resistance to antibodies and suggesting that variant-induced modulation of conformational plasticity can underlie the antigenic drift in the Omicron variant landscape and a remarkable increase in the antibody escape pattern [20]. Structural and biophysical studies of the RBD-ACE2 complexes for the Omicron variants revealed that the binding affinity of the Omicron BA.2 with ACE2 is stronger than the affinities of the BA.3 and BA.1 subvariants [21]. Thermal shift assays showed that the Omicron BA.2-RBD is more stable than that from BA.1, but the observed lower temperature for dissociation of the S-BA.2 trimer indicated that BA.2 is more dynamic than the S-BA.1 trimer [22]. The structural and functional analysis of the S-BA.2 protein trimer revealed three distinct states representing the closed 3-RBD-down conformation, a 1-RBD-up conformation and an RBD-intermediate conformation [23]. The cryo-EM structures of the S trimers for BA.1, BA.2, BA.3 and BA.4/BA.5 subvariants showed that the S-BA.1 trimer is stabilized in an open conformation while S-BA.2 exhibiting two conformational states is the least stable among BA.1, BA.2, BA.3 and BA.4/BA.5 variants due to more dynamic and less compact inter-protomer arrangements [24]. The cryo-EM conformations of the BA.2.75 S trimer in the open and closed forms, as well as structures of the open BA.2.75 S trimer complexes with ACE2, pointed to the increased structural heterogeneity of S1 regions in the BA.2.75 timer altering the interactions between the RBDs and yet exhibiting a more rigid and compact RBD structure [25]. This study also reported thermal stability of the Omicron variants at neutral pH, showing that the BA.2.75 S-trimer was the most stable, followed by BA.1, BA.2.12.1, BA.5 and BA.2 variants, while BA.2.75 also displayed 4–6-fold increased binding affinity to hACE2 compared with other Omicron variants [25]. Biophysical characterizations of the BA.2.75 variant unveiled a better balance between immune evasion and ACE2 binding, where the binding affinity to ACE2 is increased 9-fold as compared to the BA.2 variant [26,27,28].

The new recombinant variants such as BA.2.75.2, XBB.1 and XBB.1.5 display substantial growth advantages over previous Omicron variants, suggesting that the immune pressure promotes convergent evolution in which some RBD residues (R346, K356, K444, V445, G446, N450, L452, N460, F486, F490, R493 and S494) were observed in at least five independent Omicron sublineages that exhibited a high growth advantage [29,30]. The cryo-EM structures of the XBB.1 S ectodomain and the XBB.1 S-ACE2 complex revealed two dominant closed states and the RBD one-up state that becomes stabilized in the presence of the ACE2 receptor. These studies showed that despite the thermodynamic stability of the XBB.1 S protein, there is considerable plasticity even in the packed closed form of the S trimer [31]. XBB.1.5 is equally immune evasive as XBB.1 but may have a growth advantage by virtue of the higher ACE2 binding affinity owing to a single S486P mutation as F486S substitution in XBB.1 [32]. Subsequent functional studies confirmed the growth advantage and the increased transmissibility of the XBB.1.5 lineage due to the preserved strong neutralization resistance and the improved ACE2 binding affinity [33].

The cryo-EM and crystal structures of the S Omicron trimers and complexes with ACE2 and antibodies provided a staggering amount of accurate high-resolution structural information and unveiled the diversity of conformational states sampled by Omicron variants. Nonetheless, the details of the intrinsic conformational dynamics and drivers of functional adaptability mechanisms are often hidden or obscured due to the rapid interconversion of states and the complexity of binding-induced modulation of flexibility seen in the experiments. Hydrogen/deuterium-exchange mass spectrometry (HDX-MS) is a powerful approach for monitoring local protein structural dynamics under native solution conditions whereby by tracking deuterium uptake kinetics for peptide segments throughout a given protein, these approaches can inform of dynamics and residue-specific changes in the conformational dynamics induced by mutations or binding. HDX-MS studies uncovered an alternative open S trimer conformation that interconverts slowly with the known prefusion structures and can dynamically expose novel epitopes in the conserved region of the S2 trimer interface, unexpectedly suggesting the emergence of hidden cryptic epitopes in a highly conserved region of the protein [34]. Another HDX-MS study identified changes in the S dynamics for VOCs, revealing that Omicron mutations may preferentially induce closed conformations with dynamic core helices in the S2 subunit exploited in the fusion stage and that the NTD acts as a hotspot of conformational divergence driving immune evasion [35]. Furthermore, this study suggested that not only does the ACE2 binding promote the opening of the RBDs, but it can also allosterically enhance the dynamics of the S2 core to prime spike prefusion conformations for the transition to the post-fusion form [35]. Comparative HDXMS analysis of the WT, D614G, Alpha, Delta and Omicron BA.1 S variants tracked evolution of the intrinsic S dynamics by examining the progressively emerging VOC’s [35], showing that Omicron BA.1 can induce greater stabilization of the trimeric stalk interface in the S2 subunit concurrently with the increased NTD and RBD dynamics which may have direct implications for ACE2 binding, and proteolytic processing [36]. HDX-MS studies of the unbound S protein trimer also showed the highest relative exchange in the S2 subunit and helical segments of the stalk regions (central helix CH and heptad repeats HR1 and HR2), while a comparative analysis of the S protein and S-ACE2 complex suggested that ACE2 binding allosterically enhances dynamics at a distal S1/S2 cleavage site and can modulate dynamics of the stalk hinge region, thereby highlighting the role of stalk and proteolysis sites as dynamic allosteric centers in the prefusion state [37]. Another HDX-MS study of S-ACE2 binding confirmed that ACE2 can induce enhanced dynamics throughout the entire S protein, with allosteric changes being propagated to the S2 hinge region and to the central helical bundle of the S2 subunit [38]. Overall, these studies suggested that considerable conformational plasticity can be preserved in both open and stable closed S trimers, including more rigid S2 subunit, while Omicron mutations and ACE2 binding can allosterically modulate the S dynamics at distal functional regions.

Structural and functional studies also suggested that Omicron mutations can also impact the dynamics of binding epitopes and mediate the formation of transient cryptic binding pockets [39,40,41]. These binding pockets become exposed or “unmasked” when the S protein undergoes spontaneous conformational changes or interacts with other molecules and antibodies. The cryo-EM structure of the S protein with linoleic acid (LA) complex revealed an allosteric binding pocket in the S-RBD that is controlled through the opening of a gating helix and can allosterically induce stabilization of the closed S trimer and promote decreased S-ACE2 binding [39,40]. Structural studies discovered a highly conserved cryptic epitope that is buried deep inside the trimeric interface in the SD1 region of the S protein and is shared between the Omicron variants of SARS-CoV-2, including WT, BA.1, BQ.1.1, XBB, XBB.1.5 and the XBB.1.16 variants [42]. These highly conserved epitopes enable broad antibody recognition of Omicron variants, suggesting S trimer disassembly as a mechanism of the antibodies targeting trimeric interface epitopes. These studies further underscored the importance of identifying conserved cryptic epitopes as essential targets for the design of universal vaccines and broadly neutralizing antibodies [42]. The recent structural investigations unveiled another cryptic binding pocket located in the NTD region, also demonstrating that tetrapyrrole products of heme metabolism, biliverdin and bilirubin, can bind to this site with nanomolar affinity and inhibit the neutralizing activity of the NTD antibodies recognizing these epitopes [43,44,45]. By leveraging HDX-MS and antibody engineering approaches, structural studies confirmed a cryptic hinge epitope located in the S2 region between the CH and HR1 segments that plays a critical role in the spike conformational changes required for fusion of the viral envelope and target cell membrane [46]. This S2 hinge epitope is partially occluded by the S1 domain, and access to this epitope may be controlled through RBD conformational switching. These investigations also underscored a considerable diversity of potential cryptic pockets, with some of them being present only transiently while others showed conservation across different S protein states.

Computer simulation studies provided important atomistic insights into understanding the dynamics of the SARS-CoV-2 S protein and the effects of Omicron mutations on the conformational plasticity of the S protein states and their complexes with diverse binding partners. A significant number of computational and experimental studies emphasized the critical role of electrostatic interactions as a dominant thermodynamic force leading to the binding of the S-protein with the ACE2 receptor and antibodies. The cryo-EM structures of the S Omicron BA.1 trimer showed that the Omicron sites N856K, N969K and T547K can promote favorable electrostatic interactions and lead to the hydrogen bonds with D658, Q755 and S982 from neighboring subunits which can confer the enhanced stability for the Omicron S-trimer [13]. MD simulations of the RBD-ACE2 complexes for BA.1 BA.1.1, BA.2, an BA.3 Omicron subvariants were combined with a systematic mutational scanning of the RBD-ACE2 binding interfaces, revealing functional roles of the key Omicron mutational site R493, R498 and Y501 acting as binding energy hotspots, and drivers of electrostatic interactions [47]. Analysis of a large number of mutant S protein complexes revealed that high-affinity RBD mutations tend to cluster near known ACE2 recognition sites by optimizing the binding energies using electrostatic interaction forces [47]. Recent studies mapped the electrostatic potential surface of S protein and major variants to show accumulated positive charges at the ACE2-binding interface, revealing the critical role of complementary electrostatic interactions driving the enhanced affinity of the Omicron S-ACE2 complexes [48]. Extensive simulations investigated the electrostatic features of the RBD for the several S variants, their main physical–chemical properties, and binding affinities at several pH regimes revealing that the virus evolution may primarily exploit the electrostatic forces and Coulombic charge–charge interactions to make RBD more positively charged and improve the RBD-ACE2 binding affinity and transmission [49]. Computational analysis of the pH-dependent electric potential on the surface of the S proteins and pH-dependent free energies of S-ACE2 binding were experimentally verified by isoelectric focusing, showing that binding with the host receptor is conditioned by electrostatic attractions of the oppositely charged receptor and viral protein [50]. The comparison of the local electrostatic potentials of the Wu-Hu-1 and Omicron strains showed that the point mutations alter the electrostatic potential in a relatively small area on the RBD surface, resulting in a stronger S–ACE2 association at pH 5.5 [50]. The analysis of Omicron variants showed that the electrostatic potential of the Omicron variants is higher than in other VOCs, including the Delta variant, and indicated the relevance of the high electrostatic potential for virus transmission and infectivity [51,52,53]. Despite mutational differences between BA.1, BA.2 and BA.2.75 lineages, the recent study of the electrostatic potentials showed a strong positive electrostatic surface in the S-RBD regions for these variants supporting the mechanism in which the enhanced electropositive character of the RBD variants is strongly related with their enhanced capacity for the virus transmission [54]. Our analysis also suggested that the evolution of the SARS-CoV-2 virus and the emergence of recombinant variants may have signaled a certain critical plateau of electrostatic positively charged RBD distribution that is optimal to complement ACE2 electrostatic potential, and convergent mutations could have emerged to balance multiple tradeoffs rather than progressively improving the RBD-ACE2 binding affinities [55].

Glycosylation of the SARS-CoV-2 proteins and shielding of the receptor binding sites by glycans is a common mechanism of viral glycoproteins for evading the immune system, and recent studies showed that glycosylation on SARS-CoV proteins can camouflage immunogenic protein epitopes [56,57,58]. Site-specific mass spectrometric mapping and cryo-EM analysis of the S protein established the composition and functional role of 22 N-linked glycosylation sites and 4 O-linked glycosylation sites of S protein that occlude distinct regions on the protein surface [58]. A combination of antigenic screens and cryo-EM structure determination indicated that an N-glycan deletion at N234 results could reduce the population of the S trimer with ‘up’ state RBDs, while glycan deletion at position N165 results in an increase in the open states [59]. These studies showed that the glycans on the S protein can conceal immunogenic epitopes and modulate conformational landscapes, making it more challenging for the immune system to mount a robust response against the virus. MD simulations probed the effect of glycan heterogeneity on the antigenicity of the S glycoprotein where site-specific N-linked microheterogeneity is defined at 22 spike sites, thus highlighting roles for glycans in sterically masking epitopes and modulating S-ACE2 interactions [60]. Large-scale simulations of the S glycoprotein embedded in the viral membrane with a complete glycosylation profile suggested that the dynamics of the glycan shield can be coupled to allosteric conformational changes of the S protein and provide an additional layer of molecular control over spike response to the host receptor [61]. All-atom MD simulations of the glycosylated, full-length, membrane-bound ACE2 receptor probed the intrinsic dynamics of the ACE2 receptor in the context of the cell surface [62]. MD simulations demonstrated that glycosylation of the human ACE2 receptor contributes substantially to the binding of the virus with two glycosylation sites, N90 and N322, having opposite effects on the S protein binding [63]. In situ structural analysis of S protein using electron tomography and MD simulations showed that the opening of the RBD observed in recombinant S trimers can also take place on the virus surface and that the hinges at the stalk domain allow the protein to scan the host cell surface, shielded from antibodies by an extensive glycan coverage [64]. A number of computational studies employed advanced sampling techniques to gain insight into time scales of conformational changes in the S protein and attempted to characterize conformational landscapes of the S Omicron variants. Adaptive sampling simulations performed on a large scale for the viral proteome captured conformational heterogeneity of the S protein and predicted multiple cryptic epitopes, but the functional relevance and experimental validation of these predicted pockets were lacking [65]. Notably, through milliseconds of conformational sampling, simulations captured this rare event for both glycosylated and unglycosylated S and found that glycosylation only slightly increases the population of the open state [65]. Another large-scale MD simulation study using replica exchange umbrella sampling calculations with more than a thousand windows and an aggregate total of 160 μs of simulation found the glycosylated S protein has a higher barrier to opening and also energetically favors the closed form [66]. While the proper treatment of the glycan environment of the S protein is important to capture atomistic details of interactions and play a role in the stabilization/obstruction of cryptic pockets, these simulations are extremely time-consuming and are challenging when studies involve a wide range of S trimer protein structures for various variants and complexes with ACE2.

The replica-exchange molecular dynamics (MD) simulations examined conformational landscapes of the full-length S protein trimers, discovering the transition pathways, hidden functional intermediates along open-closed transition pathways and previously unknown cryptic pockets that were consistent with FRET experiments [67]. Computer simulation mapping of the S protein pockets using benzene probes reproduced several experimentally discovered allosteric sites and identified a spectrum of novel cryptic and potentially druggable pockets [68], but this analysis was largely focused on the closed prefusion states of the Wu-Hu-1 S protein and the effects of various Omicron variants and ACE2 binding on the evolution of cryptic sites remained unexplored. MD simulations of the unbound or ACE2-bound RBD conformations combined with pocket analysis and druggability prediction identified several promising druggable sites, including one located between the RBD monomers [69]. A network-based adaptation of the reversed allosteric communication approach identified allosteric hotspots and RBD binding pockets in the Omicron variant RBD-ACE2 complexes [70]. Integration of computational and experimental studies enabled the discovery and validation of the cryptic allosteric site located between subdomains of the S protein, with several compounds targeting this site showing characteristic binding and anti-virus activities [71,72]. Collectively, the recent structural and computational studies suggested that Omicron mutations and ACE2 binding have a significant effect on mediating conformational dynamics changes in the S protein, including not only the binding interface regions but also allosterically induced variations and altered plasticity at the remote conserved S2 regions.

Despite recent advances, our current structural and functional knowledge of the cryptic binding pockets and potential allosteric sites in the S protein is incomplete and lacks analysis of evolution, functional mechanisms and validity of the predicted cryptic sites. In addition, there is no clear understanding of the effects exerted by different Omicron variants, conformational states and ACE2 binding on the evolution and redistribution of the cryptic pockets in the S protein. There is also a steadily increasing interest in identifying novel druggable sites, particularly in the conserved S2 subunit, as most recently developed broad-spectrum fusion inhibitors and candidate vaccines can target the functional regions in the S2 subunit [73,74,75]. To investigate and explore these outstanding issues, in the current study, we perform a systematic comparative analysis of the conformational dynamics, allostery and cryptic binding pockets in the RBD-ACE2 complexes trimers and S trimer complexes with the ACE2 receptor. We first perform multiple microsecond MD simulations and Markov state model (MSM) analysis of the SARS-CoV-2 S RBD-ACE2 complexes for the Omicron BA.2, BA.2.75 and XBB.1 variants. Using a comparative MSM analysis, we discover variant-specific changes in conformational mobility and confirm the stability of the BA.2.75 RBD and increased mobility of the XBB.1 RBD. Using conformational ensembles of the SARS-CoV-2 Omicro S trimers, we conducted a systematic binding pocket screening and analysis of functional cryptic pockets in the BA.2, BA.2.75 and XBB.1 complexes with ACE2. The results of this study connect insights from conformational dynamics analysis, comparative mutational scanning of the S-ACE2 binding and the inter-protomer interactions with the evolution of cryptic binding sites. We show that our approach can capture and reproduce all experimentally known allosteric sites and discover networks of conserved cryptic pockets preserved in different conformational states of the Omicron variants and in the S-ACE2 complexes. The results of our study suggest that despite the general rigidity of the S2 regions in comparison with the more adaptable S1 subunit, there is still an appreciable level of conformational adaptability in the S2, resulting in a significant number of dynamic cryptic pockets. The results detailed how mutational and conformational changes in the BA.2 and BA.2.75 spike trimers can modulate the distributions and mediate networks of inter-connected conserved and variant-specific druggable allosteric pockets. The results of this study can be important for understanding mechanisms underlying functional roles of druggable cryptic pockets that can be used for both site-specific and allostery-inspired therapeutic intervention strategies targeting distinct conformational states of the SARS-CoV-2 Omicron variants.

2. Materials and Methods

2.1. All-Atom Molecular Dynamics Simulations

The crystal structures of the BA.2 RBD-ACE2 (pdb id 7XB0) [21] and BA.2.75 RBD-ACE2 complexes (pdb id 8ASY) [76] and the atomic coordinates for the structures of XBB.1 S RBD bound to ACE2 (pdb id 8IOV) [31] as well as the cryo-EM structures of the BA.2 S trimer with two human ACE2 bound (pdb id 7XO7) [22], BA.2 S trimer with three ACE2 bound (pdb id 7XO8) [22], BA.2.75 S trimer with one human ACE2 bound (pdb id 7YR2) [25] and XBB.1 S trimer with one human ACE2 bound [31] (Figure 1) were obtained from the Protein Data Bank [77]. During structure preparation stage, protein residues in the crystal structures were inspected for missing residues and protons. Hydrogen atoms and missing residues were initially added and assigned according to the WHATIF program web interface [78]. The missing loops in the studied structures of the SARS-CoV-2 S protein were reconstructed and optimized using template-based loop prediction approach ArchPRED [79]. The side chain rotamers were refined and optimized by SCWRL4 tool [80]. The protonation states for all the titratable residues of the ACE2 and RBD proteins were predicted at pH 7.0 using Propka 3.1 software and web server [81,82]. The protein structures were then optimized using atomic-level energy minimization with composite physics and knowledge-based force fields implemented in the 3Drefine method [83]. In the current study, we considered glycans that were resolved in the cryo-EM structures. The structurally resolved 2-acetamido-2-deoxy-beta-D-glucopyranose, 2-acetamido-2-deoxy-beta-D-glucopyranose-(1-4)-2-acetamido-2-deoxy-beta-D-glucopyranose and chloride ions present in the RBD-ACE2 structures were included and optimized.

We employed CHARMM36M force field [84] with the TIP3P water model [85] to conduct microsecond all-atom MD simulations for structures of the BA.2, BA.2.2.75 and XBB.1 RBD-ACE2 complexes. The protein systems were solvated in 130 Å × 85 Å × 75 Å water boxes. In each system, sodium and chloride ions were added to maintain an ionic strength of 0.1 M. After energy minimization, the systems were first heated up from 100 to 300 K with a temperature increment of 20 K per 50 picoseconds (ps). Consequently, the systems were subjected to 1.5 nanoseconds (ns) isothermal–isobaric (NPT) equilibrations at 300 K (equilibrium run), followed by 1 microsecond (µs) canonical (NVT) simulations (production run) at 300 K. While the production stage of MD simulations is often performed under the NPT ensemble to enable more accurate comparison with the experimental data, the NPT ensemble often needs significantly longer time to reach convergence because both temperature and volume fluctuate during the simulations. In this study, we performed the equilibration stage in the NPT ensemble to allow the total energy and volume of the system to fluctuate to reach conditions suitable for the target pressure, which is often employed in MD simulations to maintain a specific pressure. In this case, the production simulations could be carried out in the NVT ensemble using the equilibrium volume sampled in the equilibration NPT simulations. By using the equilibrium volume sampled during the NPT equilibration in the NVT ensemble, the desired pressure can be maintained while allowing the system to evolve in terms of its energy and temperature. Although the binding association constants for complexes are often measured under constant (atmospheric) pressure and temperature (NPT ensemble), corresponding to Gibbs free energy, the respective experimentally measured macroscopic parameters obtained in the NVT ensemble, corresponding to Helmholtz free energy, are typically very similar.

While most differences between MD simulations are often attributed to the force fields, other factors can influence the outcome, including the water model, algorithms that constrain motion, how atomic interactions are handled, and the simulation ensemble employed. The protocol used in our study aims to take advantage of the NVT ensemble for better convergence than NPT ensemble, given the same length of simulations. Using the NPT ensemble for equilibration helps in achieving a stable and realistic pressure for the system. Once the system has equilibrated in terms of both pressure and temperature, transitioning to the NVT ensemble for production allows you to maintain a consistent temperature while studying the dynamics and properties of the system. This approach helps ensure convergence and stability in the simulation by allowing the system to equilibrate under the appropriate thermodynamic conditions (NPT) before moving into the production phase (NVT). A detailed comparison of MD simulations using different force fields and ensembles (including NAMD with the CHARMM36 force field and NVT for the production phase) produced models and equilibrium ensembles that are in similar agreement with experimental results [86,87,88]. Moreover, the agreement between MD simulations and experimental observables may not necessarily constitute a validation of the conformational ensemble produced by MD simulations, as diverse conformational ensembles may produce averages that are consistent with the experiment [86]. Importantly, as the main objective of our study is the examination of the effects of different mutations on the SARS-CoV-2 S protein complexes, the employment of NPT or NVT ensemble for the production simulations is not expected to produce appreciable differences and lead to different conclusions.

Snapshots of the production run were saved every 100 ps. In all simulations, the SHAKE constraint was used to constrain bonds associated with hydrogen atoms in the solvent molecules and the proteins [89]. The nonbonding interactions within 10 Å were calculated explicitly. The Lennard-Jones interactions were smoothed out to zero at 12 Å. The long-range electrostatic interactions were calculated using the particle mesh Ewald method [90] using a fourth-order (cubic) interpolation and a cutoff of 10.0 Å. The simulations were conducted using OpenMM (version 7.6.0) software package [91]. For each system, MD simulations were conducted three times in parallel to obtain comprehensive sampling. Each individual simulation has 10,000 frames. For BA.2 RBD-ACE2, BA.2.75 RBD-ACE2 and XBB.1 RBD-ACE2 systems, MD simulations were conducted three times in parallel for each system to obtain comprehensive sampling. Each individual simulation of the studied Omicro n RBD-ACE2 complexes stored 10,000 frames; thus, we collected 30,000 frames for each of the complexes.

The root mean square deviation (RMSD) between one reference state and a trajectory of structures is often calculated to measure the dissimilarity of the trajectory conformational ensemble to the reference. This reference is frequently the first frame of the trajectory to monitor the evolution of dynamics from the initial starting point. The stable RMSD values from a reference structure are frequently used as a measure of conformational convergence.

$$\mathrm{RMSD}(\mathbf{x},{\mathbf{x}}^{\mathrm{ref}})=\sqrt{\frac{1}{n}{\sum }_{i=1}^n|{\mathbf{x}}_i-{\mathbf{x}}_i^{\mathrm{ref}}{|}^2}$$

(1)

We used MDTraj Python library [92] for the analysis, manipulation and visualization of MD simulation trajectories. MDtraj provides a wide range of functionalities to extract structural and dynamic information from MD trajectories, including a fast RMSD procedure [93] that executes QCP algorithm [94] three times faster than the original implementation. The RMSD calculations were also performed using MDAnalysis Python toolkit (www.mdanalysis.org, accessed on 10 May 2023), utilizing the fast QCP algorithm and the optimal rotation matrix for superposition of simulation frames [95].

The root-mean-square fluctuation (RMSF) of a structure is the time average of the RMSD. The RMSF is calculated according to the equation blow, where $〈{\mathbf{R}}_i〉$ is the mean atomic coordinate of the ith C_α atom and ${\mathbf{R}}_i$ is its instantaneous coordinate:

$${\rho }_i=\sqrt{⟨({\mathbf{R}}_i-⟨{\mathbf{R}}_i⟩{)}^2⟩}$$

(2)

The RMSF analysis is restricted to backbone or alpha-carbon atoms of the protein conformations for better characterization of conformational dynamics and structural changes.

2.2. Markov State Model Analysis of the Omicron RBD-ACE2 Complexes

The time-structure Independent Components Analysis (tICA) method identifies the slowest degrees of freedom and can preserve the kinetic information present in the MD trajectories by maximizing the auto-correlation function [96,97,98,99]. In contrast to principal component analysis (PCA), which finds coordinates of maximal variance, tICA finds coordinates of maximal autocorrelation at the given lag time. Therefore, the tICA approach is useful for finding the slow components in a dataset and a robust method to transform MD simulation information before clustering data for the construction of a Markov model. Using a time-series of molecular coordinates provided by an n-dimensional MD trajectory x(t) = (x₁(t), …, x_n(t))^T ∈ ℝⁿ with Cartesian coordinates (x₁, …, x_n) tICA can reduce the dimensionality of the trajectories and determine the slowest independent collective degrees of freedom onto which the projections of the initial dataset have the largest time-autocorrelation. The tICA approach identifies the slowest degrees of freedom by solving the generalized eigenvalue problem:

$$\overline{C}F=CFK$$

(3)

where K = diag(k₁, …, k_n) and F = (f₁, …, f_n) are the eigenvalue and eigenvector matrices, respectively; $C$ and $\overline{C}$ are the covariance matrix and the time-lagged covariance matrix of the coordinate vector defined as follows:

$$C=\langle {(x(t)-\langle x(t)\rangle )}^T(x(t)-\langle x(t)\rangle )\rangle$$

(4)

$$\overline{C}=\langle {\left(x\right(t)-\langle x(t)\rangle )}^T(x(t+t_0)-\langle x(t)\rangle )\rangle$$

(5)

To obtain a symmetric time-lagged covariance matrix, $\frac{1}{2}(\overline{C}+{\overline{C}}^T)$ is calculated. The latter step assumes the time reversibility of the process, which is satisfied in MD simulations. The featurization and dimensionality reduction were performed using the PyEMMA package [10]. We employed Stochastic Markov state models (MSMs) [100,101,102,103,104,105] to characterize conformational landscapes of the Omicron RBD-ACE2 complexes and describe the transitions between functional protein states. In MSM, protein dynamics is modeled as a kinetic process consisting of a series of Markovian transitions between different conformational states at discrete time intervals. A specific time interval, referred to as lag time, needs to be determined to construct transition matrix. First, k-means clustering method is conducted on projected low-dimensional space, and each simulation frame is assigned to a microstate. The transition counting is constructed based on a specific time interval lag time τ. Macrostates are kinetically clustered based on the Perron-cluster cluster analysis (PCCA++) [106] and considered to be kinetically separate equilibrium states. The transition matrix and transition probability were calculated to quantify the transition dynamics among macrostates. The corresponding transition probability from state i to state j is calculated as

$${\mathrm{P}}_{\mathrm{i}\mathrm{j}}(\mathsf{\tau })=\mathrm{P}\mathrm{r}\mathrm{o}\mathrm{b}\left({\mathrm{x}}_{\mathrm{t}}+\mathsf{\tau }\in {\mathrm{S}}_{\mathrm{j}}|{\mathrm{x}}_{\mathrm{t}}\in \mathrm{S}\mathrm{i}\right)$$

(6)

A proper lag time is required for MSM to be Markovian. The value of the lag time and the number of macrostates are selected based on the result of estimated relaxation timescale [107]. The implied timescales can be calculated using the eigenvalues (${\lambda }_i$) in the transition matrix as

$$t_i=-\frac{\tau }{\ln \left|{\lambda }_i\left(\tau \right)\right|}$$

(7)

The number of protein metastable states associated with these slow relaxation timescales can be inferred based on the convergence of implied relaxation time scale. These metastable states effectively discretize the conformational landscape. The Markov state model building was conducted using PyEMMA package (v2.5.12) [108]. Based on the transition matrix we obtain implied timescales for transitioning between various regions of phase space and use this information to determine the number of metastable states.

2.3. Coarse-Grained Dynamics Simulations of the Omicron S Trimer Complexes with ACE2

Coarse-grained Brownian dynamics (CG-BD) simulations have been conducted for the cryo-EM structures of the BA.2 S trimer with two human ACE2 bound (pdb id 7XO7) [22], BA.2 S trimer with three ACE2 molecules bound (pdb id 7XO8) [22], BA.2.75 S trimer with one human ACE2 bound (pdb id 7YR2) [25] and XBB.1 S trimer with one human ACE2 bound [31] (Figure 1). We employed ProPHet (Probing Protein Heterogeneity) approach and program [108,109,110]. BD simulations employed a high-resolution CG protein representation of the SARS-CoV-2 S Omicron trimer structures that can distinguish different residues. In this model, each amino acid is represented by one pseudo-atom at the Cα position and two pseudo-atoms for large residues. The interactions between the pseudo-atoms are treated according to the standard elastic network model (ENM) in which the pseudo-atoms within the cut-off parameter, R_c = 9 Å are joined by Gaussian springs with the identical spring constants of γ = 0.42 N m⁻¹ (0.6 kcal mol⁻¹ Å⁻². The simulations use an implicit solvent representation via the diffusion and random displacement terms and hydrodynamic interactions through the diffusion tensor using the Ermak-McCammon equation of motions and hydrodynamic interactions [111,112]. The stability of the SARS-CoV-2 S Omicron trimers was monitored in multiple simulations with different time steps and running times. We adopted Δt = 5 fs as a time step for simulations and performed 100 independent BD simulations for each system using 500,000 BD steps at a temperature of 300 K. The CG-BD conformational ensembles were also subjected to all-atom reconstruction using PULCHRA method [113] and CG2AA tool [114] to produce atomistic models of simulation trajectories.

While the proper treatment of glycan environment of the S protein is important to capture atomistic details of interactions and play a role in stabilization/obstruction of cryptic pockets, these simulations are extremely time-consuming and are challenging when studies involve a wide range of S trimer protein structures for various variants and complexes with ACE2. It is important to note that the choice to simulate the S protein with or without glycans should be based on the research objectives and the specific scientific questions being addressed. In the current study, we included glycans that were resolved in the cryo-EM structures of S Omicron trimer complexes (on both S protein and ACE2), which allowed for adequate characterization of functional dynamics, energetic changes and robust screening of potential cryptic pockets. The reconstructed atomistic structures from simulation trajectories were decorated by N-acetyl glycosamine (NAG) glycan residues and optimized. The structure of glycans at particular glycosites of the closed and open states of SARS-CoV-2 S protein was previously determined, and these glycans were incorporated in atomistic modeling of the SARS-CoV-2 mutant structures. The glycosylated microenvironment for atomistic models of the simulation trajectories was mimicked by using the structurally resolved glycan conformations for 16 out of 22 most occupied N-glycans in each protomer (N122, N165, N234, N282, N331, N343, N603, N616, N657, N709, N717, N801, N1074, N1098, N1134, N1158) as determined in the cryo-EM structures of the SARS-CoV-2 spike S trimer in the closed state (K986P/V987P,) (pdb id 6VXX) and open state (pdb id 6VYB) and the cryo-EM structure SARS-CoV-2 spike trimer (K986P/V987P) in the open state (pdb id 6VSB). The glycan-decorated atomistic models were subsequently used for the ensemble-based structural and dynamic analyses of the Omicron mutants examined in this work.

2.4. Mutational Scanning Analysis of the RBD-ACE2 Binding Interactions and Inter-Protomer Trimer Interactions

Mutational scanning analysis of the binding epitope residues for the Omicron RBD-ACE2 and S trimer-ACE2 complexes was conducted using PoPMuSiC approach [115,116,117] that is based on statistical potentials which include contributions from 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$$

(8)

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

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

(9)

In this approach, the binding interface residues are defined based on the condition that the difference between a residue’s solvent accessibility in the complex and apo-protein is at least 5%. The solvent accessibility is defined as the ratio of the solvent-accessible surface in the considered structure relative to an extended tripeptide Gly-X-Gly. Each binding epitope residue was systematically mutated using all substitutions and corresponding protein stability and binding free energy changes were computed. We computed the ensemble-averaged binding free energy changes using equilibrium samples from simulation trajectories. The binding free energy changes were computed by averaging the results over 1000 equilibrium samples that were evenly distributed along each of the three independent MD simulations for each of the studied systems. Hence, a total of 3000 protein snapshots were employed in the computations of the mutation-induced binding free energies for each of the studied systems.

In addition to mutational scanning of the RBD-ACE2 binding interfaces, the S protein residues involved in the inter-protomer contacts in the S trimer structures were also systematically mutated to evaluate the protein stability differences between BA.2, BA.2.75 and XBB.1 Omicron variants. If a free energy change between a mutant and the wild type (WT) proteins ΔΔG = ΔG (MT) − ΔG (WT) > 0, the mutation is destabilizing, while when ΔΔG < 0 the respective mutation is stabilizing. We computed the average ΔΔG values using 1000 samples from the CG-BD equilibrium ensembles.

2.5. Machine Learning Detection of Cryptic Pockets

We used two different complementary approaches for identification of the cryptic binding pockets in the conformational ensembles of the Omicron RBD-ACE2 and S timer-ACE2 complex structures. A template-free P2Rank approach is among the most efficient and fast available ML tools for prediction of ligand binding sites that combines sequence and structural data to rank potential binding sites based on their likelihood of binding a specific ligand [118,119]. P2Rank uses support vector machine (SVM), random forests (RF) and artificial neural networks (ANNs) to learn the ligandability (or druggability) of a local chemical environment that is centered on points placed on the protein’s solvent-accessible surface of the protein [118,119]. P2Rank v2.4 with default parameters was deployed to identify pockets across all of the representative states from our simulations. This method belongs to the class of pocket-centric methods capable of uncovering novel binding sites and was shown to outperform existing traditional template-free methods such as Fpocket [120], SiteHound [121], and comprehensive consensus-based tool MetaPocket 2.0 [122]. By combining eXtreme gradient boosting (XGBoost) and graph convolutional neural networks (GCNNs), a robust approach for allosteric site identification and Prediction of Allosteric Sites Server (PASSer) was developed [123,124,125]. We also employed the PASSer Learning to Rank (LTR) model, which is capable of ranking pockets in order of their likelihood to be allosteric sites [112]. Using P2Rank [118,119] and PASSer LTR [125] approaches, we identified binding pockets in the conformational ensembles and computed P2Rank-predicted residue pocket probability. The reported top binding pockets for each protein structure correspond to top-ranked consensus P2Rank/LTR predicted sites. By expanding previous approaches for pocket detection in the S proteins that were based on static structures [126,127], the employed in this study approach includes the dynamics of the S structures and monitors the occurrence and stability of identified pockets in the conformational ensembles. The pocket detection analysis is performed on MD and CG-BD-generated simulation trajectories. A total of 1000 ensemble conformations sampled at regular intervals along the trajectories were used for analysis of pocket occurrence and stability. This allows for the identification of novel pockets during dynamics and allows the monitoring of emergence, separation or fusion of pockets occurring in transient and allosteric pockets. In addition to employment of P2Ranl/LTR consensus pocket detection that considers ligandability/druggability criteria and scores, we used several additional druggability assessment tools to validate the predictions. For this additional analysis, we used DoGSiteScorer (a web server for automatic binding site prediction, analysis and druggability assessment) [128] and PockDrug tool that provides a prediction of the druggability score which, if greater than 50%, indicates a druggable pocket [129,130]. The pockets that are formed by a minimum of 8 amino acids and defined with the ligandability/druggability score of at least 2.0 were considered based on the premise that a minimal pocket size and sufficient druggability may be required to ensure the productive ligand binding with the protein target. The consensus dynamics-based analysis of pocket frequency, variability and druggability allowed the identification of promising cryptic druggable pockets in the studied SARS-CoV-2 S proteins and complexes.

3. Results

3.1. Atomistic MD Simulations and Markov State Model Analysis of the Omicron BA.2, BA.2.75 and XBB.1 RBD-ACE2 Complexes

We first conducted an analysis of the conformational dynamics and distribution of states in the Omicron BA.2, BA.2.75 and XBB.1 RBD-ACE2 complexes (Figure 1, Supplementary Materials, Figure S1) using microsecond atomistic MD simulations combined with subsequent dimensionality reduction and MSM analysis of the conformational landscapes. Here, we considerably expanded our recent simulations of the Omicron RBD-ACE2 complexes [131] by performing a systematic comparative analysis using multiple microsecond MD simulations based on the experimental structures for all studied Omicron RBD-ACE2 complexes (Figure 1), which included the recently released structure of the XBB.1 RBD-ACE2 complex [31]. The receptor-binding motif (RBM) is the main functional motif in RBD and is composed of the RBD residues 437–508 that form the binding interface between the RBD and ACE2 (Supplementary Materials, Figure S2). The recognition loop residues (residues 470–491) of the RBM form a dense interaction network with the ACE2 receptor that determines the RBD-ACE2 binding affinity.

We performed three independent 1 µs MD simulations for each of the Omicron RBD-ACE2 complexes. The RMSF profiles of three independent MD simulations for the Omicron BA.2, BA.2.75 and XBB.1 RBD-ACE2 complexes are shown in the Supplementary Materials, Figure S3. As may be expected, the individual trajectory profiles revealed a common shape and considerable similarities in the stable RBD core region while displaying variations in the flexible regions (Supplementary Materials, Figure S3). For clarity of presentation and discussion, we also combined the results from three individual simulations and reported the average RMSF profiles for the S-RBD residues in the Omicron BA.2, BA.2.75 and XBB.1 RBD-ACE2 complexes (Figure 2A).

The RMSF profiles confirmed the stability of the β-sheet RBD core region (residues 350–360, 375–380, 394–403) (Figure 2A) while revealing characteristic RMSF peaks corresponding to the flexible RBD regions (residues 360–373 and residues 380–396) (Figure 2A). In general, the RMSF profiles showed that thermal displacements were observed in atomistic simulations of the BA.2, BA.2.75 and XBB.1/XBB.1.5 RBDs are similar but also highlighted greater fluctuations of the flexible RBD regions in the XBB.1 RBD-ACE2 complex (Figure 2A). We observed some differences in the flexible RBD regions (residues 355–375, 380–400, 480–490), also pointing to the greater local mobility of the RBM residues in the XBB.1 RBD (Figure 2A). While the flexible RBD loop (residues 440–452) contains several convergent mutational sites present in the XBB.1 variant (K444, V445, G446, N450, L452), we observed no appreciable differences in its mobility between Omicron variants. It may be argued that moderately increased RBM mobility in XBB.1 may result from F486S and F490S substitutions that could weaken the RBM-ACE2 interfacial contacts and allow for greater flexibility in this region (Figure 2A,B). The RBM tip is structurally ordered (“hook-like” fold) in both BA.2 and BA.2.75 conformational samples throughout the MD trajectory, while F486S mutation in the XBB.1 RBD-ACE2 complex may induce partly disordered configurations of the RBM tip (Figure 2E).

The RMSF analysis of the ACE2 residues showed similar profiles across variants (Supplementary Materials, Figure S4). The stable ACE2 residues are in the core and the binding interface positions. The key ACE2 α-helix (residues 24–31) and a β-sheet (residue 350–356) show small RMSF values in all complexes (Supplementary Materials, Figure S4).

Another useful metric of the local conformational dynamics patterns is the distance fluctuation stability index, which evaluates the extent of fluctuations in the mean distance between each residue and all other protein residues [86,87,88,101]. The high values of distance fluctuation stability indexes are associated with structurally rigid residues as they display small fluctuations in their distances to all other residues, while small values of this index are characteristic of flexible sites [86,87,88]. The distributions predicted the local maxima for structurally stable regions in the RBD core (residues 400–420, 450–460) and the RBD binding interface region (residues 484–505) (Figure 2B). The differences in the stability index profiles are mostly in the RBM region (residues 470–491), showing greater stability for the BA.2.75 RBD, while smaller indexes are seen in the XBB.1 RBD profile in support of the greater mobility of this region (Figure 2B). The conformational dynamics profiles projected onto the RBD-ACE2 structures showed a similar and strong stabilization of the RBD core, RBD binding interface and the interfacial helices on ACE2 (Figure 2C–E). Although the conformational dynamics profiles revealed certain important differences between Omicron variants pointing to more flexible XBB.1 RBD and confirming the stability of the BA.2.75 RBD, these metrics are poorly suited to distinguish salient features of the underlying landscapes and provide a quantitative characterization of the dynamic state equilibrium.

To clarify key differences between conformational landscapes of the studied Omicron RBD-AC E2 complexes and identify the distribution of states, we used the tICA dimensionality reduction method [96,97,98,99] by projecting the results of MD simulations onto low dimensional space. Using this approach, we can reduce a multidimensional complex set of variables to a lower dimension and extract information from sampled conformations by determining the slowest degrees of freedom. The low-dimensional projection of the MD ensemble obtained from microsecond trajectories of Omicron BA.2 RBD-ACE2 complex (Figure 3A) showed a fairly broad distribution along one of the slowest degrees of freedom. In sharp contrast, the two-dimensional density for the BA.2.75 RBD-ACE2 complex is extremely narrow and localized in a single basin corresponding to the crystallographic conformation that dominates the thermodynamic equilibrium for this variant (Figure 3B). We observed a much broader density distribution for the XBB.1 RBD-ACE2 complex (Figure 3C) where variations are seen along both slowest degrees of freedom. Interestingly, although both BA.2 and XBB.1 RBD-ACE2 complexes displayed more variability as compared to structurally rigid BA.2.75 complex, their low-dimensional maps are distinct. As a result, the recombinant XBB.1 variant may induce a specific distribution of favorable conformations. Moreover, each of the Omicron variants displayed a specific low-dimensional map, as evident from the combined representation of densities for all RBD-ACE2 complexes (Figure 3D), showing only partial overlap between conformational ensembles.

Together, the low-dimensional projection analysis emphasized the key differences in the conformational dynamics of the Omicron variants, revealing the existence of multiple functional RBD conformations in the BA.2 and XBB.1 RBD-ACE2 complexes and the increased RBD plasticity induced by these two variants. In contrast, the striking convergence of three independent microsecond simulations of the BA.2.75 RBD-ACE2 complex confirmed structural stabilization of the RBD for this Omicron variant. These findings corroborate with the experimental data showing a more rigid and compact BA.2.75 RBD structure as compared to other Omicron variants [25].

Using low-dimensional characterization of the conformational ensembles, we proceeded with the MSM analysis and partitioning of the MSM macrostates (Figure 4). For clustering and MSM characterization, we employed conformational ensembles of all studied Omicron RBD-ACE2 complexes. Notably, the conformational ensemble for each complex is obtained by combining three independent 1 µs MD simulations. The entire data set used for partitioning combined multiple MD trajectories of all studied complexes, resulting in a total of nine macrostates (Figure 4). The stationary distribution of states was evaluated for BA.2 (Figure 4A), BA.2.75 (Figure 4B) and XBB.1 RBD-ACE2 complexes (Figure 4C). For the BA.2 RBD-ACE2 ensemble, macrostates 2, 3, 5, 6 and 9 contribute to the distribution, with macrostate 9 dominating ~ the conformational population (Figure 4A). This microstate corresponds to the crystallographic conformation of the complex, showing that other RBD conformations (macrostates 5, 6) showing variations in the flexible RBD loops could appreciably contribute to the equilibrium and binding with ACE2. A single peak distribution in the BA.2.75 RBD-ACE2 complex confirmed that this ensemble is overwhelmingly dominated by the native conformation (Figure 4B), suggestive of markedly rigidified RBD structure for this variant. In contrast, we detected a broader distribution of macrostates for the XBB.1 RBD-ACE2 complex (Figure 4C) where macrostates 1, 3, 4, 8 and 9 contribute to the equilibrium population. Nonetheless, the distribution is still dominated by the crystallographic state (macrostate 9), contributing 62%, and yet there are several additional macrostates present in the equilibrium, which supported the notion of the greater XBB.1 RBD plasticity (Figure 4C). By expanding to the microsecond time scale, a combination of MD simulations and MSM analysis provided evidence that despite highly similar RBD-ACE2 structures, the distributions of states for different Omicron variant complexes are considerably different.

These results revealed important differences in the conformational landscapes and distribution of functional conformations in the Omicron variants, still showing that the topology of the RBD fold and RBD-ACE2 binding interface are fully preserved across all variants. As a result, functional and binding differences between BA.2, BA.2.75 and XBB.1 variants can be determined through moderate variations in conformational flexibility that enable modulation and balancing of multiple fitness tradeoffs between immune evasion, ACE2 affinity and conformational adaptability, which may potentially be an important driver behind the evolution of Omicron variants.

The transition probabilities were determined among different macrostates for all systems, with the 5 ns lag time. The high percentage of self-conserved probability shows the stability of macrostates. The results showed significant and informative differences in the distributions and transitional probability maps between these structurally similar RBD-ACE2 complexes. The transitional probability map for the BA.2 RBD complex revealed that the macrostates 9 and 6 can interconvert, representing two major dense areas (Figure 4D), indicating that the equilibrium population for this variant is determined by these RBD conformations. Structural inspection of the microstates 9 and 6 (Supplementary Materials, Figure S5) showed significant differences in the flexible regions, particularly in the RBM tip region that adopted an alternative orientation as compared to the crystallographic conformation. Hence, BA.2 RBD maintains considerable plasticity in the complex with ACE2, suggesting that the entropic component may contribute to the binding affinity of this complex. Remarkably, the transition map for the BA.2.75 RBD-ACE2 complex showed only a single conformation that is virtually identical to the crystallographic state (Figure 4E). In sharp contrast, transition maps for XBB.1 (Figure 4F) revealed that dominant macrostate 9 can interconvert with macrostates 1 and 3, where microstate 3 is connected with other identified macrostates of the ensemble. Structural analysis of the macrostates showed a considerable remodeling of the flexible RBD region (residues 460–490), including different orientations of the RBM tip (Supplementary Materials, Figure S5). The MSM analysis also revealed that dominant macrostate 9 featured well-ordered and stable “hook-like” conformation of the RBM tip. This state becomes less favorable in XBB.1 (F486S) while the contribution of other macrostates with a partly disordered conformation of the RBM tip and more dynamic RBD loops increased. XBB.1 subvariant is a descendant of BA.2 and recombinant of BA.2.10.1 and BA.2.75 sublineages, also featuring a group of specific RBD mutations (G339H, R346T, L368I, V445P, G446S, N460K, F486S, F490S and reversed R493Q) (

Table 1. Mutational landscape of the Omicron BA.2, BA.2.75 and XBB.1 Spike Variants.

Omicron Variant	Mutational Landscape
BA.1	A67, T95I, G339D, S371L, S373P, S375F, K417N, N440K, G446S, S477N, T478K, E484A, Q493R, G496S, Q498R, N501Y, Y505H, T547K, D614G, H655Y, N679K, P681H, N764K, D796Y, N856K, Q954H, N969K, L981F
BA.2	T19I, G142D, V213G, G339D, S371F, S373P, S375F, T376A, D405N, R408S, K417N, N440K, S477N, T478K, E484A, Q493R, Q498R, N501Y, Y505H, D614G, H655Y, N679K, P681H, N764K, D796Y, Q954H, N969K
BA.2.75	T19I, G142D, K147E, W152R, F157L, I210V, V213G, G257S, G339H, S371F, S373P, S375F, T376A, D405N, R408S, K417N, N440K, G446N, N460K, S477N, T478K, E484A, Q498R, N501Y, Y505H, D614G, H655Y, N679K, P681H, N764K, D796Y, Q954H, N969K
XBB.1	T19I, V83A, G142D, Del144, H146Q, Q183E, V213E, G252V, G339H, R346T, L368I, S371F, S373P, S375F, T376A, D405N, R408S, K417N, N440K, V445P, G446S, N460K, S477N, T478K, E484A, F486S, F490S, R493Q reversal, Q498R, N501Y, Y505H, D614G, H655Y, N679K, P681H, N764K, D796Y, Q954H, N969K

). Many of these RBD mutations are known for their immune evasion functions, including R346T, G446S and F486S. Based on the structural comparison of the macrostates, it becomes apparent that a combination of F486S and F490S mutations in XBB.1 may induce increased RBD mobility in regions important for antibody binding and thereby enhance the immune evasion potential of this Omicron variant. The MSM analysis suggested that BA.2 and XBB.1 RBDs may exhibit a similar degree of conformational plasticity in the complex, while BA.2.75 RBD is characterized by structural rigidity featuring a single dominant microstate.

3.2. Coarse-Grained Simulations of the Omicron BA.2, BA.2.75 and XBB.1 S Trimer-ACE2 Complexes

CG-BD simulations of the BA.2 S trimer with two human ACE2 bound, BA.2 S trimer with three ACE2 bound, BA.2.75 S trimer with one human ACE2 bound and XBB.1 S trimer with one human ACE2 bound [31] showed the appreciable level of thermal fluctuations in the generally adaptable S1 subunit (residues 14–530), including NTD (residues 14–306) and moderate displacements of the RBD-up regions (residues 331–528) that are stabilized via binding interactions with ACE2 and a more rigid S2 subunit (residues 686–1200) (Figure 5). We observed markedly larger fluctuations in the NTD and RBD regions for the XBB.1 1RBD-up complex where both ACE2-bound RBD-up and two closed RBDs displayed appreciable mobility (Figure 5D) as compared to more stable S1 regions in the BA.2 (Figure 5A,B) and BA.2.75 trimer complexes (Figure 5C). We also observed differences in the NTD mobility, including NTD N-terminus (residues 14–20), N3 (residues 141–156) and N5 (residues 246–260) that form an antigenic supersite on the NTD. These regions showed increased RMSF values for all variants, but conformational heterogeneity in the NTD was particularly apparent in the BA.2 and XBB.1 complexes (Figure 5A,B,D). Three major fluctuation peaks were observed in the Wuhan NTD, residues 62–83, 140–158, 177–189 and 239–260, corresponding to the loop regions N2, N3, N4 and N5, respectively. The RMSF values for the ACE2-interacting RBD residues 440–456, 470–491 and 491–505 were small in the BA.2 and BA.2.75 complexes (Figure 5A–C), but these fluctuations were larger in the XBB.1 complex (Figure 5D). Conformational flexibility of the NTD and RBD in the open BA.2.75 trimer complexes is reduced, displaying smaller displacements within RMSF < 2.0–2.2 Å (Figure 5C). Hence, ACE2 binding may induce stabilization of the BA.2.75 S trimer by affecting not only the interacting RBD regions but also allosterically promoting the enhanced inter-protomer packing and rigidification of the S2 subunit. The observed dynamics pattern in the BA.2.75 S trimer complex is consistent with the experimental findings that the BA.2.75 S-trimer was the most stable, followed by BA.1, BA.2.12.1, BA.5 and BA.2 variants, also featuring the increased RBD rigidity and enhanced ACE2 binding affinity compared with other Omicron variants [25].

In the Omicron BA.2 S trimer-ACE2 complexes, the RMSF profiles showed only relatively moderate fluctuations in the C-terminal domain 1, CTD1 region (residues 528–591) but more significant variations in the C-terminal domain 2, CTD2 (residues 592–686) (Figure 5A,B, Supplementary Materials, Figure S6). Strikingly, we observed greater stability of both S1 and S2 regions in the BA.2.75 S trimer complex (Figure 5C), where NTD regions and RBD-up regions exhibited RMS displacements within 2.0–2.2 Å. The conformational dynamics profiles revealed larger deviations and greater flexibility for CTD2 (residues 592–686), particularly near the furin cleavage region (residues 670–690) (Figure 5). The interesting patterns of conformational fluctuations were observed in the functionally important S2 regions including fusion peptide proximal region (FPPR) (residues 828–853), immediately downstream of the fusion peptide, the upstream helix (UH) (residues 736–781), central helix (CH) (residues 986–1035) are particularly rigid, UH (residues 736–781), heptad repeat HR1 (residues 910–985), central helix CH (residues 986–1035), CD, connector domain (residues 1035–1068) and heptad repeat HR2, heptad repeat 2; (residues 1069–1163) (Figure 5, Supplementary Materials, Figure S6). In all complexes, with the exception of BA.2.75 complex, the RMSF displacements are significant near the S1/S2 cleavage site and in the FPPR region, while uniformly across all trimers, only small deviations were seen in structurally stable UH and CH regions (Figure 5). Consistent with strong sequence conservation, the CH region is rigid in all S trimer complexes. We found that ACE2 binding may promote stabilization of the trimeric interface stalk region of S2 (residues 899–913, 988–998, 1013–1021) (Figure 5, Supplementary Materials, Figure S6). The stalk immobilization is particularly noticeable in the BA.2 and BA.2.7 complexes (Figure 5A–C). Only moderate displacements were seen in the HR1 region (residues 910–985) for BA.2 and BA.2.75 trimers, but the flexibility in this region increases in the XBB.1 complex. In this context, it may be worth noting the results of HDX-MS experiments suggesting partial destabilization not only near the S1/S2 cleavage site but also in the HR1 region, which may be functionally required to promote fusion stage. Our results showed considerable stability of the highly conserved hinge epitope (residues 980–1010) located in the S2 region between the CH and HR1 segments that are involved in regulation of conformational changes required for fusion of the viral envelope [46].

3.3. Mutational Scanning of the Binding and Inter-Protomer Interactions in the Omicron BA.2, BA.2.75 and XBB.1 S Trimer-ACE2 Complexes

To provide a systematic residue-based mutational analysis of binding and stability in the Omicron S trimer-ACE2 complexes, we conducted mutational scanning of the RBD interface residues to evaluate ACE2 binding interactions and comprehensive scanning of the inter-protomer positions to estimate the packing interactions and differences in stability between Omicron trimers. Structure-based analysis of the Omicron RBD0ACE2 complexes showed considerable similarities, as the differences in the number of intermolecular contacts were fairly small (Supplementary Materials, Tables S1–S3). We first constructed mutational heatmaps for the RBD binding interface residues (Figure 6). Unlike previous studies, these maps are based on scanning the RBD residues in the full-length S-ACE2 complexes with a single or two ACE2 molecules present. In general, the results are in excellent agreement with the experimental deep mutagenesis scanning (DMS) data, revealing a consistent group of predominantly hydrophobic energy hotspots (Figure 6). Indeed, common energetic hotspots Y453, L455, F456, F486, N487, Y489, Y501 and H505 also emerged as critical stability and binding hotspots in the experimental DMS studies [132,133]. Mutational heatmaps of the RBD-ACE2 binding interactions for all Omicron variants showed that the majority of substitutions in these key interfacial positions could lead to a loss in the stability and binding affinity with ACE2 (Figure 6).

Hence, the conserved and stable hydrophobic RBD hotspots that are critical for RBD stability may be universally important for binding. Not surprisingly, Omicron mutations typically target more dynamic and vulnerable positions in the RBD to optimize virus fitness and balance the ACE2 binding affinity with immune escape potential. Interestingly, the mutational heatmaps for BA.2 trimer complexes with 2 ACE2 molecules (Figure 6A) and three bound ACE2 proteins (Figure 6B) showed considerable destabilization of the interfacial RBD residues where multiple positions (Y449, Y453, L455, F456, Y473, F486, N487, Y489, R493, R498, T500 and Y501) contributed significantly to the ACE2 binding. Compared to BA.1, variant BA.2 additionally contains S371F, T376A, D405N and R408S mutations (Table 1). However, these RBD positions are outside of the binding interface, and corresponding mutations may exert their effect through couplings with the interfacial sites. Overall, the structure of the S-BA.2 trimer complex with ACE2 revealed a more extensive interaction network at the binding interface and, combined with the stability of the BA.2 RBD, could contribute to the high binding affinity of the BA.2 variant. BA.2.75 variant has additional mutations in the RBD (D339H, G446S, N460K and R493Q).

The mutational heatmap showed an increase in tolerance for positions R403, Y449, S493 and R498, where mutations result in moderate changes (Figure 6C). At the same time, the key hydrophobic sites Y453, L455, F456, F486, N487, Y489 as well as Y501 are preserved as binding hotspots in the BA.2.75 complex with a single ACE2 bound molecule. Previous studies suggested that R493Q reversal could restore the lost receptor-binding affinity due to N460K mutation [134]. In the BA.2 RBD/hACE2 complex, both R493 of RBD and K31 of human ACE2 are positively charged, which could decrease their binding, but a highly favorable salt bridge was observed between R493 of RBD and E35 of ACE2. Biophysical studies argued that R493Q substitution could improve the binding affinity towards human ACE2 [135]. At the same time, structural studies of the BA.2 trimer binding complexes with ACE2 used in our analysis suggested that Q493R forming a new salt bridge with E35 of human ACE2 and Q498R forming a new salt bridge with D38 of human ACE2 are key mutations resulting in an enhanced binding to human ACE2 [22]. This may explain the emergence of Q493R mutation in the first place in the BA.1 and BA.2 Omicron variants (Table 1). Our results showed that even though mutations of Q493 in BA.2.75 typically result in a significant loss of binding affinity, Q493R substitution may not cause a dramatic decrease in binding (Figure 6C). Based on our findings, the increased binding potential of BA.2.75 may be a cumulative effect also resulting from the increased stabilization of the RBD and S1–S2 regions.

Strikingly, the mutational heatmap for the XBB.1 trimer complex revealed a smaller number of binding hotspots (Y453, F456, Y489 and Y501), while other positions, including L455, Y473, A475, G476, N477, S486 and S490, are more tolerant to substitutions that induce only small free energy changes (Figure 6D). These findings indicated that the increased flexibility of the XBB.1 trimer and binding interface combined with F486S and F490S mutations can result in decreased binding efficiency. Indeed, DMS studies have shown that F486P might enhance the affinity to ACE2 compared with F486S [103]. Moreover, the binding affinity of the BA.2.75 (K_D 1·8 nM) is much stronger than that of XBB.1 (K_D 19 nM) [32]. Mutational scanning analysis of the XBB.1 S-ACE2 binding is consistent with these experiments, revealing a more adaptable and mutation-tolerant interface, which may also explain the significant immune escape potential of this variant traded in exchange for reduced ACE2 binding.

In addition, we also performed a systematic mutational scanning of the inter-protomer interfacial residues in the Omicron S trimer complexes to evaluate differences in the trimer packing and identify important structural stability hotspots (Figure 7). Strikingly, mutational analysis revealed clear and important differences between the variants. We found that the inter-protomer interaction contacts in the BA.2 trimer are generally sensitive to mutations, revealing strong destabilization changes in the S1 and S2 regions, including peaks aligned with NTD (residues F135, W152, F168, Y170, G199, Y200, P230, G257 and CTD1 interfaces (residues 542–547) (Figure 7A). Some of these residues, W152 and G257, play a key role in NTD binding with antibodies. The analysis showed that the inter-protomer NTD and RBD positions may form small clusters, but the density of these interfacial sites is smaller in BA.2 as compared to BA.2.75 trimer (Figure 7B).

The analysis of mutational scanning showed large destabilization changes induced by mutations in the inter-protomer β-sheet segments involving residues 701–705 on one protomer and residues 780–790 from the other protomer (Figure 7A). These inter-protomer hotpots are preserved in all Omicron trimers. Interestingly, this region was identified in structure-guided mutagenesis to facilitate allosteric communication between NTD ligation and proteolytic fusion activation in the S2 region [136]. A spectrum of destabilization changes is seen in the S2 regions (residues 910–1150) (Figure 7A), reflecting the S2 inter-protomer interfaces while also showing some degree of conformational plasticity and energetic tolerance in the S2. A significantly denser pattern of inter-protomer interface sites was seen in the BA.2.75 trimer complex (Figure 7B), revealing a number of inter-protomer residues in the NTD and RBD regions. These findings indicate that the inter-protomer contacts in BA.2.75 are stronger than in both BA.2 and XBB.1 trimers, which is consistent with the increased structural stability of this Omicron variant. This may also have implications for the detection and distribution of cryptic pockets. We argue that the greater conformational plasticity of the BA.2 trimer as compared to the more densely packed BA.2.75 trimer would affect the evolution and emergence of dynamic cryptic sites. These results are consistent with the notion that Omicron open trimer structures are packed more tightly with enhanced inter-domain and inter-subunit interactions as compared to the Wu-Hu-1 WT S protein as each of the RBD-down protomers packs more tightly with its neighboring NTD [19]. In this mechanism, Omicron mutations promote a more stable open conformation. Our analysis showed that the inter-protomer stabilization progressively increased from BA.2 to BA.2.75 variant.

Strikingly, mutational scanning of the inter-protomer contacts in the XBB.1 trimer showed a distinct pattern, revealing fewer inter-protomer sites and visibly sparse distribution of these sites in the NTD and RBD (Figure 7C). Nonetheless, the distribution displayed several stability hotspots across all S regions, including NTD (residues 38, 43, 47, 135, 168, 199, 200, 225, 232, 235), RBD (residues 396, 413), CTD1 (residues 560–563), (residues 1047, 1079). These results suggested that the XBB.1 trimer may be less densely packed and has a smaller inter-protomer interface, which may contribute to the reduced stability and greater plasticity of the XBB.1 variant. It may be argued that through enhanced plasticity, the XBB.1 variant may improve virus adaptability and enhance immune evasion potential. This may affect the distribution of cryptic pockets that could be more broadly distributed across both S1 and S2 regions.

3.4. Detection of Cryptic Binding Pockets in the Conformational Ensembles of S Omicron Trimer Complexes with ACE2: The Effects of Binding and Structural Plasticity in Mediating Networks of Conserved Allosteric Sites

Using the conformational ensembles of the Omicron S RBD-ACE2 complexes an S trimer-ACE2 complexes for BA.2, BA.2.75 and XBB.1 variants, we performed a large-scale cryptic pocket detection and comparative analysis of the putative binding sites in these systems. In particular, we investigated how variant-induced dynamic changes and ACE2 binding can affect the residue-based propensities for ligand binding (often termed ligandability preferences) and the distribution of cryptic pockets. Through this analysis, we examined the role of Omicro variant evolution and the impact of recombinant Omicron sublineages from BA.2 to BA.2.75 and XBB.1 on the stability and population of cryptic binding sites. Here, we particularly focus on identifying conservation patterns and notable divergencies in the distribution of cryptic pockets for studied Omicron variants. We will investigate and validate a hypothesis that functionally important for allostery cryptic pockets formed at the inter-domain and inter-protomer interfaces in the S trimer complexes with ACE2 may be protected through evolution and shared between Omicron variants. There is an increasing interest in identifying druggable sites in the S2 subunit as most recently developed broad-spectrum fusion inhibitors, and candidate vaccines can target the conserved elements in the S2 subunit [54].

First, we compared the residue-based pocket propensity profiles (Figure 8A–C) and the predicted RBD cryptic pockets in the RBD-ACE2 complexes (Figure 8D–F). The pocket propensity profiles for the BA.2 RBD (Figure 8A) featured a major dense peak in the RBD core (residues 360–390) and two small peaks (residues 426–432 and 510–515). The pocket detection analysis revealed two highly probable pockets. One pocket is formed by residues Y365, S366, Y369, F377, V382, S383, P384, L387, N388, F392, C432 and F515, while pocket 2 (residues F342, N343, A344, V367, L368, N370, A372, F374, W436, N437, L441 and R509) is immediately adjacent to pocket 1 (Figure 8D). In general, the predicted pockets corresponded to a conserved and experimentally validated allosteric site in the RBD core where the essential free fatty acid LA binds [39,40,41]. According to the cryo-EM structure of the S complex with this small molecule, the allosteric pocket is lined up by residues F338, V341, F342, F377, F374, F392 and W436 [39], which is consistent with our predictions. Notably, in the experimental S trimer complex, LA binding occurs in a bipartite binding pocket in which the hydrophobic tail of LA molecule is bound to the RBD on one protomer, and the carboxyl headgroup of LA interacts with R408 Q409 of the neighboring RBD [39,40,41]. Interestingly, our analysis revealed that this cryptic pocket is highly probable and conserved even in the absence of the adjacent protomers. A small pocket on the other side of the BA.2-RBD is minor and much less probable. As a result, the population of druggable RBD pockets in the BA.2 complex is dominated by the experimentally known LA allosteric site. A similar pocket propensity distribution in the BA.2.75 RBD-ACE2 complex displayed a dominant dense cluster peak for RBD residues 360–390 and smaller peaks for RBD core 430–440 and ACE2 interface positions 485–490 (Figure 8B). The most probable pocket (residues 338, 339, 342, 343, 367, 368, 371, 372, 374, 436, 437, 441) is aligned precisely along the LA allosteric site on the RBD (Figure 8E). The emergence of a single conserved allosteric pocket tin the BA.2.75 RBD may partly reflect the structural rigidity of BA.2.75 and the structural conservation of the major allosteric site. While the pocket propensity distribution for XBB.1 RBD is generally similar, we noticed the appearance of multiple isolated peaks (residues 367–371, 377–379, 432–434, 513–515), which may be due to the increased conformational plasticity of XBB.1 RBD (Figure 8C). The analysis revealed three equally probable smaller pockets (pocket 1—residues 365, 368, 369, 377, 387, 432, 434, 513, 515; pocket 2—residues 343, 344, 345, 347, 371, 372, 436, 441, 509; and pocket 3—residues 363, 365, 386, 387, 390, 391, 392, 395, 515, 524, 525) (Figure 8F). The three RBD cryptic pockets are adjacent to each other and cover a much larger portion of the RBD surface when combined (Figure 6F). Although these pockets overlay with the experimental allosteric site, we found a partial fragmentation of the cryptic site, leading to the formation of several smaller pockets. Hence, our analysis suggested that variant-induced modulation of conformational plasticity may result in structural variations and expansion of the cryptic RBD site while the location and hydrophobic nature of this allosteric pocket is preserved in all complexes. We argue that this may have implications for small molecule binding in the allosteric site, likely reducing the efficiency of ligands targeting the RBD.

We then used conformational ensembles obtained from molecular simulations of the S trimer-ACE2 complexes to identify the available spectrum of potential cryptic sites. By ranking the pockets using a residue-based propensity score for ligand binding, we systematically analyzed the distribution and functional significance of these cryptic sites. The reported pocket propensity scores represent ensemble-averaged values obtained from P32Rank calculation based on the equilibrium trajectories (Figure 9). Here again, despite structurally similar structural arrangements of the open S trimer complexes, we observed significant differences between pocket propensity profiles. A broad distribution with dense clusters of residues featuring peaks in both S1 and S2 subunits was seen for the BA.2 S-ACE2 complex (Figure 9A). Notably, one could differentiate three major peaks, with one corresponding to residue clusters in the RBD core, another cluster at the CTD1 region (residues 528–591), and a pronounced peak in the S2 subunit pointing to the CH region (residues 986–1035). Smaller peaks were also observed in the UH region (residues 736–781) and HR1 (residues 910–985). These findings suggested that the BA.2 S trimer could retain considerable conformational plasticity in the ACE2-bound complex and feature a significant number of cryptic pockets in both S1 and S2 regions. Interestingly, the pocket propensity profile changes in the BA.2 S trimer complex with three ACE2 molecules bound to RBDs in the up position (Figure 9B). In this ensemble, we observed several sharp, narrow peaks in the RBD and inter-domain CTD1 regions, whereas smaller density was seen in the S2 regions. A different distribution profile was obtained for the BA.2.75 1 RBD-up trimer complex (Figure 9C), in which pronounced sharp peaks are aligned with the NTD and RBD regions while the pocket propensities of residues in the S2 subunit are markedly reduced. Accordingly, this may indicate that diverse cryptic pockets in the BA.2.75 trimer complex may be localized in the NTD and RBD, while the tightly packed S2 subunit may harbor a smaller number of conserved binding pockets. Notably, the pocket propensity profile of the XBB.1 1RBD-up trimer complex (Figure 9D) is similar to BA.2 2RBD-up trimer (Figure 9A) featuring multiple cluster peaks in the NTD, RBD, S1/S2 cleavage region and stalk region of the S2 subunit.

It is interesting to relate our findings with the HDX-MS experiments showing that ACE2 binding may result in stabilization of the CH (residues 986–1035) and CD regions (resides 1035–1068), including S2 stalk segments (regions 899–913, 988–998, 1013–1021) while induce the increased mobility in the S1/S2 cleavage site and HR1 domain [36,37]. According to our predictions, the high pocket propensity regions are generally distributed in the dynamic NTD and RBD regions as well as in the S2 subunit. For BA.2 and XBB.1 trimer complexes, we observed a dense, broad peak corresponding to HR1 and CH regions (Figure 9A,D), suggesting a potential for the emergence of cryptic pockets in these conserved S2 regions. Overall, the results revealed that cryptic binding pockets may emerge along the S trimer architecture in the S1 regions and conserved S2 segments.

By using the distribution-based ranking of the binding pockets, we characterized and classified the probable cryptic pockets in the Omicron complexes. In particular, we focused on the identification of evolutionary conserved and stable sites that may be shared among variants as well as potentially variant-specific unique pockets. The structural analysis of the top-ranked binding pockets in the S trimer complexes revealed conserved allosteric pockets shared between Omicron variants while displaying considerable differences in the composition, size and distribution in the S1 and S2 regions (Figure 10). The distribution of the top five predicted cryptic pockets in the Omicron S trimer complexes showed that they predominantly occupy the NTD and NTD-RBD regions. However, we also detected top-ranked binding pockets in the stalk region of the S2 trimer interface, which emerged in all three complexes (Figure 10). Interestingly, among the top cryptic pockets are the inter-domain binding sites seen in the BA.2 and BA.2.75 trimer complexes (Figure 10A,B) but not found in the XBB.1 complex (Figure 10C).

Overall, the predictions are consistent and expand previous studies of cryptic binding pockets in the SARS-CoV-2 S structures, including the RBD-ACE2 complexes and S trimers [68,69,126]. In agreement with these studies and experimental data, we identified several binding pockets in the RBD, including conserved and experimentally validated allosteric sites where the essential free fatty acid LA binds [39,40,41]. In addition, our pocket analysis of the S trimers captured a bipartite binding pocket formed by two protomers where LA molecule can bind in with the hydrophobic tail of LA molecule is bound to the RBD on one protomer, and the carboxyl headgroup of LA interacts with the neighboring RBD [39,40,41].

The details and amino acid composition of the top 10 predicted stable binding pockets, along with the ligandability/druggability scores and occurrence frequencies in the course simulations, are reported for the Omicron BA.2 S trimer-ACE2 complex (Table S4), Omicron BA.2.75 S trimer-ACE2 complex (Table S5) and Omicron XBB.1 S trimer-ACE2 complex (Table S6). Structural mapping of the detected pockets (Figure 10, Supplementary Materials, Figures S7–S9) revealed a general correspondence between the distribution of pockets in the original Wu-Hu-1 S trimers and Omicron variant S trimers which is consistent with structural and dynamic similarities of the S trimer structures. However, while the critical inter-protomer binding pocket at the S1–S2 hinge region is conserved across Omicron variants, we found that the composition and size of this critical site can change (Figure 10, Supplementary Materials, Figures S7–S9). Indeed, while in the BA.2 trimer complex, this pocket incudes about 14 residues (B_740 B_741 B_744 B_745 B_855 B_856 B_966 B_976 C_548 C_549 C_571 C_572 C_573 C_589) with the ligandability score of 3.7 (Supplementary Materials, Table S4) this cryptic is considerably enlarged in the BA.2.75 S trimer complex with ACE2 and consists of 22 residues (C_541 C_546 C_547 C_548 C_549 C_570 C_572 C_573 C_587 C_589 C_592 D_1000 D_740 D_741 D_744 D_745 D_855 D_856 D_966 D_976 D_977 D_978) with a considerably greater pocket score of 8.95 (Supplementary Materials, Table S5). Remarkably, this pocket is markedly transformed and engages different and fewer residues on adjacent protomers (B_733 B_735 B_768 B_771 B_772 B_775 B_861 C_312 C_314 C_596 C_613 C_666) in a more dynamic XBB.1 S trimer complex. Moreover, the newly reshaped inter-protomer pocket is less stable, with a pocket score of 1.12 (Supplementary Materials, Table S6).

As modulation of this inter-protomer interface can shift the equilibrium in the S protein and affect the population of closed and open states [2], our findings suggested that this functionally important site is stable and highly druggable in the more rigid BA.2.75 S trimer complex and can be exploited for potential therapeutic intervention, while the increased plasticity of the XBB.1 S trimer may eliminate this favorable cryptic pocket and make ligand-targeted modulation of this hinge region far more challenging. We argue that the recombinant Omicron variant may exploit the greater flexibility and adaptability of the S trimer to facilitate immune escape and also potentially reduce the landscape of feasible cryptic sites for potential targeting and modulation by small molecules.

To provide a systematic characterization of predicted pockets and identify conserved and variant-specific binding sites, we embarked on a more detailed structure-functional analysis. First, we found that the distribution of the pockets in the BA.2 complex is consistent with the conformational adaptability of the S-BA.2 trimer featuring pockets in the NTD-RBD regions, the inter-protomer S1–S2 regions and also in the S2 subunit (Figure 10A, Supplementary Materials, Figure S7 and Table S4). The cryptic pockets formed in the inter-protomer NTD-RBD regions include pocket 1 (NTD:B_115, B_132, B_167, B_168, B_170, B_230, B_231, B_232, RBD: C_357, C_393, C_394, C_518, C_520, C_521, C_523) and pocket 2 (RBD: A_357, A_359, A_360, A_393, A_394, A_520, A_521, A_523, NTD: C_115, C_130, C_168, C_230, C_231) (Supplementary Materials, Figure S7B,C). A close-up view of the predicted NBD-RBD pocket in the up protomer interacting with ACE2 highlighted the exposed nature of this pocket in the open conformation, which is not obstructed by ACE2 binding and interactions with other protomers.

Our analysis suggested that targeting this site may allosterically weaken the RBD-ACE2 interactions even in the RBD-up open form. Our analysis also predicted a cryptic pocket formed at the RBD-NTD interface between the RBD-up protomer bound to ACE2 and NTD of the adjacent protomer (Supplementary Materials, Figure S7C). In this case, RBD residues R357, T393, N394, L518 and T523 from the other side of the RBD-up protomer form a pocket with the NTD residues Q115, E132, T167, F168, P230, I231 and G232 of the neighboring protomer located near the NTD supersite region. Interestingly, both pockets engage RBD residues from the ACE2-interacting up protomer. It is possible that targeting this binding pocket can alter the binding strength of the BA.2 trimer with ACE2 and restrict NTD/RBD movements, thus limiting the conformational adaptability of the RBD-up states. We also found a cryptic binding pocket in the S2 subunit formed by all three protomers in which residues 756, 759, 994, 995, 998 and 1002 from each protomer create a well-defined pocket (Supplementary Materials, Figure S7D). This pocket is formed by conserved residues from the structurally stable UH region (residues 736–781) and CH region (residues 986–1035). Strikingly, the predicted inter-protomer S2 pocket overlaps with the highly conserved, conformational hinge epitope spanning residues 970–1006 of the SARS-2 spike at the apex of the S2 domain [46]. This conserved pocket is shared among Omicron variants and is located at the critical functional hinge region involved in the modulation of S1–S2 opening motions. The experimental studies indicated that S dynamics may impact hinge epitope accessibility [46], and our results showed that this cryptic epitope may be accessible for ligand binding in the S trimer complexes with ACE2.

The recent HDX-MS studies of ACE2-induced allosteric activation of the S protein demonstrated allosteric changes upon ACE2 binding that extend to the hinge region and the top of the central helical bundle of the S2 subunit [38]. The experiments showed that ACE2 binding allosterically perturbs and primes HR1 and regions flanking the S1/S2 cleavage site for cleavage [36,37,38]. Targeting the predicted conserved cryptic pocket in the hinge apex region may potentially interfere with ACE2-induced destabilization of the HR1 helical structure (residues 930–960) that facilitates the pre-to-post-fusion spike conformational change. This pocket is shared between BA.2 and BA.2.75 complexes but was not stable in the XBB.1 complex, suggesting some expansion and thermal breezing of the apex helices preventing three protomers from coming together to form a binding site. Hence, the HR1 region in the S2 subunit was identified in our study as a promising conserved druggable region that may have important Implications and is consistent with structural studies. Since HR1 and HR2 regions interact with each other to form the 6-helical bundle, which is for viral and cell membrane fusion, this region may be an attractive target for HR1- and HR2-derived peptides that can inhibit the fusion process [137,138].

We also found the inter-protomer cryptic pocket formed near the S1–S2 hinge region formed by residues functional residues of one protomer (740, 741, 744, 745, 855, 856, 966, 976) and CTD1 residues from the neighboring protomer (548, 549, 571, 572, 573, 589) (Supplementary Materials, Figure S7E). The importance of these inter-protomer interface residues has been firmly established in experimental studies, showing the quadruple mutant (A570L/T572I/F855Y/N856I) introducing modifications in these positions can shift the equilibrium in the S protein and affect the population of closed and open states [2]. These findings showed that potentially druggable pockets could be formed in the conserved S2 regions of the Omicron trimer complexes. It was also proposed that S2-targeted antibodies and vaccines that are developed to target rigid S2 are likely to be effective for various Omicron variants, provided that evolution would continue to favor stabilization of the S2 stalk regions. As a result, the design of small molecules for specific targeting of these pockets may potentially affect the ACE2 binding and interfere with S activation.

A considerable redistribution of the cryptic pockets was detected in the BA.2.75 complex (Figure 10B), where most of the predicted sites are in the NTD and RBD regions, while only several minor pockets were found in the peripheral S2 regions (Supplementary Materials, Figure S7). The differences in the distribution of cryptic binding pockets between BA.2 and BA.2.75 complexes are consistent with the increased stability of BA.2.75. We found that BA.2.75 S trimer complex featured well-defined and unobstructed RBD allosteric pockets in all protomers (residues 336, 338, 342, 358, 363, 364, 365, 368, 387, 392, 395, 397, 511, 513, 515, 524) (Supplementary Materials, Figure S8B). Given that the binding of small molecules in this LA-binding site is known to stabilize the closed S trimer [39,40,41], our predictions suggest that the RBD allosteric site can be targeted in the ACE2-bound Omicron complex and potentially induce allosteric destabilization of the RBD-up open state. The analysis also revealed several binding pockets in the NTD: pockets formed by residues 101, 104, 117, 119, 121, 128, 170, 172, 175, 190, 192, 194, 201, 203, 227, 228, 229, 231; and a pocket formed by residues 100, 101, 239, 240, 242, 248, 249, 250, 263, 264, 265, 65, 66, 81, 84, 94, 95, 96 (Supplementary Materials, Figure S8C). Some of these NTD residues include residues 245–264 of the NTD supersite loop. The important revelation of this analysis is the detection of several inter-protomer cryptic pockets formed near the S1–S2 hinge region. These pockets bring CTD1 residues on one protomer (positions 541, 546, 547, 548, 549, 570, 572, 573, 587, 589, 592) and S2 residues from the adjacent protomer (positions 1000, 740, 741, 744, 745, 855, 856, 966, 976, 977, 978) (Supplementary Materials, Figure S8D). A similar pocket was also detected among top-ranked sites in the BA.2 S trimer complex. Our previous studies [139,140,141,142] and experimental data [2] showed that hinge residue F592 in one protomer could form the inter-protomer cluster with K854, F855, N856 and T859 hinge sites of the adjacent protomer. This cryptic site provides a fairly small and deep pocket formed by critically important inter-protomer residues. Indeed, positions S591 and F592 have conserved structurally stable hinge sites of collective motions that can modulate both the inter-protomer and inter-domain changes [139,140,141,142]. This binding pocket is located immediately next to the ordered FPPR motif (residues 823–858) that engages in allosteric cross-talk with the RBD regions. Moreover, residues K854, Y855 and I856 are critical sites involved in the inter-protomer contacts, and changes in these positions can modulate the shifts between the open and closed sites. The emergence of viable cryptic pockets formed by these regulatory positions can provide a rational strategy for allosteric targeting and modulation of S activity. The pocket detection also revealed the presence of a conserved cryptic pocket formed between two protomers in the connector domain CD of S2, which links CH and the C-terminal HR2 (Supplementary Materials, Figure S8E and Table S5).

The sensitivity of the cryptic pockets formation to the conformational dynamics is more pronounced for the XBB.1 trimer complex (Figure 10C, Supplementary Materials, Figure S8 and Table S6). Consistent with the greater flexibility of the XBB.1 complex observed in simulations, we found a heterogeneous distribution of the pocket propensity scores and, respectively, a more diverse structural allocation of the predicted pockets (Supplementary Materials, Figure S9 and Table S6). The predicted cryptic pockets are located in the NTD, RBD, CTD1) and S2 regions. Interestingly, the composition of the pockets indicated an increased number of small pockets, as conformational variability may limit the formation of large and rigid pockets in the XBB.1 complex. The observed plasticity of the XBB.1 pockets also indicates that these pockets may readily adapt and alter, which may be evolutionary beneficial for the virus to increase immune evasion and mediate antibody resistance. We detected pockets in the NTD and RBD regions that are similar to the ones seen in other Omicron variants. The divergences in the composition of cryptic pockets for XBB.1 trimer are mainly in S2 regions where the predicted binding sites engage residues from UH, CH and CD functional segments (Supplementary Materials, Figure S9 and Table S6). Interestingly, we did not find the inter-protomer CTD1-S2 hinge pocket in the XBB.1, while this pocket is found as highly probable in more stable BA.2 and BA.2.75 trimer structures.

To summarize, our detailed analysis of the cryptic binding pockets in the Omicron complexes revealed a number of conserved and functionally important sites in the NTD and RBD regions, as well as near the inter-protomer and inter-domain hinges. Unexpectedly, the results unveiled severa1 conserved cryptic pockets located deeply in the S2 subunit at the intersection of UH, CH and HR1 regions. Importantly, our findings pointed to divergences in the distribution of cryptic sites in a more dynamic XBB.1 trimer complex, resulting in the increased number of smaller pockets in less rigid S2 subunits while preserving the composition of cryptic binding sites in the NTD and RBD regions. Overall, the predicted binding pockets captured the experimentally known allosteric sites in the RBD, NTD and S2 apex hinge regions that emerged as top-ranked pockets in the computational analysis.

4. Discussion

In this study, we performed a systematic comparative analysis of conformational dynamics, binding energetics and protein stability in the RBD-ACE2 complexes and full S trimer-ACE2 complexes for BA.2, BA.2.75 and XBB.1 variants. Based on simulations and using conformational ensembles for these systems, we also carried out a comprehensive analysis of the evolution of cryptic binding pockets in the Omicron variants. XBB.1 subvariant is a descendant of BA.2 and recombinant of BA.2.10.1 and BA.2.75 variants. Conformational dynamics and MSM analysis confirmed that a combination of F486S and F490S mutations in XBB.1 may induce increased RBD mobility in regions important for antibody binding and thereby enhance the immune evasion potential of this Omicron variant. The MSM analysis suggested that BA.2 and XBB.1 RBDs may exhibit a similar degree of conformational plasticity in the complex, while BA.2.75 RBD is characterized by pronounced structural rigidity featuring a single dominant microstate.

Among the central results of our study is evidence of progressive stabilization of the RBD regions and inter-protomer interfaces in BA.2 and BA.2.75 variants. At the same time, XBB.1 trimer may be less densely packed and has a smaller inter-protomer interface, which may contribute to the reduced stability and greater plasticity of the XBB.1 variant. It may be argued that through enhanced plasticity, the XBB.1 variant may improve virus adaptability and enhance immune evasion potential. The importance of these findings can be appreciated in the context of experimental data showing that Omicron BA.1, BA.2 and BA.2.75 open trimer structures can be characterized by progressively enhanced inter-domain interactions and improved packing of the inter-protomer interfaces, leading to more stable open states of the S protein for these variants [19]. Our analysis of the XBB.1 trimer complex showed looser packing and weakened interfacial contacts resulting from the increased conformational flexibility. As it appeared, this may alter the composition of functional pockets near inter-protomer hinges and affect the population distribution of closed and open states. Strikingly, cryo-EM studies of XBB.1 ectodomain and the XBB.1 S-ACE2 complex showed two closed states (closed-1 and closed-2) for the unbound XBB.1 trimer [31]. Moreover, the structure of the RBD one-up state, which is stable for BA.2 and BA.2.75 even in the absence of ACE2, was hardly observed in the structures of XBB.1 S protein [31]. Our analysis is consistent with these data, confirming the dynamic nature of the XBB.1 trimer, which may alter cryptic sites and remodel conserved inter-protomer hinge sites seen in BA.2 and BA.2.75 structures. The increased conformational mobility and variability of cryptic pockets in XBB.1 can also be related to HDX-MS data showing S2 adaptability in the Omicron S protein, which is propagated to HR1 and CH regions [35]. As a result, sequestering some of these S2 pockets in XBB.1 through selective targeting may rigidify specific XBB.1 state and disrupt dynamic changes priming S protein for the fusion events.

One of the important issues addressed in our study is how evolutionary changes in the Omicron variants affect the stability and distribution of cryptic sites and what is the functional role of various binding sites in the context of S protein mechanisms and activity. Mutational changes in the Omicron trimers preserved druggable NTD and RBD pockets that are ranked as the most probable across all studied systems. The conformational flexibility of the NTD regions that harbor sites of conformational and mutational divergence still retains several cryptic pockets that are shared among Omicron variants. The consistent emergence of cryptic pockets in the NTD across all Omicron complexes showed that these hidden epitopes may be available for targeting even in the presence of bound ACE2 in open trimers. These results support recent data suggesting that the NTD can serve as an adaptable antigenic surface capable of unlocking cryptic binding pockets, which may enable efficient immune escape [35]. The pocket detection analysis consistently uncovered the allosteric RBD site targeted by the LA molecule [39,40,41].

The results of our study suggested that despite the general rigidity of the S2 regions in comparison with more adaptable S1 subunit undergoing functional movements, there is also an appreciable level of conformational adaptability in the S2, resulting in a significant number of dynamic cryptic pockets in the S2. In particular, our data pointed to several conserved pockets in the HR1 and CH regions of the rigid S2 subunit, thereby indicating that an appreciable level of plasticity is present in S2 regions, giving rise to the broader accumulation of dynamic cryptic pockets. There has been a surge of interest in developing broad-spectrum fusion inhibitors, including antibodies, peptides and small molecules, as well as vaccines targeting the conserved elements in the S2 subunit such as the fusion peptide, stem helix and heptad repeats 1 and 2 (HR1-HR2) bundle. These targetable elements emerged as promising targets for the design of small molecules that can modulate various mechanisms of action, such as fusion mechanism of action and allostery-based modulators regulating the experimentally established cross-talk between fusion peptide regions and RBD regions. The results of our study revealed several conserved and druggable cryptic pockets in the regulatory hinge regions of CTD1-S2 and S2 that are involved in communication between functional elements of S2 and RBD. Targeted ligand screening in the predicted druggable sites can allow for the design of modulators of the S activity and facilitate the development of chemical probes of S functions.

5. Conclusions

In the current study, we performed a systematic comparative analysis of the conformational dynamics, allostery and cryptic binding pockets in the RBD-ACE2 complexes trimers and S trimer complexes with the ACE2 receptor. Multiple microsecond MD simulations and Markov state model (MSM) analysis of the RBD-ACE2 complexes for the Omicron BA.2, BA.2.75 and XBB.1 variants enabled a detailed analysis of conformational states and populations. Using a comparative MSM analysis, we showed an exquisite stability signature of the BA.2.75 RBD, which could be contrasted with the increased mobility of the XBB.1 RBD. Using conformational ensembles of the SARS-CoV-2 Omicro S trimers, we conducted a systematic binding pocket screening and analysis of functional cryptic pockets in the BA.2, BA.2.75 and XBB.1 complexes with ACE2. The results of this study connected insights from conformational dynamics analysis, comparative mutational scanning of the S-ACE2 binding and the inter-protomer interactions with the evolution of cryptic binding sites. We demonstrated that our approach could reproduce all experimentally known allosteric sites and identify networks of conserved cryptic pockets preserved in different conformational states of the Omicron variants and in the S-ACE2 complexes. The results of our study suggest that despite the general rigidity of the S2 regions in comparison with the more dynamic S1 subunit, there is still an appreciable level of conformational adaptability in the S2, resulting in a significant number of dynamic cryptic pockets. The determined cryptic binding pockets at the inter-protomer regions and in the functional regions of the S2 subunit, such as the HR1-HR2 bundle and stem helix region, are consistent with the role of the pocket residues in modulating conformational transitions and antibody recognition. Of particular interest is the detection of highly probable pockets in the S2 hinge regions known to be allosterically linked with the S1 movements during activation. The results detailed how mutational and conformational changes in the BA.2 and BA.2.75 spike trimers can modulate the distributions and mediate networks of inter-connected conserved and variant-specific druggable allosteric pockets. These findings can be important for understanding mechanisms underlying functional roles of druggable cryptic pockets that can be used for both site-specific and allostery-inspired therapeutic intervention strategies targeting distinct conformational states of the SARS-CoV-2 Omicron variants. This may enable the engineering of allosteric modulators that could rationally target a complex functional landscape of virus transmissibility.