Balancing Functional Tradeoffs between Protein Stability and ACE2 Binding in the SARS-CoV-2 Omicron BA.2, BA.2.75 and XBB Lineages: Dynamics-Based Network Models Reveal Epistatic Effects Modulating Compensatory Dynamic and Energetic Changes

Evolutionary and functional studies suggested that the emergence of the Omicron variants can be determined by multiple fitness trade-offs including the immune escape, binding affinity for ACE2, conformational plasticity, protein stability and allosteric modulation. In this study, we systematically characterize conformational dynamics, structural stability and binding affinities of the SARS-CoV-2 Spike Omicron complexes with the host receptor ACE2 for BA.2, BA.2.75, XBB.1 and XBB.1.5 variants. We combined multiscale molecular simulations and dynamic analysis of allosteric interactions together with the ensemble-based mutational scanning of the protein residues and network modeling of epistatic interactions. This multifaceted computational study characterized molecular mechanisms and identified energetic hotspots that can mediate the predicted increased stability and the enhanced binding affinity of the BA.2.75 and XBB.1.5 complexes. The results suggested a mechanism driven by the stability hotspots and a spatially localized group of the Omicron binding affinity centers, while allowing for functionally beneficial neutral Omicron mutations in other binding interface positions. A network-based community model for the analysis of epistatic contributions in the Omicron complexes is proposed revealing the key role of the binding hotspots R498 and Y501 in mediating community-based epistatic couplings with other Omicron sites and allowing for compensatory dynamics and binding energetic changes. The results also showed that mutations in the convergent evolutionary hotspot F486 can modulate not only local interactions but also rewire the global network of local communities in this region allowing the F486P mutation to restore both the stability and binding affinity of the XBB.1.5 variant which may explain the growth advantages over the XBB.1 variant. The results of this study are consistent with a broad range of functional studies rationalizing functional roles of the Omicron mutation sites that form a coordinated network of hotspots enabling a balance of multiple fitness tradeoffs and shaping up a complex functional landscape of virus transmissibility.


Structural Modeling and Refinement
The crystal structures of the BA.2 RBD-hACE2 (pdb id 7XB0), and BA.2.75 RBD-hACE2 complexes (pdb id 8ASY) (Supporting Figure S1) were obtained from the Protein Data Bank [88]. During the 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 [89]. The missing loops in the studied cryo-EM structures of the SARS-CoV-2 S protein were reconstructed and optimized using the template-based loop prediction approach ArchPRED [90]. The side chain rotamers were refined and optimized by the SCWRL4 tool [91]. The protein structures were then optimized using atomic-level energy minimization with composite physics and knowledge-based force fields implemented in the 3Drefine method [92,93]. The refined structural models of the XBB.1 RBD-ACE2 and XBB.1.5 RBD-ACE2 complexes were obtained with the aid of the MutaBind2 approach that utilizes molecular mechanics force fields and fast side-chain optimization algorithms via the random forest (RF) method [94,95]. MutaBind2 utilizes the FoldX approach [96,97] to introduce single or multiple point mutations on the crystal structure followed by robust side-chain optimization and multiple rounds of energy minimization using the NAMD 2.9 program [98] with CHARMM36 force field [99].

Coarse-Grained Brownian Dynamics Simulations
Coarse-grained Brownian dynamics (CG-BD) simulations have been conducted using the ProPHet (Probing Protein Heterogeneity) approach and program [100][101][102][103]. BD simulations are based on a high resolution CG protein representation where 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. 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 as described in the original pioneering studies that introduced Brownian dynamics for simulations of proteins [104,105]. 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 100,000 BD steps at a temperature of 300 K. The CG-BD conformational ensembles were also subjected to all-atom reconstruction using the PULCHRA method [106] and CG2AA tool [107] to produce atomistic models of simulation trajectories.

Molecular Dynamics Simulations
NAMD 2.13-multicore-CUDA package [98] with CHARMM36 force field [99] was employed to perform 500 ns all-atom MD simulations for each of the Omicron RBD-hACE2 complexes. The structures of the SARS-CoV-2 S-RBD complexes were prepared in Visual Molecular Dynamics (VMD 1.9.3) [108] by placing them in a TIP3P water box with 20 Å thickness from the protein. Assuming normal charge states of ionizable groups corresponding to pH = 7, sodium (Na + ) and chloride (Cl − ) counter-ions were added to achieve charge neutrality and a salt concentration of 0.15 M NaCl was maintained. All Na + and Cl − ions were placed at least 8 Å away from any protein atoms and from each other. The long-range non-bonded van der Waals interactions were computed using an atombased cutoff of 12 Å with the switching function beginning at 10 Å and reaching zero at 14 Å. The SHAKE method was used to constrain all bonds associated with hydrogen atoms. Simulations were run using a leap-frog integrator with a 2 fs integration time step. The ShakeH algorithm of NAMD was applied for water molecule constraints. The long-range electrostatic interactions were calculated using the particle mesh Ewald method [109] with a cut-off of 1.0 nm and a fourth order (cubic) interpolation. Simulations were performed under NPT ensemble with Langevin thermostat and Nosé-Hoover Langevin piston at 310 K and 1 atm. The damping coefficient (gamma) of the Langevin thermostat was 1/ps. The Langevin piston Nosé-Hoover method in NAMD is a combination of the Nose-Hoover constant pressure method [110] with piston fluctuation control implemented using Langevin dynamics [111,112]. Energy minimization was conducted using the steepest descent method for 100,000 steps. All atoms of the complex were first restrained at their crystal structure positions with a force constant of 10 Kcal mol −1 Å −2 . Equilibration was done in steps by gradually increasing the system temperature in steps of 20 K starting from 10 K until 310 K and at each step 1 ns equilibration was done keeping a restraint of 10 Kcal mol −1 Å −2 on the protein C α atoms. After the restraints on the protein atoms were removed, the system was equilibrated for an additional 10 ns. An NPT production simulation was run on equilibrated structures for 500 ns keeping the temperature at 310 K and constant pressure (1 atm).

Distance Fluctuations Stability and Communication Analysis
We employed distance fluctuation analysis of the simulation trajectories to compute residue-based rigidity/flexibility profiles. The fluctuations of the mean distance between each pseudo-atom belonging to a given amino acid and the pseudo-atoms belonging to the remaining protein residues were computed. The fluctuations of the mean distance between a given residue and all other residues in the ensemble were converted into distance fluctuation stability indexes that measure the energy cost of the residue deformation during simulations [100][101][102][103]. The distance fluctuation stability index for each residue is calculated by averaging the distances between the residues over the simulation trajectory using the following expression: d ij is the instantaneous distance between residue i and residue j. k B is the Boltzmann constant, T = 300 K, denotes an average taken over the MD simulation trajectory and d i = d ij j * is the average distance from residue i to all other atoms j in the protein (the sum over j * implies the exclusion of the atoms that belong to the residue i). The distances between residue i and residue j are calculated for each conformation along MD trajectories and the mean values of the inter-residue distances are obtained from averaging over the complete ensemble derived from MD simulations.
The interactions between the C α atom of residue i and the C α atom of the neighboring residues i − 1 and i + 1 are excluded in the calculation since the corresponding distances are constant. The inverse of these fluctuations yields an effective force constant k i that describes the ease of moving an atom with respect to the protein structure.

Binding Free Energy Computations: Mutational Scanning and Sensitivity Analysis
The binding free energies were initially computed for the Omicron RBD-hACE2 complexes and were performed for the crystal structures and the refined structural models using a contact-based predictor of binding affinity Prodigy [113][114][115][116]. We conducted a mutational scanning analysis of the binding epitope residues for the SARS-CoV-2 S RBD-ACE2 complexes. Each binding epitope residue was systematically mutated using all substitutions and corresponding protein stability and binding free energy changes were computed. BeAtMuSiC approach [117][118][119] was employed that is based on statistical potentials describing the pairwise inter-residue distances, backbone torsion angles and Viruses 2023, 15, 1143 7 of 36 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 the 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: The change of the binding energy due to a mutation was calculated then as the following: We leveraged rapid calculations based on statistical potentials to compute the ensembleaveraged binding free energy changes using equilibrium samples from simulation trajectories. The binding free energy changes were obtained by averaging the results over 1000 and 10,000 equilibrium samples for each of the studied systems.

Dynamic Network Analysis
A graph-based representation of protein structures [120,121] is used to represent residues as network nodes and the inter-residue edges to describe non-covalent residue interactions. The network edges that define residue connectivity are based on non-covalent interactions between residue side-chains. The residue interaction networks were constructed by incorporating the topology-based residue connectivity MD-generated maps of residues cross-correlations [122] and coevolutionary couplings between residues measured by the mutual information scores [123]. The edge lengths in the network are obtained using the generalized correlation coefficients associated with the dynamic correlation and mutual information shared by each pair of residues. The length (i.e., weight) of the edge that connects nodes i and j are defined as the element of a matrix measuring the generalized correlation coefficient between residue fluctuations in structural and coevolutionary dimensions. Network edges were weighted for residue pairs with in at least one independent simulation. The matrix of communication distances is obtained using generalized correlation between composite variables describing both dynamic positions of residues and coevolutionary mutual information between residues. Residue Interaction Network Generator (RING) program [124] was employed for generation of the residue interaction networks using the conformational ensemble where edges have an associated weight reflecting the frequency in which the interaction present in the conformational ensemble. The residue interaction network files in xml format were obtained for all structures using RING v3.0 webserver [124]. Network graph calculations were performed using the python package NetworkX [125]. Using the constructed protein structure networks, we computed the residue-based betweenness parameter. The short path betweenness of residue i is defined to be the sum of the fraction of shortest paths between all pairs of residues that pass through residue i: where g jk denotes the number of shortest geodesics paths connecting j and k, and g jk (i) is the number of shortest paths between residues j and k passing through the node n i . Residues with high occurrence in the shortest paths connecting all residue pairs have a higher betweenness values. For each node n, the betweenness value is normalized by the number of node pairs excluding n given as (N − 1)(N − 2)/2, where N is the total number of nodes in the connected component that node n belongs to. The normalized short path betweenness of residue i can be expressed as follows: g jk is the number of shortest paths between residues j and k; g jk (i) is the fraction of these shortest paths that pass through residue i.

