key: cord-0740829-umn2k6fy
authors: Verkhivker, Gennady
title: Conformational Flexibility and Local Frustration in the Functional States of the SARS-CoV-2 Spike B.1.1.7 and B.1.351 Variants: Mutation-Induced Allosteric Modulation Mechanism of Functional Dynamics and Protein Stability
date: 2022-01-31
journal: Int J Mol Sci
DOI: 10.3390/ijms23031646
sha: 0b17eb73510a03f423eb3092fce52400ad122e58
doc_id: 740829
cord_uid: umn2k6fy

Structural and functional studies of the SARS-CoV-2 spike proteins have recently determined distinct functional states of the B.1.1.7 and B.1.351 spike variants, providing a molecular framework for understanding the mechanisms that link the effect of mutations with the enhanced virus infectivity and transmissibility. A detailed dynamic and energetic analysis of these variants was undertaken in the present work to quantify the effects of different mutations on functional conformational changes and stability of the SARS-CoV-2 spike protein. We employed the efficient and accurate coarse-grained (CG) simulations of multiple functional states of the D614G mutant, B.1.1.7 and B.1.351 spike variants to characterize conformational dynamics of the SARS-CoV-2 spike proteins and identify dynamic signatures of the functional regions that regulate transitions between the closed and open forms. By combining molecular simulations with full atomistic reconstruction of the trajectories and the ensemble-based mutational frustration analysis, we characterized how the intrinsic flexibility of specific spike regions can control functional conformational changes required for binding with the host-cell receptor. Using the residue-based mutational scanning of protein stability, we determined protein stability hotspots and identified potential energetic drivers favoring the receptor-accessible open spike states for the B.1.1.7 and B.1.351 spike variants. The results suggested that modulation of the energetic frustration at the inter-protomer interfaces can serve as a mechanism for allosteric couplings between mutational sites and the inter-protomer hinges of functional motions. The proposed mechanism of mutation-induced energetic frustration may result in greater adaptability and the emergence of multiple conformational states in the open form. This study suggested that SARS-CoV-2 B.1.1.7 and B.1.351 variants may leverage the intrinsic plasticity of functional regions in the spike protein for mutation-induced modulation of protein dynamics and allosteric regulation to control binding with the host cell receptor.

The coronavirus SARS-CoV-2 infection is transmitted when the viral spike (S) glycoprotein binds to the angiotensin-converting enzyme 2 (ACE2) host receptor, leading to the entry of S protein into host cells and membrane fusion [1, 2] . The full-length SARS-CoV-2 S protein consists of amino (N)-terminal S1 subunit and carboxyl (C)-terminal S2 subunit, where 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 CTD1 and CTD2 ( Figure S1 ). Structural and biochemical studies established that the mechanism of virus infection may involve conformational transitions between distinct functional forms of the SARS-CoV-2 S protein in which the RBDs continuously switch between "down" and "up" positions [3] [4] [5] [6] [7] [8] [9] [10] . The cryo-EM structure of the SARS-CoV-2 S trimer revealed a spectrum of closed states that included a structurally rigid closed form and more dynamic closed states preceding a transition to the fully open S conformation [5] . Protein engineering and structural studies also showed that specific disulfide bonds and proline mutations can modulate stability of the SARS-CoV-2 S trimer [6] and lead to the thermodynamic shifts between the closed and open forms [7] [8] [9] . Dynamic structural changes that accompany SARS-CoV-2 S binding with the ACE2 host receptor were described in cryo-EM experiments, showing a cascade of conformational transitions from a compact closed form, weakened after furin cleavage, to the partially open states, and, subsequently, the ACE2-bound open form thus priming the S protein for fusion [10] . The cryo-EM studies and biophysical studies of SARS-CoV-2 S trimer examined conformational flexibility and distribution of S trimers in situ on the virion surface, showing that spontaneous conformational changes and population shifts between different functional states are maintained in different biological environments and can coexist with receptor-induced structural adaptation, reflecting the intrinsic conformational landscapes of the SARS-CoV-2 S proteins [11, 12] .

The emergence of variants of concern (VOCs) with the enhanced transmissibility and infectivity profile, including D614G variant [13] [14] [15] [16] , B.1.1.7 (alpha) [17] [18] [19] [20] , B.1.351 (beta) [21, 22] , B.1.1.28/P.1 (gamma) [23] and B.1.1.427/B.1.429 (epsilon) variants [24, 25] have attracted enormous attention in the scientific community and a considerable variety of the proposed mechanisms explaining functional observations from structural and biochemical perspectives. S-B.1.1.7 variant of the SARS-CoV-2 has 9 of the 17 mutations (∆69-70 deletion, ∆144 deletion, N501Y, A570D, D614G, P681H, T716I, S982A, D1118H) in the S protein, featuring N501Y mutation, which can increase binding affinity with ACE2 while eliciting immune escape from RBD-targeting antibodies [17] [18] [19] [20] . S-B.1.351 variant is characterized by 21 mutations with 9 mutations (L18F, D80A, D215G, R246I, K417N, E484K, N501Y, D614G, and A701V) in the spike protein, of which three mutations (K417N, E484K, and N501Y) are located in the RBD and increase the binding affinity for the ACE receptors and induce significant immune escape [21, 22] . The cryo-EM studies provided atomistic characterization of the closed and open SARS-CoV-2 S trimers for these variants, suggesting subtle structural and functional effects of mutations that can modulate dynamics and stability of the closed and open forms, increase binding to the human receptor ACE2, and confer resistance to neutralizing antibodies [26] [27] [28] [29] . The cryo-EM structures of the S-B.1.1.7 variant in the apo form identified four distinct conformational substates, all describing one or two RBDs in the open conformation without populating the fully closed state [30] .

