Molecular Simulations and Network Modeling Reveal an Allosteric Signaling in the SARS-CoV-2 Spike Proteins.

The development of computational strategies for the quantitative characterization of the functional mechanisms of SARS-CoV-2 spike proteins is of paramount importance in efforts to accelerate the discovery of novel therapeutic agents and vaccines combating the COVID-19 pandemic. Structural and biophysical studies have recently characterized the conformational landscapes of the SARS-CoV-2 spike glycoproteins in the prefusion form, revealing a spectrum of stable and more dynamic states. By employing molecular simulations and network modeling approaches, this study systematically examined functional dynamics and identified the regulatory centers of allosteric interactions for distinct functional states of the wild-type and mutant variants of the SARS-CoV-2 prefusion spike trimer. This study presents evidence that the SARS-CoV-2 spike protein can function as an allosteric regulatory engine that fluctuates between dynamically distinct functional states. Perturbation-based modeling of the interaction networks revealed a key role of the cross-talk between the effector hotspots in the receptor binding domain and the fusion peptide proximal region of the SARS-CoV-2 spike protein. The results have shown that the allosteric hotspots of the interaction networks in the SARS-CoV-2 spike protein can control the dynamic switching between functional conformational states that are associated with virus entry to the host receptor. This study offers a useful and novel perspective on the underlying mechanisms of the SARS-CoV-2 spike protein through the lens of allosteric signaling as a regulatory apparatus of virus transmission that could open up opportunities for targeted allosteric drug discovery against SARS-CoV-2 proteins and contribute to the rapid response to the current and potential future pandemic scenarios.

Introduction

The coronavirus SARS-CoV-2 is associated with acute respiratory distress syndrome[1,2] and is similar to the severe acute respiratory syndrome (SARS) and Middle East respiratory syndrome (MERS) viruses.[3] The genomic sequences of SARS-CoV-2 showed a high level of sequence similarity between the SARS-CoV-2, SARS, and MERS proteins involved in the replication cycle.[4−6] SARS-CoV-2 has four main structural proteins: spike (S) glycoprotein, small envelope (E) glycoprotein, membrane (M) glycoprotein, and nucleocapsid (N) protein, along with several accessory proteins.[7−9] Recent studies have identified that SARS-CoV-2 uses angiotensin-converting enzyme 2 (ACE2)[10−12] and the cellular protease transmembrane protease serine 2 (TMPRSS2) as cell entry receptors.[10] SARS spike proteins are trimeric fusion proteins with two main domains, an amino (N)-terminal S1 subunit and a carboxyl (C)-terminal S2 subunit. SARS-CoV-2 binds to the ACE2 receptor on the surface of the host cell using binding of the S1 region of the virus spike (S) protein followed by the fusion of the viral and cellular membranes mediated by the S2 subunit of the S protein.[13] The interaction between S and ACE2 results in fusion of the viral and cell membranes, releasing the viral genome and allowing for viral replication and host cell takeover. The crystal structure of the receptor-binding domain (RBD) of SARS-CoV-2 S protein bound to human ACE2 showed a common binding mode shared with SARS-CoV-S protein, and yet the SARS-CoV-2-RBD binding affinity is four times stronger.[14,15] Binding to the ACE2 receptor is a critical initial step for SARS-CoV to enter into target cells, and the increased infectiousness of SARS-CoV-2 is believed to arise from subtle structural and dynamic alterations in the RBD regions of the SARS-CoV2-S protein, which can increase the binding affinity for ACE2 as compared with SARS-CoV-S. Strikingly, however, the full-length SARS-CoV and SARS-CoV-2 S proteins bind ACE2 with comparable affinities, which was attributed to the less-exposed and semioptimal orientation of the SARS-CoV-2-RBD regions to interactions with ACE2 within the full-length protein.[16] The high-resolution crystal structures of the RBD of the spike protein for SARS-CoV[17−19] and SARS-CoV-2 bound to the cell receptor ACE2[20,21] have identified residues and two major hotspots in the RBD that may be important for the stronger binding of SARS-CoV-2 with ACE2 binding. Despite the sequence and structure similarities between these complexes, a number of residues in the highly variable receptor binding motif (RBM) region are different between SARS-CoV-2 RBD and SARS-CoV RBD.[17−21] It was argued that although the binding affinity changes are important, these differences may not be sufficient to adequately explain the unusual transmissibility of the SARS-CoV-2 virus to the host receptor protein.[20,21]

The full-length SARS-CoV-2 spike (S) protein includes the N-terminal domain (NTD), RBD, C-terminal domain 1 (CTD1), C-terminal domain 2 (CTD2), S1/S2 cleavage site, S2′ cleavage site, fusion peptide (FP), fusion peptide proximal region (FPPR), heptad repeat 1 (HR1), central helix region (CH), connector domain (CD), heptad repeat 2 (HR2), transmembrane anchor, and cytoplasmic tail (CT) (Figure ). The molecular basis of the virus transmission implies that the process starts with the binding of the spike protein to ACE2, which leads to a conformational change in the S1 and S2 domains of the spike protein and the exposure of the “fusion peptide” of the S2 domain, mediating the fusion of the viral and host cell membranes. The cryo-EM structures of MERS-CoV and SARS-CoV spike glycoproteins revealed that the spike protein trimer can exist in a dynamic equilibrium between a metastable closed (“RBD-down”) prefusion conformation that undergoes large structural rearrangements and hinge-based functional motions to adapt a receptor-accessible open (“RBD-up”) state in which the S protein can fuse the viral membrane with the host cell membrane.[22−24]

Figure 1

Domain and structural organization of the full-length SARS-CoV-2 spike (S) protein. (A) Schematic representation of domain organization and residue range for the full-length SARS-CoV-2 spike (S) protein. The subunits S1 and S2 include the N-terminal domain (NTD), receptor-binding domain (RBD), C-terminal domain 1(CTD1), C-terminal domain 2 (CTD2), S1/S2 cleavage site (S1/S2), S2′ cleavage site (S2′), fusion peptide (FP), fusion peptide proximal region (FPPR), heptad repeat 1 (HR1), central helix region (CH), connector domain (CD), heptad repeat 2 (HR2), transmembrane anchor (TM), and cytoplasmic tail (CT). (B) Structure of the wild-type full-length SARS-CoV-2 spike (S) protein trimer in the locked closed (“RBD-down”) prefusion conformation (PDB ID 6XR8). The structure of the S protein is shown as ribbons, and structural components of the trimer are annotated and colored according to the color scheme shown in panel A. The annotated structural components include the NTD, RBD, FP, FPPR, HR1, and CH. (C) Structure of the SARS-CoV-2 spike (S) protein trimer in the open (“RBD-up”) prefusion conformation state (PDB ID 6VYB) and the cryo-EM structure of the SARS-CoV-2 spike trimer (K986P/V987P) in the open state (PDB ID 6VSB). (D) Monomer form of the SARS-CoV-2 spike (S) protein trimer in the prefusion conformation state. Structural components of the trimer are annotated and colored according to the color scheme shown in panel A.

The cryo-EM structures of the SARS-CoV-2 spike ectodomain trimer in two distinct conformations detailed the conformational changes of the SARS-CoV-2 S trimer and characterized the receptor-accessible epitopes, confirming that SARS-CoV-2 and SARS-CoV bind with similar affinities to human ACE2.[25] The recent cryo-EM structure of the SARS-CoV-2 trimeric spike in the prefusion conformation reconstructed an asymmetrical trimer with a single RBD in the optimal “up” conformation in which SARS-CoV-2 trimer binds ∼10 times stronger to the ACE2 host cell receptor as compared with SARS-CoV protein.[26] The cryo-EM structures of the full-length human ACE2 in the presence of the neutral amino acid transporter B0AT1 with or without the RBD of the SARS-CoV-2 unveiled a complex structural arrangement with the concurrent binding of two S trimers to an ACE2 homodimer.[27] The two cryo-EM structures of the full-length S protein in the prefusion and postfusion conformations revealed the structural basis of the spontaneous transition to the postfusion state, discovering for the first time an ordered conformation of the ∼25-residue FPPR segment that is disordered in other structures and demonstrating that RBDs of the prefusion trimer can be locked in the closed form, arguably to allow long-range allosteric couplings to emerge between the stabilized RBD and FPPR regions.[28]