Network-Based Mutational Profiling of Allosteric Residue Potentials and Epistasis
Through mutation-based perturbations of protein residues we compute dynamic couplings of residues and changes in the average short path length (ASPL) averaged over all possible modifications in a given position. The change of ASPL upon mutational changes of each node is reminiscent of the calculation of residue centralities by systematically removing nodes from the network.
where i is a given site, j is a mutation and · · · denotes averaging over mutations. ∆L node i (j) describes the change of ASPL upon mutation j in a residue node i. ∆L i is the average change of ASPL triggered by mutational changes in position i.
The Z-score is then calculated for each node as follows: ∆L is the change of the ASPL under mutational scanning averaged over all protein residues in the S-RBD and σ is the corresponding standard deviation. The ensembleaveraged Z-scores ASPL changes are computed from a network analysis of the conformational ensembles using 10,000 snapshots of the simulation trajectory. Through this approach, we evaluate the effect of mutations in the RBD residues on long-range allosteric couplings with the other residues in the RBD-ACE2 complex. We used a measurement based on the Jensen-Shannon divergence (JS) for measuring the similarity between the two distributions of mutation-induced ASPL changes in the Omicron variants relative to the original Wu-Hu-1 strain. Given two distributions, p and q, both with g categories, the Kullback-Leibler (KL) divergence is defined as follows: Given two distributions, p and q, both with g categories, the JS divergence is defined as follows: JS(p, q) = 0.5 KL p p+q 2 + 0.5 KL q p+q 2 (10)

Network-Based Community Decomposition Analysis
The analysis of the interaction networks was done using network parameters such as cliques and communities. The Girvan-Newman algorithm [126] is used to identify local communities. In this approach, edge centrality (also termed edge betweenness) is defined as the ratio of all the shortest paths passing through a particular edge to the total number of shortest paths in the network. The method employs an iterative elimination of edges with the highest number of the shortest paths that go through them. By eliminating edges, the network breaks down into smaller communities. The algorithm starts with one vertex, calculates edge weights for paths going through that vertex, and then repeats it for every vertex in the graph and sums the weights for every edge. However, in complex and dynamic protein structure networks it is often that number of edges could have the same highest edge betweenness. An improvement of the Girvan-Newman method was implemented, and the algorithmic details of this modified scheme were given in our recent studies [82]. In this modification of Girvan-Newman method, instead of a single highest edge betweenness removal, all highest betweenness edges are removed at each step of the protocol. This modification makes community structure determination invariant to the labeling of the nodes in the graph and leads to a more stable solution. The modified algorithm proceeds through the following steps: (a) Calculate edge betweenness for every edge in the graph; (b) Remove all edges with the highest edge betweenness within a given threshold; (c) Recalculate edge betweenness for remaining edges; (d) Repeat steps b-d until the graph is empty. By eliminating edges, the network breaks down into smaller communities.

Network Clique-Based Model of Epistatic Interactions
The k-cliques are complete sub graphs of size k in which each node is connected to every other node. In our application, a k-clique is defined as a set of k nodes that are represented by the protein residues in which each node is connected to all the other nodes. A k-clique community is determined by the Clique Percolation Method [127] as a subgraph containing k-cliques that can be reached from each other through a series of adjacent kcliques. We have used a community definition according to which in a k-clique community two k-cliques share k − 1 or k − 2 nodes. Computation of the network parameters was performed using the Clique Percolation Method as implemented in the CFinder program [128]. Given the chosen interaction cutoff I min we typically obtain communities formed as a union of k = 3 and k = 4 cliques. The interaction cliques were considered to be dynamically stable if these interaction networks remained to be intact in more than 75% of the ensemble conformations.
In the dynamic network model, it is assumed that allosteric interactions and longrange communications can be propagated through stable interaction networks in which the key network hubs serve as mediators of allosteric couplings. We assume that residues that belong to the same clique during simulations would have stronger dynamic and energetic couplings leading to synchronization and potentially epistatic effects. To examine the epistatic effect of a mutational site, we compared changes in the k-clique community distributions induced by single and double mutations and calculated the probability by which the two mutational sites belong to the same interfacial 3-clique [129].
We computed the proportion P ab of snapshots in the ensemble in which the two mutational sites (a, b) belong to the interfacial 3-clique: C ab (i) = 1 if (a, b) belong to the same 3-clique. P ab measures the probability that two sites (a, b) are kept in some 3-clique due to either direct or indirect interactions. The closer P ab is to 1, the more likely a and b tend to have a tight connection and potential local epistasis. To further investigate the effect of mutations on the 3-clique probability, we compared changes in P ab after single and double mutations. If double mutations have a greater effect on P ab than single mutations, there may be an epistatic effect between the two sites. To quantify the degree of epistasis, we calculated the ratio of P ab after double mutations to P ab after single mutations. A ratio value of greater than 1 indicates the presence of epistasis between the two sites. If the probability of two sites belonging to the same 3-clique increases after double mutations, it would indicate that there is an epistatic effect between the two sites.

