Comparative Perturbation-Based Modeling of the SARS-CoV-2 Spike Protein Binding with Host Receptor and Neutralizing Antibodies: Structurally Adaptable Allosteric Communication Hotspots Define Spike Sites Targeted by Global Circulating Mutations.

In this study, we used an integrative computational approach to examine molecular mechanisms and determine functional signatures underlying the role of functional residues in the SARS-CoV-2 spike protein that are targeted by novel mutational variants and antibody-escaping mutations. Atomistic simulations and functional dynamics analysis are combined with alanine scanning and mutational sensitivity profiling of the SARS-CoV-2 spike protein complexes with the ACE2 host receptor and the REGN-COV2 antibody cocktail(REG10987+REG10933). Using alanine scanning and mutational sensitivity analysis, we have shown that K417, E484, and N501 residues correspond to key interacting centers with a significant degree of structural and energetic plasticity that allow mutants in these positions to afford the improved binding affinity with ACE2. Through perturbation-based network modeling and community analysis of the SARS-CoV-2 spike protein complexes with ACE2, we demonstrate that E406, N439, K417, and N501 residues serve as effector centers of allosteric interactions and anchor major intermolecular communities that mediate long-range communication in the complexes. The results provide support to a model according to which mutational variants and antibody-escaping mutations constrained by the requirements for host receptor binding and preservation of stability may preferentially select structurally plastic and energetically adaptable allosteric centers to differentially modulate collective motions and allosteric interactions in the complexes with the ACE2 enzyme and REGN-COV2 antibody combination. This study suggests that the SARS-CoV-2 spike protein may function as a versatile and functionally adaptable allosteric machine that exploits the plasticity of allosteric regulatory centers to fine-tune response to antibody binding without compromising the activity of the spike protein.

The coronavirus disease 2019 (COVID-19) pandemic associated with the severe acute respiratory syndrome (SARS)[1−5] has been at the focal point of biomedical research. SARS-CoV-2 infection is transmitted when the viral spike (S) glycoprotein binds to the host cell receptor, leading to the entry of S protein into host cells and membrane fusion.[6−8] The full-length SARS-CoV-2 S protein consists of two main domains, amino (N)-terminal S1 subunit and carboxyl (C)-terminal S2 subunit. The subunit S1 is involved in the interactions with the host receptor and includes an N-terminal domain (NTD), the receptor-binding domain (RBD), and two structurally conserved subdomains (SD1 and SD2). Structural and biochemical studies have shown that the mechanism of virus infection may involve spontaneous conformational transformations of the SARS-CoV-2 S protein between a spectrum of closed and receptor-accessible open forms, where RBD continuously switches between “down” and “up” positions where the latter can promote binding with the host receptor ACE2.[9−11] The crystal structures of the S-RBD in the complexes with human ACE2 enzyme revealed structurally conserved binding mode shared by the SARS-CoV and SARS-CoV-2 proteins in which an extensive interaction network is formed by the receptor binding motif (RBM) of the RBD region.[12−16] These studies established that binding of the SARS-CoV-RBD to the ACE2 receptor can be a critical initial step for virus entry into target cells. The rapidly growing body of cryo-EM structures of the SARS-CoV-2 S proteins detailed distinct conformational arrangements of S protein trimers in the prefusion form that are manifested by a dynamic equilibrium between the closed (“RBD-down”) and the receptor-accessible open (“RBD-up”) form required for the S protein fusion to the viral membrane.[17−34] The cryo-EM characterization of the SARS-CoV-2 S trimer demonstrated that S protein may populate a spectrum of closed states by fluctuating between structurally rigid locked-closed form and more dynamic, closed states preceding a transition to the fully open S conformation.[26] Conformational dynamics of SARS-CoV-2 trimeric spike glycoprotein in complex with receptor ACE2 suggesting considerable conformational heterogeneity of ACE2-RBD and continuous swing motions of ACE2-RBD in the context of SARS-CoV-2 S trimer. According to these experiments, the associated ACE2-RBD is relatively dynamic, showing three major conformational states with different angle of ACE2-RBD to the surface of the S trimer.[33]

Cryo-EM structural studies also mapped a mechanism of conformational events associated with ACE2 binding, showing that the compact closed form of the SARS-CoV-2 S protein becomes weakened after furin cleavage between the S1 and S2 domains, leading to the increased population of partially open states and followed by ACE2 recognition that can accelerate the transformation to a fully open and ACE2-bound form priming the protein for fusion activation.[34] These studies confirmed a general mechanism of population shifts between different functional states of the SARS-CoV-2 S trimers, suggesting that RBD epitopes can become stochastically exposed to the interactions with the host receptor ACE2.

The biochemical and functional studies using protein engineering and deep mutagenesis have quantified binding mechanisms of SARS-CoV-2 interactions with the host receptor.[35,36] Deep mutational scanning of SARS-CoV-2 RBD revealed constraints on folding and ACE2 binding showing that many mutations of the RBD residues can be well tolerated with respect to both folding and binding. A surprisingly large number of amino acid modifications could even improve ACE2 binding, including important binding interface positions that enhance RBD expression (V367F and G502D) or enhance ACE2 affinity (N501F, N501T, and Q498Y).[35] This comprehensive mutational scanning of the SARS-CoV-2 RBD residues suggested the evolutionary potential for compensation of deleterious mutations in the ACE2 interface reminiscent of multistep escape pathways and highlighted the energetic plasticity of the SARS-CoV-2 interaction networks in which mutations may enhance binding affinity, thus providing a roadmap for quantifying map immune-escape mutations. Using deep mutagenesis, it was also demonstrated that human ACE2 is only suboptimal for binding of the SARS-CoV-2 RBD as ACE2 variants near the interface can result in improved binding and simultaneously enhance folding stability.[36] Mutational landscape analysis showed a significant number of ACE2 mutations at the interface that enhance RBD binding, and the molecular basis for affinity enhancement can be rationalized from the structural analysis.[36] Functional studies characterized the key amino acid residues of the RBD for binding with human ACE2 and neutralizing antibodies, revealing two groups of amino acid residues to modulate binding, where the SARS-CoV-2 RBD mutations N439/R426, L452/K439, T470/N457, E484/P470, Q498/Y484, and N501/T487 can result in the enhanced binding affinity for ACE2.[37] Additionally, A475 and F486 in the SARS-CoV-2 RBD were identified as the key residues for the recognition of both their common functional receptor ACE2 and neutralizing antibodies, suggesting structural and energetic plasticity of the RBM residues involved in ACE2 recognition may induce mutational escape from the neutralizing antibodies targeting the RBD regions.

The rapidly growing structural studies of SARS-CoV-2 antibodies (Abs) have delineated molecular mechanisms underlying binding competition with the ACE2 host receptor, showing that combinations of Abs can provide broad and efficient cross-neutralization effects through synergistic targeting of conserved and variable SARS-CoV-2 RBD epitope.[38−49] Structural studies confirmed that the SARS-CoV-2 S protein can feature distinct antigenic sites, and some specific Abs may allosterically inhibit the ACE2 receptor binding without directly interfering with ACE2 recognition.[44] The SARS-CoV-2 Abs are divided into several main classes, of which class 1 and class 2 antibodies target epitopes that overlap with the ACE2 binding site.[50,51] The structural studies revealed binding epitopes and binding mechanisms for a number of newly reported SARS-CoV-2 Abs targeting RBD regions and competing with ACE2 include B38 and H14 Abs,[52] P2B-2F6,[53] CA1 and CB6,[54] CC12.1 and CC12.3,[55] C105,[56] and BD-23 Ab.[57] The crystal structure of a neutralizing Ab CR3022 in the complex with the SARS-CoV-2 S-RBD revealed binding to a highly conserved epitope that is located away from the ACE2-binding site but could only be accessed when two RBDs adopt the “up” conformation.[42] Subsequent structural and surface plasmon resonance studies confirmed that CR3022 binds the RBD of SARS-CoV-2, displaying strong neutralization by allosterically perturbing the interactions between the RBD regions and ACE2 receptor.[43] The crystal structure of an RBD-EY6A Fab complex identified the highly conserved epitope located away from the ACE2 binding site, showing that EY6A can compete with CR3022 by targeting residues that are important for stabilizing the prefusion S conformation.[58]

The B.1.1.7 variant of the SARS-CoV-2, a descendant of the D614G lineage, has originated in the UK and spread to 62 countries showing the increased transmissibility. Eight of the 17 mutations observed in this variant are accumulated in the S protein, featuring most prominently N501Y mutation that can increase binding affinity with ACE2 while eliciting immune escape and reduced neutralization of RBD-targeting Abs.[59−61] A new SARS-CoV-2 lineage (501Y.V2) first detected in South Africa is characterized by 21 mutations with 8 lineage-defining mutations in the S protein, including three at important RBD residues (K417N, E484K, and N501Y) that have functional significance and often induce significant immune escape.[62,63] Finally, the recently discovered new lineage, named P.1 (descendent of B.1.1.28), was observed in December in Brazil and contains a constellation of lineage defining mutations, including several mutations of known biological importance such as E484K, K417T, and N501Y mutations.[64,65]

Functional mapping of mutations in the SARS-CoV-2 S-RBD that escape antibody binding using deep mutational scanning showed that the escape mutations cluster on several surfaces of the RBD and have large effects on antibody escape while a negligible negative impact on ACE2 binding and RBD folding.[66] This illuminating study demonstrated that escape sites from antibodies can be constrained with respect to their effects on expression of properly folded RBD and ACE2 binding, suggesting that escape-resistant antibody cocktails can compete for binding to the same RBD region but have different escape mutations, which limit the virus ability to acquire novel sites of immune escape in the RBD without compromising its binding to ACE2.[66] Comprehensive mapping of mutations in the SARS-CoV-2 RBD that affect recognition by polyclonal human serum antibodies revealed that mutations in E484 site tend to have the largest effect on antibody binding to the RBD,[67] and various functional neutralization assay experiments indicated that E484 modifications can reduce the neutralization potency by some antibodies by >10-fold.[67−69] These studies also indicated that K417N and N501Y mutants can escape neutralization by some monoclonal antibodies but typically only modestly affected binding.[67,70] At the same time, mutations in the epitope centered around E484 position (G485, F486, F490) or in the 443–455 loop (K444, V445, L455, F456 sites) strongly affected neutralization for several Abs.[67−71] Functional mapping of the SARS-CoV-2 RBD residues that affect the binding of the REGN-COV2 cocktail showed that REGN10933 and REGN10987 are escaped by different mutations as mutation at F486 escaped neutralization only by REGN10933, whereas mutations at K444 escaped neutralization only by REGN10987, while E406W escaped both individual REGN-COV2 antibodies.[70] This study confirmed that escape mutations at Q493, Q498, and N501 sites may enhance binding affinity with ACE2 and that escape mutations can also emerge in positions distant from the immediate proximity of the binding epitope, highlighting structural and energetic plasticity of the RBD regions and potential allosteric-based mechanism of immune escape.[70] The REGN-COV2 cocktail (REG10987+REG10933) demonstrated significant potential in preventing mutational escape,[72] several other antibody cocktails such as COV2–2130+COV2–2196,[73] BD-368–2+BD-629,[74] and B38+H4[52] displayed promising neutralization activities. Analysis of the molecular determinants and mechanisms of mutational escape showed that SARS-CoV-2 virus rapidly escapes from individual antibodies but does not easily escape from the cocktail due to stronger evolutionary constraints on RBD-ACE2 interaction and RBD protein folding.[75] According to this study, the key RBD positions critical for the escape of antibody combinations include K444, which is an important epitope residue for CoV2–06, P2B-2F6, S309 and REG10987 Abs, as well as E484/F486 sites that are central for binding of CoV2–14 and REG10933. Functional analysis validated that mutations of these residues are responsible for viral escape from the individual Abs and, in combination with other currently circulating variants (N501Y, K417N, E484K), may induce the reduced neutralization by the antibody cocktails.[75] The SARS-CoV-2 501Y.V2 lineage that includes one cluster in NTD with four substitutions and a deletion (L18F, D80A, D215G, Δ242–244, and R246I), and another cluster of three substitutions in RBD (K417N, E484K, and N501Y) can confer neutralization escape from SARS-CoV-2 directed monoclonal antibodies and significantly increased neutralization resistance from individuals previously infected with SARS-CoV-2 virus.[76] Moreover, three class 1 antibodies (CA1, CB6, and CC12.1) that target the ACE2-binding RBM region showed a complete lack of binding for the 501Y.V2 variant, suggesting that mutations in the RBD and NTD clusters may amplify the mutational escape from RBD-targeted Abs.