Importantly, these variants of the SARS-CoV-2 share D614G mutation, which is also individually linked with the enhanced infectivity, following the evidence of the mutation enrichment via epidemiological surveillance [13] [14] [15] [16] . The initial structural studies showed that the D614G mutation can act by shifting the population of the SARS-CoV-2 S trimer from the closed form (53% of the equilibrium) in the native spike protein to a widely-open topology of the "up" protomers in the D614G mutant, with 36% of the population adopting a single open protomer, 39% with two open protomers and 20% with all three protomers in the open conformation [31] . These SARS-CoV-2 S ectodomain structures included two mutational regions in which the S1/S2 cleavage site residues 682-RRAR-685) are replaced by GSAS motif and two stabilizing prolines inserted at positions K986 and V987 to prevent conformational change to the post-fusion state [31] . The S-GSAS construct showed similar structural, antigenic, and stability behavior as the S-GSAS/PP construct, which included the K986P and V987P mutations. The cryo-EM structures of the S-GSAS/D614G mutant revealed the increased population of the 1-RBD-up open form as compared to the closed state in the S-GSAS/D614 structure [32] . The electron microscopy analysis confirmed the higher percentage (84%) of the 1-up RBD conformation in the S-G614 protein [33] . Functional studies showed that the S-G614 mutant exhibited a greater infectivity than the S-D614 protein, which was attributed to the greater stability of the S-G614 mutant and led to the reduced S1 subdomain shedding [34] . The increased stability of the D614G mutant was inferred from the recent cryo-EM structures of a full-length unmodified S-G614 trimer, which can reversibly adopt an all-down closed state and 1 RBD-up open conformation [35] .

A significant heterogeneity and plasticity of the conformational landscapes for the SARS-CoV-2 S proteins was recently unveiled in the large-scale structural investigation that reported 38 cryo-EM structures of the B.1.351, B.1.1.28/P.1, B.1.617.2 (Delta), and B.1.617.1 (Kappa) variants, featuring a broad spectrum of the RBD-up conformations intrinsically present in the apo form [36] . In particular, S-B.1.351 variant displayed a highly populated fully open state with all three RBDs in the up conformation in addition to significant fractions of the partially open states (1 RBD-up or 2 RBD-up) [36] . Importantly, among all the examined variants only S-B.1.351 exhibited a fully open 3 RBD-up conformation that resembles the structure of the fully open S-G614 variant. It was suggested that the increased conformational diversity and a greater variability of the open states can allow SARS-CoV-2 S variants to adapt their responses and allow for escape immunity from common vaccines and different classes of monoclonal antibodies. Biophysical studies using microscale thermophoresis [37] and surface plasmon resonance (SPR) experiments [38] established that SARS-CoV-2 B.1.1.7 and B.1.351 spike variants bind human ACE2 with increased affinity. These studies suggested that the N501Y and S477N mutations enhance transmission primarily by enhancing binding, and that the K417N/T mutations facilitate immune escape, and the E484K mutation enhances binding and immune escape [38] . By engineering SARS-CoV-2 S proteins with different combinations of RBD variant mutations, cryo-EM structural analysis and biochemical assays systematically dissected the mutational effect on the ACE2 affinity and antibody evasion [39] . This study showed that N501Y, E484K, and L452R mutations can modulate, either directly or allosterically, the increased ACE2 binding affinity, while E484K, L452R, and K417N/T mutations tend to primarily compromise productive antibody binding and induce immune escape. Collectively, the rapidly growing body of structural and biochemical studies broadened our understanding of the conformational diversity and adaptability of the open states for the SARS-CoV-2 S variants, suggesting that coronaviruses have a broad potential to exploit synergistic effects of mutations to readily incur structural changes and modulate binding affinity with the ACE2 receptor and antibodies without a substantial loss of function [40] . The detection of common mutational changes such as D614G, E484K, N501Y and K417N that are shared among major circulating variants B.1.1.7, B.1.351, and B.1.1.28/P.1 indicated that these positions can be particularly critical for modulation of the SARS-CoV-2 S protein responses and may induce immunity escape from vaccines and different classes of monoclonal antibodies [41] [42] [43] [44] . Together, these studies unveiled a complex balance among mutational variants in which individual modifications often act cooperatively to enhance or reduce protein stability, modulate binding to the ACE2 receptor, and confer immunity resistance to neutralizing antibodies [45] .