Atomistic MD Simulations Reveal Common and Distinct Signatures of Conformational Dynamics and Interaction Patterns in the ACE2 Complexes with the Omicron RBD Variants
To examine the dynamic signatures of the Omicron variants BA.2, BA.2.75, XBB.1 and XBB.1.5 (Table 1) we conducted multiscale simulations of the RBD-ACE2 complexes (Supporting Figure S1) that included multiple independent CG-BD simulations followed by atomistic reconstruction of the trajectories as well as all-atom MD simulations (Table 2). Through dynamics analysis, we probed the intrinsic conformational dynamics and identified differences in the RBD stability for the BA.2, BA.2.75 and XBB subvariants. Using atomistic simulation analysis, we also examined a hypothesis about whether some Omicron mutations may exert their effect on the stability and binding through long-range allosteric effects and epistatic relationships.   (Table 1). Importantly, some of these RBD mutations are known for their immune evasion functions, including R346T, G446S and F486S [130,131]. XBB.1.5 is essentially identical to XBB.1 with a critical single RBD modification F486P mutation (Supporting Figure S1). R493Q reversed mutation is present in XBB.1 and XBB.1.5 subvariants as well as in the BA.4/BA.5 variants. The binding surface patch of Omicron mutations centered around the key binding hotpots R498, Y501 and H505 and becomes broadened in the BA.2.75 (Supporting Figure S1D) and XBB.1.5 variants (Supporting Figure S1F). Notably, N440K is structurally disconnected from other Omicron sites, and it is involved in direct intermolecular contact with the ACE2. A peripheral N460K mutational site in BA.2.75 and XBB.1.5 is located away from the binding interface.
Conformational dynamics profiles obtained from CG-BD and MD simulations were similar and revealed several important trends. Here, for clarity of presentations, we focused on all-atom MD simulations and analyzed the root mean square fluctuations (RMSF) distributions for the RBD and ACE2 residues ( Figure 1A,B). In the analysis, we used 500 ns productive trajectories and a total of 100,00 stored equilibrium conformational samples. The RBD has two subdomains, where the flexible RBM with a concave surface is involved in direct interaction contacts with hACE2. The second subdomain is a five-stranded antiparallel β-sheet core region that functions as a stable core of the RBD. The conformational mobility distributions for the Omicron RBD complexes displayed several deep local minima corresponding to the RBD core residue cluster (residues 396-403) and the interfacial RBD positions involved in the contacts with the hACE2 receptor (residues 440-456 and 490-505 of the binding interface) ( Figure 1A). The observed structural stability of the RBD core regions was also seen in our earlier simulation studies of the RBD Wu-Hu-1 and Omicron complexes [86,87], further confirming that these segments remain mostly rigid across all examined RBD complexes, with hACE2. Noteworthy, the most stable RBD positions included several important hydrophobic stability centers F400, I402, F490, Y453, L455, A475, and Y489 ( Figure 1A). Some of these hydrophobic RBD positions (Y453, L455, and Y489) are also involved in the favorable interfacial contacts with hACE2 and correspond to a stable conserved region of the RBD-hACE2 interface. The RMSF profiles revealed signs of greater stability for the BA.2.75 RBD as compared to other variants, featuring small thermal fluctuations not only for the ACE2-interacting sites but also exhibiting moderate displacements in the flexible RBD regions (residues 355-375 and 380-400) ( Figure 1A). Despite similar dynamic profiles for all Omicron-hACE2 complexes, we noticed that stable RBD core regions (residues 400-420, 430-450) exhibited even smaller fluctuations in the BA.2.75 and XBB.1.5 complexes ( Figure 1A), suggesting the increased RBD stability for these variants. This may be a relevant contributing factor to the stronger binding affinities seen for BA.2.75 and XBB.1.5 RBDs. The RMSF profile for the XBB.1.5 RBD is characterized by several deep minima corresponding to stabilized regions in the RBD core and particularly the key ACE2 binding interface cluster (residues 485-505) ( Figure 1A). Furthermore, mapping of the Omicron XBB.1.5 mutational sites onto the conformational flexibility profiles highlighted the increased stabilization of P486, S490, R493Q, R498, Y501 and H505 residues that become virtually immobilized in their interfacial positions ( Figure 1A). Although this critical binding hotspot region is stable for all variants mainly due to strong interactions with ACE2, our findings indicated that the corresponding RBD positions become more rigid in BA.2.75 and XBB.1.5 ( Figure 1A). This is generally consistent with the stronger binding affinities of these variants. The increased stabilization of the core RBD regions is accompanied by the moderately increased mobility localized around specific residues including RBM mutational sites N477 and K478. However, the majority of Omicron mutational sites in the studied RBD-ACE2 complexes maintained only a moderate degree of mobility with the exception of more flexible sites L368I, F371, N477 and K478 ( Figure 1A).
It is worth noting that K440 and K460 mutational sites appeared to experience fairly moderate fluctuations, indicating that their intrinsic dynamic propensities may be partly curtailed in the BA.2.75 and XBB.1.5 RBD complexes to ensure greater stability of the RBD. The RBM residues that provide the contact interface with ACE2 also displayed relatively smaller movements and are largely stabilized in the complexes ( Figure 1A).
A generally mobile RBM tip is anchored by the F486 position and may experience appreciable flexibility but manifested in different ways depending on the Omicron variant (Supporting Figure S2). Mutations in F486 are of particular interest as F486V, F486I, F486S, F486P have been seen in other variants and arguably represent a convergent evolutionary hotspot shared by the recent wave of Omicron subvariants to optimize tradeoffs between binding affinity to ACE2 and immune evasion. According to the DMS experiments, among the most common F486 mutations (F486V/I/S/L/A/P), F486P imposes the lowest cost in RBD affinity loss and has the largest increase in RBD expression [58,59]. According to the escape calculator, the F486 position is also one of the major hotspots for escaping neutralization by antibodies [132]. The flexible RBM loops (residues 473-487) appeared to be significantly constrained in the BA.2.75 ensemble but remained more dynamic in XBB. 1.5. ing RBD positions become more rigid in BA.2.75 and XBB.1.5 ( Figure 1A). This is generally consistent with the stronger binding affinities of these variants. The increased stabilization of the core RBD regions is accompanied by the moderately increased mobility localized around specific residues including RBM mutational sites N477 and K478. However, the majority of Omicron mutational sites in the studied RBD-ACE2 complexes maintained only a moderate degree of mobility with the exception of more flexible sites L368I, F371, N477 and K478 ( Figure 1A). It is worth noting that K440 and K460 mutational sites appeared to experience fairly moderate fluctuations, indicating that their intrinsic dynamic propensities may be partly curtailed in the BA.2.75 and XBB.1.5 RBD complexes to ensure greater stability of the RBD. The RBM residues that provide the contact interface with ACE2 also displayed relatively smaller movements and are largely stabilized in the complexes ( Figure 1A).
A generally mobile RBM tip is anchored by the F486 position and may experience appreciable flexibility but manifested in different ways depending on the Omicron variant To further examine the differences in the dynamics of RBD for the Omicron variants, we employed time-structure independent component analysis (t-ICA) as the dimensionality reduction tool for the analysis of MD trajectories (Supporting Figure S2). In this method, the slowest-relaxing degree of freedom is determined by solving a generalized eigenvalue problem a CF = CFK C is the covariance matrix, C is the time-lagged covariance matrix with a certain lag time ∆t, where · · · denotes the time average, K and F are eigenvalue and eigenvector matrices, respectively. The dimensionality reduction analysis of conformational ensembles revealed that the RBM tip in the BA.2.75 (Supporting Figure S2) and XBB.1.5 RBD-ACE2 complexes (Supporting Figure S2) can be described as a "hook-like" ordered conformation and is similar to the crystal structure. Our analysis showed that in the BA.2.75 RBD-ACE2 complex, the distribution of equilibrium conformations is completely dominated by the ordered RBM conformation (~90% occupancy). While the RBM tip remained in the folded conformation in the XBB.1.5 RBD complex, the RBM loops experienced greater mobility as compared to structurally stable BA.2.75 RBD complex. In the BA.2 variant, the population of the ordered RBM tip conformations drops considerably (62% occupancy) and more flexible, partly disordered RBM tip conformations contribute measurably to the equilibrium (Supporting Figure S2). Interestingly, in the XBB.1 complex, the RBD tip becomes even more dynamic and circulates between a variety of partly disordered conformations (Supporting Figure S2). The increased variability of the RBM conformations and the ensemble of XBB.1 RBD-ACE2 conformations are reflected in the markedly increased RMSFs for this system ( Figure 1A).
Interestingly, a distal allosteric loop (residues 358-376) showed an appreciable mobility in the XBB.1.5 RBD as compared to BA.2 and BA.2.75 complexes, which may be potentially linked with the functional requirements for modulation of long-range couplings with the remote ACE2-binding interface residues ( Figure 1A). The RMSF profiles suggested that even though the distribution of rigid and flexible RBD regions is preserved and shared across all RBD complexes, the extent of rigidity and mobility in these regions may be modulated to elicit the increased rigidity of several stable RBD regions and maintaining local mobility of the flexible sites. The conformational dynamics profile of the ACE2 receptor showed a similar and strong stabilization of the interfacial helices on ACE2, indicating that dynamics signatures of the bound hACE2 receptor remain largely conserved across all Omicron RBD complexes ( Figure 1B). Importantly, we observed that the increased structural stability of the BA.2.75 and XBB.1.5 RBDs is accompanied by commensurate rigidification of the ACE2 flexible regions ( Figure 1B). Moreover, the conformational dynamics profiles showed that the induced stability of the ACE2 interfacial regions (residues 350-395) becomes amplified in the BA.2.75 and XBB.1.5 complexes.
Using conformational ensembles, we computed the fluctuations of the mean distance between each residue and all other protein residues that were converted into distance fluctuation stability indexes to measure the energetics of the residue deformations ( Figure 1C). The high values of distance fluctuation indexes correspond to more rigid residues as they display small fluctuations in their distances to all other residues, while small values of this index would point to more flexible sites that experience larger deviations of their inter-residue distances. A comparative analysis of the residue-based distance fluctuation profiles revealed several dominant and common peaks reflecting the similarity of the topological and dynamical features of the RBD-ACE2 complexes. The distributions showed that the local maxima for all RBD-ACE2 complexes are aligned with structurally stable and predominantly hydrophobic regions in the RBD core (residues 400-406, 418-421, 453-456) as well as key binding interface clusters (residues 495-505) that include binding hotspots R498 and Y501 ( Figure 1C). Among RBD positions associated with the high distance fluctuations stability indexes are F400, I402, Y421, Y453, L455, F456, Y473, A475, and Y489 ( Figure 1C). Despite a strikingly similar shape of the distributions for all Omicron RBD variants, which reflected the conserved partition of stable and flexible regions, the larger peaks were seen for BA.2.75 and XBB.1.5 RBD distributions ( Figure 1C). This implies that the RBD core regions and ACE2-binding interface positions become progressively rigidified in the BA.2.75 and XBB.1.5 variants, suggesting the improved RBD stability and further enhancement of the RBD-ACE2 binding interfaces for these variants. Importantly, common stability hotspots Y449, Y473, L455, F456, and Y489 are constrained by the requirements to maintain RBD stability and binding with the ACE2 host receptor, and therefore may be limited in evolving antibody escaping variants.
By highlighting BA.2 and XBB.1.5 mutational sites on the distributions, it may be noticed that the majority of these RBD positions are characterized by low-to-moderate stability indexes, indicating a tendency of Omicron mutations to target conformationally adaptable regions in the RBD. Interestingly, XBB.1.5 mutational positions N440K, V445P, G446S, N460K, F486P and F490S displayed moderate distance fluctuation indexes, which may indicate presence of local mobility in these regions ( Figure 1C). On the contrary, the important RBD binding interface centers R498, Y501 and H505 featured high stability indexes, reflecting a considerable rigidification of these residues due to strong interactions with ACE2. According to our previous studies [78][79][80][81][82][83] residues with the high distance fluctuation indexes often serve not only as structurally rigid centers but also as allosteric regulatory sites that control signal communication. Hence, Omicron sites R498, Y501 and H505 shared by all variants could function as key stability centers, binding hotspots as well as allosteric mediators of long-range communications in the RBD-ACE2 complexes. In this context, it may be instructive to relate our findings to the experimentally observed critical role of these hotspots in providing compensatory epistatic interactions with other Omicron sites [65][66][67] which allow us to rescue sufficiently strong ACE2 binding and offset the effects of destabilizing immune escape mutations. The observed "segregation" pattern of the Omicron sites portioning into a local cluster of structurally stable binding centers and a broadly distributed group of more flexible sites may reflect evolutionary requirements for energetically tolerant sites of immune evasion. The stability hotspots Y449, Y473, L455, F456, and Y489 featuring high indexes are constrained by the requirements for the RBD folding and binding with the ACE2 host receptor, and therefore may be limited in evolving antibody escaping mutations. At the same time, Omicron variant mutations in more dynamic sites may induce marginal destabilizing effects that are compensated for by localized clusters of structurally stable binding affinity hotspots. This interpretation is consistent with the notion that the acquisition of functionally balanced substitutions to optimize multiple fitness tradeoffs between immune evasion, high ACE2 affinity and sufficient conformational adaptability may potentially be a common strategy of SARS-CoV-2 evolution executed by Omicron subvariants.
Structural mapping of the conformational dynamics profiles further highlighted similarities between complexes, while showing some long-range stabilization of the BA.2.75 and BA.2 subvariants ( Figure 1D-F). Overall, the analysis suggested that the largest increase in stability of the RBD-hACE2 complexes may be induced by Omicron BA.2.75 mutations. The dynamic signatures of the Omicron RBD subvariants suggested that BA.2.75 and XBB.1.5 Omicron mutations may induce a cumulative effect that is manifested in small rearrangements at the intermolecular interface and moderate strengthening of the distal ACE2 regions ( Figure 1D-F). Structural maps also highlighted a progressive rigidification of the RBD-ACE2 binding interface regions in these complexes as well as the improved stability of the RBD core regions. We argue that the improved RBD stability in BA.2.75 and XBB.1.5 complexes with ACE2 may be one of the factors linked with the experimentally observed enhancement in their binding affinity with ACE2. Overall, the dynamic analysis of the RBD-ACE2 complexes suggested complementary roles of the Omicron mutation sites that may form a network of allosterically connected stable and more dynamic functional centers to enable modulation of structural stability, binding and long-range signaling.