Structural and functional studies showed that the activity of mRNA vaccine-elicited antibodies to SARS-CoV-2 and circulating variant encoding E484K or N501Y or the K417N/E484K/N501Y combination can be reduced by a small but significant margin, suggesting that these mutations in individuals with compromised immunity may erode the effectiveness of vaccine-elicited immunity.[77] Importantly, it was also shown that neutralization by 14 of the 17 most potent tested mAbs can be partly reduced or even abolished by either K417N, or E484K, or N501Y mutations. Another latest study reported the preserved neutralization of N501Y, Δ69/70 + N501Y + D614G and E484K + N501Y + D614G viruses by BNT162b2 vaccine-elicited human sera.[77] Consistent with other studies, it was shown that the neutralization against the virus with three mutations from the SA variant (E484K + N501Y + D614G) was slightly lower than the neutralization against the N501Y virus and the virus with three UK mutations (Δ69/70 + N501Y + D614G), but these differences were relatively small.[78] SARS-CoV-2-S pseudoviruses bearing either the reference strain or the B.1.1.7 lineage spike protein with sera of 40 participants who were vaccinated with the mRNA-based vaccine BNT162b2 showed largely preserved neutralization, indicating that the B.1.1.7 lineage will not escape BNT162b2-mediated protection.[79] The recent data demonstrate reduced but still significant neutralization against the full B.1.351 variant following mRNA-1273 vaccination.[80] New SARS-CoV-2 variants that resist neutralizing antibodies are now emerging in low frequencies in circulating SARS-CoV-2 populations. In particular, recent reports presented evidence of circulating SARS-CoV-2 spike N439 K variants evading antibody-mediated immunity, particularly N439 K mutation that confers resistance against several neutralizing monoclonal antibodies and reduces the activity of mRNA vaccine-elicited antibodies.[81] Computational modeling and molecular dynamics (MD) simulations have been instrumental in predicting conformational and energetic determinants of SARS-CoV-2 mechanisms and the binding affinity and selectivity with the host receptor ACE2.[82−93] Molecular mechanisms of the SARS-CoV-2 binding with ACE2 enzyme were analyzed in our recent study using coevolution and conformational dynamics.[94] Using protein contact networks and perturbation response scanning based on elastic network models, we recently discovered the existence of allosteric sites on the SARS-CoV- 2 spike protein.[95] By using molecular simulations and network modeling, we recently presented the first evidence that the SARS-CoV-2 spike protein can function as an allosteric regulatory engine that fluctuates between dynamically distinct functional states.[96] Coarse-grained normal-mode analyses combined with Markov model and computation of transition probabilities characterized the dynamics of the S protein and the effects of mutational variants D614G and N501Y on protein dynamics and energetics.[97] Using time-independent component analysis (tICA) and protein networks, another computational study identified the hotspot residues that may exhibit long-distance coupling with the RBD opening, showing that some mutations may allosterically affect the stability of the RBD regions.[98] Molecular simulations reveal that N501Y mutation increases ACE2 binding affinity and may impact the collective dynamics of the ACE2-RBD complex while mutations K417N and E484K reduce the ACE2-binding affinity by abolishing the interfacial salt bridges.[99] The growing body of computational modeling studies investigating dynamics and molecular mechanisms of SARS-CoV-2 mutational variants produced inconsistent results that propose different mechanisms. The development of a more unified view and a working theoretical model that can explain the diverse experimental observations is an important area of current efforts in the field.

In this study, we used an integrative computational approach to examine molecular mechanisms underlying the functional role of K417, N439, E484, and N501 positions targeted by novel mutational variants in the SARS-CoV-2 S protein. We combined coarse-grained (CG) simulations and atomistic reconstruction of dynamics trajectories with dynamic fluctuation communication analysis, mutational sensitivity analysis, and network community modeling to examine complexes of the SARS-CoV-2 S-RBD and dissociated S1 domain of the S protein formed with the ACE2 host receptor. Using distance fluctuations communication analysis and functional dynamics analysis, we determine and compare the distribution of regulatory centers in the RBD complexes with ACE2 and REGN-COV2 antibody cocktail (REG10987+REG10933).

Using alanine scanning and mutational sensitivity analysis, we show that K417, E484, and N501 residues correspond to key interacting centers with a significant degree of structural and energetic plasticity that allow mutants in these positions to afford the improved binding affinity with ACE2. Through perturbation-based network modeling and community analysis of the SARS-CoV-2 RBD complexes with ACE2, we demonstrate that E406, N439, K417, and N501 residues serve as effector centers of allosteric interactions and anchor major intermolecular communities that mediate long-range communication in the complexes. The results of the comparative network analysis with antibody complexes show that mutations in these positions can alter structural arrangements with antibodies and compromise their neutralization effects. These results suggest that antibody-escaping mutations target allosteric mediating hotspots with sufficient plasticity and adaptability to modulate and improve binding and allosteric signaling functions with the host receptor activity while reducing the efficiency of antibody recognition and long-range communications. This analysis suggests that the SARS-CoV-2 S protein may function as a versatile and functionally adaptable allosteric machine that exploits the plasticity of allosteric regulatory centers to generate escape mutants that fine-tune response to antibody binding without compromising the activity of the spike protein.

Materials and Methods

Coarse-Grained Molecular Simulations

Coarse-grained (CG) models are computationally effective approaches for simulations of large systems over long time scales. In this study, CG-CABS model[100−104] was used for simulations of the crystal structures of the SARS-CoV-2 RBD complex with ACE2 (PDB id 6M0J)[15] and complexes formed by the dissociated S1 domain of SARS-CoV-2 Spike bound to ACE2 (PDB id 7A91, 7A92)[34] (Figure A–C). We also simulated the cryo-EM structure of the SARS-CoV-2 RBD complex with the Fab fragments of two neutralizing antibodies, REGN10933 and REGN10987 (PDB id 6XDG)[105] (Figure D–F). In this high-resolution model, the amino acid residues are represented by Cα, Cβ, the center of mass of side chains and another pseudoatom placed in the center of the Cα-Cα pseudobond. In this model, the amino acid residues are represented by Cα, Cβ, the center of mass of side chains and the center of the Cα-Cα pseudobond. The CABS-flex approach implemented as a Python 2.7 object-oriented standalone package was used in this study to integrate a high-resolution coarse-grained model with robust and efficient conformational sampling proven to accurately recapitulate all-atom MD simulation trajectories of proteins on a long time scale.[104] Conformational sampling in the CABS-flex approach is conducted with the aid of Monte Carlo replica-exchange dynamics and involves local moves of individual amino acids in the protein structure and global moves of small fragments.[100−102]

Figure 1

Crystal structures of the SARS-CoV-2 RBD and S1 domain complexes with ACE enzyme and REGN-COV2 antibody cocktail. (A) Structural overview of the SARS-CoV-2 RBD complex with ACE2 (PDB id 6M0J). The SARS-CoV RBD is shown in cyan ribbons, and the RBM region is in blue ribbons. The subdomain I of human ACE2 is shown in red ribbons, and subdomain II is shown in green ribbons. The structure of ACE2 consists of the N-terminus subdomain I (residues 19–102, 290–397, and 417–430) and C-terminus subdomain II (residues 103–289, 398–416, and 431–615) that form the opposite sides of the active site cleft. (B) The crystal structure of the dissociated S1 domain form in the complex with ACE2 (PDB id 7A91). S1-RBD is in cyan ribbons, and ACE2 is in green ribbons. (C) The crystal structure of the fully dissociated S1 domain in the complex with ACE2 (PDB id 7A92). The S1 domain of the SARS-CoV-2 S protein is in cyan ribbons, and ACE2 is in green ribbons. (D) The cryo-EM structure of the SARS-CoV-2 RBD in the complex with REGN10933/REGN10987 antibody cocktail. The RBD region is shown by the green surface. REGN10933 Fab fragment is shown in ribbons with the heavy chain in cyan and the light chain in blue. REGN10987 is in ribbons with the heavy chain in orange and the light chain in purple. The positions of functional residues targeted by mutational variants and antibody-escaping mutations are E406, K417, E484, and N501 and are annotated and highlighted as black patches on the RBD surface. (E) A close-up of the SARS-CoV-2 RBD interactions with REGN10933. The RBD is shown on the green surface. REGN10933 Fab fragment is shown in ribbons with heavy chain in cyan and light chain in blue. The REGN10933 antibody epitope on RBD is highlighted in cyan patches on the surface. The positions of E406, K417, E484, F486, N501 are shown as black surface patches on the RBD. (F) A close-up of the SARS-CoV-2 RBD interface with REGN10987. The red patches correspond to the REGN10987 epitope. The positions of E406, K417, N439, E484, F486, N501 are shown as black surface patches on the RBD.

The default settings were applied in which soft native-like restraints are imposed only on pairs of residues fulfilling the following conditions: the distance between their Cα atoms was smaller than 8 Å, and both residues belong to the same secondary structure elements. A total of 1,000 independent CG-CABS simulations were performed for each of the studied systems. In each simulation, the total number of cycles was set to 10,000, and the number of cycles between trajectory frames was 100. MODELER-based reconstruction of simulation trajectories to the all-atom representation provided by the CABS-flex package was employed to produce atomistic models of the equilibrium ensembles for studied systems.

Structure Preparation and Analysis