Here, we employed computer simulations of the S-G614, S-B.1.1.7 and S-B.1.351 trimers in the closed and open states to systematically characterize the effect of mutational variants on conformational dynamics, and protein stability, as well as allosteric changes driven by collective motions and controlled by the regulatory hotspot centers. Computer simulations and protein modeling played an important role in shaping up our understanding of the dynamics and function of SARS-CoV-2 glycoproteins [46] [47] [48] [49] [50] [51] [52] [53] [54] [55] [56] . The latest extensive molecular dynamics (MD) simulations and free-energy landscape mapping studies of the SARS-CoV-2 S proteins and mutants detailed conformational changes and diversity of ensembles, further supporting the notion of enhanced functional and structural plasticity of S proteins [57] [58] [59] [60] [61] [62] [63] . Our recent studies offered evidence that the SARS-CoV-2 spike protein can function as an allosterically regulated machine that exploits plasticity of allosteric hotspots to fine-tune response to antibody binding [64] [65] [66] [67] [68] [69] . These studies showed that examining allosteric behavior of the SARS-CoV-2 S proteins may be useful to uncover functional mechanisms and rationalize the growing body of diverse experimental data. In particular, using an atomistic-based model of signal transmission in the SARS-CoV-2 S proteins, we determined that the D614G mutation can exert its effect through allosterically induced changes on stability and communications in the residue interaction networks [70, 71] . Previous computational studies also identified residues that exhibit long-distance couplings with the RBD opening that included sites harboring functional mutations D614G and A570D shared across SARS-CoV-2 S protein variants [59] . The free-energy landscapes of the S protein derived from MD simulations, together with nudged elastic pathway optimization mapping of the RBD opening, revealed a specific transient allosteric pocket at the hinge region, which is located near D614 position, influences RBD dynamics [60] . The energy analysis of the S-D614 and S-G614 proteins showed that local interactions near D614 position may be energetically frustrated and become stabilized in the S-G614 mutant through strengthening the inter-protomer association between S1 and S2 regions [63] .