Computational Mutational Scanning of the RBD Residues Identifies Structural Stability and Binding Affinity Hotspots in the RBD-ACE2 Complexes: Revealing Complementary Functional Effects of Omicron Mutations
Structural analysis of the RBD binding interface epitopes and the Omicron mutational sites ( Figure 2) highlighted considerable similarities in the topography of the binding epitopes and "expansionary" character of the XBB.1 and XBB.1.5 variant positions ( Figure 2E-H). The central binding interface cluster anchored by R498 and Y501 hotspots becomes further consolidated with the addition of T346, P445 and S446 mutations ( Figure 2G,H). This key binding interface region includes sites of Omicron mutations R498, Y501 and H505 that provide the bulk of the ACE2 binding affinity. A visual inspec-tion of these structural maps pointed to a denser Omicron mutational "shield" for XBB.1 ( Figure 2E,F) and XBB.1.5 ( Figure 2G,H) that encircles the core of the RBD epitope and partially overlaps with its peripheral areas. The majority of the Omicron mutations, with the exception of the R498/Y501/H505 cluster, tend to emerge near borders of the binding epitope causing relatively moderate changes in RBD stability and binding while targeting positions to induce a broad antibody escape ( Figure 2). The analysis of the intermolecular contacts in the RBD-ACE2 complexes with the cutoff for the atom contact distance of 5.5 Å and cutoff for salt bridges at 3.5 Å [115,116] yielded an overall similar number of the RBD residues forming interaction contacts with ACE2. The interaction atom pairs for the RBD-ACE2 complexes are listed in Tables S1-S4.   Among instructive observations of this structural analysis was that mutational positions N440K and N460K are not involved in the intermolecular contacts with ACE2. It is evident that structure-based considerations alone are not sufficient to dissect the functional role of these Omicron mutations that in addition to their immune evading potential may have subtle effects on the RBD stability and ACE2 binding via long-range allosteric interactions. Structure-based binding free energy analysis using a contact-based Prodigy predictor of binding affinity [115,116] revealed a similar number of the interaction contacts mediated by charged residues as well as appreciable contributions of polar-nonpolar and  Table 3). The observed differences in the number of intermolecular contacts were fairly small and resulted in the predicted binding free energies that favored BA.2 and XBB.1.5 variants, showing a moderate loss in the binding affinity for XBB.1 (Table 3). Surprisingly, these simple computations correctly predicted a moderate decrease in the binding affinity of XBB.1 as compared to parental BA.2 [54]. Indeed, these experimental studies showed that the ACE2 binding affinities of XBB and XBB.1 exhibited a modest drop relative to that of BA.2 with K D of 2.00 and 2.06 nM respectively compared to 0.95 nM for BA.2 [54]. A small loss in binding can be attributed to the F486S mutation that could remove favorable hydrophobic interactions similar to what was observed for BA.4/5 variants where F486V and R493Q induced respectively a modest impairment and restoration of binding [58]. At the same time, structure-based binding affinity predictions failed to recognize the experimentally observed strongest binding of BA.2.75 immediately followed by XBB.1.5 variant [56]. Nonetheless, structure-based assessments of the intermolecular interface contacts alone may be insufficient to fully capture subtle functional and binding affinity differences between the RBD-ACE2 complexes for different Omicron subvariants. Using the conformational ensembles obtained from MD simulations of the RBD-ACE2 complexes for BA.2, BA.2.75, XBB.1 and XBB.1.5 variants, we performed a systematic mutational scanning of the RBD residues in these complexes. In silico mutational scanning was done using the BeAtMuSiC approach [117][118][119]. We enhanced this approach by averaging the binding free energy changes over the 10,000 equally distributed sample conformations from the equilibrium ensembles. The reported binding free energy ∆∆G changes were evaluated by averaging the results of computations over 1000 samples from MD simulation trajectories. The resulting "deep" mutational scanning heatmaps are reported for the RBD binding interface residues that make stable contacts with ACE2 in the course of simulations.
To establish a relevance and validity of the in silico mutational scanning, we compared the results of the DMS experiments for the BA.1 and BA.2 RBD residues [58,59] with the computed mutational changes in the protein stability and binding for the BA.1 RBD-hACE2 ( Figure 3A) and BA.2 RBD-hACE2 complexes ( Figure 3B). A statistically significant correlation between the DMS experiments and mutational scanning data was observed, also highlighting the expected dispersion of the distributions. It is worth noting that the computed free energy changes reflected mutation-induced effects on both residue stability and binding interactions. Our findings are consistent with other simulation-based studies that showed correspondence between mutation-induced computed changes in the RBD stability and the experimental protein expression profiles [133]. It could be noticed that the computational predictions of destabilizing changes were often larger than the experimentally observed values. Nonetheless, the scatter plots showed a fairly appreciable correspondence between the predicted and experimental free energy differences, particularly for large destabilizing changes with ∆∆G > 2.0 kcal/mol ( Figure 3A,B). This ensures an adequate identification of the major stabilization and binding affinity hotspots where mutations cause pronounced energetic changes. We also compared the DMS free energies for BA.2 RBD with the computationally predicted changes in XBB.1 and XBB.1.5 RBD complexes with ACE2 ( Figure 3C,D). To provide a systematic comparison, we constructed mutational heatmaps for the RBD binding interface residues (Figure 4). Consistent with the DMS experiments [58,59], the strongest binding energy hotspots in BA.2 and BA.2.75 complexes corresponded to hydrophobic residues Y453, F456, Y473, Y489 and Y501 that play a decisive role in binding for all Omicron complexes ( Figure 4A,B). Mutational heatmaps illustrated that the majority of substitutions in these key interfacial positions can lead to a considerable loss in the stability and binding affinity with ACE2. This analysis is also consistent with our previous studies, suggesting that these conserved hydrophobic RBD residues may be universally important for binding across all Omicron variants and act as stabilizing sites of the RBD stability and binding affinity [86,87]. The common energetic hotspots Y453, F456, Y489 and Y501 also emerged as critical stability and binding hotspots in the experimental DMS studies [58,59]. In addition, mutational scanning of the RBD residues F486, N487 and H505 showed appreciable and consistent destabilization changes, indicating that these residues are important energetic centers of stability and binding (Figure 4). The computed heatmaps are consistent with the DMS experiments in which G446S and F486V mutations Figure 3. The scatter plots of the DMS-derived binding free energy changes for the RBD residues and computational mutational scanning of the RBD residues to estimate mutational effects on ACE2 binding. The effect on ACE2 receptor-binding affinity (∆log10 K D ) of every single amino-acid mutation in SARS-CoV-2 RBD was experimentally determined by high-throughput titration assays using DMS platform [58,59]. The results of computational mutational scanning of the RBD residues were averaged over conformational ensembles obtained from all-atom MD simulations. The scatter plot of the experimental and computed binding free energy changes from mutational scanning of the RBD residues in the Omicron RBD BA.1 complex, pdb id 7WBP (A) and RBD BA.2 complex, pdb id 7XB0 (B). The scatter plot of the experimental free energy changes from mutational scanning of the RBD residues in the RBD BA.2-ACE2 complex and computed binding free energy changes in XBB.1 RBD-ACE2 (C). The scatter plot of the experimental free energy changes from mutational scanning of the RBD residues in the RBD BA.2-ACE2 complex and computed binding free energy changes in XBB.1.5 RBD-ACE2 (D). The data points are shown in light-brown colored circles. The correlation coefficients between of the DMS-derived binding free energy changes for the RBD residues and computational mutational scanning of the RBD residues are shown on panels (A-D) respectively.
To analyze the contributions of the RBD residues to protein stability, we also utilized a simplified SWOTein predictor which identifies the residue contributions to the global folding free energy through three types of database-derived statistical potentials that include inter-residue distances, backbone torsion angles and solvent accessibility [134,135]. Despite its simplicity, this approach considers key contributions to the folding free energy associated with the enthalpic components, hydrophobic interactions and entropic estimations. The stability strengths and weaknesses are identified as residues that upon mutation result in strong destabilization or strong stabilization (Supporting Figure S3). Consistent with the dynamics analysis, among commonly shared stability centers are residues F400, I402, Y421, Y453, L455, F456, Y473, A475, and Y489. Although the residue stability profiles are similar for all variants, we observed favorable stability values for a number of XBB.1.5 mutational sites including V445P, G446S, N460K, S477N, T478K, E484A, F486P, F490S, R493Q, Q498R, N501Y, Y505H (Supporting Figure S3D). A strong stability peak associated with the F486P position in XBB.1.5 variant indicated that convergent mutations in this position can modulate RBD stability and affect ACE2 binding.
To provide a systematic comparison, we constructed mutational heatmaps for the RBD binding interface residues (Figure 4). Consistent with the DMS experiments [58,59], the strongest binding energy hotspots in BA.2 and BA.2.75 complexes corresponded to hydrophobic residues Y453, F456, Y473, Y489 and Y501 that play a decisive role in binding for all Omicron complexes ( Figure 4A,B). Mutational heatmaps illustrated that the majority of substitutions in these key interfacial positions can lead to a considerable loss in the stability and binding affinity with ACE2. This analysis is also consistent with our previous studies, suggesting that these conserved hydrophobic RBD residues may be universally important for binding across all Omicron variants and act as stabilizing sites of the RBD stability and binding affinity [86,87]. The common energetic hotspots Y453, F456, Y489 and Y501 also emerged as critical stability and binding hotspots in the experimental DMS studies [58,59]. In addition, mutational scanning of the RBD residues F486, N487 and H505 showed appreciable and consistent destabilization changes, indicating that these residues are important energetic centers of stability and binding ( Figure 4). The computed heatmaps are consistent with the DMS experiments in which G446S and F486V mutations decrease the ACE2 affinity of BA.2 (by −0.1 & −0.5 log10 K D , respectively), while R493Q buffers these mutations by slightly increasing ACE2 affinity [58,59]. In the context of comparative analysis, it is particularly relevant to notice that the F486 position, which is highly favorable for the RBD stability and binding is mutated to S486 in XBB.1 ( Figure 3C) and P486 in XBB.1.5 (Figure 4). The predicted binding free energy changes showed that the P486 position of XBB.1.5 ( Figure 4D) is less tolerant to modifications and more energetically favorable for the RBD stability and binding as compared to S486 in XBB.1 ( Figure 4C). Importantly, mutational scanning revealed an appreciably destabilizing free energy change ∆∆G = 0.78 kcal/mol for P486S modification in the XBB.1.5 structure, while the reversed S486P in XBB.1 yielded a modest favorable change with ∆∆G = −0.2 kcal/mol ( Figure 4C,D). These results agree with the experiments, indicating that F486P mutation in XBB.1.5 may rescue the loss of binding affinity in XBB.1. Our results provided supporting evidence to the notion that the key functional difference between XBB.1.5 and its immediate parent XBB.1 is that XBB.1.5 has traded the costly F486S mutation for a more favorable F486P mutation thus enhancing both RBD stability and ACE2 binding. We argue that a synergistic effect of the restored binding and the improved RBD stability may favor transmissibility and the observed surge of the XBB.1.5 variant. Indeed, according to functional experiments, F486P imposes the lowest cost in RBD affinity loss and has the largest increase in RBD expression [58,59].
The reversion of Q493R occurred early in BA.2 evolution and while R493Q is not a major antigenic mutation it can arguably enable both F486V (in BA.4/5) and G446S (in BA.2.75) [58,59]. According to mutational scanning results, R493Q mutation is favorable for ACE2 binding with ∆∆G = −0.41 kcal/mol in XBB.1 and with ∆∆G = −0.24 kcal/mol in XBB.1.5 ( Figure 4C,D). This is consistent with the experimental evidence that F486S/P mutations may have been selected to promote immune escape and are buffered by the favorable binding induced by the reversed R493Q change [58,59]. A similar effect was observed for F490S that is marginally unfavorable for binding while the reversed S490F modification would improve binding with ∆∆G = −0.37 kcal/mol ( Figure 4). Our data support a mechanism suggesting that the XBB lineage may have evolved to evade immune suppression and outcompete other Omicron subvariants through mutations of F486 which is a hotspot for establishing protective immunity against the virus. According to this mechanism, to provide a better tradeoff between virus fitness requirements XBB.1.5 acquired F486P which partly restored loss in ACE2 binding by acting cooperatively with R493Q, R498 and Y501 mutations [136].
ing evidence to the notion that the key functional difference between XBB.1.5 and its immediate parent XBB.1 is that XBB.1.5 has traded the costly F486S mutation for a more favorable F486P mutation thus enhancing both RBD stability and ACE2 binding. We argue that a synergistic effect of the restored binding and the improved RBD stability may favor transmissibility and the observed surge of the XBB.1.5 variant. Indeed, according to functional experiments, F486P imposes the lowest cost in RBD affinity loss and has the largest increase in RBD expression [58,59]. The reversion of Q493R occurred early in BA.2 evolution and while R493Q is not a major antigenic mutation it can arguably enable both F486V (in BA.4/5) and G446S (in BA.2.75) [58,59]. According to mutational scanning results, R493Q mutation is favorable for ACE2 binding with ΔΔG = −0.41 kcal/mol in XBB.1 and with ΔΔG = −0.24 kcal/mol in  Figure S4) displayed clear segregation between commonly shared binding hotspot sites (R493/Q493, R498, Y501, H505) and remaining mutational positions showing a notable tolerance to modifications and allowing for evolution in some of these sites without sacrificing the RBD stability or ACE2 binding. Notably, for XBB.1.5 RBD, we observed an additional small group of Omicron sites (I368, N417, P445 and P486) that are more sensitive to modifications, but other Omicron positions are tolerant to substitutions producing a small effect on the RBD stability and binding (Supporting Figure S4D). Overall, analysis of the mutational heatmaps for the RBD residues showed the presence of several structural stability centers (Y453, F456, Y489) in the RBD and binding affinity hotspots (Q493, R498, T500, N501Y), while the remaining sites can potentially tolerate functionally beneficial destabilizing mutations.
A more detailed residue-based binding free energy profile focused on the Omicron mutational changes in the RBD-ACE2 complexes for BA.2, BA.2.75, XBB.1 and XBB.1.5 variants ( Figure 5). The profiles revealed a similar trend across all studied variants. In the BA.2 variant, Omicron mutations G339D, S371F, S373P, S375F, T376A, D405N, R408S, N440K, S477N, T478K, E484A are essentially stability-neutral, having only a marginal effect on binding ( Figure 5A). K417N may induce a moderate change in ACE2 binding while promoting the increased neutralization escape potential of the Omicron variant from antibodies. At the same time Q498R, N501Y and Y505H seem to incur more significant stabilizing changes. Notably, and consistent with the experiments, N501Y modification showed a markedly larger stabilizing change in the ACE2 binding. The dominant stabilizing effect of N501Y is even more significant for BA.2.75 variant with R493Q, Q498R and Y505H also contributing appreciably to the enhanced binding affinity ( Figure 5B). While the trend of most mutations being largely neutral for binding persists for all examined Omicron variants, we observed that for RBD XBB.1.5 the appreciably favorable binding free energy changes are associated with F486P, R493Q reversed mutation, Q498R, N501Y and Y505H mutations ( Figure 5D). In addition, the other XBB.1.5 mutations are largely neutral, highlighting the tolerance of the respective positions to modifications. Consistent with the experimental data [65,67], these profiles further exemplified that functionally beneficial effects of destabilizing mutations (such as immune escaping potential and conformational adaptability) may be contingent on compensatory effects provided by a binding hotspot cluster centered on the R498/Y501 positions. These findings highlighted the important feature of a mechanism that may be characterized by the stability hotspots that ensure sufficient RBD stability and a spatially localized group of key binding affinity centers, while allowing for functionally beneficial but binding neutral (or moderately destabilizing) mutations in other positions to balance tradeoffs between immune evasion and ACE2 binding. The revealed patterns are reminiscent of direct evolution studies showing that enhanced protein stability in key sites can promote broader evolvability and expand a range of beneficial mutations while retaining the stability of the protein fold [137,138]. The related studies further elaborated that epistatic interactions between protein sites are mediated by stability, and that stabilizing mutations are often pre-requisites for adaptive destabilizing substitutions [139]. The mutational heatmaps for Omicron mutational sites are generally consistent with this mechanism, revealing a small group of shared stabilizing RBD positions that protect the stability and binding affinity allowing for substantial evolvability at the remaining tolerant sites.
We argue that by expanding the range of stability/binding affinity hotspots and recruiting F486P/R493Q positions to maintain stability and restore binding affinity, XBB.1.5 variant can optimize both ACE2 binding and immune escaping profiles to achieve superior viral fitness. These findings may be relevant in the context of epistasis in which nonlinear couplings between mutations may determine balance and tradeoffs between stability, evolvability and functions. We suggest that a balance between protein stability requirements and ACE2 binding affinity can promote the evolvability of XBB variants by These findings highlighted the important feature of a mechanism that may be characterized by the stability hotspots that ensure sufficient RBD stability and a spatially localized group of key binding affinity centers, while allowing for functionally beneficial but binding neutral (or moderately destabilizing) mutations in other positions to balance tradeoffs between immune evasion and ACE2 binding. The revealed patterns are reminiscent of direct evolution studies showing that enhanced protein stability in key sites can promote broader evolvability and expand a range of beneficial mutations while retaining the stability of the protein fold [137,138]. The related studies further elaborated that epistatic interactions between protein sites are mediated by stability, and that stabilizing mutations are often pre-requisites for adaptive destabilizing substitutions [139]. The mutational heatmaps for Omicron mutational sites are generally consistent with this mechanism, revealing a small group of shared stabilizing RBD positions that protect the stability and binding affinity allowing for substantial evolvability at the remaining tolerant sites.
We argue that by expanding the range of stability/binding affinity hotspots and recruiting F486P/R493Q positions to maintain stability and restore binding affinity, XBB.1.5 variant can optimize both ACE2 binding and immune escaping profiles to achieve superior viral fitness. These findings may be relevant in the context of epistasis in which nonlinear couplings between mutations may determine balance and tradeoffs between stability, evolvability and functions. We suggest that a balance between protein stability requirements and ACE2 binding affinity can promote the evolvability of XBB variants by tolerating mutations in positions that could confer beneficial phenotypes. Based on this assertion, we also propose that Omicron mutational effects are mediated by protein stability and that the individually destabilizing RBD mutations may be counterbalanced via allostery and epistasis by stabilizing and affinity-enhanced mutations.