All structures were obtained from the Protein Data Bank.[106,107] 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.[108,109] The structures were further preprocessed through the Protein Preparation Wizard (Schrödinger, LLC, New York, NY) and included the check of bond order, assignment and adjustment of ionization states, formation of disulfide bonds, removal of crystallographic water molecules and cofactors, capping of the termini, assignment of partial charges, and addition of possible missing atoms and side chains that were not assigned in the initial processing with the WHATIF program. The missing loops in the studied crystal structures of the dissociated S1 domain complexes with ACE2 (residues 556–573, 618–632) were reconstructed and optimized using template-based loop prediction approaches ModLoop,[110] ArchPRED server[111] and further confirmed by FALC (Fragment Assembly and Loop Closure) program.[112] The side-chain rotamers were refined and optimized by SCWRL4 tool.[113] The shielding of the receptor binding sites by glycans is an important common feature of viral glycoproteins, and glycosylation on SARS-CoV proteins can camouflage immunogenic protein epitopes.[114,115] The atomistic structures from simulation trajectories of the dissociated S1 domain complex with ACE2 (PDB id 7A92) were elaborated by adding N-acetyl glycosamine (NAG) glycan residues and optimized. The glycosylated microenvironment for atomistic models of the simulation trajectories was mimicked by using the structurally resolved glycan conformations for most occupied N-glycans as determined in the cryo-EM structures of the SARS-CoV-2 spike S trimer in the closed state (K986P/V987P,) (PDB id 6VXX) and open state (PDB id 6VYB).

Functional Dynamics and Collective Motions Analysis

We performed principal component analysis (PCA) of reconstructed trajectories derived from CABS-CG simulations using the CARMA package[116] and also determined the essential slow mode profiles using elastic network models (ENM) analysis.[117] Two elastic network models: Gaussian network model (GNM)[117,118] and Anisotropic network model (ANM) approaches[119] were used to compute the amplitudes of isotropic thermal motions and directionality of anisotropic motions. The functional dynamics analysis was conducted using the GNM in which protein structure is reduced to a network of N residue nodes identified by Cα atoms and the fluctuations of each node are assumed to be isotropic and Gaussian. Conformational mobility profiles in the essential space of low frequency modes were obtained using ANM server[119] and DynOmics server.[120]

Local Structural Parameters: Relative Solvent Accessibility

We have computed the relative solvent accessibility parameter (RSA) that is defined as the ratio of the absolute solvent-accessible surface area (SASA) of that residue observed in a given structure and the maximum attainable value of the solvent-exposed surface area for this residue.[121] According to this model, residues are considered to be solvent-exposed if the ratio value exceeds 50% and to be buried if the ratio is less than 20%. Analytical SASA is estimated computationally using analytical equations and their first and second derivatives and was computed using web server GetArea.[121]

Mutational Sensitivity Analysis and Alanine Scanning

To compute protein stability and binding free energy changes in the SARS-CoV-2 RBD structures upon complex formation with ACE2 receptor and REGN-COV2 antibody cocktail, we conducted a systematic alanine scanning of protein residues in the SARS-CoV-2 RBD and S1 domain. In addition, a complete mutational sensitivity analysis was done for binding free energy hotspots and residues E406, K417, N439, K444, E484, F486, and N501 targeted by widely circulating and antibody-escaping mutations. Alanine scanning and mutational sensitivity profiling of protein residues were performed using the BeAtMuSiC approach.[122,123] If a free energy change between a mutant and the wild type (WT) proteins ΔΔG= ΔG (MT)-ΔG (WT) > 0, the mutation is destabilizing, while when ΔΔG < 0 the respective mutation is stabilizing. The BeAtMuSiC approach is based on statistical potentials describing the pairwise inter-residue distances, backbone torsion angles, and solvent accessibilities and considers the effect of the mutation on the strength of the interactions at the interface and on the overall stability of the complex.[122,123] The reported protein stability and binding free energy changes are based on the ensemble averages of BeAtMuSiC values using equilibrium samples from reconstructed simulation trajectories.

Perturbation Response Scanning

Perturbation Response Scanning (PRS) approach[124,125] was used to estimate residue response to external forces applied systematically to each residue in the protein system. This approach has successfully identified hotspot residues driving allosteric mechanisms in single protein domains and large multidomain assemblies.[126−131] The implementation of this approach follows the protocol originally proposed by Bahar and colleagues[126,127] and was described in detail in our previous studies.[96] In brief, through monitoring the response to forces on the protein residues, the PRS approach can quantify allosteric couplings and determine the protein response in functional movements. In this approach, it 3N × 3N Hessian matrix whose elements represent second derivatives of the potential at the local minimum connect the perturbation forces to the residue displacements. The 3N-dimensional vector Δ of node displacements in response to 3N-dimensional perturbation force follows Hooke’s law . A perturbation force is applied to one residue at a time, and the response of the protein system is measured by the displacement vector ΔR(i) = H that is then translated into N × N PRS matrix. The second derivatives matrix is obtained from simulation trajectories for each protein structure, with residues represented by Cα atoms and the deviation of each residue from an average structure was calculated by ΔR(t) = R(t) – ⟨R(t)⟩, and corresponding covariance matrix C was then calculated by ΔRΔRT. We sequentially perturbed each residue in the SARS-CoV-2 spike structures by applying a total of 250 random forces to each residue to mimic a sphere of randomly selected directions. The displacement changes, Δ is a 3N-dimensional vector describing the linear response of the protein and deformation of all the residues.