It was conjectured that allosteric cross-talk between the RBD and the FPPR regions, which may be mediated by the CTD1 segment, can regulate allosteric functional transformations between closed and host receptor-accessible open RBD forms, as the presence of a structured FPPR can restrict movements and lock down the RBD segment, whereas a disordered FPPR segment may release constraints on the RBD movements and promote the transition to the open “up” form.[28] A series of structural and biophysical studies employed protein engineering approaches to generate a prefusion-stabilized SARS-CoV-2 S protein by targeting regions that dramatically change between the prefusion and postfusion conformations as well as regions involved in the transitions between the open and closed prefusion states.[29] Through the design of disulfide bonds, proline mutations to curtail the mobility of flexible loops, and the introduction of hydrophobic modifications to fill empty voids in the prefusion structure, this study discovered a SARS-CoV-2 S variant with several proline substitutions that displayed ∼10 times higher expression than the original construct.[29] The cryo-EM structure of this engineered mutant retained the topology and molecular details of the prefusion trimer conformation.[29] By using structure-based vector analysis of spike glycoprotein structures, two soluble ectodomain constructs were designed for the SARS-CoV-2 S protein, in which the mobile RBD segment is either is locked in the “down” position or adopts “up” state conformations more readily than the wild-type S protein.[30] Through targeted mutagenesis design, negative-stain electron microscopy (NSEM), and cryo-EM structure determination, this study revealed a remarkable variability in the S trimer and the differential stabilization of the mutant ectodomain constructs, showing that the conformational equilibrium of the SARS-CoV-2 S trimer can be manipulated using modifications of contact regions between the RBD and S2 (S383C/D985C) and the RBD and NTD (D398L/S514L/E516L) and the interfaces between the SD1 and S2 (N866I/A570L).[30]

This pioneering study demonstrated that whereas the overall topological arrangements of the virus spike trimers are conserved, subtle differences in the interdomain contacts could play a major role in modulating the surface antigenicity and the equilibrium partitioning between the closed and open states. The cryo-EM analysis of the engineered SARS-CoV-2 S ectodomain trimer with several intermolecular disulfide bonds demonstrated that additional stabilization factors can effectively lock the trimer in the closed state without altering its architecture, enhancing the SARS-CoV-2 S resistance to proteolysis.[31] Another illuminating biophysical investigation combined targeted design and mutagenesis with structure determination of thermostable SARS-CoV-2 spike trimers, showing the subtle heterogeneity of the SARS-CoV-2 landscape concealed by the conserved trimer topology in which disulfide-bonded S-protein trimers using pairs (S383, D385) and (G413, V987) can be trapped in the structurally distinct locked and closed states of the prefusion form with only ∼20% of the population sampling one open state.[32] In the locked form, the S1 domains undergo a movement in which RBDs move closer to the three-fold axis and form a tightly packed trimer. At the same time, it was observed that for SARS-CoV-2 variants with proline substitutions K986P and V987C, ∼20% of the trimers are in the locked or closed state, and ∼80% of molecules adopt the open “up” trimer form. This study suggested that D427 and D428 in the RBD interact with the positively charged K986, which electrostatically stabilizes the closed form. The loss of these electrostatic contacts could act as a potential regulatory switch of the spike dynamic equilibrium by promoting a large allosteric change and shifting the population toward the open trimer.[32] The cryo-EM structures of the postfusion SARS-CoV S2 trimer showed dramatic conformational changes of the SARS-CoV S2 machinery during membrane fusion, where the HR1 and the CH helical segments can form a long continuous stem helix pointing toward the target membrane in the postfusion trimer.[33]

The shielding of receptor binding sites by glycans is a common feature of viral glycoproteins, as observed on SARS-CoV proteins, where glycosylation can camouflage immunogenic protein epitopes.[34,35] Site-specific glycan mass spectrometric analysis of the SARS-CoV-2 spike further enabled mapping of the glycan-processing states,[36] revealing the structural positions of 22 glycans on the SARS-CoV-2 S protein and highlighting vulnerabilities in the RBD regions on the exposed surfaces of the trimer. The generated computational models of glycoforms of the SARS-CoV-2 spike protein and molecular dynamics (MD) simulations determined the extent to which glycan heterogeneity can affect the antigenicity of the S glycoprotein.[37] This computational study detailed at the molecular level that the protein surface is extensively shielded from antibody recognition by glycans, with the exception of the RBD regions that present points of vulnerabilities for the glycan shield.[37] The development and all-atom modeling of a fully glycosylated full-length SARS-CoV-2 S protein in a viral membrane were recently reported, representing a significant technical breakthrough in enabling rigorous simulations of the SARS-CoV-2 S trimer structures in a glycosylated environment (Figure S1).[38] MD simulations have been instrumental in predicting the shapes and motions of glycans for glycoproteins.[39−41] The unprecedented level of atomistic detail of the all-atom MD simulations of the full-length SARS-CoV-2 S glycoprotein with a complete glycosylation profile established that the dynamics of the glycan shield can be allosterically coupled to the conformational changes, dictating the heterogeneity of the conformational landscape and the response to the host receptor.[42] MD simulations of the SARS-CoV-2 spike glycoprotein identified the changes in the molecular properties due to conformational flexibility.[43] Network-based computational methods were developed to identify allosteric sites on the SARS-CoV-2 spike protein.[44] Several computational studies examined SARS-CoV-2 interactions with the ACE2 enzyme.[45−47] MD simulations also elucidated the determinants of the higher affinity of SARS-CoV-2 with ACE2, highlighting the role of the RBM motif[46] as well as role of the hydrophobic interactions and the elaborate hydrogen-bonding network in the SARS-CoV-2-RBD interface.[47]

Computational investigations have led to important methodological advances and an improved atomistic understanding of allosteric regulation in proteins.[48−51] The growing number of the SARS-CoV-2 S protein structures in different conformational states suggests that allosteric mechanisms may regulate functional transitions and the acquisition of host-accessible conformations for the SARS-CoV-2 S prefusion trimer. The current study introduces a hypothesis that the SARS-CoV-2 spike protein can function as an allosteric regulatory engine that fluctuates between dynamically distinct functional states and is controlled by allosteric switch centers that may determine the signal transmission and the regulation of virus entry. Coarse-grained (CG) simulations are employed together with perturbation-response scanning (PRS) of allosteric interactions and hierarchical network modeling to identify the salient signatures of allosteric interaction networks in the distinct conformational states of the SARS-CoV-2 spike trimers. The analysis of functional dynamics and the perturbation-response modeling identify the regulatory hotspots that drive transformations between locked, closed, and open SARS-CoV-2 spike forms. The results show that the stabilized locked form of the prefusion trimer is characterized by a broad allosteric network with key regulatory centers located in the RBD, FPPR, and HR regions that dictate the dynamic switching between conformational states of the SARS-CoV-2 spike trimer. This study offers a useful complementary perspective on the mechanisms underlying the conformational and dynamic changes in the SARS-CoV-2 spike trimer through the lens of functional dynamics and allosteric communications that enable us to identify the key residues and interactions responsible for the control of signal transmission and the regulation of virus entry to the host receptor.

Materials and Methods

Sequence Conservation Analysis

Sequence conservation for SARS spike proteins was estimated using the ConSurf approach[52−56] by computing the residue-based conservation score profiles that measure evolutionary conservation. The low-score values are associated with the most conserved position in the protein. Multiple sequence alignment (MSA) was obtained using the MAFFT approach,[57] and homologues were obtained from UNIREF90.[58,59] MSA profiles were also obtained from the Pfam database of protein families (P59594, SPIKE_CVHSA, and CoV_S1_C, PF19209).[60−62]

Coarse-Grained Simulations and Elastic Network Models

CG models are computationally effective approaches that leverage a topology-based framework for the characterization of the protein structure and dynamics, enabling simulations of large systems over long time scales. In this study, the CABS CG model is used for simulations of the cryo-EM structures of the SARS-CoV-2 spike trimer in distinct conformational states (Table S1). This approach efficiently combines a CG model and Monte Carlo local moves, allowing the accurate reproduction of conformational ensembles of protein structures from all-atom MD simulations.[63−67] The sampling scheme of the CABS model employed for this study is based on replica-exchange simulations. The number of cycles during the CABS-CG simulations was set to 10 000, with the number of cycles between frames set to 100. Multiple CG simulations were performed for the distinct functional states of the SARS-CoV-2 S spike protein trimer (Figure ). A total of 1 000 000 samples was generated for each system, and the total number of saved models in the trajectory used for the analysis was 10 000.