Dynamic-Based Network Modeling and Community Decomposition Analysis of the Omicron RBD-ACE2 Complexes Detail Allosteric Role of the Binding Energy Hotspots
Functional studies [65][66][67] suggested that the evolutionary potential of the Omicron variants could be enhanced through epistatic interactions between variant mutations and the broader mutational repertoire available to the S proteins. In particular, it was suggested that weak epistasis in the Wu-Hu-1 original strain may become much stronger in newly emerged Omicron variants as a potential virus mechanism to counteract structural limitations and stability constraints for continued evolution [67]. Here, we employed the ensemble-based modeling of the residue interaction networks utilizing a graph-based description [120,121] in which both dynamics [122] and coevolutionary couplings between protein residues [123] determine the strength of the interaction links. As dynamic couplings between RBD interface residues can be determined in simulations, we propose that strongly coupled residue positions may communicate and affect their ACE2 binding interactions via epistatic relationships.
First, we explored the network modeling of the RBD-ACE2 complexes to explore potential epistatic relationships between RBD residues by quantifying mutational effects on allosteric communications and identifying critical dynamically coupled mediating centers. In this model, by probing dynamic coupling between different positions and allosteric interactions, we evaluate how the network can be altered and rewired through coupled substitutions, which may manifest their epistatic relationships. The mutation-induced changes in the short path length of connecting every pair of residues in the network were calculated from MD trajectories. By introducing single mutations in each residue and double mutations we evaluated the effect of mutations on dynamic allosteric couplings with the other residues in the RBD-ACE2 system. This topological network parameter is used as a proxy to probe potential epistatic couplings between RBD sites and their allosteric potential. By identifying RBD residues where single and double mutations induce a synergistic effect and cause a significant increase in the network's short path length, we locate dynamically coupled mediating hotspots that may be involved in epistasis driving spike stability and binding ( Figure 6).
The distribution of mutation-induced ASPL network changes in the RBD-ACE2 complexes revealed several distinct residue clusters that are characterized by significant peaks (Figure 6). This implies that mutations (or deletions) in these positions can on average lead to significantly increased network path length and therefore have a considerable impact on the fidelity of long-range signaling in the RBD-ACE2 complexes. The relative importance of the major peaks associated with N501Y. Q498R and Y505H positions are amplified as compared to the background original strain, suggesting a synergistic effect induced by these Omicron mutations. Moreover, these sites may be also critical for the efficiency of signal transmission in the complexes ( Figure 6). Interestingly, we observed an appreciable increase in the density of epistatic sites for the RBD residues 501-506, 491-496, 449-453 which is particularly strong in the BA.2.75 RBD-ACE2 complex displaying also greater RBD stability and enhanced AC2 binding ( Figure 6B). The results are with the experimental evidence, showing that epistatic couplings are primarily mediated by affinity-enhancing mutations Q498R and N501Y [59,65,67]. The emergence of clusters of peaks corresponding to structurally proximal mediating sites exhibiting epistatic relationships supported the notion that the epistatic effect is stronger for mediating sites that are close to each other. We performed community decomposition and focused on intrinsically stable RBD modules in which residues are densely interconnected through coupled interactions and dynamic correlations (Figure 7). Here, we will refer to RBD communities as small local interacting modules of inter-connected and dynamically coupled residues within the RBD domain. The distribution of local communities obtained from decomposition is optimal given the constructed residue interaction network. To characterize the "minimalist" community partitioning in the RBD, we varied the criteria for the residue interconnectivity during the network construction and applied the community decomposition protocol to the resulting networks. The reported community partition corresponds to the converged solution with the modules consisting of at least three RBD residues that are dynamically and coevolutionary coupled and can switch their conformational states cooperatively. Assuming that the distribution and number of the intrinsically stable modules provide an estimate of the RBD stability, structural mapping revealed an appreciable and progressive increase in the RBD communities from the BA.2-ACE2 complex ( Figure 7A) to BA.2.75 ( Figure 7B) and XBB.1.5 ( Figure 7C). In particular, the analysis showed a broadly distributed and dense network of inter-connected stable communities in BA.2.75 and XBB.1.5 complexes. We examined the involvement of the Omicron mutational sites in the RBD communities which could shed some light on their role in mediating the RBD stability. Although the key binding role of R498, Y501 and H505 sites is well-established, their role in mediating intrinsic RBD stability is less obvious. In fact, in the BA.2 complex none of the Omicron sites participates in the RBD communities, and structural mapping illustrated the peripheral location of these sites relative to the community distribution ( Figure  7A). We performed community decomposition and focused on intrinsically stable RBD modules in which residues are densely interconnected through coupled interactions and dynamic correlations (Figure 7). Here, we will refer to RBD communities as small local interacting modules of inter-connected and dynamically coupled residues within the RBD domain. The distribution of local communities obtained from decomposition is optimal given the constructed residue interaction network. To characterize the "minimalist" community partitioning in the RBD, we varied the criteria for the residue interconnectivity during the network construction and applied the community decomposition protocol to the resulting networks. The reported community partition corresponds to the converged solution with the modules consisting of at least three RBD residues that are dynamically and coevolutionary coupled and can switch their conformational states cooperatively. Assuming that the distribution and number of the intrinsically stable modules provide an estimate of the RBD stability, structural mapping revealed an appreciable and progressive increase in the RBD communities from the BA.2-ACE2 complex ( Figure 7A) to BA.2.75 ( Figure 7B) and XBB.1.5 ( Figure 7C). In particular, the analysis showed a broadly distributed and dense network of inter-connected stable communities in BA.2.75 and XBB.1.5 complexes. We examined the involvement of the Omicron mutational sites in the RBD communities which could shed some light on their role in mediating the RBD stability. Although the key binding role of R498, Y501 and H505 sites is well-established, their role in mediating intrinsic RBD stability is less obvious. In fact, in the BA.2 complex none of the Omicron sites participates in the RBD communities, and structural mapping illustrated the peripheral location of these sites relative to the community distribution ( Figure 7A). At the same time, in the XBB.1.5 complex, N417 and I368 mutational sites are involved in stable RBD clusters ( Figure 7C). N417 is linked in one community with Y453-V350-W353-I402-I418-N422-Y423-Y495-D398-F400 which includes some of the most stable RBD positions F400, I402, I418 and Y453. This community becomes fragmented in a less stable XBB.1 variant. In addition, while L368 is involved in local community F338-F342-L368 in BA.2, the mutational site I368 helps to consolidate stable clusters in the XBB.1.5 RBD-ACE2 complex that include a group of hydrophobic RBD core residues I358, V524, F338, F342, Y365, I368, V395, V513, F515, I434, F392, F374 and F377 residues ( Figure 7C).
The network modularity analysis and the obtained optimal decomposition of the RBD residue interaction networks highlighted the differences between the intrinsic RBD communities in each of the lineages. The results emphasized that Omicron mutations play a relatively peripheral role in mediating the network modular organization and the RBD stability. These findings reinforced the notion that Omicron positions display significant plasticity and tolerance to substitutions without causing appreciable deleterious effects on the RBD stability. Another important implication of the network community analysis is a clear trend showing the increased density and spatial "expansion" of the stable RBD communities in the BA.2.75 and especially XBB.1.5 complexes (Figure 7). These results support the notion that the improved viral fitness of the XBB.1.5 variant may arise from the enhanced RBD stability and stronger ACE2 binding while maintaining the favorable immune escaping profile by utilizing the functional benefits of conformationally adaptive Omicron sites. These results suggested that differences between Omicron strains can be manifested in variations of the intrinsic RBD stability in the complexes and attributed to changes in the distributions of the RBD communities.
The network modularity analysis and the obtained optimal decomposition of the RBD residue interaction networks highlighted the differences between the intrinsic RBD communities in each of the lineages. The results emphasized that Omicron mutations play a relatively peripheral role in mediating the network modular organization and the RBD stability. These findings reinforced the notion that Omicron positions display significant plasticity and tolerance to substitutions without causing appreciable deleterious effects on the RBD stability. Another important implication of the network community analysis is a clear trend showing the increased density and spatial "expansion" of the stable RBD communities in the BA.2.75 and especially XBB.1.5 complexes (Figure 7). These results support the notion that the improved viral fitness of the XBB.1.5 variant may arise from the enhanced RBD stability and stronger ACE2 binding while maintaining the favorable immune escaping profile by utilizing the functional benefits of conformationally adaptive Omicron sites. These results suggested that differences between Omicron strains can be manifested in variations of the intrinsic RBD stability in the complexes and attributed to changes in the distributions of the RBD communities.
Having established that community analysis can provide a robust network-based metric for assessment of the intrinsic RBD stabilities, we then investigated the effect of Omicron variants on the RBD-ACE2 interfacial communities that can reflect the binding strength in the RBD-ACE2 complexes (Figure 8). We examine the differences between the interfacial RBD-ACE2 communities in the complexes. This analysis provided a stronger connection between the differences in the RBD communities and the changes in the ACE2 complexes as they are observed in the different lineages. We have used this analysis to explore why these distinctions are arising and how they are contributing to the overall differences between the lineages. Omicron variants on the RBD-ACE2 interfacial communities that can reflect the binding strength in the RBD-ACE2 complexes (Figure 8). We examine the differences between the interfacial RBD-ACE2 communities in the complexes. This analysis provided a stronger connection between the differences in the RBD communities and the changes in the ACE2 complexes as they are observed in the different lineages. We have used this analysis to explore why these distinctions are arising and how they are contributing to the overall differences between the lineages. Indeed, due to the structural similarity of the binding interfaces, all RBD-ACE2 complexes share a number of major common communities that persist throughout simulations. These interfacial communities are present in the BA.2 RBD-ACE2 complex and include D38-Y49-R498, Q24-Y83-N487, Y489-F456-K31, N330-D355-T500, D355-T500-Y41, R493-H34-Y453, Y41-K353-Y501, K353-Y501-H505, Y41-Y501-R498 and Y489-K31-R493 ( Figure 8A). In the BA.2.75 RBD-ACE2 complex we detected several additional stable modules (Q24-G476-N487, Q24-Y83-F486, S19-Q24-N477) ( Figure 8B). A further expansion of the key interfacial communities was seen in the XBB.1.5 complex, including Y489-K31-F456, Q325-N439-Q506, N330-D355-T500-Y41-L45 and Y41-K352-Y501-H505-R498 modules ( Figure 8C).
The analysis showed the increased density of major interfacial communities in the BA.2.75 and XBB.1.5 variants as compared to the BA.2 variant which is consistent with the experimentally observed stronger binding of the BA.2.75 and XBB.1.5 variants. Importantly, the commonly shared RBD-ACE2 interfacial communities are primarily mediated by R498, Y501 and Y489 hotspots making them indispensable hubs of the interfacial network and controlling long-range interactions across the binding interface. These networks of local communities the R498 and Y501 binding hotspots may enhance binding strength due to epistatic couplings with other interfacial sites.