Using the residue displacements upon multiple external force perturbations, we compute the magnitude of the response of residue k as ⟨∥ΔR(∥2⟩ averaged over multiple perturbation forces F(, yielding the ikth element of the N × N PRS matrix.[126,127] The average effect of the perturbed effector site i on all other residues is computed by averaging over all sensors (receivers) residues j and can be expressed as ⟨(Δ)2⟩effector. The effector profile determines the global influence of a given residue node on the perturbations in other protein residues and can be used as a proxy for detecting allosteric regulatory hotspots in the interaction networks. In turn, the jth column of the PRS matrix describes the sensitivity profile of sensor residue j in response to perturbations of all residues, and its average is denoted as ⟨(Δ)2⟩sensor. The sensor profile measures the ability of residue j to serve as a receiver of dynamic changes in the system.

Protein Structure Networks and Community Analysis

A graph-based representation of protein structures[132,133] is used to represent residues as network nodes and the inter-residue edges to describe noncovalent residue interactions. The details of graph construction using residue interaction cutoff strength (Imin) were outlined in our previous studies.[96] The network edges that define residue connectivity are based on noncovalent interactions between residue side chains that define the interaction strength I according to the following expression used in the original studies:[132,133]where n is the number of distinct atom pairs between the side chains of amino acid residues i and j that lie within a distance of 4.5 Å. N and N are the normalization factors for residues i and j. We constructed the residue interaction networks using both dynamic correlations[134] and coevolutionary residue couplings[135] that yield robust network signatures of long-range couplings and communications. The details of this model were described in our previous studies.[135−137] More specifically, the edges in the residue interaction network are then weighted based on dynamic residue correlations and coevolutionary couplings measured by the mutual information scores. 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) w = −log[ (, )] of the edge that connects nodes i and j is defined as the element of a matrix measuring the generalized correlation coefficient (, ) as between residue fluctuations in structural and coevolutionary dimensions. Network edges were weighted for residue pairs with (, ) > 0.5 in at least one independent simulation as was described in our initial study.[135] The matrix of communication distances is obtained using the generalized correlation between composite variables describing both dynamic positions of residues and coevolutionary mutual information between residues. As a result, the weighted graph model defines a residue interaction network that favors a global flow of information through edges between residues associated with dynamics correlations and coevolutionary dependencies. To characterize allosteric couplings of the protein residues and account for cumulative effect of dynamic and coevolutionary correlations, we employed the generalized correlation coefficient first proposed by Lange and Grubmüller.[138] The g_correlation tool in the Gromacs 3.3 package was used that allows computation of both linear or nonlinear generalized correlation coefficients.[139] The protocol was previously introduced and detailed in our earlier study,[135] showing that the generalized correlation coefficient based on dynamic and coevolutionary couplings provided a robust metric for detecting the cross-correlation between protein residues. A similar strategy for analysis of allosteric motions and interactions was successfully undertaken and improved in series of illuminating studies by Palermo and colleagues[140−143] where the introduced generalized correlation (GC) matrix proved to be a sensitive and accurate method for detecting the interdependence of spatially distant residues, providing a reliable and reproducible measure of how much the motion of one residue is dependent on the fluctuations of another spatially separated residue.

The RING program[144,145] was also employed for the initial generation and analysis of residue interaction networks. The ensemble of shortest paths is determined from the matrix of communication distances by the Floyd-Warshall algorithm.[146] Network graph calculations were performed using the python package NetworkX.[147] The betweenness of residue i is defined as the sum of the fraction of shortest paths between all pairs of residues that pass through residue i:where g denotes the number of shortest geodesics paths connecting j and k, and g (i) is the number of shortest paths between residues j and k passing through the node n.

The Girvan–Newman algorithm[148−150] is used to identify local communities. In this approach, edge centrality (also termed as 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 Girvan–Newman method was implemented, and the algorithmic details of this modified scheme were given in our recent studies.[151,152] Briefly, 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 highest edge betweenness within a given threshold; (c) recalculate edge betweenness for remaining edges; (d) repeat steps b–d until the graph is empty.

The residue betweenness is then used to rank the most influential nodes in the network and communities. For defining community leader nodes, we follow the Leader-Follower algorithm, in which a community is defined as a clique and is characterized by the presence of a leader and at least one “loyal follower”.[153] Community leaders are defined as nodes that (a) are connected not only to members of the local community but also have neighbors outside of the community; and (b) whose distance to other nodes in the network is less than the neighbors in their respective communities. These nodes could either directly link different communities or are connected to isolated bridging nodes between communities. A loyal follower in a community is defined as a residue node that only has neighbors within this single community. To characterize global bridges from a community structure, we introduce community bridgeness metric similar to Rao-Stirling index.[154−156] this parameter uses as input a prior categorization of the nodes into distinct communities:where the sum is over communities J (different from the community of node i, denoted as I), δ is equal to 1 if there is a link between node i and community J and 0 otherwise. l corresponds to the effective distance between community I and community J as measured by the inverse of the number of links between them. Nodes that are only linked to nodes of their own community, i.e., loyal follower nodes have G(i) = 0, while community leader nodes involved in bridging two (or more) communities have a positive value of the index. All topological measures were computed using the python module the python package NetworkX[147] and Cytoscape platform for network analysis.[157]

Results and Discussion

Conformational Dynamics Profiles of the SARS-CoV-2 S RBD and Dissociated S1 Domain Binding with ACE2: Balancing Structural Rigidity and Plasticity at Binding Interfaces

Multiple CG-CABS simulations of the SARS-CoV-2 RBD and S1 domain complexes with the ACE2 host receptor were performed to analyze similarities and differences in th conformational dynamics profiles of the RBD regions and specifically binding interface residues (Figure ). We combined CG-CABS simulations with the atomistic reconstruction of simulation trajectories to characterize regions of structural stability and plasticity at the ACE2 binding interfaces and determine the effect of the complete S1 domain in the complex on flexibility of the interacting residues and functional RBD regions. To our knowledge, this is the first comparative computational analysis of S1 domain binding with ACE2 initiated to understand structural and dynamic rearrangements of the S1 domain to form a stable monomeric complex with ACE2. The conserved core of SARS-CoV-RBD consists of five antiparallel β strands with three connecting α-helices (Figures and 2). The central β strands (residues 354–363, 389–405, 423–436) in SARS-CoV-2 RBD are stable and, as expected, only small thermal fluctuations were observed in these regions (Figure A). The antiparallel β-sheets (β5 and β6) (residues 451–454 and 492–495 in SARS-CoV-RBD) that anchor the RBM region to the central core also displayed a significant stabilization in the complex with ACE2. The small α-helical segments of the RBD (residues 349–353, 405–410, and 416–423) also displayed significant stability in simulations. These regions become even more rigidified in the complex formed by the dissociated S1 domain (Figure B). Of special interest were the ACE2 induced changes in the RBM region (residues 437–508) and particularly in a stretch of residues 471–503 involved in multiple contacts with the ACE2 receptor. The overall similar RBM profiles in all systems were seen, but the greater stabilization of the interfacial residues in the ACE2 complex was observed for the bound S1 domain structure (Figure A,B). The interfacial loop residues 436–455 containing an important motif 444-KVGGNYNY-451 displayed significantly reduced fluctuations. Among residues that experience a more pronounced stabilization in the complex are K417, G446, Y449, Y453, L455, F456, Y473, A475, and G476 positions in the middle segment of the RBM (Figure C-E). The analysis of the intermolecular contacts in the SARS-CoV-2 RBD and S1 complexes with ACE2 (Supporting Information, Tables S1–S3) indicated a significant number of RBM interactions formed by the middle segment of the interface (K417, Y453, L455, F456, and Q493) with K31 and E35 of ACE2 which may explain a pronounced stabilization of these positions in simulations. Indeed, K417 forms contacts with D30 and H34 hotspot ACE2 residues, L455 is involved in interactions with K31 and D30, and F456 forms stabilizing contacts with D30, K31, and T27 ACE2 hotspots (Supporting Information, Tables S1–S3). Another group of residues in the RBM ridge involved in multiple binding contacts included E484 and F486 sites that interact with K31, Q24, M82, and L79 of ACE2 (Supporting Information, Tables S1–S3). Notably, this analysis showed that the interacting RBD motif 495-YGFQPTNG-502 is involved in the most persistent interaction contacts with ACE2 that is exemplified by stabilization of Y489, F490, Q493, Y495, G496, Q498, T500, and Y505 residues (Figure ). These RBD residues form the largest number of contacts with ACE2 (Supporting Information, Tables S1–S3) and experienced the most significant stabilization in the complex (Figure A,B). Importantly, not only these residues showed markedly reduced fluctuations but also the large interfacial stretch of residues across the entire binding interface (residues 486-FNCYFPLQSYGFQ-498) including key Q493 and Q498 interacting sites exhibited even a stronger stabilization in the complex formed by the dissociated S1 domain (Figure B).

Figure 2

CABS-GG conformational dynamics of the SARS-CoV-2 S-RBD and S1 domain complexes with ACE2. (A) The root-mean-square fluctuations (RMSF) profiles from simulations of the structures of the SARS-CoV-2 S-RBD complex with ACE2, PDB id 6M0J (in orange lines), and S1 domain complex with ACE2, PDB id 7A91 (in maroon lines). (B) The RMSF profiles from simulations of the structures of the SARS-CoV-2 S-RBD complex with ACE2, PDB id 6M0J (in orange lines), and complete S1 domain complex with ACE2, PDB id 7A92 (in maroon lines). The S-RBD sites L455, F456, S459, Q474, A475, F486, F490, Q493, and P499, whose mutations can abolish binding affinity with ACE2, are shown in yellow filled circles. The S-RBD sites N439, L452, E484/P470, Q498, and N501, whose mutations could enhance binding affinity for ACE2, are shown in blue filled circles. The S-RBD sites E406, N439, K417, E484, and N501 targeted by novel circulating mutations and antibody-escaping mutations are highlighted in red filled squares. (C) Structural mapping of the conformational dynamics profiles in the SARS-CoV-2 S–RBD complex with ACE2 (PDB id6M0J). A ribbon-based protein representation is used with coloring (blue-to-red) according to the protein residue motilities (from more rigid–blue regions to more flexible–red regions). (D, E) Structural mapping of the conformational dynamics profiles in the S1 domain complexes with ACE2 (PDB id 7A91 and 7A92, respectively). The stability profile for protein residues is shown as in panel C using a coloring spectrum from blue to red to highlight changes from rigid to more flexible regions.

Structural maps of the conformational dynamics profiles highlighted these observations showing a more uniform and broad stabilization of the RBD regions in the complex formed by the dissociated S1 domain (Figure C-E). Although the overall fluctuation profile of the RBM residues remained largely unchanged, we noticed small but important differences pointing to the greater stability of the 495-YGFQPTNG-502 loop in the ACE2 complex with the dissociated S1 domain. These observations are consistent with the latest structural studies showing that ACE2 binding can induce disassembly of the SARS-CoV-2 S trimer and promote the formation of a stable dissociated monomeric S1 complex with ACE2 receptor. At the same time, a modest mobility of the RBM positions N439, L452, T470, E484, Q498, and N501 was seen in the S1-ACE2 complex (Figure B). Importantly, several of these positions, N49, E484, and N501, correspond to sites that confer mutational variants with the increased binding to ACE2 and elevated level of transmission and infectivity. The observed partial flexibility of this SARS-CoV-2 RBM motif in the complex formed by the dissociated S1 domain may allow for tolerance and adaptability of these sites to specific modifications resulting in the improved binding affinity with ACE2.

We also report the relative solvent accessibility (RSA) ratio in the SARS-CoV-2 RBD and S1 domain complexes with ACE2 that were obtained by averaging the SASA computations over the simulation trajectories (Supporting Information, Figure S1). The antiparallel β-sheet regions in the SARS-CoV-2 (residues 451–454 and 492–495) are deeply buried at the interface. The key RBM residues in the central segment (K417, L456, F456, Y473, F490, and Q493) also showed small RSA values, indicating that these positions are buried in the ACE2 complex. Of particular interest were the average RSA values for functional sites E406, K417, N439, E484, and N501 that are the central focus of our investigation. We found that E406, K417, and N439 showed moderate RSA values (∼20–30%) indicating that these positions could retain a certain degree of plasticity in the RBD and S1 domain complexes with ACE2 (Supporting Information, Figure S1). The more extreme cases were exemplified by E484 that maintains significant solvent exposure (RSA ∼ 65%) in the ACE2 complexes, while N501 is largely buried with very small RSA values (Supporting Information, Figure S1). However, conformational dynamics profiles indicated that N501 may still maintain some level of plasticity in the ACE2 complex.

To compare the differences in the local flexibility with experimental functional data, we specifically analyzed a group of SARS-CoV-2 RBD residues L455, F456, S459, Q474, A475, F486, F490, Q493, and P499 whose mutations to their SARS-CoV RBD counterpart positions resulted in the abolished binding affinity.[37] It could be noticed that, in the S1 domain complex with ACE2, these residues become appreciably more stable than residues from another group (N439, L452, T470, E484, Q498, N501) that are more susceptible to affinity-improving mutations. Functional studies showed that N439/R426, L452/K439, T470/N457, E484/P470, Q498/Y484, and N501/T487 modifications of these SARS-CoV-2 RBD residues to their respective position in SARS-CoV-RBD can, in fact, result in the enhanced binding affinity for ACE2.[37] Our analysis allowed to capture these subtle differences showing that this group of RBD residues may experience larger fluctuations (Figure A,B). Interestingly, these differences become more evident only in the ACE2 complex with the dissociated S1 domain, suggesting that the partial redistribution of mobility in the S1-ACE2 complex could provide more room for structural adaptation of N439, E484, and N501 positions (Figure B). Hence, a moderate level of residual fluctuations can be preserved even when RBM residues are involved in strong stabilizing contacts with ACE2. Although E484 interacts with the K31 interaction hotspot residue of hACE2, this residue retains a more significant degree of mobility and plasticity in the RBM region which may be associated with the mutational variability and emergence of the E484K variant that can improve binding affinity with the host receptor. Interestingly, we found that escape mutations and variants improving binding affinity with the ACE receptor may emerge in sites that are moderately flexible in the S1-ACE2 complex. This suggested that, although some of these positions such as K417 and N501 are involved in multiple contacts with ACE2, there should be substantial energetic plasticity in the interaction network. According to our findings, there may be more room for tolerant modifications of N439 and E484 positions, while the potential for favorable mutations at K417 and N501 sites could be more limited.

To summarize, the central finding of this analysis is that the conformational dynamics profile for the SARS-CoV-2 S-RBD residues remained largely conserved among ACE2 complexes of the S1-RBD and fully dissociated S1 domain. Another important observation is a consistent trend for moderate residual mobility of RBD residues whose mutations may often lead to the enhanced binding with ACE2. Conformational dynamics analysis also indicated that RBM residues targeted by novel mutational variants may be adaptable and display a range of flexibility, from more dynamic positions at N439 and E486 to more constrained K417 and N501 residues.

Essential Dynamics of the SARS-CoV-2 S Complexes with ACE2 and the REGN-COV2 Cocktail Unveils Regulatory Roles of Functional Sites Targeted by Mutational Variants

To characterize collective motions and determine the distribution of hinge regions in the SARS-CoV-2 S-RBD and SARS-CoV-2 S1 domain complexes with ACE2 (Figure ) and REGN-COV2 antibody combination (Figure ), we performed PCA of trajectories derived from CABS-CG simulations and also determined the essential slow modes using ENM analysis. The reported functional dynamics profiles were averaged over the first three major low-frequency modes. For comparison of ACE2-induced changes in the functional dynamics, we first characterized the slow mode profiles for the unbound forms of the RBD and S1 domain extracted from the crystal structures of the complexes (Figure A–C). For the unbound RBD structure, the local minima associated with local hinge points corresponded to W353, F374, F400, L452, R466, Q493, and V510 residues. Some of these residues, L452 and Q493, are involved in the interactions in the complex, and their immobilized hinge position in the unbound form may be important to induce the optimal intermolecular association with ACE2. Interestingly, none of the functional positions targeted by novel mutational variants that promote infectivity and antibody resistance (E409, K417, N439, E484, and N501) corresponded to hinge positions in the unbound RBD form. In fact, E484 residue is located in the moving region of the unbound RBD structure in slow modes (Figure A). The unbound form of the dissociated S1 domain featured local hinge positions in residues F518 and V539, C538, and F592 (Figure B,C). Strikingly, these findings are consistent with our recent analysis of collective dynamics in the open and closed forms of the SARS-CoV-2 S trimer structures showing that these residues correspond to major regulatory centers of functional motions and coordinate global displacements of the S1 regions with respect to more rigid S2 subunit in distinct timer states.[158]

Figure 3

Functional dynamics of the SARS-CoV-2 S-RBD and S1 domain complexes with ACE2. The mean square displacements in functional motions are averaged over the three lowest frequency modes. The essential mobility profiles for the unbound forms of the SARS-CoV-2 S-RBD (A) and dissociated forms of the S1 domain (B, C). (D) The slow mode profile of the bound SARS-CoV-2 S-RBD structure in the complex with ACE2 (PDB id 6M0J). (E, F) The slow mode profiles of the dissociated S1 domain complexed with ACE2. The essential mobility profiles are shown in maroon lines, and positions of key functional residues E406, K417, N439, E484, and N501 are highlighted by filled green circles. Structural maps of the essential mobility profiles for the SARS-CoV-2 S-RBD complex with ACE2 (G) and complexes formed by the dissociated S1 domain with ACE2 (H, I). Structural maps of collective dynamics are derived from fluctuations driven by the slowest three modes. The color gradient from blue to red indicates the decreasing structural stability (or increasing conformational mobility) of protein residues. The key functional residues E406, K417, N439, E484, and N501 are shown in spheres colored according to the level of mobility in the low frequency slow modes.

Figure 4

Functional dynamics of the SARS-CoV-2 S-RBD complex with REGN-COV2 antibody cocktail (REGN10933 and REGN10987). The mean square displacements in functional motions are averaged over the three lowest frequency modes. (A) The slow mode profile of the bound SARS-CoV-2 S-RBD structure in the complex with REGN-COV2 antibody combination (PDB id 6XDG). The essential mobility profiles are shown in maroon lines, and positions of E406, K417, N439, E484, and N501 are highlighted by filled green circles. (B) Structural map of the essential mobility profiles for the SARS-CoV-2 S-RBD complex with REGN-COV2 antibody cocktail. The color gradient from blue to red indicates the decreasing structural stability (or increasing conformational mobility) of protein residues. The key functional residues E406, K417, N439, E484, and N501 are annotated and shown in colored spheres according to the level of their respective mobility in the low frequency slow modes.

Consistent with these studies, the functional movements of RBDs can be determined by the main hinge centers located near F318, S591, F592, and V539 residues. These results highlighted the conserved nature of global hinges in the dissociated S1 domain and in the SARS-CoV-2 S complete trimer. Notably, the RBD residues targeted by novel mutational variants are located in the flexible moving regions of the unbound S1 domain. The distribution of hinge sites is altered in the SARS-CoV-2 RBD complex with ACE2, and strikingly the sites subjected to circulating variants and escape mutations often coincided with hinge clusters anchored by E406, K417, N439, and N501 residues (Figure D). E484 is located near F486 that is aligned with another local hinge position, while N501 together with T500 and Y505 can form a dominant hinge center in the complex. Hence, a group of functional residues that include N439, F486, T500, N501, and Y505 sites may form a network of regulatory centers coordinating global movements of the RBD and ACE2 molecules. These findings are in line with studies of conformational dynamics of the SARS-CoV-2 trimeric spike glycoprotein in complex with ACE2 revealing these positions may be involved in the regulation of continuous swing motions of ACE2-RBD relative to the SARS-CoV-2 S trimer.[33] The slow mode profiles obtained for complexes of the dissociated S1 domain with ACE2 further highlighted the role of these functional residues in collective motions (Figure E,F). The entire RBD region and structurally conserved C-terminal domain 1, CTD1 (residues 528–591) become largely immobilized in the S1 domain complex with ACE2 complex, while C-terminal domain 2, CTD2 region (residues 592–686) can undergo large movements (Figure F). By zooming on the RBD regions, one could see that local minima and corresponding hinge sites are almost precisely aligned with residues E406, K417, N501, and Y505. Hence, the RBD residues targeted by mutant variants may play an important role in coordinating the relative orientation and approach angle of the ACE2 receptor in the complex and, consequently, affect recognition and signal transmission in the functional complex. Structural maps of functional dynamics profiles illustrated these findings showing that the RBD regions that include functional positions E406, N439, and N501 are aligned with immobilized in slow motions hinge centers, while K417 and E484 residues could be less constrained during collective movements (Figure G–I). The central finding of this analysis is the unique role that positions targeted by novel variants (N439, N501) and escape mutations (E406) could play in concerted functional movements of the S-RBD and S1 domain when bound to ACE2. On the basis of the results, we argue that the functional role of these sites in controlling global motions and long-range interactions in ACE2 complexes could be an important reason for mediating escape from antibody binding while maintaining and enhancing binding with the host receptor.

To compare the slow mode profiles of the SARS-CoV-2 S1/RBD complexes with ACE2 and neutralizing antibodies, we performed ENM-based modeling of slow mode profiles in the S-RBD complex with REGN-COV2 cocktail of two antibodies REGN-10933 and REGN10987 (Figure ). This analysis showed that sites E406 and K417 corresponded to local hinge positions, while N439 and K444/G446 residues are now aligned with the dominant hinge center of the SARS-CoV-2 RBD complex with the REGN-COV2 cocktail (Figure A).

Importantly, antibody binding altered the dynamic role of residues E484 and N501 that become aligned with moving regions in the collective dynamics of the complex (Figure A). Structural mapping of the essential profiles further illustrated this point, showing that functionally immobilized in collective motions hinge centers are localized near E406, N439, and K444 sites, while E484 and N501 positions could undergo some movements in the complex (Figure B). In this context it is particularly interesting to compare our observations with functional studies showing that K417 and F486 are sites of escape from RERGN10933, while mutations in K444 and G446 escape neutralization by REGN10987 and E406 is a unique site susceptible to mutations escaping both antibodies.[70] In line with these experiments, we found that E406 and K444 positions may correspond to the antibody-specific unique hinge centers of collective motions that control relative orientation and rigid body movements of REGN10933 and REGN10987 molecules (Figure B). This is in some contrast to SARS-CoV-2 RBD complexes with ACE2 in which N439 and N501 form the major hinge center of functional dynamics. As a result, it is possible that mutations in K444 and E406 positions may perturb not only local interactions with antibody molecules but alter the global collective movements and long-range communication, which may be sufficient to trigger mutational escape from antibody binding. At the same time, these mutations could only moderately change the SARS-CoV-2 RBD local interactions with ACE2 without affecting the collective movements in the complex.

To summarize, this analysis suggested that mutational variants and escape mutations may preferentially target specific positions involved in regulation and coordination of functional dynamics motions and allosteric changes in the SARS-CoV-2 complexes with ACE2 and REGN-COV2 antibody cocktail.

Mutational Sensitivity Profiling of the SARS-CoV-2 RBD Binding Interfaces Reveals Energetic Plasticity in Sites Susceptible to Circulating Mutational Variants and Antibody Escaping Modifications

We first performed a systematic alanine scanning of the SARS-CoV RBD S protein residues (Figure A) and residues from the dissociated S1 domain in the complexes with the ACE2 host receptor (Figure B,C). Using the equilibrium ensembles obtained from simulation trajectories, we evaluated the average cumulative mutational effect of alanine substitutions on protein stability and binding affinity with the host receptor. The alanine scanning of the SARS-CoV-2 RBD residues highlighted a significant destabilization effect caused by mutations of G446, Y453, L455, F456, F486, Y489, Y495, T500, and Y505 residues (Figure A). These residues also corresponded to the binding free energy hotspots in the complexes formed by the dissociated S1 domain with ACE2 (Figure B,C). In particular, large destabilization effects were observed upon mutations of Y453, L455, F456, Y489, and F490 residues in the S1-ACE2 complexes. Notably, the largest destabilization changes were produced by alanine mutations of F456 and Y489 residues, displaying clear and pronounced peaks of the profile and pointing to these positions as key binding affinity hotspots in the S1-ACE2 complexes. (Figure B,C). Several key binding energy hotspot sites (Y453, Y489, and Y505) are conserved between SARS-CoV and SARS-CoV-2 proteins and are located in the central segment of the interface. A detailed analysis of the intermolecular contacts in the SARS-CoV-2 RBD and S1 complexes with ACE2 aided in understanding the binding energy preferences of RBD residues (Supporting Information, Tables S1–S3). This analysis is particularly instructive by considering contact distributions with two virus-binding hotspots on ACE2 formed by interacting residues K31 and E35 as well as K353 and D38.

Figure 5

Alanine scanning of the RBD residues in the SARS-CoV-2 S-RBD and S1 domain complexes with ACE2 and REGN-COV2 antibody cocktail. (A) The binding free energy changes upon alanine mutations for the RBD residues in the SARS-CoV-2 S-RBD complex with ACE2 (PDB id 6M0J). (B) The binding free energy changes upon alanine mutations for the S1-RBD residues in the S1 domain complex with ACE2 (PDB id 7A91). (C) The binding free energy changes upon alanine mutations for the dissociated S1 domain residues in the complex with ACE2 (PDB id 7A92). (D) The binding free energy changes upon alanine mutations for the SARS-CoV-2 S-RBD residues in the complex with the REGN-COV2 cocktail (PDB id 6XDG). The binding energy changes for the protein residues are shown in maroon bars. The binding interface residues are depicted in orange filled circles and functional residues K417, E484, and N501 targeted by mutational variants are highlighted in magenta filled circles.

In particular, Y489 residue makes numerous favorable contacts with multiple ACE2 residues K31, F28, Y873, L79, and T27 residues, while another hotspot position, F456, forms interactions with T27, D30, and K31 positions on ACE2 (Supporting Information, Tables S1–S3). L455 residue of the RBD makes favorable contacts with the key hotspots on ACE2 K31, H34, and D30 while Q493 contacts K31, H34, and E35 positions. The hydrophobic residue F486 forms multiple interactions with M82, L79, Y83, Q24, while another hydrophobic RBD site F490 interacts with K31 and Y473 with T27 (Supporting Information, Tables S1–S3). Overall, the alanine scanning highlighted the importance of the SARS-CoV-2 RBD interactions formed by the middle segment of the RBM interface (K417, Y453, L455, F456, Y489, and Q493) as mutations of these residues resulted in a significant loss of binding affinity (Figure A-C). These results are consistent with recent functional studies, indicating that mutations of the SARS-CoV-2 RBD residues (L455/Y442, F456/L443, F486/L472, F490/W476, Q493/N479) result in a significant reduction of their binding affinity with ACE2.[37]

The alanine scanning of the SARS-CoV-2 RBD residues in the complex with REGN-COV2 antibody cocktail revealed large destabilization effects and strong peaks at positions K444, V445, F456, F486, and Y489 (Supporting Information, Table S4, Figure D). The largest free energy changes exceeding 3 kcal/mol were observed for alanine modifications of F456 and F486 residues which is consistent with the prominent role these residues play in eliciting antibody-escaping mutations.[70] We also specifically highlighted binding free energy changes caused by alanine modifications in sites targeted by circulating and antibody-escaping variants E406, K417, N439, E484, and N501. The results showed that alanine substitutions in these positions induced only minor destabilization changes in the RBD-ACE2 and S1-ACE2 complexes, and these values were particularly small when mutations were introduced in E484 and N501 positions (Figure A–C). A slightly different pattern was seen in the RBD-REGN-COV2 complex, where alanine modifications in K417, N439 and E484 residues led to appreciable >1.0 kcal/mol binding free energy loss, while mutations in E406 and N501 positions produced only a small destabilization effect. These patterns indicated that antibody binding may induce changes in the binding interactions and distribution of the binding energy hotspots.

To further quantify the effects of the binding energy hotspots, we followed up with a complete mutational sensitivity analysis of these RBD positions (Figure ). The profiling showed that all mutations in Y453, F456, F486, and Y505 positions were highly destabilizing (Figure A–D), while, for L455 and F490 positions, the loss of binding affinity was only moderately destabilizing for the majority of substitutions (Figure E,F). The destabilization pattern observed for all modifications of F456 and F486 residues is consistent with the functional experiments[37] highlighting the importance of these positions in binding affinity. We also conducted mutational sensitivity scanning of K417, E484, and N501 residues that are targeted by circulating variants in the UK (B.1.1.7/501Y.V1), South Africa (501Y.V2), and Brazil (B1.1.28/501.V3) lineages[59−65] as well as profiling of E406 and N439 sites that are of considerable interest due newly emerging circulating variants and antibody-escaping mutations (Figure ). Of special interest was the analysis of the protein stability and binding free energy changes incurred by N501Y mutation that is prominently featured in the UK B.1.1.7 variant and mutations K417N, E484K that together with N501Y modifications are central to the increased transmission and infectivity effects seen in the South Africa (501Y.V2) and Brazil (B1.1.28/501.V3) lineages. In addition, we considered other important sites of newly circulating mutations N439 and the unique position E406 giving rise to unique antibody-escaping mutations.[70,81]

Figure 6

Mutational sensitivity analysis of binding free energy hotspots in the SARS-CoV-2 S-RBD complex with ACE2 (PDB id 6M0J). (A) Mutational sensitivity scanning of the Y453 residue. (B) Mutational sensitivity scanning of the F456 residue. (C) Mutational sensitivity scanning of the F486 residue. (D) Mutational sensitivity scanning of the Y505 residue. (E) Mutational sensitivity scanning of the L455 residue. (F) Mutational sensitivity scanning of the F490 residue. The protein stability changes are shown in maroon filled bars.

Figure 7

Mutational sensitivity analysis of functional RBD residues targeted by novel mutational variants and antibody-escaping mutations in the SARS-CoV-2 S-RBD complex with ACE2 (PDB id 6M0J). (A) Mutational sensitivity scanning of the E406 residue. (B) Mutational sensitivity scanning of the K417 residue. (C) Mutational sensitivity scanning of the N439 residue. (D) Mutational sensitivity scanning of the E484 residue. (E) Mutational sensitivity scanning of the N501 residue. The protein stability changes are shown in maroon filled bars.

Consistent with the deep mutational scanning experiments,[35,66] we found that E406, N439, and E484 sites are energetically adaptable and can effectively tolerate different mutations without incurring significant changes in protein stability and binding affinity (Figure A,C,D). Somewhat larger but still relatively tolerable were binding free energy changes induced by mutations in K417 and N501 positions (Figure B,E). K417 is a unique ACE2-interacting residue that forms favorable contacts with central residues of the ACE2 interface H34 and D30 (Supporting Information, Table S1–S3). However, an appreciable energetic plasticity could be seen in mutational sensitivity profiling of K417 residue (Figure B). Although K417 mutations to alanine or glycine produced fairly significant destabilization changes, K417N and K417D mutations led to only small perturbations (∼0.4–0.5 kcal/mol). Indeed, deep mutational scanning suggests that the K417N mutation has minimal impact on binding affinity with ACE2.[35]

The results also predicted the marginal improvement in the binding free energy mediated by E484K mutation, and only a very modest increase in the binding affinity upon K417N modification (Figure D). The experimental studies indicated that the E484K mutation may induce a moderate improvement in binding affinity and showed that other single mutations of E484 may only slightly compromise spike folding stability and binding affinity for ACE2.[35,70] According to our analysis, several hydrophobic substitutions in this position (E484I, E484 V, E484F, E484W, and E484P) may, in fact, lead to the moderately improved affinity, while other mutations appeared to produce only marginal destabilization (Figure D). These results indicated significant plasticity of this important RBD position that is relatively exposed and may favor hydrophobic residues in this position to improve both stability and binding. The mutational sensitivity profiling at the N501 position is consistent with deep mutational scanning experiments[35,66] reproducing the improvements in binding mediated by N501F and N501Y mutations (Figure E). Indeed, deep mutation scanning showed that N501F, N501T, and N501Y mutations may lead to moderate enhancement of binding with ACE2, while N501D is an affinity-decreasing mutation.[35] We observed only a small destabilization effect for N501T and more significant destabilization upon N501D and N501A/G mutations (Figure E). Importantly, these results supported the notion that N501Y mutational variant could be beneficial for ACE2 binding, while escaping neutralizing antibodies targeting the same region.

We also performed a mutational sensitivity analysis of the key functional positions in the SARS-CoV-2 RBD complex the REGN-COV2 (Figure A–D). The results revealed moderate changes upon mutations at E406, E484, N501, and N439 positions. Interestingly, although E406W escaped both individual REGN-COV2 antibodies,[70] our results indicate that this mutation would not drastically perturb the RBD region and significantly affect the binding interactions with the REGN-COV2 antibody cocktail (Figure A). Hence, the antibody-escaping effect of E406W substitution may not be trivially linked to the local interaction effects. In this context, given the results of functional dynamics analysis, it is tempting to argue that mutations at the E406 position may instead alter collective movements and compromise long-range allosteric couplings in the SARS-CoV-2 RBD. This may ultimately affect the neutralization activity of the REGN-COV2 antibody combination. At the same time, a wide range of modifications at K417, K444, and F486 sites resulted in significant destabilization changes and loss of the binding affinity (Figure E–G). These results are in excellent agreement with functional mapping of the SARS-CoV-2 RBD residues that affect the binding of the REGN-COV2 cocktail showing that F486 mutations are predominant for escaping neutralization by REGN10933 and mutations at K444 evade binding of REGN10987 antibody.[70]

Figure 8

Mutational sensitivity analysis of functional RBD residues in the SARS-CoV-2 S-RBD complex with REGN-COV2 antibody cocktail (REGN10933 and REGN10987). (A) Mutational sensitivity scanning of the E406 residue. (B) Mutational sensitivity scanning of the E484 residue. (C) Mutational sensitivity scanning of the N501 residue. (D) Mutational sensitivity scanning of the N439 residue. (E) Mutational sensitivity scanning of the K417 residue. (F) Mutational sensitivity scanning of the K444 residue. (G) Mutational sensitivity scanning of the F486 residue. The protein stability changes are shown in maroon filled bars.

Consistent with the functional analysis of the immune-selected mutational landscape in the S protein, we found that a wide spectrum of K444 modifications induced a significant loss of binding free energy (Figure E) including K444E and K444N mutations that showed a broad-range resistance against multiple antibodies.[69] The large destabilization changes caused by F486 mutations can be contrasted to fairly small changes incurred by E484 mutations, indicating that the E484 site is characterized by a sufficient level of structural plasticity and energetic adaptability to readily accommodate mutations in complexes with ACE2 and REGN-COV2 cocktail. These findings may explain why single-site mutations of these residues can only slightly change binding affinity for ACE2 and folding stability, while double-site mutations of proximal E484 and F486 can significantly weaken the fitness of the SARS-CoV-2 RBD region and binding.[66]

To summarize, our results pointed to several interesting trends. First, mutational sensitivity profiling of the conserved hydrophobic binding energy hotspots Y453, L455, F456, F486, and Y505 consistently yielded large destabilization changes affecting folding stability and binding to ACE2 receptor, making these positions unlikely candidates for antibody escaping mutations as even small modifications in these positions could have a severely detrimental effect on the spike activity. Second, we found that SARS-CoV-2 binding affinity could be strongly influenced by the virus-binding hotspot K31 and H34 in the middle of the interface through an extensive interaction network with K417, Y453, L455, F456, and Q493 residues. Finally, mutational analysis of K417, E484, and N501 positions implicated in new mutational strains and antibody-escaping changes showed that these residues correspond to important interacting centers with a significant degree of structural and energetic plasticity. Indeed, N501Y, E484K, and K417N mutations can result in an improved or only slightly decreased affinity with ACE2. These results suggest a hypothesis that antibody-escaping mutations target residues with sufficient plasticity and adaptability to preserve a sufficient spike activity while having a more detrimental effect on antibody recognition. These findings are particularly interesting in light of recent functional studies[66] showing that escape mutations target a subset of sites in the antibody-RBD interfaces corresponding to binding energy hotspots. Importantly, these experiments suggested that escape mutations are consistently those that have significant deleterious effects on antibody binding but little negative impact on ACE2 binding and RBD folds. On the basis of our findings, we argue that escape mutations constrained by the requirements for ACE2 binding and preservation of RBD stability may preferentially select structurally plastic and energetically adaptable allosteric centers at the key interfacial regions to compromise antibody recognition through modulation of global motions and allosteric interactions in the complex.

Perturbation Response Scanning Reveals Structurally Adaptable Allosteric Effector Hotspots in Sites Targeted by Global Circulating Mutations

Allosteric molecular events involve a complex interplay of thermodynamic and dynamic changes taken place on a long-time scale that are difficult to directly observe and simulate. Perturbation-based computational approaches based on linear response theory allow for sequential perturbation of the SARS-CoV-2 S protein residues by applying random forces to each residue while monitoring the protein response. By using this approach, we can examine and quantify long-range couplings between the sites of local perturbations and response.

Using the PRS method,[124−127] we quantified the allosteric effect of each residue in the SARS-CoV-2 complexes. The effector profiles estimate the propensities of a given residue to influence dynamic changes in other residues and can be applied to identify regulatory hotspots of allosteric interactions as the local maxima along the profile. First, we computed the residue-based effector response profiles for the SARS-CoV-2 RBD complex with ACE2 (Figure A) and the complexes formed by the dissociated S1 domain with ACE2 (Figure B,C). By comparing the PRS profiles in the ACE2 complexes with SARS-CoV-2 S1/RBD and REGN-COV2 antibody cocktail, we determined the distribution of regulatory allosteric centers and highlighted a potential role of sites targeted by global circulating mutations.

Figure 9

PRS effector profiles and distance fluctuation communication indexes for the SARS-CoV-2 S-RBD and S1 domain complexes with ACE2. (A) The PRS effector distribution profile for the SARS-CoV-2 S-RBD complex with ACE2 (PDB id 6M0J). (B) The PRS effector profile for the S1-RBD complex with ACE2 (PDB id 7A91). (C) The PRS effector profile for the dissociated S1 domain complex with ACE2 (PDB id 7A92). The positions of functional residues E406, K417, N439, E484, and N501 are indicated by filled green circles. (D) The distance fluctuation communication index profile for the SARS-CoV-2 S-RBD complex with ACE2 (PDB id 6M0J). (E) The distance fluctuation communication index profile for the S1-RBD complex with ACE2 (PDB id 7A91). (F) The distance fluctuation communication index profile for the dissociated S1 domain complex with ACE2 (PDB id 7A92). The positions of functional residues E406, K417, N439, E484, and N501 are indicated by filled orange circles.

Strikingly, the effector profile of the SARS-CoV-2 RBD complex with ACE2 featured two major peaks corresponding to residues E606 and T500/N501, indicating that E406 and N501 positions are aligned with the regulatory centers that may control allosteric communications in the complex (Figure A). Several other notable peaks corresponded to W353, K417, N439, and L452 residues. Hence, all known positions targeted by novel circulating mutational variants with the exception of E484 corresponded to the effector peaks and are involved in the coordination of allosteric communications in the SARS-CoV-2 RBD complex with ACE2. Moreover, the effector profiles indicated that these regulatory sites may function in a coordinated manner and maintain an allosteric cross-talk to control signal transmission “traffic” and long-range interactions in the RBD-ACE2 complex. Interestingly, several of these effector centers L452, N439, and N501 were among SARS-CoV-2 RBD residues whose mutations to the SARS-CoV RBD counterparts N439/R426, L452/K439, and N501/T487 enhanced the binding affinity.[37] The prominent role of these residues as regulatory effector centers becomes even more apparent in the S1 domain complexes with ACE2 (Figure B,C). It is evident that E406, K417, N439 and especially N501 positions corresponded to sharp peaks of the effector profile. This implies that these sites may be collectively responsible for the coordination of long-range communication in the system.

The central result of this analysis is that circulating and escape mutations appeared to target residues corresponding to structurally and energetically adaptable regulatory control points that can tolerate individual mutations and often enhance ACE2 binding, while at the same time allowing for coordinated modulation of allosteric communications. We suggest that allosteric signaling in the SARS-CoV-2 RBD complex with ACE2 is adaptable where a mutation of a regulatory control point can be functionally compensated through energetic rebalancing of structurally plastic allosteric hotspots.

Using a protein mechanics-based approach,[159] we also employed distance fluctuations analysis of the conformational ensembles to further probe allosteric communication preferences of the RBD residues in the SARS-CoV-2 RBD and S1 complexes with ACE2. The residue-based distance fluctuation communication indexes measure the energy cost of the dynamic residue deformations and could serve as a robust metric for the assessment of allosteric propensities of protein residues. In this model, dynamically correlated residues whose effective distances fluctuate with low or moderate intensity are expected to communicate with higher efficiency than the residues that experience large fluctuations. Notably, structurally stable and densely interconnected residues as well as moderately flexible residues that serve as a source or sink of allosteric signals could feature a high value of these indexes.

The distance fluctuation profile of the SARS-CoV-2 RBD and S1 domain complexes with ACE2 showed a small but important redistribution of major peaks, pointing to sites E406, W436, N439, N501, and Y505 (Figure D). Notably, this group of residues is featured prominently among peaks of the profile when the entire S1 domain monomer forms complex with ACE2 (Figure E,F). We also noticed that the overall shape and distribution of the peaks are similar between the PRS effector profiles and distance fluctuation communication index profiles. Notably, E406, N439, and N501 sites were featured as recurring peaks in both distributions, strengthening the proposed notion that positions targeted by the emerging mutational variants can cooperate and play a central role in the regulation of long-range couplings and allosteric communications in the complexes with the ACE2 host receptor. Hence, the distance fluctuation profiling and analysis of communication indexes provide important supporting evidence to the PRS modeling, suggesting that structurally stable positions and potential allosteric hotspot residues only partially overlap, and allosteric hubs may exhibit a certain degree of structural plasticity and energetic adaptability to enable a balance between binding and signaling function.

To understand a potential role of the E484 residue, it is instructive to analyze the PRS sensor profile (Supporting Information, Figure S2). A comparison between sensor profiles obtained for the unbound and bound forms of the SARS-CoV-2 RBD showed that the E484 position is aligned with the dominant peak of the sensor profile in the unbound form (Supporting Information, Figure S2A). Interestingly, in the complex with ACE2, this site also corresponded to a major sensor peak at the binding interface (Supporting Information, Figure S2B). Structural mapping of sensor profiles in the unbound and bound RBD forms illustrated these observations, pointing to the role of E484 residue as a major receiver site of allosteric signaling in the RBD-ACE2 complex (Supporting Information, Figure S2C,D). Hence, the PRS analysis of the RBD-ACE2 and S1-ACE2 complexes demonstrated that, while E406, K417, N439, and N501 are aligned with dominant effector positions representing the source and regulatory points of allosteric signaling, E484 corresponded to a major sensor/receiver site that may absorb signal information. Collectively, these sites may represent key nodes of the allosteric interaction network in the functional ACE2-bound complexes and determine the robustness and efficiency of signal transmission.

We also computed the PRS effector profiles for the SARS-RBD complex with REGN-COV2 antibody combination (Figure ). The effector profile revealed some redistribution of peaks, featuring V401/E406, N439, K444/G446, and G496 positions as major effector centers (Figure A). At the same time, residues E484/F486 and N501 were aligned with the local sensor peaks (Figure B). These results could provide a feasible rationale for the critical role of K444 and F486 positions in escaping antibody combinations. Indeed, K444 is a central epitope residue for REG10987, while F486 residue is fundamental for the recognition of the REG10933 antibody.[70] Our findings also indicated that E406 and K444 are the dominant effector centers in the RBD complex with REGN-COV2 (Figure A) and may be functionally important not only for binding affinity but also for mediating signaling and long-range communications in the complex. In the context of the perturbation-based PRS model, this implies that single mutations at these positions could affect collective movements and allosteric couplings between many residues in the system and potentially compromise the functional activity of the REGN-COV2 cocktail. Interestingly, other positions targeted by antibody-escaping mutants E484 and F486 are the major sensor sites (Figure B). On the basis of these observations, we suggest that allosteric control of the RBD-REGN COV2 complex is provided through a cross-talk between major effector sites (E406, K444) and receiver sites (E484 and F486).

Figure 10

PRS profiles for the SARS-CoV-2 S-RBD and S1 domain complexes with REGN-COV2 antibody cocktail. (A) The PRS effector profile for the SARS-CoV-2 S-RBD complex with REGN-COV2 antibodies (REGN10933 and REGN10987) (PDB id 6XDG). (B) The PRS sensor profile for the SARS-CoV-2 S-RBD complex with REGN-COV2 (PDB id 6XDG). The positions of functional residues E406, K417, N439, K444, E484, and N501 are indicated by filled green circles.

To summarize, perturbation-based modeling of the SARS-CoV-2 RBD complexes suggested that functional residues targeted by global circulating variants and antibody-escaping mutants could form a network of structurally adaptable allosteric hotspots that collectively coordinate allosteric interactions in the system. These results bear some significance and support the latest illuminating study suggesting a model functional plasticity and evolutionary adaptation of allosteric regulation.[160] This function-centric model of allostery revealed remarkable functional plasticity of allosteric switches allowing modulate and restore regulatory activity through mutational combinations or ligand interactions. Our results similarly suggested that functional plasticity and cross-talk of allosteric control points in the SARS-CoV-2 RBD region can allow for differential modulation of recognition and long-range communication with ACE2 and antibodies.

Network Modeling Reveals that Sites Targeted by Circulating Mutations are Mediating Anchors of the Intermolecular Communities with ACE2 and REGN-COV2 Antibodies

Structure-based network approaches have offered a powerful conceptual formalism and an array of robust computational tools for describing allosteric interactions.[161] Network-centric models of protein structure and dynamics used in the current study provide a complementary perspective to the physics-based analysis of conformational dynamics landscapes and allow for quantitative analysis of allosteric changes in which conformational landscapes of protein systems can be remodeled by various perturbations such as mutations, ligand binding, or interactions with other proteins.[162] Another emerging concept central to computational models of allostery is the identification of regulatory control points that mediate long-range communications and allosteric pathways between conformational switch centers.[161,162] Our current understanding of communication pathways in proteins is based on the ensemble-based statistical model that often invokes community-based methods for modeling ensembles of intermodular pathways and analysis of the modular organization of protein structure networks.[135−137] Using community decomposition, the residue interaction networks can be divided into local interaction modules in which residues are densely interconnected and highly correlated during simulations, while different communities are connected through long-range couplings. A community-based model of allosteric communications is based on the notion that groups of residues that form local interacting communities are correlated and switch their conformational states cooperatively. In this model, allosteric communications can be transmitted through a chain of stable local modules connected via intercommunity bridges.[135−137] In the present study, we leveraged the results of community decomposition and used edge betweenness in the global interaction network as a proxy for modeling of allosteric communication pathways and assessment of antibody-induced modulation of allosteric interactions in the SARS-CoV-2 S proteins. Using this computational framework, we set out to explore allosteric communications and identify functional centers in the spike protein that propagate cooperative structural changes through modular community structure.

Using this network-centric description of residue interactions, we compared the organization of stable local communities in the SARS-CoV-2 RBD complexes with ACE2 (Figure ) and REGN-COV2 cocktail (Figure ). In the RBD-ACE2 and S1-ACE2 complexes, we found a deeply interconnected community organization (Figure ) where stable modules in the RBD core are tightly linked with the interfacial clusters. A number of stable intramolecular communities in the RBD core contribute to the stability of the RBD regions. Some of these communities are formed by hydrophobic core residues including F342–V511–F374–W436-F347-R509, W353–F400–Y423, as well as E406-Q409-I418 and N439–443-P499 centered around functional residues E406 and N439 (Figure ). Of particular interest and importance was a more detailed comparative analysis of the intermolecular communities. In the RBD-ACE2 and S1-ACE2 complexes, these modules are integrated around key anchor residues K417, F456, Y489, N501, and Y505. The largest and most stable community in which each node is strongly linked with each other is centered on N501 (K353-D38-Y41-Q498–N501-Y505) and engaged ACE2 hotspots K353, D38, and Y41 (Figure A). This interfacial module anchored by N501 allows for persistent interactions by N501, Y505, and Q498 RBD residues. Moreover, this community may be instrumental for signal transmission between RBD and ACE2 molecules, highlighting the important role of the N501 position. By introducing the N501Y mutation, we rebuilt the residue interaction network and performed community decomposition, which revealed the preservation of this major community. This observation supported our energetic analysis indicating structural plasticity and stability of the key intermolecular communities in the RBD-ACE2 complex.

Figure 11

Structural maps of the intra- and intermolecular local communities formed in the SARS-CoV-2 S-RBD complex with ACE2. (Central panel) Molecular topography of major communities in the SARS-CoV-2 S-RBD complex with ACE2. SARS-CoV-2 RBD is shown in red ribbons and ACE2 in cyan ribbons. The RBD residues involved in local communities are shown in red spheres, and ACE2 residues involved in the intermolecular communities only are shown in cyan spheres. (Left panel A) A close-up of the major intermolecular community that is anchored by N501 and includes Y505, Q498 residues of RBD, and K353, and Y41 hotspot residues of ACE2 are highlighted in red and cyan spheres, respectively, and annotated. (Right panel B). A close-up of another dominant interfacial community anchored by F456 is shown. The community includes K417, F456, Y489, and E484 residues of RBD and K31/D30 hotspot residues of ACE2. The community residues are shown in spheres and annotated. Notably, these two major intermolecular communities that mediate communications and stability of the interface include several key functional sites, K417, F456, E484, and N501, targeted by mutational variants and antibody-escaping mutations.

Figure 12

Structural maps of local communities in the SARS-CoV-2 S-RBD complex with REGN-COV2 cocktail of antibodies. (A) A general community overview of the SARS-CoV-2 S-RBD complex with REGN-COV2. The RBD is in red ribbons, and interfacial communities are shown in red spheres. (B) A close-up of the intermolecular communities formed by the S-RBD. The communities at the interface with REGN10933 include F486(RBD)-E484(RBD)-Y489(RBD)-R100(REGN10933)-Y50(REGN10933)-Y59(REGN10933)-L94(REGN10933) and K417(RBD)-E406(RBD)-Q409(RBD)-Y32(REGN10933)-T102(REGN10933). The interfacial community with another antibody in the complex is W99(REGN10987)-W47(REGN10987)-G446(RBD)-V445(RBD)-K444(RBD)-Y59(REGN10987) that is anchored by K444, V445, and G446 residues connected with Y59 and W99 of REGN10987.

Another group of interfacial communities is anchored by K417 and F456 positions that couple modules Y489–F456–K31 and D30–K417-F456 (Figure B). Interestingly, these intermolecular communities are directly coupled through K417 with the intramolecular module I402–I418-E406-Y495-Q409 centered on the E406 residue. Hence, the community organization revealed strong interconnectivity between key functional sites E406, K417, F456, and N501 that integrate the residue interaction network and enable allosteric couplings between RBD and ACE2 molecules. The important revelation of this analysis is that only a fraction of the RBD residues anchor the intermolecular community organization and mediate long-range interactions in the RBD-ACE2 complex. Furthermore, the binding free energy hotspots are not necessarily involved in community-mediating functions. Instead, a group of structurally plastic allosteric centers such as E406, K417, F456, and N501 plays key roles in integrating local communities into a robust and adaptable global network that can mediate signal transmission and communication between SARS-CoV-2 RBD and ACE2.

The REGN-COV2 antibody binding induced a partial and yet significant reorganization of the interface communities (Figure A). We found that key RBD residues F486, Y489, and K417 anchor major intermolecular communities with REGN10933 including R100–Y50–F486-W47-L94, Y32-T102-K417, and Y33-T52-Y489 (Figure B). Interestingly, K417 and another site of escape mutations E406 are interconnected in the local RBD community E406-Q409–K417-I418. E484 is involved in the formation of the intermolecular contacts with T52, Y53, T57, Y59 residues of the heavy chain of REGN10933 and could bridge several interfacial communities anchored by F486 and Y489 residues. Among major interfacial communities formed by the RBD with REGN10987 is the W99–W47–V445–K444-Y59 module, in which V445 and K444 play a key role in mediating intermolecular communication. Another notable community is anchored by T500, which engages N439, P499, K444 residues of the RBD and W99 of the heavy chain of REGN10987 (Figure B). These findings are consistent with the experimental deep mutational screening showing that K444 is a critical RBD residue for REGN10987, while F486 and Y489 sites are essential in inducing the activity of REGN10933. According to these experiments, mutations F486I and Y489H escape REGN10933, and E484A/F486I combinations evade REGN10933. Interestingly, it was also established that mutations at F486 escaped neutralization only by REGN10933, whereas mutations of N439 and K444 escaped neutralization only by REGN10987.

E406 site is at the center of the largest intramolecular community in the SARS-CoV-2 RBD (R403–Y453–V350-I418–N422–Y423-Y495-F497-E406-Q409) that connects the intermolecular interfaces with the RBD core (Supporting Information, Figure S3). This site could anchor large communities in the core with the interfacial regions facing both antibodies. On the basis of these observations, we suggested that unique escape mutations in the E406 position may be largely determined by its allosteric mediating role in interconnecting functional regions of the RBD. Given the fact that E406 is only involved in several contacts with T28 of REGN10933, the strong mutational escape effect may be mainly driven by long-range allosteric effects and attributed to the strategic position of this residue in the global network. Noticeably, E406 is closely connected with the Y453 residue, where another REGN10933 specific escape mutant Y453F was detected.[70] Another key member of this community is F497 that effectively bridges local intermolecular modules interacting with REGN10933 and REGN10987 antibodies (Supporting Information, Figure S3). Recent studies indicated that SARS-CoV-2 neighbor residues G496 and F497 are critical for the RBD–ACE2 interaction and that F497 may play an important role in enhancing the RBD–ACE2 interaction for SARS-Cov-2 RBD. Our network analysis quantifies the structural hypothesis offered in the experimental study according to which E406 escape mutation may affect recognition by REGN10987 through cascading effect onto adjacent structural elements across the RBD and propagating changes through aromatic residues Y453, Y495, and F497.[70] Importantly, network analysis revealed that these hydrophobic residues belong to the tightly packed stable community anchored by the E406 residue. Owing to the modular interconnectivity where each of these residues is connected with every community neighbor, it is likely that E406W mutation may simultaneously perturb multiple contacts and alter couplings between these residues, thus adversely affecting the fidelity of allosteric communication with REGN10987. Hence, although this residue makes no persistent contact with either of the interacting antibodies, its position in the largest community anchored by E406 could be important for the integrity of the network organization in the complex with antibodies.

The network analysis showed that not only are these residues involved in favorable contacts with these antibodies but also they define key regulatory nodes that mediate stability and connectivity of the intermolecular communities and may be responsible for the control of signal transmission between SARS-CoV-2 RBD and REGN10933/REGN10987 antibodies. In agreement with functional data, the network community analysis singled out N439, K444, F486, and Y489 sites as allosteric network hubs where mutations would result in weakening of the entire interface and compromise the efficiency of allosteric interactions.

Spotting Allosteric Functional Hotspots in the SARS-CoV-2 S Proteins through Integrative Computational Approaches

This integrative computational investigation combined molecular simulations and functional dynamics analysis with mutational energetic profiling of the SARS-CoV-2 S protein binding and network-based community modeling to delineate specific allosteric signatures of functionally important residues that are subjected to novel circulating variants. The proposed computational framework is based on the notion that dynamic characterization of functional allosteric states and atomistic reconstruction of conformational energy landscapes enables to unravel mechanisms of binding-induced modulation of SARS-CoV-2 spike protein activity. Although all-atom MD simulations with the explicit and detailed characterization of the glycosylation shield can provide a truly comprehensive and rigorous assessment of the conformational landscape of the SARS-CoV-2 S proteins, such simulations remain to be extremely computationally demanding due to the size of the system. The adopted approach in our study represented an opportunistic strategy in which the robustness and simplicity of high-resolution CG simulations to adequately probe conformational space were combined with atomistic reconstruction of the simulation trajectories. The global dynamic and network features of the spike protein landscapes are assumed to be largely determined by the underlying topology of the SARS-CoV-2 S structures and can be captured through a proposed combination of CABS-CG dynamics and subsequent atomistic refinement of trajectories with structurally resolved glycans. This simulation approach enabled efficient conformational sampling and topological characterization of the energy landscapes, hinge regions, and regulatory switch centers that may control binding-induced modulation of conformational changes in the SARS-CoV-2 S structures.

Using the conformational energy landscapes of the SARS-CoV-2 S proteins as the conceptual core of our strategy, we conducted the ensemble-based mutational scanning of spike protein residues to assess local binding propensities and affinities with specific antibodies as well as perturbation-response scanning in which global allosteric propensities of spike residues are profiled to determine functional allosteric hotspots of the SARS-CoV-2 S regulation. These approaches mimic deep mutational scanning experiments by averaging over equilibrium samples and provide valuable complementary information about binding energy hotspots that drive local changes in the binding affinity and global allosteric hotspots that are involved in the coordination of allosteric communication and regulatory control of SARS-CoV-2 S binding with antibodies Through mutational and perturbation-based scanning we independently profiled binding and allosteric propensities of the SARS-CoV-2 S proteins, showing structural and energetic plasticity of the important allosteric hotspots that corresponded to the sites of the emerging global circulating mutations. The important revelation of this dual profiling of the SARS-CoV-2 spike protein is the fact that allosteric regulatory centers may not necessarily correspond to the binding free energy hotspots, suggesting that antibody-escaping mutational variants would likely to emerge in dynamically adaptable regions important for allostery as this “strategy” my enable virus to preserve the stability of the spike fold and retain favorable interactions with the host receptor while compromising protein response and allosteric signaling to antibodies. By combining the dynamics and ensemble-based scanning of the SARS-CoV- S binding with modeling of the residue interaction networks and community analysis, our integrated strategy linked allosteric properties of functional spike regions with their mediating role in global interaction networks. Using this computational framework, we found that efficient allosteric communications in SARS-CoV-2 spike proteins could be controlled by allosteric functional centers that bridge local stable communities and serve as key “stepping stones” of allosteric communications and long-range interactions in the SARS-CoV-2 S complexes with antibodies. Collectively, the application of synergistic computational approaches allowed for quantitative characterization of molecular mechanisms of SARS-CoV-2 S proteins, suggesting that functionally important spike sites subjected to global circulating mutations may emerge due to evolutionary adaptation in structurally plastic and energetically tolerable positions that play a unique and critical role as global mediators of allosteric interactions and communications of SARS-CoV-2 S protein binding. We argue that evolutionary adaptation in the SARS-CoV-2 S proteins may operate on functionally critical allosteric sites, such as sites of global circulating mutations, by using only minimal perturbations while achieving global protein response through allosteric signaling propagated through conserved interaction network scaffolds.

Conclusions

Comparative modeling of the SARS-CoV-2 S complexes with ACE2 and REGN-COV2 antibody combination revealed several important trends and characterized the unique allosteric-centric signatures of functional spike residues. Conformational dynamics and analysis of global motions in the SARS-CoV-2 S complexes demonstrated that the RBD residues targeted by novel mutational variants may be structurally adaptable and play a central role in the regulation of collective movements and long-range couplings. Through mutational sensitivity analysis of the SARS-CoV-2 RBD residues, we accurately reproduced the binding affinity changes for N501Y, E484K, and K417N mutations and found that these important interacting centers are characterized by a significant degree of structural and energetic plasticity. Using PRS analysis and community modeling of the SARS-CoV-2 S complexes, we characterized signatures of functional residues implicated in novel variants and demonstrated that these sites form a network of energetically adaptable regulatory centers modulating long-range communication of the SARS-CoV-2 S-RBD regions with ACE2 and antibodies. The results of this study demonstrated that the SARS-CoV-2 S protein may function as a versatile and functionally adaptable allosteric machine that exploits the plasticity of allosteric regulatory centers to fine-tune the response to antibody binding without compromising the activity of the spike protein. By examining molecular mechanisms of the SARS-CoV-2 S protein binding with antibodies through the lens of dynamic allosteric analysis and network modeling, we identified the functional centers in the spike protein that can be exploited for manipulating conformational landscapes of the SARS-CoV-2 S proteins, design of allosteric modulators, and potential engineering of ligand-specific regulatory responses ,which may be useful for tailoring new therapeutic interventions to combat potential antibody-escaping resistance.