In the current study, we utilized molecular simulations of the studied SARS-CoV-2 S variants to perform the ensemble-based local frustration analysis of the S protein in different functional forms, and systematically characterize mutational adaptability of the spike residues. Conformational dynamics analysis and characterization of collective motions in the SARS-CoV-2 S variants are combined to identify the regulatory hotspots that control allosteric structural changes between the closed and opened states. We also perform a systematic mutational scanning of the SARS-CoV-2 S residues that are modified in the S protein variants in order to determine protein stability hotspots and identify potential energetic drivers underlying thermodynamic preferences of the B. Figure 1A ,C,E). The S1 subunit, which includes NTD (residues 14-306), RBD (residues 331-528), CTD1 (residues 528-591), and CTD2 (residues 592-685), showed a greater mobility level as compared to a more rigid S2 subunit (residues 686-1162). Common to these structures, 630-loop (residues 617 to 644) and FPPR segment (residues 828 to 853) are ordered in the RBD-down form, while the 630-loop becomes disordered in the 1-up protomer of the open state ( Figure S1 ). (Figure 2A ). Of particular interest are the large thermal fluctuations of the NTD regions and the receptor binding motif (RBM) (residues 424-494) where one of the key mutational sites E484 showed a considerable variability. In some contrast, the conformational mobility of these positions in the S-G614 and S-B.1.1.7 closed states may be more moderate. The enhanced mobility of this functional RBD region is important to enable spontaneous transitions between the closed and open forms required for engagement with the ACE2 receptor. It is worth noting that the conformational dynamics profiles also revealed the increased thermal "breezing" that is spread across both S1 and S2 subunits in the S-B.1.1.7 and S-B.1.351 closed states (Figure 2A ), providing support for the notion that these mutational variants may perturb and destabilize the closed form, thereby shifting the equilibrium to a more dynamic open form [34] [35] [36] (Table S1 ). The convergence of CG simulations was assessed by using the average Spearman's correlation coefficient (r s ) for residue mobility between different simulations based on the respective RMSF profiles [73] ( Table S1 ). The reported average correlation coefficients for residue mobility between different CG simulations correspond to the mean values obtained from pairwise comparisons of 100 independent trajectories for each studied system. The average statistics showed a relatively high average Spearman's correlation for S-G614 trajectories (r s = 0.835 for the closed state and r s = 0.804 for the open state) (Table S1) (Table S1 ). The global similarity of conformations generated in CG simulations was obtained by computing the average root mean square deviation (RMSD) between the snapshots in the trajectories (Table S1) . The position of the RBM motif is indicated by an arrow. The S1 subunit (residues 14-685) and S2 subunit (residues 686-1163). The S1 domains include NTD (residues 14-306), RBD (residues 331-528), CTD1 (residues 528-591), CTD2 (residues 592-685). The Figure 3 ). Of particular notice is a significant softening of the closed S-G614 and S-B.1.351 conformations. In these closed states, we detected the increased fluctuations in both S1 and S2 subunits ( Figure 3A ,E). For the S-B.1.1.7 closed trimer, the fluctuations of the S2 regions were smaller and a more significant flexibility was observed for the NTD regions ( Figure 3C ). These results suggest that the increased preferences of the S protein variants towards 1 RBD-up open conformation may be partly caused by the increased mobility and destabilization of the closed forms. One of the key observations of the conformational dynamics analysis is a steady stabilization trend in the open forms of the S protein variants ( Figure 3B ,D,F). In particular, we found that S-G614 and S-B.1.1.7 conformations displayed a broad stabilization in both S1 and S2 subunits but pointed to plasticity at the inter-protomer interfaces particularly near D614G site. Noticeably, the S2 subunit regions displayed only small fluctuations and are considerably more stable than more dynamic S1 domains ( Figure 3B [36] . Accordingly, our findings agree with the experimental studies indicating that the greater variability of the open states can allow S-B.1.351 variant to readily adapt a spectrum of specific RBD-up conformations, thus providing a plausible mechanism for enhancing binding with the ACE2 receptor while modulating the immunity response.

By highlighting the sites of mutational variants, we gained some additional insight into mutation-induced modulation of conformational mobility in the S proteins. In particular, it is instructive to examine dynamic changes in the inter-protomer interface regions near G614 ( Figure 3 ). This key mutational site shared by all studied variants experienced appreciable fluctuations in the 3 RBD-down states of the S-G614 ( Figure 3A To characterize the ensemble-based structural environment of the mutational sites in the S proteins, we evaluated relative solvent accessibility (RSA) of protein residues. A residue-specific local RSA measure is defined as the ratio of the observed solventaccessible surface area for a residue to the expected unfolded state value for that amino acid type [74, 75] . The RSA values can be used as proxy for predicting intrinsic residue flexibility and the extent of conformational changes that may be induced by mutational variants. According to this model, residues are considered to be completely solvent-exposed if the ratio value exceeds 50%, to be buried if the ratio is less than 20%, and to be designated as partially buried for residues with 20% < RSA < 50%. We monitored the RSA values for the S protein residues averaged over simulation trajectories (Figure 4 ). For the S-G614 ensemble of closed conformations, only E484 residue is solvent-exposed, while K417 and N501 are buried ( Figure 4A ). The functional position A570 is totally buried and for the mutational variant position G614 the RSA value of~19.8% is right at the border line between buried (RSA < 20 %) and partially buried sites (20% < RSA < 50%), reflecting some level of plasticity at the inter-protomer interface near the G614 l site ( Figure 4A ). In the S-G614 ensemble of open states, residues K417 (RSA = 27.7%) and N501 (RSA = 34.6%) become partially buried, while E484 (RSA = 65.4%) is solvent-exposed and the inter-protomer A570 position remains totally buried (RSA = 2.2%) ( Figure 4B) . Notably, the mutational site G614 becomes only partially buried (RSA = 35%) to the proposed scale, which reflected conformational flexibility of the RBD-up protomers and widening of the inter-protomer interfaces near the mutational site G614 ( Figure 4B ). Strikingly, in both closed and open states of the S-B.1.1.7 variant, A570D and S982A positions remain completely buried, featuring RSA < 10% ( Figure 4C,D) . At the same time, the functional RBD sites K417, E484, and N501Y become solvent-exposed, particularly E484 and N501Y positions with RSA > 70% ( Figure 4C,D) . Interestingly, while the G614 position is solvent-exposed in the closed state of S-B. (panels A,B) . The S-B.1.1.7 mutational sites are N501Y, A570D, D614G, P681H, T716I, S982A, and D1118H (panels C,D). The S-B.1.351 mutational sites are L18F, D80A, D215G, R246I, K417N, E484K, N501Y, D614G, and A701V (panels E,F). The structures are in sphere-based representation and are rendered using UCSF ChimeraX [72] with the rigidity-to-flexibility sliding scale colored from blue (highly rigid) to red (highly flexible). The S1 and S2 subunits are indicated by arrows. The NTD(protomer B)-RBD(protomer A) interfaces and inter-protomer interfaces near G614 site are also indicated by arrows, and annotated. These revelations may be particularly important given that the low-frequency soft modes are cooperative and robust as they are determined by the protein fold topology. Furthermore, allosterically-driven conformational changes typically evolve along slow modes of motions intrinsically accessible to folded structures [77, 78] . As a result, allosteric responses in proteins can be efficiently triggered when mutations modulate protein propensities for energetically favorable movement along the slow modes. The observed correspondence between mutational sites and regulatory hinge points driving collective motions indicated that the virus may target these functional positions to modulate global dynamic response and promote allosteric transitions to the open state required for efficient binding with the host-cell receptor. E,F) . The structures are in sphere-based representation and are rendered using UCSF ChimeraX [72] with the rigidity-to-flexibility sliding scale colored blue (highly rigid) to red (highly flexible). The inter-protomer interfaces and hinge near G614 site are indicated by arrows, and annotated.

We employed the conformational ensembles of the S-G614, S-B.1.1.7 and S-B.1.351 trimers to estimate the average local frustration profiles of the S residues and quantify the relationship between structural plasticity and mutational frustration. This analysis is based on scanning of the conformational ensembles by local frustratometer, which computes the local frustration index based on the contribution of a residue or residue pairs to the energy in a given conformation as compared to what it would contribute to decoy conformations [79] [80] [81] [82] [83] . The distribution of local mutational frustration showed a low relative density of highly frustrated population for residues targeted by mutations in the S trimer mutants (Figure 7 ). This is particularly apparent for the S-G614 conformations where K417, E484, N501, A570 and G614 sites featured a very small density of highly frustrated positions ( Figure 7A ). In some contrast, sites N501 and T716 in the S-B.1.1.7 conformations may be highly frustrated, while N501 and E484 may also feature an appreciable relative density of highly frustrated positions ( Figure 7D ). The role of frustration in the RBD regions is particularly interesting in light of recent evidence that the disordered or highly flexible regions can be critically important for mediating allostery and binding to multiple protein partners [84] [85] [86] [87] [88] . Notably, the population of highly frustrated conformations is relatively minor for all S conformations (Figure 7A Figure 7A ). The relative density of neutral frustration is pronounced for all mutational sites, especially for the N501 and G614 positions, showing no appreciable differences between the closed and open states ( Figure 7B) . Notably, the hinge site A570 showed an appreciable density of minimal frustration, which is consistent with the regulatory role of this position in coordinating global allosteric changes. In general, sites of mutational variants located in the flexible RBD regions showed a propensity for neutral mutational frustration.

There are several revealing differences in the local frustration distributions of the S-B.1.1.7 ( Figure 7D Combined with the collective dynamics analysis, these observations suggest that moderately-to-minimally frustrated A570 and A570D hinge positions may allow for conformational switches at the inter-protomer interfaces that drive functional transitions. This is consistent with the experimental studies showing that A570D and D614G hinge positions in the S-B.1.1.7 protein form a molecular switch that incurs allosteric structural changes that can enhance the RBD motions [30] . In this S-B.1.1.7-specific switch, A570D forms the inter-protomer salt bridges with K854 and K964 that compensate for the loss of the salt bridge between D614 and K854 due to the D614G mutation [30] .

In the S-B.1.351 conformations, we observed that sites targeted by mutations in the exposed RBD and NTD regions also featured neutral local frustration, suggesting a moderate mutational adaptability in these positions ( Figure 7G conformations could allow for some mutational adaptability in the hinge position, which retains its regulatory role in the S variants, as confirmed by the structural and functional studies [30] . The generally prevailing pattern of frustration neutrality for sites targeted by mutations across all studied variants and significant contributions of high and neutral frustration density in the RBD sites E484 and N501 are important findings of this analysis, since local frustration is often an important mediator of allosteric changes. The relatively high frustration in these positions may allow for a discrete set of configurations involving local motions of the frustrated residues. Combined with modulation of the inter-protomer hinge regions by A570D and D614G positions, the local frustration in the RBD mutational sites could drive dynamical transitions between closed and open states, accompanied by local adjustments of the RBD residues.

The relationship between frustration in proteins and their function has been explored in illuminating studies by Wolynes [89] . The protein folding landscape theory established that while minimally frustrated interactions may have evolved in proteins to enable strongly funneled landscapes, a number of functional regions and specific protein sites could have been selected to be frustrated to allow for modulation of global motions and binding adaptability with various interaction partners. In this context, our results suggest that the detected neutral frustration patterns could allow for suboptimal inter-protomer interactions and adaptable substates to efficiently control access and binding affinity with the host-cell receptor ACE2.

Using the equilibrium conformational ensembles, we performed mutational scanning of the spike protein residues and mutational sensitivity analysis of the S-B.1.1.7 ( Figure 8 ) and S-B.1.351 mutant variants in the closed and open states (Figure 9 ). The protein stability ∆∆G changes were computed by averaging the results of computations over 1000 samples obtained from simulation trajectories. We first examined the pattern of free-energy changes for the S-B.1.1.7 closed and open states (Figure 8) . Interestingly, mutational scanning of A570D and D614G produced relatively few destabilization changes in the open state ( Figure 8B ). Hence, structural stability constraints tend to protect the key positions that are responsible for modulation of the inter-protomer interactions. This is consistent with the experimental structural studies [26] showing that A570D and D614G favor the acquisition of the open S-B.1.1.7 state. In some contrast, mutations in the functional RBD positions (K417, E484) and modifications of mutational variant sites T716I, S982A and D1118H resulted in minor energy changes in both closed ( Figure 8A ) and open S-B.1.1.7 conformations ( Figure 8B ). These results are also consistent with the experiments [26] showing that N501Y, T716I and D1118H mutations induce only minimal local conformational changes without affecting stability of the S protein. Interestingly, S982A substitution in the S-B.1.1.7 conformations abolished hydrogen bonding between central helices of the S2 domain and the CTD1 region [26] . As a result, mutations of this position in the S-B.1.1.7 states resulted in moderate and often stabilizing changes, indicating a considerable mutational tolerance and conformational plasticity in this region (Figure 8 ).

The mutational cartography analysis revealed that A570D and D614G mutational positions in the S-B.1.1.7 conformations are the most sensitive to modifications that result in more significant destabilizing changes ( Figure 8) . Importantly, the effect of mutations in these positions on the protein stability is greater in the open state ( Figure 8B ). The S-B.1.1.7 open state featured a salt bridge involving interactions of A570D with K854 of the other protomer [70] . In the closed S-B.1.1.7 states A570D can form the inter-protomer interactions with K964 and N856 that together comprise an important hinge cluster. Functional dynamics analysis of slow modes confirmed the role of A570D as a potential regulatory switch that controls RBD movements. Although modifications of A570D are generally destabilizing, the range of free-energy changes associated with this position suggested a moderate level of residual energetic frustration and suboptimal interactions ( Figure 9B ) and allow for more optimal interactions of the RBD-up protomer with the host receptor. In both states of S-B.1.351 variant, we observed appreciable destabilization changes induced by modifications in the D614G position ( Figure 9 ). This is consistent with the increased density of minimal frustration for this position in the open S-B.1.351 state. Mutational site A701V is located in the surface-exposed region of S2 and caused minimal structural changes [26] . Consistent with these experimental observations, modifications in this position have a minor effect on the protein stability. Importantly, the results show a significant difference between S-B.1.1.7 and S-B.1.351 structures in the mutational profiles for K417/E484/N501 positions (Figures 8 and 9 ). We found that these K417N, E484K and N501Y sites can afford greater mutational tolerance and conformational plasticity in the S-B.1.351 conformations. This suggests that the level of mutational and conformational plasticity can progressively increase from S-G614 to S-B.1.1.7 and S-B.1.351 variants.

The results also suggest an interesting interplay between mutation-induced protein stability and local frustration patterns. It is interesting to examine these relationships in the context of allosteric regulation models centered on the energetic frustration that could emerge at the inter-domain interfaces [79] [80] [81] [82] [83] . In these models, activation or repression functions may be encoded in the conformational ensemble and reveal through frustrationbased allosteric regulation of the inter-domain interactions. The energetic scanning and local frustration analyses revealed the inherent mutational plasticity of the sites targeted by circulating variants. Mutational positions that are involved in hinge clusters are characterized by some degree of energetic frustration, thereby allowing for allosteric couplings and modulation of the RBD motions leading to the greater diversity of RBD-exposed conformations. These results suggested that S protein variants may uniquely exploit the intrinsic conformational and mutational plasticity of the S proteins, which is broadly distributed and is characteristic of not only RBD regions but also present in the inter-protomer and inter-domain regions. We argue that S-B.1.1.7 and S-B.1.351 variants may leverage this plasticity to adopt a mechanism of frustration-based allosteric modulation at the interprotomer interfaces to differentially control binding with the host cell receptor ACE2 and interacting proteins.

The All structures were obtained from the Protein Data Bank [90, 91] . During structure preparation stage, protein residues in the crystal structures were inspected for missing residues and protons. Hydrogen atoms and missing residues were initially added and assigned according to the WHATIF program web interface [92, 93] . The structures were further pre-processed through the Protein Preparation Wizard (Schrödinger, LLC, New York, NY, USA) and included the check of bond order, assignment and adjustment of ionization states, formation of disulphide bonds, removal of crystallographic water molecules and co-factors, capping of the termini, assignment of partial charges, and addition of possible missing atoms and side chains that were not assigned in the initial processing with the WHATIF program. The missing loops in the studied cryo-EM structures of the SARS-CoV-2 S protein were reconstructed and optimized using template-based loop prediction approaches ModLoop [94] , ArchPRED server [95] and further confirmed by FALC (Fragment Assembly and Loop Closure) program [96] . The side chain rotamers were refined and optimized by SCWRL4 tool [97] . The conformational ensembles were also subjected to all-atom reconstruction using PULCHRA method [98] and CG2AA tool [99] to produce atomistic models of simulation trajectories. The protein structures were then optimized using atomic-level energy minimization with composite physics and knowledge-based force fields, as implemented in the 3Drefine method [100] . The atomistic structures from simulation trajectories were further 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 22 most occupied N-glycans [101, 102] in each as determined in the cryo-EM structures of the SARS-CoV-2 spike S trimer in the closed state (K986P/V987P,) (pdb id 6VXX) and open state (pdb id 6VYB), and the cryo-EM structure SARS-CoV-2 spike trimer (K986P/V987P) in the open state (pdb id 6VSB).

CABS-flex approach was employed in CG simulations that efficiently combine a highresolution interaction model and efficient search protocol shown to recapitulate all-atom MD simulation trajectories and dynamic profiles of large biomolecules on a long timescale [103] [104] [105] [106] [107] [108] . In this model, the amino acid residues are represented by Cα, Cβ, the center of mass of side chains and another pseudo atom placed in the center of the Cα-Cα pseudo-bond. The CABS-flex approach implemented as a Python 2.7 object-oriented standalone package [107] was employed in this study. 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 [103] [104] [105] . CABS-flex protocol, in which long series of random local Monte Carlo changes are performed, can accurately recapitulate all-atom long-time dynamics [73] . The default settings were used, 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. The CABS-flex default distance restraints moderately penalize the position of restrained residues if their distance differs from the distance in the original cryo-EM structure, becoming more than 1 Å. In these settings loop regions are fully unrestrained.

A total of 100 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. The convergence of CG-CABS simulations was ensured by performing multiple long-series Monte Carlo runs followed by statistical analysis of the RMSF profiles obtained in each of 100 long trajectories. In this analysis, the average Spearman's correlation coefficients were computed between the RMSF profiles for individual trajectories [73] . The reported average correlation coefficients represent the mean values obtained from pairwise comparisons of 100 independent trajectories for each of the studied structures of the S-G614, S-B.1.1.7 and S-B.1.351 variants (Table S1) . MODELLER-based reconstruction of simulation trajectories to all-atom representation provided by the CABS-flex package was employed to produce atomistic models of the equilibrium ensembles for studied systems [107] . Principal component analysis (PCA) of reconstructed atomistic trajectories was derived from CABS-CG simulations using the CARMA package [109] .

A systematic mutational scanning and sensitivity analysis of the functional RBD residues (K417, E484, N501), S-B.1.1.7 mutational sites (N501Y, A570D, P681H, T716I, S982A, and D1118H) and S-B.1.351 mutational sites (L18F, D80A, D215G, K417N, E484K, N501Y, D614G, and A701V).

Each position was systematically mutated using all 19 substitutions and corresponding protein stability changes were computed. The BeAtMuSiC approach was employed, which 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 [110] [111] [112] . If a free-energy change occurs between a mutant and the wild-type (WT) proteins ∆∆G = ∆G (MT)-∆G (WT) > 0, the mutation is destabilizing, while when ∆∆G < 0 the respective mutation is stabilizing. We computed the ensemble-averaged free energy changes using 1,000 equilibrium samples from simulation trajectories of the S-G614, S-B. 

In this work, we combined molecular simulations and collective dynamics analysis with the ensemble-based frustration analysis to characterize conformational plasticity and functional adaptability We found that mutations in the S variants may promote movements of the RBD regions, which may increase in the S-B.1.1.7 and especially S-B.1.351 conformations. The local frustration analysis shows a prevailing pattern of frustration neutrality for sites targeted by mutations across all studied variant and dynamic contributions of high and neutral frustration in the RBD sites E484 and N501 positions. A strong preference of the mutational sites towards neutrally frustrated environment may allow for suboptimal inter-protomer interactions and regulatory control of SARS-CoV-2 S binding with the ACE2 host cell receptor. In agreement with structural studies, the ensemble-based energetic scanning of protein stability for mutated positions in the S protein variants revealed that level of conformational and mutational plasticity can progressively increase from S-G614 to the S-B.1.1.7 and S-B.1.351 variants. We found that K417N, E484K and N501Y sites are relatively tolerant to modifications in the S-B.1.351 closed and open conformations. The results also suggest an interplay between mutation-induced protein stability, local frustration and allosteric modulation of the S protein dynamics. By combining conformational dynamics and local frustration analysis with the ensemble-based mutational scanning, our analysis demonstrated that A570D and D614G mutational sites emerge as key inter-protomer hinges that control functional movements and allosteric conformational changes in the S-B.1.1.7 variant ( Figure 6C,D) . The results suggest that mutational variants target key residues of the spike protein that coordinate an allosteric cross-talk between rigid hinge sites (A570D, D614G) and more flexible RBD mutational sites (K417N, E484K, and N501Y ). We argue that using a frustration-driven allosteric mechanism, mutational variants can impose regulatory control over functional movements and conformational adaptability of the RBD regions, thus modulating the accessibility and binding affinity with the ACE2 receptor. The Acknowledgments: This work was supported by the Kay Family Foundation Grant A20-0032. The author thanks Schmid College of Science and Technology, Chapman University for providing computing resources at the Keck Center for Science and Engineering.

The authors declare that the research was conducted in the absence of any commercial or financial relationship that could be construed as a potential conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results. 

Characterization of the receptor-binding domain (RBD) of 2019 novel coronavirus: Implication for development of RBD protein as a viral attachment inhibitor and vaccine

Structural and functional basis of SARS-CoV-2 entry by using human ACE2

Structure, Function, and Antigenicity of the SARS-CoV-2 Spike Glycoprotein

Cryo-EM structure of the 2019-nCoV spike in the prefusion conformation

Distinct conformational states of SARS-CoV-2 spike protein

Structure-based design of prefusion-stabilized SARS-CoV-2 spikes

Controlling the SARS-CoV-2 spike glycoprotein conformation

Structure-guided covalent stabilization of coronavirus spike glycoprotein trimers in the closed conformation

A thermostable, closed SARS-CoV-2 spike protein trimer

Receptor binding and priming of the spike protein of SARS-CoV-2 for membrane fusion

In situ structural analysis of SARS-CoV-2 spike reveals flexibility mediated by three hinges

Real-time conformational dynamics of SARS-CoV-2 spikes on virus particles

Tracking changes in SARS-CoV-2 Spike: Evidence that D614G increases infectivity of the COVID-19 virus

Spike mutation D614G alters SARS-CoV-2 fitness

SARS-CoV-2 D614G variant exhibits efficient replication ex vivo and transmission in vivo

Functional importance of the D614G mutation in the SARS-CoV-2 spike protein

First detection of SARS-CoV-2 spike protein N501 mutation in Italy in August

Increased mortality in community-tested cases of SARS-CoV-2 lineage B.1.1.7

Neutralization of SARS-CoV-2 lineage B.1.1.7 pseudovirus by BNT162b2 vaccine-elicited human sera

Detection of a SARS-CoV-2 variant of concern in South Africa

Sixteen novel lineages of SARS-CoV-2 in South Africa

Rapid Increase in SARS-CoV-2 P.1 Lineage Leading to Codominance with B.1.1.7 Lineage

Rapid characterization of spike variants via mammalian cell surface display

Broad sarbecovirus neutralization by a human monoclonal antibody

Structural basis for enhanced infectivity and immune evasion of SARS-CoV-2 variants

Effect of natural mutations of SARS-CoV-2 on spike structure, conformation, and antigenicity

SARS-CoV-2 immune evasion by the B.1.427/B.1.429 variant of concern

Structural and functional ramifications of antigenic drift in recent SARS-CoV-2 variants

Effect of SARS-CoV-2 B.1.1.7 mutations on spike protein structure and function

Structural and functional analysis of the D614G SARS-CoV-2 spike protein variant

D614G Mutation Alters SARS-CoV-2 Spike Conformation and Enhances Protease Cleavage at the S1/S2 Junction

D614G Spike Mutation Increases SARS CoV-2 Susceptibility to Neutralization

SARS-CoV-2 spike-protein D614G mutation increases virion spike density and infectivity

Structural impact on SARS-CoV-2 spike protein by D614G substitution

Structure-activity relationships of B.1.617 and other SARS-CoV-2 spike variants

SARS-CoV-2 B.1.1.7 and B.1.351 spike variants bind human ACE2 with increased affinity

Effects of common mutations in the SARS-CoV-2 Spike RBD and its ligand, the human ACE2 receptor on binding affinity and kinetics

Structural analysis of receptor binding domain mutations in SARS-CoV-2 variants of concern that modulate ACE2 and antibody binding

The emerging plasticity of SARS-CoV-2

Cross-Reactive Neutralizing Antibody Responses Elicited by SARS-CoV-2 501Y.V2 (B.1.351)

Sensitivity of infectious SARS-CoV-2 B.1.1.7 and B.1.351 variants to neutralizing antibodies

Genomic characterization of a novel SARS-CoV-2 lineage from Rio de Janeiro, Brazil

In vivo monoclonal antibody efficacy against SARS-CoV-2 variant strains

Probing Affinity, Avidity, Anticooperativity, and Competition in Antibody and Receptor Binding to the SARS-CoV-2 Spike by Single Particle Mass Analyses

Conformational transition of SARS-CoV-2 spike glycoprotein between its closed and open states

Enhanced receptor binding of SARS-CoV-2 through networks of hydrogen-bonding and hydrophobic interactions

Developing a Fully Glycosylated Full-Length SARS-CoV-2 Spike Protein Model in a Viral Membrane

Beyond Shielding: The Roles of Glycans in the SARS-CoV-2 Spike Protein

A multiscale coarse-grained model of the SARS-CoV-2 virion

Computational epitope map of SARS-CoV-2 spike protein

The SARS-CoV-2 exerts a distinctive strategy for interacting with the ACE2 human receptor

Critical Sequence Hotspots for Binding of Novel Coronavirus to Angiotensin Converter Enzyme as Evaluated by Molecular Simulations

Critical Interactions between the SARS-CoV-2 Spike Glycoprotein and the Human ACE2 Receptor

Structural Impact of Mutation D614G in SARS-CoV-2 Spike Protein: Enhanced Infectivity and Therapeutic Opportunity

Enhanced binding of the N501Y-mutated SARS-CoV-2 spike protein to the human ACE2 receptor: Insights from molecular dynamics simulations

SARS-CoV-2 simulations go exascale to predict dramatic spike opening and cryptic pockets across the proteome

A rigorous framework for detecting SARS-CoV-2 spike protein mutational ensemble from genomic and structural features

Distant residues modulate conformational opening in SARS-CoV-2 spike protein

Free Energy Landscapes from SARS-CoV-2 Spike Glycoprotein Simulations Suggest that RBD Opening can be Modulated via Interactions in an Allosteric Pocket

The SARS-CoV-2 Spike variant D614G favors an open conformational state

Modelling conformational state dynamics and its role on infection for SARS-CoV-2 Spike protein variants

D614G substitution at the hinge region enhances the stability of trimeric SARS-CoV-2 spike protein

Coevolution, dynamics and allostery conspire in shaping cooperative binding and signal transmission of the SARS-CoV-2 spike protein with human angiotensin-converting enzyme 2

Molecular simulations and network modeling reveal an allosteric signaling in the SARS-CoV-2 spike proteins

Dynamic Network Modeling of Allosteric Interactions and Communication Pathways in the SARS-CoV-2 Spike Trimer Mutants: Differential Modulation of Conformational Landscapes and Signal Transmission via Cascades of Regulatory Switches

Integrated Biophysical Modeling of the SARS-CoV-2 Spike Protein Binding and Allosteric Interactions with Antibodies

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

Dynamic Profiling of Binding and Allosteric Propensities of the SARS-CoV-2 Spike Protein with Different Classes of Antibodies: Mutational and Perturbation-Based Scanning Reveals the Allosteric Duality of Functionally Adaptable Hotspots

Computational analysis of protein stability and allosteric interaction networks in distinct conformational forms of the SARS-CoV-2 spike D614G mutant: Reconciling functional mechanisms through allosteric model of spike regulation

Landscape-Based Mutational Sensitivity Cartography and Network Community Analysis of the SARS-CoV-2 Spike Protein Structures: Quantifying Functional Effects of the Circulating D614G Variant

Structure visualization for researchers, educators, and developers

Consistent view of protein fluctuations from all-atom molecular dynamics and coarse-grained dynamics with knowledge-based force-field

Relative solvent accessible surface area predicts protein conformational changes upon binding

Exact and efficient analytical calculation of the accessible surface areas and their gradients for macromolecules

Conservation and prediction of solvent accessibility in protein families

Adaptability of protein structures to enable functional interactions and evolutionary implications

Intrinsic dynamics is evolutionarily optimized to enable allosteric behavior

Protein frustratometer: A tool to localize energetic frustration in protein molecules

Protein Frustratometer 2: A tool to localize energetic frustration in protein molecules, now with electrostatics

Frustration in Fuzzy Protein Complexes Leads to Interaction Versatility

Fuzziness and Frustration in the Energy Landscape of Protein Folding, Function, and Assembly

Localization of Energetic Frustration in Proteins

Structural and energetic basis of allostery

The ensemble nature of allostery

Ensemble allosteric model: Energetic frustration within the intrinsically disordered glucocorticoid receptor

Genetically tunable frustration controls allostery in an intrinsically disordered transcription factor

Functional plasticity and evolutionary adaptation of allosteric regulation

Evolution, energy landscapes and the paradoxes of protein folding

The Protein Data Bank

The RCSB protein data bank: Integrative view of protein, gene and 3D structural information

Positioning hydrogen atoms by optimizing hydrogen-bond networks in protein structures

WIWS: A protein structure bioinformatics web service collection

Automated modeling of loops in protein structures

ArchPRED: A template based loop structure prediction server

The FALC-Loop web server for protein loop modeling

Improved prediction of protein side-chain conformations with SCWRL4

Fast procedure for reconstruction of full-atom protein models from reduced representations

CG2AA: Backmapping protein coarse-grained structures

An interactive web server for efficient protein structure refinement

Vulnerabilities in coronavirus glycan shields despite extensive glycosylation

Site-specific glycan analysis of the SARS-CoV-2 spike

Characterization of protein-folding pathways by reduced-space modeling

Coarse-grained protein models and their applications

Modeling of protein structural flexibility and large-scale dynamics: Coarse-grained simulations and elastic network models

Modeling of disordered protein structures using monte carlo simulations and knowledge-based statistical force fields

CABS-flex standalone: A simulation environment for fast modeling of protein flexibility

Protocols for fast simulations of protein structure flexibility using CABS-Flex and SURPASS

A fully automated task-oriented interface for the analysis of molecular dynamics trajectories

BeAtMuSiC: Prediction of changes in protein-protein binding affinity on mutations

A new generation of statistical potentials for proteins

Fast and accurate predictions of protein stability changes upon mutations using statistical potentials and neural networks: PoPMuSiC-2.0