A Clique-Based Network Model of Epistatic Interactions Reveals Key Mediators of Binding and Cooperativity in the RBD-ACE2 Complexes
To characterize and rationalize the experimentally observed epistatic effects of the Omicron mutations [65][66][67], we further explored the network analysis and proposed a simple clique-based network model for describing the non-additive effects of the RBD residues. Using equilibrium ensembles and dynamic network modeling of the original RBD-ACE2 complexes for BA.2, BA.2.75, XBB.1 and XBB.1.5 variants, we used mutational scanning to perturb modular network organization represented by a chain of inter-connected sable 3-cliques. Specifically, we calculated the probability by which the two mutational sites belong to the same interfacial 3-clique. For this, we generated an ensemble of 10,000 protein conformations from MD simulations of the studied RBD-ACE2 complexes. By systematically using double mutational changes of the Omicron positions over the course of the MD simulation trajectory for the original RBD-ACE2 complexes we attribute mutational sites that belong to the same 3-clique to have local non-additive effects, while the effects of specific mutations on changes of the entire distribution and the total number of 3-cliques at the RBD-ACE2 interface will be attributed to long-range epistatic relationships.
In network terms, a k-clique is defined a complete sub-graph of k nodes in which each pair of nodes is connected by an edge, reflecting strong mutual interactions and dynamic coupling between every node in the clique with all other nodes that belong to the same clique. A collection of all interconnected k-cliques in a given network defines a k-clique community. If the mutational sites are arranged in a 3-clique structure, all three sites are connected to each other. As a result, when one site is mutated, it will have a greater effect on the stability of the complex because the other two sites will also be affected. This is in contrast to a situation where the sites are not arranged in a 3-clique structure, in which case a mutation in one site would only have a limited effect on the other two sites. Therefore, the presence of a 3-clique structure can be used as a predictor of potential local non-additive effects. At the same time, if Omicron mutations induce global changes in the distribution of 3-cliques along the RBD-ACE2 binding interface, these Omicron mutations may have long-range effect on other interfacial positions that can be propagated via allosterically coupled network of the interfacial 3-cliques. In other words, mutational positions that provide an indispensable stabilizing anchor to multiple intermolecular RBD-ACE2 cliques are assumed to have a stronger non-additive effect on binding.
Using topological network analysis, we reported the overall distributions and composition of stable 3-cliques in the interfacial network that are induced by the key Omicron mutations by Q498R, N501 and Y505H known to be involved in strong compensatory epistasis ( Figure 9) [67]. Our results suggested that non-additive effects may depend on both structural topology and dynamic couplings between mutation sites, which can be quantitatively determined using the distribution of the inter-connected 3-cliques. In the Omicron RBD BA.2 complex, we found a significant accumulation of stable interfacial 3-cliques most of which were directly anchored by Y501 and H505 mutations (K353-Y501-H505, D38-G396-Y501, Y41-K353-Y501, K353-Y501-H505, T500-D355-Y401 and Y41-Y501-R498) ( Figure 9A). We also found that the interfacial 3-cliques anchored by Y501 are the most stable and persist throughout the course of simulations. In network terms, the involvement of N501Y in multiple 3-cliques implies that this mutational site not only enhances ACE2 binding but is strongly dynamically coupled with several RBD positions (Y449, G496, R498 and H505) and allows for the strengthening of their binding interactions thus amplifying the effect of Omicron mutations. The network distribution revealed that R498 and Y501 sites can promote a larger number of stable 3cliques at the central interfacial patch including D38-R498-Y501, R493-H34-Y453, R493-K31-Y489, H34-K31-R493, and R493-K31-F456 cliques ( Figure 9A). Although most of these 3-cliques can be formed independently of N501Y and Q498R mutations, it appeared that an additional D38-R498-Y501 clique can be formed while the remaining cliques along the central patch become more persistent in simulations. These results highlighted the role of the R493 position in anchoring multiple interaction clusters with the ACE2 residues, also indicating some level of dynamic coupling with Y453 and Y489 residues. The conserved Y453 and Y489 positions are involved in favorable hydrophobic interactions and dynamic coupling with the interactions mediated by R493 suggests some level of cooperativity and synchronicity between these contributions. At the same time, our analysis showed that the interfacial cliques in the distal flexible region of the interface (Q24-Y83-N487, Q24-Y83-F486, f486-Y83-N487, Y489-K31-F456, T27-Y473-Y489) are persistent independently of the presence of R498 and Y501 mutations ( Figure 9A). Notably, we observed that F486 and N487 act as local mediators of stable cliques in this region.
Overall, while the topography of the RBD-ACE2 binding interface plays a fundamental role in determining the distribution and composition of the intermolecular network cliques, the dynamic residue couplings can modulate the strength and number of stable cliques through epistatic relationships. These observations agree with the functional studies showing that epistatic shifts in the RBD are primarily driven by the Y501 site [65]. According to these experiments, the largest epistatic shift in mutational effects is In network terms, the involvement of N501Y in multiple 3-cliques implies that this mutational site not only enhances ACE2 binding but is strongly dynamically coupled with several RBD positions (Y449, G496, R498 and H505) and allows for the strengthening of their binding interactions thus amplifying the effect of Omicron mutations. The network distribution revealed that R498 and Y501 sites can promote a larger number of stable 3-cliques at the central interfacial patch including D38-R498-Y501, R493-H34-Y453, R493-K31-Y489, H34-K31-R493, and R493-K31-F456 cliques ( Figure 9A). Although most of these 3-cliques can be formed independently of N501Y and Q498R mutations, it appeared that an additional D38-R498-Y501 clique can be formed while the remaining cliques along the central patch become more persistent in simulations. These results highlighted the role of the R493 position in anchoring multiple interaction clusters with the ACE2 residues, also indicating some level of dynamic coupling with Y453 and Y489 residues. The conserved Y453 and Y489 positions are involved in favorable hydrophobic interactions and dynamic coupling with the interactions mediated by R493 suggests some level of cooperativity and synchronicity between these contributions. At the same time, our analysis showed that the interfacial cliques in the distal flexible region of the interface (Q24-Y83-N487, Q24-Y83-F486, f486-Y83-N487, Y489-K31-F456, T27-Y473-Y489) are persistent independently of the presence of R498 and Y501 mutations ( Figure 9A). Notably, we observed that F486 and N487 act as local mediators of stable cliques in this region.
Overall, while the topography of the RBD-ACE2 binding interface plays a fundamental role in determining the distribution and composition of the intermolecular network cliques, the dynamic residue couplings can modulate the strength and number of stable cliques through epistatic relationships. These observations agree with the functional studies showing that epistatic shifts in the RBD are primarily driven by the Y501 site [65]. According to these experiments, the largest epistatic shift in mutational effects is associated with non-additive contribution of Q498R and N501Y, followed by less pronounced epistatic shifts at sites 491-496, 505-506, and 446-449 [65].
A considerably larger number of stable interfacial cliques are mediated by BA.2.75 Omicron mutations ( Figure 9B), showing the increased number of cliques in the key binding region that are anchored by Y449, G496, R498, and Y501 mutations. Interestingly, we found that the BA.2.75 variant may induce the formation of additional stable cliques (Q42-R498-S446, L45-R498-V445, S446-Q42-Y449, V445-L45-T500) in which Y449 and R498 sites engage specific Omicron mutational positions V445 and S446 ( Figure 9B). The denser network of inter-connected cliques in this region suggested that R498 and Y501 mutations may promote further strengthening of the binding interactions in the BA.2.75 variant amplifying the individually moderate effect of V445 and S446 mutations. We also noticed the increased number of 3-cliques in a more dynamic region where several additional cliques are mediated by F486, N487 and Y489 positions thus suggesting a partial stabilization of the interfacial interactions in this region ( Figure 9B). However, this consolidation of the binding interface in BA.2.75 occurred independently of the R498 and Y501 mutational sites and is mainly determined by reduced mobility in F486 and N487 positions. A similarly significant number of the interfacial cliques are mediated by Y449, G496, R498 and Y501 residues in the critical binding region for the XBB.1 ( Figure 9C) and XBB.1.5 complexes ( Figure 9D). Importantly, network modeling highlighted a noticeable drop in the stable 3-cliques formed in the flexible interface region for the XBB.1 variant, revealing that F486S mutation may not only compromise the strength of local binding interactions but also weaken the network of RBD interfacial contacts in this region ( Figure 9C). Strikingly, the network of 3-cliques in this region is fully restored and further enhanced in the XBB.1.5 variant as F486P mutation can reduce the flexibility in the interfacial interactions, likely providing an allosteric contribution to the improved binding affinity of XBB.1.5 ( Figure 9D).
The findings of this analysis also suggested that the extent of epistatic contributions in the Omicron RBD BA.2.75 and XBB.1.5 complexes may be stronger than in the other subvariants. It is worth emphasizing that the stabilizing cliques in the critical binding region are anchored by the hotspots R498 and Y501 which also emerged as potential epistatic hubs that can mediate strong dynamic and energetic couplings with other RBD residues. We also found that Y449, G496, and T505 are dynamically strongly coupled with Y501 through a persistent network of the interfacial cliques, confirming that these residues could form a second group of important epistatic centers. These observations are in agreement with functional studies which showed that the major sites exhibiting epistatic shifts in the presence of Y501 include G446, Y449, G496, Y499, and H505 residues [65,67]. Our results also indicated the presence of dynamic couplings between G446/Y449, R498/Y501 and R493 that may be important in mediating broad epistatic shifts which is consistent with the experimental data [67]. Hence, a clique-based network model can identify highly correlated and potentially non-additive mutational sites in the Omicron RBD complexes and distinguish them from other mutational sites that are less likely to experience epistatic shifts. These results provide a plausible rationale for the experimentally observed epistatic relationships in which mutations G446S, Q493R, and G496S individually reduce ACE2 binding but via strong epistasis with the pair R498/Y501 these losses can be fully compensated [67]. According to our findings, R498/Y501 can promote the formation of an extensive stable network of inter-connected 3-cliques that enables us to rescue the weaker binding potential of other mutational sites by amplifying their binding contributions. The non-additive effects are ensured by a chain of linked 3-cliques in which each pair of nodes/residues is connected by an edge, indicating a strong and intense mutual interaction among amino-acids on these nodes.
Another interesting finding of the network analysis is the role of F486, F486S and F486P mutations in BA.2.75, XBB.1 and XBB.1.5 variants respectively on modulating the density of the interfacial 3-cliques in the flexible interfacial region. Our analysis suggested that mutations of F486 can radically alter the strength and density of the interfacial cliques in this region and these F486-mediated changes are independent of the R498/Y501 mutations ( Figure 9). We found that F486S can appreciably reduce both local and global binding interactions in this region, while F486P can remarkably restore binding not only via the improved local packing but also by promoting stabilization of 3-cliques formed by RBD residues P486 N487 and ACE2 residues Q24 and Y83 ( Figure 9D). These findings are consistent with the experimentally established role of F486 as a critical evolutionary hotspot in which mutations at residue F486, such as F486V, F486I, and F486S, have been recurring among prior Omicron subvariants.
To summarize, the network-based community analysis provided additional insight into the mutational scanning data showing that mutational changes in these three positions may be coupled and lead to negative epistatic effects. However, importantly, these effects are secondary as the major non-additive contributions are likely to arise from the presence of a large number of stable cliques mediated by the Y501 position. Our findings suggested that the extent of non-additive contributions to the binding affinity may be greater for the Omicron BA.2.75 and XBB.1.5 complexes that displayed the strongest binding affinity among the examined Omicron subvariants.