Figure 2

Cryo-EM structures of the SARS-CoV-2 spike (S) protein trimer in the locked and closed states of the prefusion form used in this study. Upper Panel: The following structures of the SARS-CoV-2 spike (S) protein trimer in the locked form were used in the simulations and modeling: the cryo-EM structure of the wild-type full-length S protein in the locked prefusion conformation (PDB ID 6XR8),[28] the cryo-EM structure of the disulfide-stabilized SARS-CoV-2 spike trimer (S383C, D985C, K986P, V987P, single R S1/S2 cleavage site) in the locked state (PDB ID 6ZOZ),[32] and the cryo-EM structure of the SARS-CoV-2 spike trimer in the locked state (K986P, V987P, single R S1/S2 cleavage site) (PDB ID 6ZP2).[32] Lower Panel: The following structures of the SARS-CoV-2 spike (S) protein trimer in the closed form were used for the simulations and modeling: the cryo-EM structure of the SARS-CoV-2 spike trimer in the closed state (K986P/V987P) (PDB ID 6VXX),[25] the cryo-EM structure of the disulfide-stabilized SARS-CoV-2 spike trimer (S383C, D985C, K986P, V987P, single R S1/S2 cleavage site) in the closed state (PDB ID 6ZOY),[32] and the cryo-EM structure of the SARS-CoV-2 spike trimer in the closed state (K986P, V987P, single R S1/S2 cleavage site) (PDB ID 6ZP1).[32] The following structures of the SARS-CoV-2 spike (S) protein trimer in the open form were subjected to simulations and modeling: the cryo-EM structure of the SARS-CoV-2 spike trimer ectodomain (K986P/V987P) in the open state (PDB ID 6VYB)[25] and the cryo-EM structure of the SARS-CoV-2 spike trimer (K986P/V987P) in the open state (PDB ID 6VSB).[26]