Discussion
The results of this study provided molecular rationale and support to the experimental evidence that the acquisition of functionally balanced substitutions that optimize multiple fitness tradeoffs between immune evasion, high ACE2 affinity and sufficient conformational adaptability might be a common strategy of the virus evolution and serve as a primary driving force behind the emergence of new Omicron subvariants. Functional studies suggested that the evolutionary paths for significant improvements in the binding affinity of the Omicron RBD variants with hACE2 are relatively narrow and require balancing between various fitness tradeoffs of preserving RBD stability, maintaining binding to ACE2, and allowing for immune evasion [65][66][67]. These factors may limit the "evolutionary opportunities" for the virus to adapt new mutations that markedly improve ACE2 binding affinity without compromising immune evasion and stability. As a result, it led to a growing realization that evolutionary pressure invokes a complex interplay of thermodynamic factors to "designate" a privileged group of Omicron mutational hotspots that drive binding affinity with the ACE2, while allowing other Omicron sites to readily evolve immune escape capabilities with minor destabilizing liabilities. Moreover, some studies proposed that immune evasion may be a primary driver of Omicron evolution that sacrifices some ACE2 affinity enhancement substitutions to optimize immune-escaping mutations [36][37][38][39]. By examining forces driving the accelerated emergence of RBD mutations it was suggested that the immune pressure on the RBD becomes increasingly focused and promotes convergent evolution on the same sites including R346, K444, V445, G446, N450, L452, N460, F486, F490, R493, and S494 most of which are antibody-evasive [62]. These findings indicated that Omicron subvariants may evolve to accumulate convergent escape mutations while protecting and maintaining mutations that enable sufficient ACE2-binding capability. Analysis of the convergent evolution provided a useful summary of the observed tuning of the ACE2 binding affinity seen in new Omicron sub-lineages [63]. In particular, XBB.1 features E484A inherited from the BA.2 parent but further mutated into A484T in the child XBB.1.3 while XBB.1.5 adopted F486P mutation (F486S in XBB.1), and BA.2.75.2 inherited F486S from BA.2.75 but further mutated into F486L in the child CA.4 [63]. Remarkably, convergent evolution seen in these examples allowed for improved or neutral binding affinity changes and immune escape. Several lines of evidence indicated that the observed coordination of evolution at different sites is largely due to epistatic, rather than random selection of mutations [140].
Epistatic interactions between mutations add substantial complexity to their adaptive landscapes and are believed to play an important role in the virus evolution. The proposed in our study a network-based model for the analysis of non-additive contributions of the RBD residues indicated that some convergent Omicron mutations such as G446S (XBB.1.5) can display epistatic relationships with the major stability and binding affinity hotspots which may allow for the observed broad antibody resistance induced by these mutations [68]. Using atomistic simulations, the ensemble-based mutational scanning of binding/stability and network-based approaches, we showed that the binding affinity hotspots R498 and Y501 serve as central mediators of the interfacial communities in the RBD-ACE2 complexes. As a result, epistatic couplings mediated by R498 and Y501 hotspots with the RBD residues 491-496, 446-449 can be exemplified at the network level as these positions are involved in strong independent interactions within stable network cliques. This may allow for moderate negative effects on ACE2 binding in various Omicron immune evasion sites to be mitigated by strong compensatory epistasis exhibited by Y501.
An important lesson of this study is the spatially localized dependence of mutational effects on preexisting mutations R498 and Y501. According to the experimental data, while the effect of V445P and G446S mutations in XBB.1/XBB.1.5 can be compensated through epistatic couplings with R498/Y501, a single mutation F486S in XBB.1 results in an appreciable loss of binding affinity and could not be offset by the presence of the background R498/Y501 pair. Our results suggested that the restored binding strength mediated by F486P mutation in XBB.1.5 variant may arise from a spatially localized redistribution of local network cliques in the flexible interface region which is located on the other side of the binding interface from the Y501 position. In global epistasis, the fitness effect of a particular mutation can be determined by the fitness of its genetic background. According to presented evidence, pairs of Omicron substitutions with strong epistasis tend to be spatially proximal, and form localized stable modules allowing for compensatory energetic changes. Previous studies noted similar patterns of epistatic substitutions co-occurring spatially more frequently than expected if the substitutions had occurred randomly [141]. Statistical sequence-based landscape analysis based on the direct coupling analysis (DCA) incorporated pairwise epistatic terms which allowed us to capture local evolutionary constraints specific to the SARS-CoV-2 sequence background and identify K417, N440, E484, Q493, Q498, and N501 as sites of mutational enrichment [142]. Other studies have also found longer-range signals of co-occurring substitutions where allostery and epistasis conspire to facilitate the evolution of new functions through coordinated mutations at distal sites [143]. The results of this study are consistent with the idea that protein stability in key sites can promote evolvability via epistasis and enhance tolerance to destabilizing mutations that often contribute to immune escape [137][138][139]. During the evolution of a virus, a mutation that helps the virus evade the human immune system might only be tolerated if the virus has acquired other mutations beforehand. This type of mutational interaction between stability hotspots and evolvability sites would constrain the evolution of the virus, since its capacity to take advantage of the second mutation depends on the first mutation having already occurred. In general. the interactions between Omicron mutational sites may be controlled by spatially localized compensatory epistatic relationships in which key binding hotspots can rescue binding affinities to offset the effects of destabilizing immune escape mutations.