The structures of the SARS-CoV-2 spike proteins were obtained from the Protein Data Bank (PDB).[68,69] The recently reported fully glycosylated SARS-CoV-2 spike trimer atomistic models of the closed form (PDB ID 6VXX) and open forms (PDB ID 6VSB) (residues 1–1146)[38] were obtained from the CHARMM-GUI COVID-19 archive (http://www.charmm-gui.org/docs/archive/covid19) (Figure S1). The structure preparation process of other cryo-EM structures leveraged the developed full atomistic models for the closed and open states.[38] The missing residues were initially modeled based on these structures. Additionally, and for comparison, the missing loops in the cryo-EM structures were also reconstructed using template-based loop prediction approaches ModLoop,[70] ArchPRED server,[71] and DaReUS-Loop.[72]

The generated CABS-CG conformational ensembles were also subjected to all-atom reconstruction using the PULCHRA method[73] and the CG2AA tool.[74] The protein structures were then optimized using atomic-level energy minimization using the 3Drefine method.[75] The principal modes of protein motions are extracted from CG approaches and also using elastic network models (ENMs). In this work, the collective motions of the protein structures were also modeled using the Gaussian network model (GNM)[76,77] and the anisotropic network model (ANM) approaches[78] that are the most widely used ENM-based methods, computing the amplitudes of isotropic thermal motions and the directionality of anisotropic motions. The algorithmic details of this approach were discussed at length in related studies of molecular chaperones.[79] The essential mobility profiles were obtained using the DynOmics server[80] and the ANM server.[78]

Perturbation-Response Scanning

The PRS approach[81,82] estimates the residue response to external forces and was successfully used to locate allosteric hotspots and regulatory centers in various protein systems.[83−89] The PRS approach can also be efficiently combined with the dynamic analysis of residue interaction networks to identify and characterize allosteric hotspots and pathways of allosteric communications.[85−91] By monitoring the response to forces on the protein residues, the PRS approach can quantify allosteric couplings and link the directionality of the inserted force to the protein response in functional movements.[90,91] In the PRS approach, the 3N × 3N Hessian matrix of the second derivatives of the potential at the local minimum is computed. The 3N-dimensional vector Δ of node displacements in response to the 3N-dimensional perturbation force is evaluated according to Hooke’s law = *Δ. A perturbation force is applied to each residue one at a time, and the response of the protein is estimated using the displacement vector Δ() = –1. The deviations of each residue in the protein from the mean structure of the equilibrium ensemble and the covariance matrix ΔΔ were calculated.

Each residue in the SARS-CoV-2 spike structures is perturbed one at a time by applying random forces distributed along a sphere of randomly selected directions. Using the residue displacements obtained from external force perturbations, the magnitude of the response of residue k is computed as ⟨|Δ(|2⟩ averaged over multiple perturbation forces , yielding the ikth element of the N × N PRS matrix. A measure of the response of residue k is the magnitude ⟨|Δ(|2⟩ of the kth block of Δ averaged over multiple F(, expressed as the ikth element of the N × N PRS matrix, SPRS. The effector profile ⟨(Δ)2⟩effector represents the average effect of the perturbed effector site i on all other residues and is computed by averaging over all sensor residues j. Similarly, the jth column of the PRS matrix S̅PRS represents the sensitivity profile of sensor residue j in response to perturbations of all residues, and its average is denoted as ⟨(Δ)2⟩sensor.

Modeling and Community Analysis of Residue Interaction Networks

The residue interaction networks are represented as graphs with residues as network nodes and the inter-residue edges describing the inter-residue interactions.[92−94] The algorithmic details of constructing the residue interaction networks were presented and discussed in the previous studies.[95−97] The edges in the residue interaction networks are weighted by dynamic cross-correlations[94] and coevolutionary couplings.[97] The ensemble of shortest paths was determined by the Floyd–Warshall algorithm.[98] Network graph calculations were performed using the Python package NetworkX.[99] The betweenness of residue i is computed as the sum of the fraction of shortest paths between all pairs of residues that pass through residue iwhere 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 modified version of the Girvan–Newman algorithm[100−102] is used to identify local communities. The community centrality measure was computed on a hierarchical meta-network whose nodes are the local communities using the following expression[103,104]N is the number of communities to which node i belongs, and S is the Jaccard similarity coefficient between communities j and k. The network parameters were also evaluated using the ModuLand program in the Cytoscape package.[105−107]

Results and Discussion

Sequence Analysis of the SARS-CoV-2 Spike Glycoproteins Reveals the Sharp Partition of Evolutionary Conservation between S1 and S2 Subunits

To understand the interplay between evolutionary, topological, and dynamic signatures of the SARS-CoV-2 spike glycoprotein trimer, the sequence conservation profiles are computed to identify regions of high conservation and variability (Figure ). The sequence conservation ConSurf approach[52−56] was used to determine the residue-based conservation score profiles in which the low score is associated with the most conserved position in the protein. Consistent with previous studies,[108−110] this analysis showed that residues in the S1 domain are considerably less conserved than the S2 fusion domain. In particular, the NTD regions and the exposed RBM segment that forms a contact interface with the host receptor displayed a considerable variability (Figure A). In some contrast, the core region in SARS-CoV-2 RBD (residues 331–438) exhibited strong conservation, exemplified in this model by the low Consurf scores (Figure A). The RBM motif includes highly variable residues 472–486 that are not present in S proteins of coronavirus isolated in bats.[108] Among the variable residues are G482, Y484, F486, F490 N501, E471, E472, E468, and K444, which are consolidated near the flexible ridge loop subjected to structural changes upon the binding of SARS-CoV-2 RBD to the ACE2 receptor (Figure A). The structural studies of the SARS-CoV spike protein in the complex with ACE2 suggested that variable positions L472, N479, and T487 that are substituted for F486, Q493, and N501 in the SARS-CoV-2-RBD could be critical for binding selectivity.[20,21] Structural mapping of highly variable sites in the SARS-CoV-2 trimer highlighted the density of variable positions in the NTD and RBD regions (Figure A).

Figure 3

Sequence conservation analysis of the SARS-CoV-2 spike (S) glycoprotein. (A) Schematic representation of the domain organization and residue range for the full-length SARS-CoV-2 spike (S) protein. The subunits S1 and S2 include the N-terminal domain (NTD), receptor-binding domain (RBD), C-terminal domain 1 (CTD1), C-terminal domain 2 (CTD2), S1/S2 cleavage site (S1/S2), S2′ cleavage site (S2′), fusion peptide (FP), fusion peptide proximal region (FPPR), heptad repeat 1 (HR1), central helix region (CH), connector domain (CD), heptad repeat 2 (HR2), transmembrane anchor (TM), and cytoplasmic tail (CT). Annotation is as in Figure . (B) Normalized ConSurf conservation scores for the SARS-CoV-2 spike (S) glycoprotein projected onto the structure of the wild-type full-length SARS-CoV-2 spike (S) protein trimer in the locked closed (“RBD-down”) prefusion conformation (PDB ID 6XR8).[28] The ConSurf profiles are shown as colored bars highlighting the conservation in different segments of the S protein. In particular, NTD residues are shown as maroon bars and RBD residues are shown as orange bars. The low negative ConSurf scores correspond to highly conserved sites, and high positive scores depict highly variable positions. (C) Structural mapping of variable positions onto the structure of the wild-type full-length SARS-CoV-2 spike (S) protein trimer (PDB ID 6XR8). The structure is shown as ribbons with reduced transparency, and the three protomers (A, B, and C) are colored green, red, and blue, respectively. The highly variable positions are shown as spheres colored according to the corresponding protomer (green, red, and blue). A high density of variable sites can be seen in the NTD and RBD regions.

The furin cleavage site at the S1/S2 boundary is in a surface-exposed and disordered loop, and a small cluster of variable positions in this region is consistent with its functional role (Figure ). Although the evolutionary analysis identified important variable RBM regions implicated in binding, a mix of conserved and variable positions in the RBD suggested that the interplay of stability and flexibility may enable dynamic adaptation of the spike protein to the host receptor and can be important for rendering the binding selectivity of the SARS-CoV-2 glycoprotein.

The S2 domain showed considerably higher conservation, which was pronounced for the FP region (812-SFIEDLFNKVTLADAGF-829), the HR1 region (residues 910–985), the CH regions (residues 986–1035), and the CD residues (Figure ). The high sequence conservation of the FP region is consistent with its functional role as a viral fusion motif that is composed of both small and bulky hydrophobic residues. The FPPR segment (residues 828–853) that is immediately downstream of the FP region also featured a high degree of conservation (Figure A). Interestingly, this functional segment was disordered in the SARS-CoV-2 ectodomain structures but adapted a well-defined helix-turn-helix motif in recently reported cryo-EM structures of the wild type and K986P/V987P mutants of the SARS-CoV-2 S trimer.[28]

Given the amino acid sequence conservation of the S2 subunit (Figure ), it is possible that coronaviruses can operate through similar regulatory hotspots and regions to mediate virus–host membrane fusion in a unified manner. Interestingly, the cascade of tectonic conformational changes leading to the formation of the postfusion form implies a mechanism in which RBD binding to the host receptor and cleavage of the S2′ site can accelerate the release of the S1 subunit and refolding of helical regions in the remaining S2 trimer.[33,111] These dramatic rearrangements involve highly conserved HR1 and CH motifs that become exposed and vulnerable to the external intervention after dissociation of the S1 subunits.

To summarize, the striking partition of conservation and variability in the SARS-CoV-2 spike proteins between S1 and S2 subunits can be also embedded in the topology and drive the functional dynamics of the spike trimer in the prefusion form. In particular, this analysis suggested that conserved and structurally rigid elements of the S2 subunit can control and mediate dynamic changes on the landscape of the prefusion trimer, whereas highly variable NTD and RBD regions provide structural adaptation to external binding partners and host receptors.

Molecular Simulations of the SARS-CoV-2 Spike Trimers Characterize Salient Dynamic Signatures and the Redistribution of Mobility in Distinct Functional States

Ideally, all-atom MD simulations of the full-length SARS-CoV-2 S structures in different functional states with a complete glycosylation shield profile and a rigorous solvent description are required for the rigorous assessment of atomistic fluctuations and interaction details. The recent pioneering computational studies of the SARS-CoV-2 S protein with a complete and explicit glycosylation shield description[42] underscored the technical challenges of performing these simulations, which are time-consuming even with dedicated supercomputer resources. Accordingly, these simulations would require significant time and enormous resources for the execution of multiple SARS-CoV-2 S structures in fully glycosylated environments. The objective of this investigation is to provide a useful complementary analysis of topology-based dynamic signatures that are characteristic of the major functional states of the full-length SARS-CoV-2 S trimer by employing a large number of efficient CABS-CG simulations and subsequent atomistic reconstruction. Using these simulation approaches applied to full atomistic high-resolution structures of the SARS-CoV-2 S prefusion trimer in multiple functional states, a detailed comparative analysis of the conformational dynamics profiles was performed, and salient features of the dynamic conformational landscapes for the SARS-CoV-2 S trimer were analyzed (Figure ). Whereas the dynamic fluctuations and variations of individual residues may be partly affected due to the absence of the explicit glycosylation shield, these molecular simulations can adequately reproduce global topological and dynamic patterns, allowing for a comparative analysis and enabling key differences between the dynamics of functional states to be pinpointed.

Figure 4

CABS-CG conformational dynamics of the SARS-CoV-2 spike (S) protein trimer in the locked (A), closed (B), and open states (C) of the prefusion form. (A) Root-mean-square fluctuation (RMSF) profiles obtained from CABS-CG simulations of the cryo-EM structure of the disulfide-stabilized SARS-CoV-2 S trimer (S383C, D985C, K986P, V987P, single R S1/S2 cleavage site) in the locked state (PDB ID 6ZOZ).[32] (B) RMSF profiles from simulations of the cryo-EM structure of the SARS-CoV-2 S trimer in the closed state (K986P/V987P) (PDB ID 6VXX).[25] The complete atomistic model of the full-length SARS-CoV-2 S protein developed in the study[38] was used in the simulations. (C) RMSF profiles from simulations of the cryo-EM structure of the SARS-CoV-2 S trimer ectodomain (K986P/V987P) in the open state (PDB ID 6VYB).[25] The profiles for protomer chains A, B, and C are shown as green, red, and blue lines, respectively. The conformational dynamics profile of the spike monomer is shown as an orange line.

CABS-CG simulations of SARS-CoV-2 spike trimer in the locked, closed, and open forms of the perfusion state revealed important dynamic signatures of these states (Figure ). Simulations of single monomers for each of the SARS-CoV-2 spike structures were also performed. A comparative analysis of the dynamic profiles generally showed very moderate thermal fluctuations of both the S1 and S2 subunits in the locked state (Figure A). Whereas the NTD regions and the furin cleavage site (residues 680–688) showed larger displacements, the mobility of the RBD residues was markedly suppressed in this form as compared with the more dynamic closed and open prefusion states (Figure A–C). Another important dynamic signature of the locked state is the thermal stabilization of the FP region and the structurally ordered FPPR motif (residues 828–853).

It was conjectured that CTD1 (residues 529–591) is a structural relay between RBD and FPPR (residues 828–853) that can sense the displacement on either side and communicate the signal from and to the fusion peptide. These segments of the S1 subunit are largely stable in the locked state of the S trimer (Figure A). This is in sharp contrast with the dynamic profiles of the closed state (Figure B) and open state (Figure C), where this segment becomes highly flexible and can be partially disordered. In addition, it was observed that both closed and open states featured significantly larger thermal fluctuations in the NTD and RBD regions, in particular, of the exposed RBM segment involved in the recognition of the host receptor. As expected, CG simulations of the open form of the trimer revealed the large fluctuations in the “up” monomer that become exposed for interactions with the host receptor (Figure C). Domains HR1, CH, and CD close to the viral transmembrane exhibited the least movement during MD simulations.

Of particular interest was a comparison with multiple simulations of the monomeric forms. Somewhat unexpectedly, the RBD fluctuations were also reduced in the monomer of the locked SARS-CoV-2 spike trimer (Figure A). Although this analysis is based on CABS-CG simulations that may underestimate the extent of the mobility of the isolated monomer, the comparison with other states clearly indicated that the monomeric form becomes highly flexible in the closed and open forms, with the RBD undergoing large conformational changes between the “down” and “up” forms (Figure ). These simulations captured the fundamental topological features of the distinct conformational states, suggesting that despite the structural similarities, the dynamics of the locked and closed states could be quite different.

Structural analysis suggested that the absence of an ordered FPPR segment in the closed states of the SARS-CoV-2 spike trimer may be linked to the higher mobility of the RBD and the emergence of the RBD-up conformation.[32] In fact, the “up” form of the RBD was not detected in the locked states of the trimer with the ordered FPPR segment. These simulations are consistent with this notion, showing that the dynamics of the FPPR and RBD regions may be allosterically coupled in all states of the SARS-CoV-2 spike trimer. The results also suggested that the RBD segment may be intrinsically predisposed to conformational selection between the “up” and “down” forms that can be modulated and stimulated (or suppressed) through the long-range couplings with the FP and FPPR motifs. Moreover, such redistribution of thermal fluctuations is more pronounced in SARS-CoV-2-RBD, which is coupled to a stronger stabilization of the entire binding interface. These findings imply that the binding mechanism of virus entry into the host receptor may exhibit signs of dynamically driven allostery, which is typically exemplified by the lack of structural changes between the unbound and bound forms, coupled to an exchange of conformational mobility between local protein regions.[99−102]

Functional Dynamics of the SARS-CoV-2 Spike Trimer Reveals the Migration of the Hinge Sites and the Release of Dynamic “Breaks” on the RBD Motions in the Closed State

To characterize collective motions and determine the distribution of hinge regions in the SARS-CoV-2 spike trimers, principal component analysis (PCA) of all-atom reconstructed trajectories was undertaken along with ENM analysis of the essential slow modes.[80] Previous studies indicated that CABS-CG MC simulations and ENM computations can provide similar dynamic profiles that closely reproduce the experimental dynamic data for a wide range of protein structures.[65] The analysis of collective motions revealed functionally important domain motions and hinge regions in the SARS-CoV-2 spike trimer structures (Figure ). The simplified CG models used in the analysis did not affect the accuracy of the results and significantly reduced the computational cost so that the collective motions can be analyzed in great detail for all available SARS-CoV-2 trimer structures in the locked, closed, and open forms. By leveraging the efficiency of these methods, the essential profiles for the slowest modes were computed, and the functional dynamics profiles averaged over the first three major low-frequency modes were analyzed in detail (Figure ). As expected within the realm of CG models, the overall shape of slow mode profiles was generally preserved in structurally similar locked and closed states but was largely altered in the open form of the trimer. Nonetheless, a detailed analysis of the hinge regions that correspond to local minima along the slow mode profiles revealed subtle important changes between the locked and closed states of the prefusion trimer form (Figure , Table S1).

Figure 5

Functional dynamics and analysis of collective motions in the locked, closed, and open states of the SARS-CoV-2 spike trimer prefusion form. (A) Mean-square fluctuations averaged over the three lowest frequency modes for the locked state of the cryo-EM structure of the disulfide-stabilized SARS-CoV-2 S trimer (PDB ID 6ZOZ).[32] (B) Essential mobility profiles averaged over the three lowest frequency modes for the cryo-EM structure of the SARS-CoV-2 S trimer in the closed state (PDB ID 6VXX).[25] (C) Essential mobility profiles along the three slowest modes for the cryo-EM structure of the SARS-CoV-2 S trimer in the open state (PDB ID 6VSB).[25] The profiles for protomer chains A, B, and C are shown as green, red and blue lines, respectively. Structural maps of the essential mobility profiles for the locked state of the SARS-CoV-2 S prefusion trimer (D), closed state (E), and open state (F). The color gradient from blue to red indicates the decreasing structural rigidity and increasing flexibility. Notably, in the locked state of the SARS-CoV-2 spike trimer, the RBD residues are generally aligned with shallow minima positions, suggesting that large functional movements of RBD may be largely restricted and suppressed in the locked trimer form (Figure A). Several other local hinge points were precisely aligned with the RBD residues K386/L387 and D427/D428, and these hinge sites are unique for the locked form of the trimer (Figure A). The latter finding is of special interest and significance, as the biophysical studies showed that the stabilizing interactions of D427 and D428 with K986 are vital for the thermodynamic stabilization of the closed trimer and the formation of the rigid locked state.[32] The presented results indicated that these RBD residues could act as important functional hinge sites that regulate and restrict conformational transformations of RBD in the locked state.

The RBD regions of the S1 subunit are linked through two antiparallel hinge linkers to the NTD and CTD2 and allow the transition between closed and open conformations. Uniquely for the locked state, the hinge regions that can regulate the interdomain movements between RBD and NTD are fairly broadly distributed, forming shallow local minima aligned with a helical linker (residues 298–306) and extending through a flexible loop to another hinge site corresponding to a β-strand segment (residues 324-ESIV-327) that is connected to the RBD region (Figure A, Table S1). Interestingly, the recent protein engineering studies demonstrated that the dynamic equilibrium of the SARS-CoV-2 spike trimer can be manipulated using modifications of contact regions between the RBD and NTD through a triple D398L/S514L/E516L mutant[30] which can now be explained by the fact that D398 belongs to an important hinge site, and mutations in this regulatory switch position could alter the global equilibrium and distribution of states in the prefusion form of the trimer.

Together with the experimental evidence,[32] this suggested that D427/D428 and K986 may play the role of the specific regulatory switch of the spike dynamic equilibrium, where mutations or dynamic changes in these positions could trigger an allosteric conformational change and population shift between the closed and open trimers. In addition, the hinge sites in the locked state also corresponded to residues 570–572, 750–755, and 986–990 (Figure A). Similarly, modifications of hinge points A570 and T572 near interfaces between SD1 and S2 and SD1 and S2 domains can change the population of locked, closed, and open forms.[30] These positions were selected in the mutagenesis engineering studies to produce variants A570L, T572I, F855Y, and N856I of the SARS-CoV-2 S trimer that destabilized the closed state rather than stabilize the open state.[30]

Despite obvious structural similarities of the locked and closed trimer states, the essential mobility profile showed important changes reflecting a more flexible nature of the metastable closed form (Figure B). The most significant difference can be seen in the S1 subunit with only a single well-defined hinge site observed at the border of the NTD and RBD regions and aligned with the β-sheet motif (residues 315–320) (Figures B and 6A,B). Interestingly, the key hinge site near D427/D428 residues considerably weakened and corresponds to a very minor local minimum along the profile (Figure B), indicating that dynamic “brakes” on the RBD motions can be largely released in the closed state that becomes prone to large conformational transformations. Moreover, the RBD regions were mostly aligned with the local maxima of the slow mode profile, indicating that these regions in the closed form can experience functional movements between “down” and “up”’ forms as the dynamic constraints on the RBD repositioning are removed in a more flexible closed state (Figure B). Other hinge sites were preserved in the closed trimer, showing the presence of local minima for residues 570–572 and 856–862 (Figure B). These positions may be involved in regulating functional motions near the interdomain regions and overlapped with positions used in protein design experiments (N866I/A570L, A570L/T572I/F855Y/N856I) to manipulate the dynamic equilibrium and promote a population shift of the SARS-CoV-2 spike protein.[30] Structural projection of the essential mobility profiles onto the SARS-CoV-2 S trimer structures in the locked state (Figure D) and closed state (Figure E) illustrated the migration of the immobilized regions away from the RBD during the transition to the dynamic closed form. The observed changes in the distribution of rigid and flexible regions in slow modes reflected the release of dynamic breaks on functional motions of the RBD in the closed state (Figure D,E). In the open form, the main hinge sites corresponded to residues 42–44 in the NTD and β-sheet motif (residues 315–320) (Figure C).

Figure 6

Structural maps of the hinge regions in the closed and open states of the SARS-CoV-2 spike trimer prefusion form. (A) Projection of hinge sites (shown as yellow spheres) onto the cryo-EM structure of the SARS-CoV-2 S trimer in the closed state (PDB ID 6VXX).[25] The SARS-CoV-2 S trimer is shown as ribbons with protomers A, B, and C in green, red, and blue, respectively. (B) Close-up map of hinge sites presenting a single protomer from the closed prefusion state. The protomer is shown as red ribbons, and hinge sites are highlighted with yellow spheres and annotated. (C) Projection of hinge sites (shown as yellow spheres) onto the cryo-EM structure of SARS-CoV-2 S trimer in the open state (PDB ID 6VSB).[25] (D) Close-up map of hinge sites presenting a single protomer undergoing the transition to the “up” prefusion state. The protomer is shown as red ribbons, and hinge sites are highlighted with yellow spheres and annotated.

These regulatory positions of collective dynamics are associated with concerted functional movements of the NTD and RBD regions around the respective hinges. Notably, in the monomer undergoing a transition to the “up” conformation, the RBD regions are aligned with maxima along the profiles reflecting this large change (Figures C and 6C,D). The collective motions of the other two monomers corresponded to synchronous changes in the NTD, whereas the RBD regions remained in their locked closed positions (Figure C), suggesting that NTD adjustments may provide necessary room for the upward movement of the RBD in the single protomer of the trimer. Structural mapping of the slow mode profiles in the open state (Figure F) showed redistribution in the density of moving regions in which the RBD of the upward moving monomer cooperates with functional motions of the NTDs in the other two monomers. Overall, the results suggested that collective dynamics of the spike trimer may be driven by coordinated movements of the NTD and RBD regions around the more dynamically restricted S2 subunit.

The positions of 22 glycosylation sites of the S protein that constitute the glycan shield of the SARS-CoV-2 spike trimer[36] were aligned onto structural maps of functional motions in slow modes for the locked, closed, and open states (Figure S2). It is evident that the glycan shield can be especially effective at enhancing the protection of the stable locked state, as the majority of the moving regions tend to be localized in the NTD regions. In some contrast, in a more dynamic closed state of the prefusion S trimer, the glycan shield may display potential vulnerabilities in the RBDs, exposing the evolutionary variable RBD regions that undergo collective movements (Figure S2). As a result of the reduced protection, the more flexible closed state of the spike trimer could become susceptible to large-scale transformations and readily undergo allosteric structural changes to the open form, in which one RBD protomer is exposed for binding with the host-cell receptor.

Perturbation-Response Scanning Identifies Regulatory Hotspots of Allosteric Interactions in Different Conformational States of the SARS-CoV-2 Spike Trimer

Using the PRS method,[81−91] the allosteric effect of each residue in the protein structures in response to external perturbation was quantified and examined. PRS analysis produced the residue-based effector and sensor response profiles in different functional states of the SARS-CoV-2 S trimer (Figures and 8). The effector profiles assess the propensities of a given residue to influence dynamic changes in other residues and are often applied to identify the regulatory hotspots of allosteric interactions as the local maxima along the profile (Figure ), whereas the sensor/receiver profiles measure the ability of residues to serve as transmitters of allosteric perturbations and allow for the determination of the regions undergoing large structural changes during allosteric transformations. The peaks of the effector profiles may point to the allosteric sites that control the allosteric signal transmission and conformational transitions between the distinct functional states of the SARS-CoV-2 S protein. First, the effector profiles of the SARS-CoV-2 S homotrimer were analyzed in the locked closed form (Figure A), the closed form (Figure B), and the open receptor-accessible form (Figure C). In the locked form, the effector profiles featured three major peaks corresponding to the RBD residues 426-LPDDF-429, residues 1026–1040 near the homotrimer cavity, residues 883–893 and 986/987 (K986P/V987P), and also residues 850–855 from the FPPR region (residues 833–855) (Figure A). The peaks in the RBD region (D427/D428) and the FPPR motif appeared to be uniquely characteristic of the effector profile in the locked state (Figure A), as these peaks and the effector density in the RBD regions were considerably weakened in both the closed (Figure B) and open forms (Figure C, Table S1).

Figure 7

PRS effector profiles in the locked, closed, and open states of the SARS-CoV-2 spike trimer prefusion form. (A) PRS effector profile for the SARS-CoV-2 S trimer in the locked state (PDB ID 6ZOZ).[32] (B) Effector profile for the SARS-CoV-2 S trimer in the closed state (PDB ID 6VXX).[25] (C) Effector profiles for the SARS-CoV-2 S trimer in the open state (PDB ID 6VSB).[25] The profiles for protomer chains A, B, and C are shown as green, red, and blue lines, respectively. Structural maps of the PRS effector profiles for the locked state of the SARS-CoV-2 S prefusion trimer (D), closed state (E), and open state (F). The color gradient from blue to red indicates the increasing effector propensities. The clusters of residues with the high allosteric potential corresponding to the peaks of the effector profile are shown as red spheres.

Figure 8

Detailed structural analysis of PRS effector profiles in the locked state of the SARS-CoV-2 spike trimer prefusion form. (A) Structural map of the PRS effector profiles for the wild-type full-length SARS-CoV-2 spike (S) protein trimer in the locked closed (“RBD-down”) prefusion conformation (PDB ID 6XR8).[28] The color gradient from blue to red indicates the increasing effector propensities. The key effector centers and corresponding interaction clusters are shown as spheres colored according to the effector propensity. (B) Close-up of the interprotomer interaction cluster formed by effector hotspots D427, D428, and K986 residues. The residues are shown as sticks colored by atom type. (C) Close-up of the interprotomer interaction cluster formed by the effector center K854 with D614.

Interestingly, the effector peak in the RBD region is precisely aligned with the hinge residues D427 and D428 that form the interprotomer electrostatic interactions with K986. Although K986 and V987 are mutated to the neutral proline in the locked structures of the SARS-CoV-2 trimer, this interacting cluster retains its structural integrity and is aligned with the unique allosteric hotspot in the RBD (Figure A). Another important signature of allosteric interactions in the locked trimer form is the emergence of the effector peak (residues 850–855) in the FPPR region, which is N-terminal to the FP and S2′ site. The key structural difference of the locked trimer form is the ordering of the FPPR segment into a well-defined motif, whereas in the majority of SARS-CoV-2 spike trimer structures of the closed and open forms, this region is poorly defined and highly flexible. In addition, the effector sites in FPPR are adjacent to K854, which forms a salt bridge with residue D614 that is known to be critical for the stability of the closed form (Table S1).

It may be argued that the regulatory center in the RBD (D427/D428) can communicate with the effector center in the FPPR region to cooperatively control allosteric couplings and interactions in the locked closed state of the SARS-CoV-2 spike trimer. Interestingly, some of the predicted effector sites are located in close proximity to residues that were engineered to produce a prefusion-stabilized SARS-CoV-2 spike protein that is more stable than the original construct.[29] For example, the effector residues 886–890 include S884 and A893 positions that were covalently linked by cysteine substitutions S884C/A893C and T791C/A879C, leading to the increased thermal stability of the closed form by improving the intermonomer packing.[29] Notably, some of these sites also coincided with the major hinge points of collective motions, suggesting that these positions may be involved in the regulation of allosteric interactions in the SARS-CoV-2 spike trimer (Table S1).

The effector profiles were also computed for other SARS-CoV-2 trimer structures in the locked state, including the full-length, wild-type SARS-CoV-2 spike trimer with the structurally determined FPPR region (Figure S3). The effector profiles were similar, highlighting the unique peak in the RBD region that is aligned with the D427/D2428 residues and indicating the critical role of this allosteric hotspot in signal transmission.

In the closed form of the SARS-CoV-2 spike trimer, the major effector peaks corresponded to residues 567–569 and 756–758 in the S2 CH repeat region, residues 886–890, as well as residues 1038–1040 in the CH/CD regions (Figure B, Figure S4). These characteristic peaks are shared in other closed trimer conformations (Figure S4), revealing an important redistribution of the effector sites and indicating the reduced allosteric potential of the RBD regions that in the closed form function mainly as sensors of allosteric signals. This redistribution of the effector allosteric centers and migration of the effector density is evidenced from structural mapping (Figure D,E). Notably, this dynamic switching of control points occurs in the absence of significant structural transformations between the locked and closed forms, suggesting that dynamically driven allostery[112−115] can govern the equilibrium between the locked and closed states.

In the open trimer conformation, the shape of the effector profile (Figure C) displayed the further migration of the effector density away from the RBD regions (Figure F) and exhibited density in the CTD1 region (residues 529–591) that is implicated as a structural relay between RBD and FPPR. Nonetheless, several clusters of conserved effector peaks were retained in the S2 subunit (residues 756–758, 887–891, and 1035–1045 in the S2 domain) (Figure C, Table S1). Structural mapping of effector hotspot clusters highlighted key findings, showing the evolution of the effector density in the trimer, from high and broad density reaching the RBD regions in the locked closed state to a narrower density in the closed and open forms, in which allosteric couplings between FPPR and RBD regions are weakened (Figure D–F). These results are consistent with the notion that the locked conformation can represent an early intermediate state before opening of the RBD and binding to the host receptor.[32] The key finding of the PRS effector analysis is the emergence of a broad and dense allosteric network with multiple regulatory centers in the locked state and a cross-talk between the effector hotspots in the RBD and FPPR regions (Figure ). Through these regulatory control points, allosteric signaling dictates the stable, down-regulated arrangement of the RBD regions. Interestingly, the allosteric interaction network may be weakened, become smaller, and diffuse in a more dynamic metastable closed form, where the effector centers regulate allosteric transitions. The allosteric couplings between effector hotspots from the conserved interdomain regions of the S2 subunit with FPPR (residues 850–855) and RBD regulatory centers (D427/D428) may be central to the coordination of allosteric interactions in the trimer, so that dynamic changes in these regions can determine the extent and directionality of the RBD motions. These results also provide an additional insight into the functional role of the intermonomer interactions of D427/D428 with K986 and K854 contacts with D614 to promote the stability of the closed state (Figure ). The results are consistent with protein engineering studies,[29−32] suggesting that these predicted effector centers may function as regulatory switch points of the dynamic equilibrium in the SARS-CoV-2 spike trimer. The structural rigidity and the strategic position of the effector centers as mediating hubs in the interaction networks suggest that mutations or ligand-based targeting of these positions can affect the dynamic equilibrium and allosteric conformational transitions between the functional states of the SARS-CoV-2 spike protein, which may be relevant for the discovery of allosteric modulators.

Sensor Profiling of SARS-CoV-2 Spike Trimer Structures: Signal Transmitters of Allosteric Changes Are Localized in the RBD Regions and Expose the Vulnerabilities of Glycosylation Shielding

The computed sensor profiles quantified the propensity of residues to serve as carriers or transmitters of long-range allosteric interactions, revealing the key regions that participate in this function. The major clusters of the sensor density in the locked form of the trimer are aligned with a small cluster of RBM residues 478–491 that are involved in the recognition of the host receptor ACE2 (Figure A). According to the results, the allosteric control over RBM residues may be exerted through regulatory centers D427/D428 that are allosterically coupled to the FPPR effector sites. In the locked state, the RBD sensor residues maintain their down-regulated positions, in which the RBD region is shielded from access to the host receptor (Figure A,D). In addition, the sensor sites are also consolidated near the FP motif (residues 812–829), with the peaks aligned with residues L821, L822, and F823 (L803, L804, and F805 in SARS-CoV) (Figure A, Table S1). Mutagenesis studies and infectivity assays showed the critical role of these residues in membrane fusion, as alanine mutants in these positions produced only limited fusion events in comparison with wild-type trypsin-treated cells.[116] Only small groups of RBM and FP residues form two clusters of sensor residues in the locked state of the trimer. Interestingly, these clusters of sensor residues are immediately proximal to their respective effector controllers (D427/D428 in the RBD and 850–855 residues of the FPPR).

Figure 9

PRS sensor profiles in the locked, closed, and open states of the SARS-CoV-2 spike trimer prefusion form. (A) Sensor profile for the SARS-CoV-2 S trimer in the locked state (PDB ID 6ZOZ).[32] (B) Sensor profile for the SARS-CoV-2 S trimer in the closed state (PDB ID 6VXX).[25] (C) Sensor profile for the SARS-CoV-2 S trimer in the open state (PDB ID 6VSB).[25] The profiles for protomer chains A, B, and C are shown as green, red, and blue lines, respectively. Structural maps of the PRS sensor profiles for the locked state of the SARS-CoV-2 S prefusion trimer (D), closed state (E), and open state (F). The color gradient from blue to red indicates the increasing sensor (receiver) propensities. The clusters of residues with the high sensor potential corresponding to the peaks of the effector profile are shown as red spheres.

It can be suggested that in the locked state these effectors could mediate the local mobility in the narrowly localized sensor regions, as allosteric signaling dictates only minor conformational changes of the rigid trimer form (Figure A,D). The sensor profiles were also computed for other SARS-CoV-2 trimer structures in the locked state, including the full-length, wild-type SARS-CoV-2 spike trimer with the structurally determined FPPR region (PDB ID 6XR8) (Figure S5). This comparison demonstrated the strong similarity of the sensor distributions in all locked states, where most sensor residues are localized in the flexible NTD regions. Hence, allosteric communications in the stable locked state may be primarily directed from the dominant effector region in the S2 subunit through the effector switch in the RBD (D427/D428) to the flexible NTD residues acting as primary receivers of the allosteric signal and enabling thermal breezing without altering the structural stability of the locked state.

A notable change can be observed in the closed form, where the major cluster of sensor sites is now localized in the flexible RBD regions (Figure B,E), and the second cluster in the RBD region becomes much broader and covers a wide range of exposed residues in the RBD region (residues 350–505). As a result of this expansion in the distribution of sensor sites, the entire RBD region becomes susceptible to large conformational changes that are governed by the effector centers in the conserved regions of the S2 subunit (Figure B, Table S1). In the closed form, the furin cleavage site (residues 682–685) also overlapped with the increased density of the sensor sites. Together, the effector and sensor sites can form a dynamic allosteric network that drives signal transmission in the prefusion trimer. The sensor profile of the open form of the trimer featured a different allocation of sensor sites in the monomers (Figure C,F). Most of the sensor sites in the monomer undergoing “up” transition are consolidated in the flexible RBD region (Figure C). The entire RBM motif in the “up” monomer becomes enriched with sensor sites transmitting the allosteric signal in the open trimer. At the same time, the sensor profiles of the other two monomers revealed that the density localized mainly in the NTD regions.

The PRS analysis showed that sensor positions often overlap with evolutionary variable sites in the flexible RBD region that undergo mutations between SARS-CoV and SARS-CoV-2 glycoproteins. The partitioning between effector and sensor sites in the SARS-CoV-2 S trimer also reflects the evolutionary diversity elevated at the NTD and RBD regions and the prevalent conservation in the regions in the S2 domain, including the FP, FPPR, HR1, and the CH domains, which are subject to greater functional constraints. These findings are consistent with the properties of sensor sites in protein systems that are often tolerant of mutational changes and tend to be located in the mobile loops and exposed binding interfaces.[117] Computational studies of molecular chaperones conducted in my laboratory[89] and other groups[83−87] demonstrated that sensors are frequently found close to recognition sites and protein regions undergoing allosteric changes to adapt to binding and to accommodate diverse interactions.

Mechanistic details of the binding mechanism with the host receptors result from a complex interplay of thermodynamic forces determined by the intrinsic dynamics of the spike trimer coupled to the variability and adaptation of glycosylation shields.[42] To access the role of glycosylation sites in shielding sensor/receiver sites in the trimer, the results of the site-specific glycan analysis of the SARS-CoV-2 spike protein were utilized to project the determined 22 glycosylation sites on structural maps of allosteric sensor profiles of the locked closed, closed, and open forms (Figure S6). This analysis showed a considerable overlap of the glycosylation sites dispersed across both the S1 and S2 subunits, with the regions of increased sensor density in the locked trimer state (Figure S6). Strikingly, some of the important glycosylation sites N122, N149, N74, N165, N234, and N343 correspond precisely to the distinct and sharp peaks of the sensor profile in the locked state. The analysis indicates that N-linked glycans can shield the most sensor-sensitive regions of the SARS-CoV-2 in the NTD regions of the locked state (Figure S6). In the closed flexible state, an increasing density of sensor residues in the RBD regions that are not properly shielded may create vulnerabilities, leading to the reduced protection against conformational changes and ultimately promoting the transition to the open form. These findings are consistent with the pioneering studies by Amaro and colleagues,[42] who performed a series of unprecedented atomistic MD simulations of the SARS-CoV-2 spike trimers in a fully glycosylated environment, showing that the N-glycans can function to protect the glycoprotein but also to allow for a coordinated response to allosteric signaling and to promote the “down–up” RBD transition through observed shield vulnerabilities in the closed state.

Hierarchical Analysis of the Residue Interaction Networks: Betweenness Centrality and Community Centrality Profiling Unveil the Mediators of Allosteric Interactions

Using a graph-based representation of protein structures,[92−94] the residue interaction networks were constructed, in which the inter-residue edges were weighted using the residue cross-correlations obtained from simulations[94] and coevolutionary residue correlations.[97] Using the ensemble-averaged model of the residue interaction networks, residue betweenness (or residue centrality) is computed and is used to identify the key mediating centers of the allosteric interaction networks. The high-centrality positions in the interaction networks were aligned with many effector centers in the PRS profiles and have high allosteric propensities (Figure ). The important result of this analysis is the emergence of several clusters of high-centrality sites in the locked prefusion state that mediate allosteric interactions in the global interaction network (Figure A). The first dense cluster of mediating residues is located in the RBD regions, with major centrality peaks corresponding to residues Y369, L387, and V382. The second cluster corresponds to residues in the CTD1 (residues 529–591), which is believed to function as an allosteric relay between the RBD and FPPR regions by communicating the signal from and to the FP. Another cluster was found near the HR1 motif and CH regions (Figures A and 11A). In the prefusion state, the HR1 (residues 910–985) and CH regions (residues 986–1035) form separate helices, and the HR1 helices and the CH helix are arranged in an antiparallel orientation. The high-centrality sites in the HR1 region are aligned with residues I980, R983, L984, K/P986, V/P987, and R1000 (Figure B). Interestingly, some of these sites act as regulatory switches of the dynamic equilibrium of the SARS-CoV-2 spike trimer (K986, V987), as proline mutations of these residues can dramatically reduce the fraction of the closed state and favor the open state of the trimer.[32] A high density of high-centrality sites in this important region of the S2 subunit is present and is shared in distinct conformational states of the SARS-CoV-2 spike protein (Figure A–C). Notably, the number of high-centrality peaks was reduced in the more dynamic closed form, particularly indicating the reduced density of global mediating centers in the RBD regions, which may be a consequence of the increased flexibility of the RBD regions (Figure B). These results are consistent with the PRS analysis in which the RBD regions in the dynamic closed form are mainly aligned with sensors of allosteric signals. In network terms, this implies that the allosteric interaction network in the closed form may be weakened and become more diffuse, reflecting the reduced stability and showing signs of the intermediate conformation preceding the transition to the open form. The distribution of the high-centrality sites was further altered in the open state, where the high density of mediating centers was detected near the RBD regions in the monomer undergoing “up” transition (Figure C).

Figure 10

Residue-based betweenness centrality profiles in the locked, closed, and open states of the SARS-CoV-2 spike trimer prefusion form. The residue-based centrality values are computed by averaging the results over 1000 representative equilibrium samples of the SARS-CoV-2 spike trimer structure from CABS-CG simulations and atomistic reconstruction of trajectories. (A) Betweenness centrality profile for the SARS-CoV-2 S trimer in the locked state (PDB ID 6ZOZ).[32] (B) Betweenness centrality profile for the SARS-CoV-2 S trimer in the closed state (PDB ID 6VXX).[25] (C) Betweenness centrality profile for the SARS-CoV-2 S trimer in the open state (PDB ID 6VYB). The profiles for protomer chains A, B, and C are shown as green, red, and blue lines, respectively.

Figure 11

Structural map of high-centrality clusters for the locked state of the SARS-CoV-2 S prefusion trimer. (A) Projection of high-centrality clusters (shown in red spheres) onto the cryo-EM structure of SARS-CoV-2 spike trimer in the locked state (PDB ID 6ZOZ).[32] The S trimer structure is in ribbon representation with reduced transparency. The protomers A, B, and C are green, red, and blue, respectively. (B) Close-up of the S2 subunit with high-centrality sites (in red spheres) mapped in the HR1 (residues 910–985) and CH regions (residues 986–1035). The NTD is green, RBD is black, HR1 motif is cyan, CH is red, FPPR is blue, and upstream helices (UH) motif is magenta.

Importantly, the cluster of high-centrality positions in HR1/CH regions retained its mediating role in the closed form of the trimer, suggesting that these regions play a key role during allosteric conformational changes from the locked to the open state. It is worth noting that global rearrangements in the HR1/CH regions harboring the proposed regulatory sites act as a nucleus during the formation of the postfusion state.[33] The important result of this analysis was a revelation that key regulatory sites that could dictate dynamic switching between conformational states of the SARS-CoV-2 spike trimer overlap with the high-centrality mediators of the interaction network that are also predicted to be effector centers with high allosteric potential (Figures and 11). To further examine the network properties, the hierarchical multilayer community decomposition was performed, and the community centrality metric[118] of the residue interaction networks was computed for distinct functional states of the SARS-CoV-2 spike trimer (Figure ). The community centrality profile revealed that network hubs that could bridge local modules in the locked state are located at the borders of the RBD regions (Figures A and 13A). The denser distribution of bridging centers in a stable locked state of the SARS-CoV-2 spike trimer (Figure A) is indicative of a broad and stable allosteric interaction network. The community centrality hubs in the allosteric network are located near the experimentally known regulatory switches of dynamic equilibrium in the SARS-CoV-2 spike trimer. The structural mapping of the community centrality sites in the locked state (Figure A) highlighted a network of switch points that could enable allosteric couplings between spatially separated communities in the interaction network. The results also indicated that the number of local communities and community centrality sites can be reduced in the more dynamic intermediate closed state (Figure B) and open state of the prefusion trimer (Figure C).

Figure 12

Community analysis and community centrality profiles in the locked, closed, and open states of the SARS-CoV-2 spike trimer prefusion form. The community centrality values are computed using community decomposition of the residue interaction networks and by averaging the results over 1000 representative equilibrium samples of the SARS-CoV-2 spike trimer. (A) Community centrality profile for the SARS-CoV-2 S trimer in the locked state (PDB ID 6ZOZ).[32] (B) Community centrality profile for the SARS-CoV-2 S trimer in the closed state (PDB ID 6VXX).[25] (C) Community centrality profile for the SARS-CoV-2 S trimer in the open state (PDB ID 6VYB). The profiles for protomer chains A, B, and C are shown as green, red, and blue lines, respectively.

Figure 13

Structural map of community centrality clusters for the locked state of the SARS-CoV-2 S prefusion trimer. (A) Projection of community centrality clusters onto the cryo-EM structure of SARS-CoV-2 spike trimer in the locked state (PDB ID 6ZOZ).[32] (B) Mapping of community centrality clusters onto the SARS-CoV-2 S trimer in the closed state (PDB ID 6VXX).[25] (C) Projection of community centrality clusters on the SARS-CoV-2 S trimer in the open state (PDB ID 6VYB). The profiles for protomer chains A, B, and C are shown as green, red, and blue lines, respectively. The S trimer structure is in ribbon representation with reduced transparency. The community centrality clusters are depicted as spheres colored by their respective protomers.

Hence, the hierarchical community analysis confirmed a distinct organization of the interaction networks in functionally different locked and closed prefusion states. The regulatory sites that could dictate dynamic switching between the conformational states of the SARS-CoV-2 spike trimer overlap with the high-centrality sites and effector hotspots with high allosteric potential.

Conclusions

This study systematically examined the functional mechanisms of the SARS-CoV-2 spike prefusion trimer through the lens of allosteric regulation and signal transmission. Functional dynamics analysis and PRS characterized the hinges of conformational transitions and identified the regulatory hotspots that control the signal transmission and allosteric interactions in distinct states of the SARS-CoV-2 spike trimer. The results revealed the migration of the hinge sites during the transition from the stable locked to the more dynamic closed trimer state, leading to the release of dynamic breaks on the RBD motions. Through perturbation-based network analysis, it was demonstrated that the stabilized locked form of the prefusion trimer is characterized by a broad allosteric network with key regulatory centers located in the RBD, FPPR, and HR regions that dictate the dynamic switching between conformational states of the SARS-CoV-2 spike trimer. The PRS analysis identified key regulatory hotspots in the SARS-CoV-2 spike, showing that these allosteric centers are precisely aligned with known regulatory switches in the RBD and FPPR regions. The results are consistent with protein engineering studies and suggest that the predicted effector centers may function as regulatory switch points of the dynamic equilibrium in the SARS-CoV-2 spike trimer. The discovered effector centers can alter the distribution of functional states and conformational transitions between the closed and the host-accessible open states of the SARS-CoV-2 spike trimer, thus presenting opportunities for rational allosteric drug targeting of the SARS-CoV-2 spike proteins.