Conclusions
In this study, we systematically examined conformational dynamics, stability and binding of the Omicron RBD BA.2, BA.2.75, XBB.1 and XBB.1.5 RBD complexes with ACE2 using multiscale molecular simulations, in silico mutational scanning of the RBD residues and network-based community analysis of allosteric communications and epistatic interactions. Using a multiscale simulation approach and conformational landscapes derived from all-atom MD simulations, we found a progressive rigidification of the RBD-ACE2 binding interface regions in the BA.2.75 and XBB.1.5 variants as well as the improved stability of the RBD core regions. The results of simulations showed that the improved RBD stability in BA.2.75 and XBB.1.5 complexes with ACE2 may be one of the factors linked with the experimentally observed enhancement in their binding affinity with ACE2. By using distance fluctuations analysis of rigidity/flexibility in the complexes, we found that the important RBD binding interface centers R498, Y501 and H505 featured high stability indexes, reflecting a considerable rigidification of these residues due to strong interactions with ACE2. Our results suggested that these RBD residues function as key stability centers, binding hotspots as well as allosteric mediators of long-range communications in all RBD-ACE2 complexes. A systematic mutational scanning of the RBD residues in the complexes with ACE2 identified a conserved group of protein stability centers and binding affinity hotspots that determine the binding thermodynamics. Our data provided support to the emerging mechanism that the XBB lineage may have evolved to evade immune suppression and outcompete other Omicron subvariants through mutations of F486 which is a hotspot for establishing protective immunity against the virus. Consistent with the experimental data, our results revealed that functionally beneficial effects of destabilizing mutations may be contingent on compensatory effects provided by a binding hotspot cluster centered on the R498/Y501 positions. These findings highlighted the important feature of a mechanism in which key binding affinity hotspots R498 and Y501 enable compensatory epistatic interaction with other binding-neutral Omicron positions to balance tradeoffs between immune evasion, protein stability and ACE2 binding. To characterize and rationalize the experimentally observed epistatic effects of the Omicron mutations we explored the network analysis and proposed a clique-based network model for describing non-additive effects of the RBD residues. We found that the network analysis and community-based assessment of dynamic and energetic couplings can identify highly inter-dependent mutational sites in the Omicron RBD complexes that experience epistatic shifts. These results provided a plausible rationale for the experimentally observed epistatic relationships in which effects of Omicron mutations that can reduce ACE2 binding but are important for immune escape are compensated via strong epistatic interactions with the binding affinity hotspots R498 and Y501. The results of this study suggested distinct and yet complementary roles of the Omicron mutation sites forming a coordinated network of hotspots that enable efficient modulation and balance of multiple fitness tradeoffs including structural stability, host receptor binding, immune evasion and conformational adaptability which create a complex functional landscape of virus transmissibility.