key: cord-0773338-kmcnmb52 authors: Verkhivker, Gennady M.; Di Paola, Luisa title: 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 date: 2021-01-15 journal: J Phys Chem B DOI: 10.1021/acs.jpcb.0c10637 sha: ef7b3555639aac8f9a04cf075ce264d2c9181a7f doc_id: 773338 cord_uid: kmcnmb52 [Image: see text] The rapidly growing body of structural and biochemical studies of the SARS-CoV-2 spike glycoprotein has revealed a variety of distinct functional states with radically different arrangements of the receptor-binding domain, highlighting a remarkable function-driven conformational plasticity and adaptability of the spike proteins. In this study, we examined molecular mechanisms underlying conformational and dynamic changes in the SARS-CoV-2 spike mutant trimers through the lens of dynamic analysis of allosteric interaction networks and atomistic modeling of signal transmission. Using an integrated approach that combined coarse-grained molecular simulations, protein stability analysis, and perturbation-based modeling of residue interaction networks, we examined how mutations in the regulatory regions of the SARS-CoV-2 spike protein can differentially affect dynamics and allosteric signaling in distinct functional states. The results of this study revealed key functional regions and regulatory centers that govern collective dynamics, allosteric interactions, and control signal transmission in the SARS-CoV-2 spike proteins. We found that the experimentally confirmed regulatory hotspots that dictate dynamic switching between conformational states of the SARS-CoV-2 spike protein correspond to the key hinge sites and global mediating centers of the allosteric interaction networks. The results of this study provide a novel insight into allosteric regulatory mechanisms of SARS-CoV-2 spike proteins showing that mutations at the key regulatory positions can differentially modulate distribution of states and determine topography of signal communication pathways operating through state-specific cascades of control switch points. This analysis provides a plausible strategy for allosteric probing of the conformational equilibrium and therapeutic intervention by targeting specific hotspots of allosteric interactions and communications in the SARS-CoV-2 spike proteins. The coronavirus disease 2019 (COVID-19) pandemic associated with the severe acute respiratory syndrome (SARS) has emerged as a global international health crisis spread over the world with serious challenges for global economy, peace, and security. 1, 2 Genomic studies have shown that sequences of the novel coronavirus SARS-CoV-2 virus displayed high sequence similarity with SARS and Middle East respiratory syndrome (MERS) viral proteins involved in the replication cycle. 3−6 SARS-CoV-2 infection is transmitted when one of the main structural proteins, the viral spike (S) glycoprotein, 7−9 binds to the host cell receptor angiotensinconverting enzyme 2 (ACE2) 10−12 followed by entry of the spike protein into host cells via membrane fusion. SARS spike (S) proteins are trimeric fusion proteins, with two main domains, an amino (N)-terminal S1 subunit and a carboxyl (C)-terminal S2 subunit. The full-length SARS-CoV-2 S protein consists of two main domains, an amino (N)-terminal S1 subunit and a carboxyl (C)-terminal S2 subunit. The subunit S1 is involved in interactions with the host receptor and includes an N-terminal domain (NTD), the receptorbinding domain, RBD), and two structurally conserved subdomains (SD1 and SD2). Structural and biochemical studies of SARS coronavirus established that the key interactions between SARS-CoV/SARS-CoV-2 proteins and the host receptor ACE2 enzyme are mediated by the receptor binding motif (RBM) of the RBD, suggesting that these interaction mechanisms may regulate both the cross-species and human-to-human infection transmissions. 14, 15 The S1 subunit undergoes spontaneous conformational transformations between ensembles of closed and open receptoraccessible forms, switching between a standing up RBD position amenable for ACE2 binding and a closed-down position eliciting immune evasion. The S2 subunit is highly evolutionary conserved fusion machinery that contains an Nterminal hydrophobic fusion peptide (FP), fusion peptide proximal region (FPPR), heptad repeat 1 (HR1), central helix region (CH), connector domain (CD), heptad repeat 2 (HR2), transmembrane domain (TM), and cytoplasmic tail (CT). The S1 domains are situated above the S2 subunit, covering and protecting the fusion apparatus. Upon proteolytic activation at the S1/S2 boundary and transition to the fully open form of the S trimer, the S1 subunit can release the dynamic breaks, promote dissociation from S2 followed by a cascade of massive structural rearrangements of the S2 subunit to mediate the fusion of the viral and cellular membranes and host cell takeover. 13−16 The cryo-EM structures of MERS-CoV and SARS-CoV spike glycoproteins revealed that S protein trimers can exist in a dynamic equilibrium between a closed ("RBD-down") prefusion conformation that undergoes large structural rearrangements and hinge-based functional motions to adopt a receptor-accessible open ("RBD-up") state in which the S protein can fuse the viral membrane. 17−26 The early studies showed that engineered coronavirus spikes through S-2P mutations (K986P and V987P in SARS-CoV-2) are resistant to conformational changes induced by receptor recognition and can prevent the transition from the prefusion to postfusion states. 27 These early studies of SARS-CoV S proteins showed that the majority of the S-2P proteins were in the single RBDup conformation (58%), with double-up and triple-up S conformations populating a smaller fraction of the equilibrium (39% and 3% respectively), indicating that S-2P mutations can destabilize the closed form but still maintain stability of the prefusion trimer. The structures of the full-length S protein in the prefusion and postfusion conformations provided an unprecedented look and insight into tectonic rearrangements during transition to the postfusion state. 28−30 Using mass spectrometry and single-particle EM, the molecular events associated with refolding of the prefusion S glycoprotein to the postfusion form were mapped, revealing a strong similarity between postfusion organization of coronaviruses S and a conserved mechanism of mediating membrane fusion entry of these viruses. 28 The cryo-EM structure of the postfusion SARS-CoV trimer showed dramatic conformational changes of the SARS-CoV S2 machinery during membrane fusion, where the HR1 and the CH helical segments, which are arranged into an antiparallel orientation in the prefusion state, can form a long continuous stem helix pointing toward the target membrane in the postfusion trimer. 29 The rapidly growing body of cryo-EM structures of the SARS-CoV-2 S protein in distinct functional conformations detailed structural changes in the trimer that are manifested in large movements of the RBD regions, highlighting a remarkable function-driven conformational plasticity and adaptability of the spike proteins. 31−34 The recent cryo-EM characterization of the SARS-CoV-2 S trimer demonstrated that S protein may not only fluctuate between the dynamic open and closed states but also exhibit subtle dynamic variability by adopting structurally rigid locked-closed form that revealed the ordered FPPR segment and cross-talk between stabilized RBD and FPPR regions. 34 Structural and biophysical studies employed protein engineering to generate prefusion-stabilized SARS-CoV-2 S variants by introducing disulfide bonds and proline mutations to modulate the stability of the S2 subunit and intersubunit boundaries and consequently prevent refolding changes that accompany acquisition of the postfusion spike state. 35 Cryo-EM structures of SARS-CoV-2 S variant with several proline substitutions (F817P, A892P, A899P, A942P) featured two distinct conformationsone with a single RBD in the up conformation and the other with two RBDs in the up conformation. These results suggested that modulating stability of S2 domains in the SARS-CoV-2 spike trimer may have a strong allosteric effect on dynamic couplings between S1 and S2 regions, altering the equilibrium and leading to a significant population of partially open and open states. 35 By combining targeted mutagenesis and cryo-EM structure determination, a recent biophysical study revealed a differential stabilization of the mutant constructs in the closed and open spike forms. 36 While modifications of contact regions between the RBD and S2 domains via S383C/D985C mutations leads to the thermodynamic prevalence of the closed-down conformation, the quadruple mutant (A570L/T572I/F855Y/ N856I) introducing modifications at the interprotomer contacts between SD1 and S2 can shift the equilibrium away from the closed-down state, yielding a significant fraction of spike conformations in the open receptor-accessible form and showing the enhanced binding with ACE2 receptor. 36 This pioneering study indicated that the key interdomain contacts may act as potential global dynamic switches and constitute regulatory centers responsible for modulation of large scale changes and the equilibrium partitioning between the closed and open states. Protein engineering and cryo-EM studies of a prefusion-stabilized SARS-CoV-2 S ectodomain trimer using the interprotomer disulfide bonds (S383C/D985C, G413C/ P987C, T385C/T415C) between RBD and S2 regions can lock the trimer in the closed state enhancing the SARS-CoV-2 S resistance to proteolysis. 37 Targeted design and structure determination of thermostable SARS-CoV-2 spike trimers unmasked the heterogeneity of the SARS-CoV-2 landscape showing that disulfide-bonded S-protein trimer variants imposing stabilization in strategically located interprotomer positions (S383-D385) and (G413-V987) can lead to a predominant population of the prefusion closed states with only ∼20% population sampling the open state. 38 In contrast, proline substitutions K986P and V987P switch this equilibrium toward the receptor-accessible open form with only ∼20% of the trimers in the closed state and ∼80% of molecules adopting the open form. 38 SARS-CoV-2 S mutants with the enhanced infectivity profile 39−41 have been a primary focus of the recent structural and biochemical studies, revealing that D614G mutation can radically shift the population of the SARS-CoV-2 S trimer toward open states. 42 The cryo-EM and tomography tools were employed to image the intact SARS-CoV-2 virions and determine the highresolution structure, conformational flexibility and distribution of S trimers in situ on the virion surface. 43 These studies confirmed that the opening of the RBD observed in recombinant S trimers can also take place on the virus surface, and that the structure of SARS-CoV-2 S trimers on the virion surface is very similar to the closed prefusion form of the trimeric ectodomain stabilized by double proline mutations. Importantly, these studies confirmed that a general mechanism of spontaneous conformational changes and population shifts between different functional states can be maintained in different biological environments, reflecting the intrinsic properties of conformational landscapes for SARS-CoV-2 S trimers and suggesting that RBD epitopes become stochastically exposed to interactions with ACE2 and with antibodies. 43 Single-molecule fluorescence (Forster) resonance energy transfer (smFRET) studies of SARS-CoV-2 S trimer on virus particles revealed at least four distinct conformational states for the SARS-CoV-2 S trimer, also detailing a sequence of conformational transitions from the closed state to the receptor-bound open state through one necessary intermediate in which all three RBD domains in the closed conformations are oriented toward the viral particle membrane. 44 This study showed that the host receptor ACE2 activates spike proteins on the virus by shaping the conformational landscape toward stabilizing the "RBD-up" conformation, suggesting that mechanisms of conformational selection and receptor-induced structural adaptation can be intertwined and showing that allosteric stabilization of the closed S protein by antibodies can be effective strategy to achieve SARS-CoV-2 neutralization 44 Cryo-EM structural studies also mapped a mechanism of conformational events associated with ACE2 binding, showing the compact closed form of the SARS-CoV-2 S protein becomes weakened after furin cleavage between the S1 and S2 domains, leading to the increased population of partially open states, while ACE2 binding to an open RBD can further accelerate transformation to a fully open and ACE2-bound form, priming the protein for fusion activation. 45 Glycosylation of the SARS-CoV-2 proteins and shielding of the receptor binding sites by glycans is a common mechanism of viral glycoproteins for evading the immune system and recent studies showed that glycosylation on SARS-CoV proteins can camouflage immunogenic protein epitopes. 46−48 The site-specific analyses of N-linked glycosylation on soluble SARS, MERS, and SARS-CoV S glycoproteins revealed an extensive heterogeneity and identified hotspot clusters suggesting that N-linked glycan modifications of SARS-CoV S proteins may not constitute an effective shield and exhibit numerous vulnerabilities in comparison to the glycan shields of other viruses. 47 This study highlighted that boundaries between buried and exposed surface regions may often define potential sites of glycan shield vulnerability. 47 Site-specific mass spectrometric mapping of the glycan-processing states and cryo-EM analysis of the SARS-CoV-2 S protein showed that SARS-COV-2 is highly glycosylated and established the composition and functional role of 22 N-linked glycosylation sites and 4 O-linked glycosylation sites of S protein that occlude distinct regions on the protein surface. 48 This study confirmed that the glycan shield of SARS-CoV-2 S proteins may exhibit various vulnerabilities in the RBD regions, particularly highlighting the role of proximal glycosylation sites N165, N234, and N343 in shielding of the receptor binding sites especially when the receptor binding domain is in the closed-down conformation. 48 Using a combination of antigenic screens and high-resolution cryo-EM structure determination, it was shown that an N-glycan deletion at N234 results in a dramatically reduced population of the SARS-CoV-2 S trimer with "up" state RBDs, while glycan deletion at position N165 results in a the increase of the open states. 49 This study demonstrated that these glycans can act not only as a elements of the protective shield to antibodymeditated immunity but are also capable of modulating the conformational landscape of the SARS-CoV-2 S protein, thereby providing another layer of regulatory switch points to control shifts between different functional states. Computational modeling and molecular dynamics (MD) simulations have been instrumental in predicting shapes and motions of glycans for glycoproteins. 50−52 The developed structural models of glycoforms of the SARS-CoV-2 spike protein and MD simulations probed the effect of glycan heterogeneity on the antigenicity of the S glycoprotein showing that different glycan forms can shield ∼40% of the protein surface on the S glycoprotein from antibody recognition, while the most accessible antigenic surface on the S protein and regions of vulnerabilities correspond to the ACE2 binding domain, where glycan shielding or mutational changes are unable to evade host immune response. 53, 54 The development of a fully glycosylated full-length SARS-CoV-2 S protein in a viral membrane enabled atomistic simulations of the SARS-CoV-2 S trimer structures in a properly glycosylated environment. 55 MD simulations of the SARS-CoV-2 S glycoprotein embedded in the viral membrane with a complete glycosylation profile reported a great level of details about conformational dynamics of the open and closed structures, suggesting that the shape and dynamics of the glycan shield can be coupled to allosteric conformational changes of the S protein and provide an additional layer of molecular control over spike response to the host receptor. 56 The development of a "bottom-up" coarse-grained model of the SARS-CoV-2 virion was recently presented in which cryo-EM, X-ray crystallography, and computational predictions were integrated to build a multiscale model of structural SARS-CoV-2 proteins and a complete virion model. 57, 58 Another large-scale MD study of a 4.1 million atom system containing a patch of viral membrane with four full-length, fully glycosylated S proteins was reported, where structural dynamics and the analysis of steric accessibility were used for epitope mapping, revealing signatures of the flexible glycan coat that can shield a larger surface area than can be derived from static structures. 59 MD simulations of the SARS-CoV-2 spike glycoprotein identified differences in flexibility of functional regions that may be important for modulating the equilibrium changes and binding to the ACE2 host receptor. 60, 61 Our recent study examined the molecular basis of the SARS-CoV-2 binding with the ACE2 enzyme through the lens of coevolution and conformational dynamics that conspire to drive cooperative binding interactions and signal transmission. 62 Using protein contact networks and perturbation response scanning based on elastic network models, we recently presented an integrated approach for detecting allosteric sites on the SARS-CoV-2 spike protein, showing convergence to a single allosteric region. 63 Molecular simulations and network modeling approaches were also used on our most recent investigation to present evidence that the SARS-CoV-2 spike protein can function as an allosteric The Journal of Physical Chemistry B pubs.acs.org/JPCB Article regulatory engine that fluctuates between dynamically distinct functional states. 64 Allosteric molecular events can involve complex cascades of thermodynamic and dynamic changes that occur on different spatial and temporal scales. The thermodynamic-centric energy landscape concepts and conformational selection models of allosteric regulation have gained a considerable prominence based on the notion that statistical ensembles of preexisting conformational states and communication pathways are intrinsic to a given protein and provide molecular machinery for allosteric modulation of the underlying landscapes induced by external perturbations, ligand binding, and mutations. 65−69 These models can also link the allosteric mechanisms of proteins to their functional role in regulatory processes and signaling events. 69 In this study, we expanded our comparative analysis of mechanisms of allosteric signal transmission for the SARS-CoV-2 S states in the closed-a down form, partially open (1up) and open (2-up) conformations determined for a double cysteine mutant S383C/D985C (rS 2d, RBD to S2 double mutant) quadruple mutant A570L T572I F855Y N856I (u1S2q, subdomain 1 to S2 quadruple mutant). 36 Using an integrated approach that combined coarse-grained (CG) simulations, protein stability analysis, and multifaceted dynamic modeling of residue interaction networks, we examined how mutations can differentially affect dynamics, stability and allosteric signaling in distinct functional states if the SARS-CoV-2 S protein. A detailed atomistic network analysis and modeling of allosteric communications in the SARS-CoV-2 S trimer mutants revealed key regulatory centers and ensemble of pathways that mediate allosteric interactions and dictate conformational changes, thereby governing mechanisms of signaling and binding. The results of this study provide a novel insight into the allosteric regulatory mechanisms of SARS-CoV-2 S proteins showing that mutations at the key regulatory positions can differentially modulate distribution of states and shape up signal transmission pathways using cascades of control switch points. This analysis provides a plausible strategy for therapeutic intervention and allosterically driven modulation of the conformational equilibrium by targeting specific hotspots of allosteric interactions and communications in the SARS-CoV-2 proteins. Coarse-grained (CG) models are computationally effective approaches that leverage topology-based framework for characterization of protein structure and dynamics, enabling simulations of large systems over long time scales. In this study, CG-CABS model was used for simulations of the cryo-EM structures of the SARS-CoV-2 spike trimer mutants in distinct conformational states. In this high-resolution model, the amino acid residues are represented by Cα, Cβ, the center of mass of side chains and another pseudoatom placed in the center of the Cα−Cα pseudobond. The sampling scheme of the CABS model used in our study is based on Monte Carlo replicaexchange dynamics and involves a sequence of local moves of individual amino acids in the protein structure as well as moves of small fragments. 70−72 the solvent effect is incorporated implicitly through statistical potentials of the CABS force-field. We employed CABS-flex approach that efficiently combines a high-resolution coarse-grained model and efficient search protocol capable of accurately reproducing all-atom MD simulation trajectories and dynamic profiles of large biomolecules on a long time scale. 70−74 CABS-flex standalone package dynamics implemented as a Python 2.7 object-oriented package was used for fast simulations of protein structures. 74−76 This approach also allowed for MODELER-based reconstruction of generated models and simulation trajectories to all-atom representation. 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. 75, 76 The loop regions are completely flexible in CABS-CG simulations and regions of secondary structure are restrained only locally by applying soft distance restraints to residues if their distance deviates from the distance in the crystal structure by more than 1 Å. 75, 76 A total of 1000 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 structures of the SARS-CoV-2 spike (S) protein trimer were used in simulations and modeling (Table 1, Figure 1 ): the cryo-EM structure of the full-length S protein mutant rS 2d S383C/D985C in the closed prefusion conformation (pdb id 6X29), 36 the cryo-EM structure of the full-length S u1S2q quadruple mutant (A507L T572I F855Y N856I) in the all-down closed prefusion conformation (pdb id 6X29), 36 the cryo-EM structure of the full-length S u1S2q quadruple mutant (A507L T572I F855Y N856I) in the partially open (1-up) conformation (pdb id 6X2A), 36 and the cryo-EM structure of the full-length S u1S2q quadruple mutant (A507L T572I F855Y N856I) in the open (2-up) conformation (pdb id 6X2B). 36 All structures were obtained from the Protein Data Bank. 77, 78 Structure preparation process included checking and correcting for errors, adding missing residues and atoms, filling valences with hydrogens. First, protein residues in all studied structures were systematically inspected for missing residues and protons. Hydrogen atoms and missing residues were initially added and assigned according to the WHATIF program web interface. 79, 80 The structures were then preprocessed through the Protein Preparation Wizard (Schrodinger, LLC, New York, NY) and included a systematic checking of bond order, assignment, and adjustment of ionization states, checking of disulfide bonds, removal of crystallographic water molecules and cofactors, capping of the termini, assignment of partial charges, and addition of missing atoms and side chains that were not assigned in the initial processing with the WHATIF program. The cryo-EM structures of the SARS-CoV-2 S proteins contained a number 81 and ArchPRED server 82 and further confirmed by the FALC (fragment assembly and loop closure) program. 83 The side chain rotamers were refined and optimized by the SCWRL4 tool. 84 The CG conformational ensembles were also subjected to all-atom reconstruction using CABS-flex standalone package dynamics and checked using PULCHRA method 85 and CG2AA tool 86 to produce atomistic models of simulation trajectories. The protein structures were then optimized using atomic-level energy minimization with a composite physics and knowledge-based force fields as implemented in the 3Drefine method. 87 In all-atom MD simulations, the glycans are typically added to the protein structure and a complete model of the glycoprotein is then subjected to simulations. 56−59 In our study, we employed CABS-CG model for dynamics 36 The structure is shown in ribbons with protomers colored in green, cyan and blue, respectively. The down positions of the RBDs are annotated. A close-up of the mutational site S383C/D985C at the RBD to S2 boundaries. C383 and C985 are shown in sticks and colored according to their respective protomer. (B) Cryo-EM structure of the full-length S u1S2q quadruple mutant (A507L T572I F855Y N856I) in the all-down closed prefusion conformation (pdb id 6X29). 36 The protomers are colored as in part A. A close-up of the mutational site at the boundary between subdomain 1 to S2 is shown with L570, I572, Y855, and I856 residues at the interprotomer interfaces. The residues are in sticks and colored according to their protomer color scheme. (C) Cryo-EM structure of the full-length S u1S2q quadruple mutant (A507L T572I F855Y N856I) in the partially open (1-up) prefusion conformation (pdb id 6X2A). 36 (D) The cryo-EM structure of the full-length S u1S2q quadruple mutant (A507L T572I F855Y N856I) in the open (2-up) conformation (pdb id 6X2B). 36 A close-up of the mutational site at the boundary between subdomain 1 to S2 is shown as in part B. The Journal of Physical Chemistry B pubs.acs.org/JPCB Article simulations followed by all-atom reconstruction of trajectories. In addition, the reconstructed atomistic structures from simulation trajectories were decorated by N-acetyl glycosamine (NAG) glycan residues and optimized. The structure of glycans at particular glycosites of the closed and open states of SARS-CoV-2 S protein was previously determined and these glycans were incorporated in atomistic modeling of the SARS-CoV-2 mutant structures. The glycosylated microenvironment for atomistic models of the simulation trajectories was mimicked by using the structurally resolved glycan conformations for 16 out of 22 most occupied N-glycans in each protomer (N122, N165, N234, N282, N331, N343, N603, N616, N657, N709, N717, N801, N1074, N1098, N1134, and N1158) as determined in the cryo-EM structures of the SARS-CoV-2 spike S trimer in the closed state (K986P/V987P,) (pdb id 6VXX) and open state (pdb id 6VYB), 31 and the cryo-EM structure SARS-CoV-2 spike trimer (K986P/V987P) in the open state (pdb id 6VSB). 32 The glycan-decorated atomistic models were subsequently used for the ensemblebased structural, energetic and network analyses of the SARS-CoV-2 mutants examined in this work. We performed principal component analysis (PCA) of reconstructed trajectories derived from CABS-CG simulations using the CARMA package 88 and also determined the essential dynamics profiles in slow modes using elastic network models (ENM) analysis. 89 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. 71, 72 Two elastic network models: Gaussian network model (GNM) 90, 91 and anisotropic network model (ANM) approaches 92 were adopted to compute the amplitudes of isotropic thermal motions and directionality of anisotropic motions. The functional dynamics analysis was conducted using the GNM in which protein structure is reduced to a network of N residue nodes identified by C α atoms, and the fluctuations of each node are assumed to be isotropic and Gaussian. The topology of the protein structure is described by N × N Kirchhoff matrix of inter-residue contacts Γ, where the off-diagonal elements are −1, if the nodes are within a cutoff distance r c , and zero otherwise. In the GNM approach 90,91 the interaction potential for a protein of N residues was: where R ij 0 and R ij were the equilibrium and instantaneous distance between residues i and j, and Γ was the N × N Kirchhoff matrix: Thus, the square fluctuations were: The normal modes were extracted by eigenvalue decomposition of the matrix Γ=UΛU T , U being the orthogonal matrix whose k th column u k was the k th mode eigenvector. Λ was the diagonal matrix of eigenvalues λ k . In ANMapproach 92 the interaction potential for a protein of N residues was: The 3N × 3N Hessian matrix, H, determines the systems dynamics. where X i , Y i , and Z i represented the Cartesian components of the ith residue and V was the potential energy of the system. We selected r c = 13 Å. ANMs provide information about the amplitudes, as well as the direction of residue fluctuations. Conformational mobility profiles in the essential space of low frequency modes were obtained using DynOmics server 91 and ANM server. 92 Protein Stability Analysis. To compute protein stability changes in the SARS-CoV-2 trimer mutants, we conducted a systematic alanine scanning of protein residues in the SARS-CoV-2 trimer mutants. A BeAtMuSiC approach was employed that 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. 93− 95 We leveraged rapid calculations based on statistical potentials to compute the ensembleaveraged binding free energy changes using equilibrium samples from MD trajectories. The binding free energy changes were computed by averaging the results over 1000 equilibrium samples for each of the studied systems. Protein Contact Networks and Network Clustering. The protein contact network is a network whose nodes are the protein residues and links are active contacts between residues in the protein structure. The protein contact network is an undirected, unweighted graph; it is built on the basis of the distance matrix d, whose generic element d ij records the Euclidean distance between the ith and the jth residue (measured between the corresponding α carbons). A detailed description of network construction and significance of network descriptors is presented in our previous studies. 96 −98 The active network links are defined using a range of contacts between 4 and 8 Å. The description of the network is given by the following adjacency matrix: where d ij is the distance between the residues. The node degree describes the number of links each residue has with other residues, defined as We previously demonstrated that spectral network clustering targets functional modules in proteins. 88 The network clustering is based on the spectral decomposition of the network Laplacian L defined as where D is the degree matrix, a diagonal matrix whose diagonal is the degree vector, and A is the adjacency matrix. The eigenvector corresponding to the second minor eigenvalue of the Laplacian (Fiedler's vector) assigns nodes (residues) to two different clusters on the basis of the sign of its corresponding components. The method is then binary and hierarchical. Once the network is divided into a given number of clusters (powers of 2), we define the participation coefficient, defined as where k i is the overall node degree, while k si is the node degree including only links with nodes (residues) that belong to their own cluster. The participation coefficient P describes the propensity of residue nodes to participate into intercluster communication. Functional regions in proteins have been identified on the basis of values of non null participation coefficients, pointing to an active role of corresponding residues in signal transmission between different clusters (functional domains). We highlighted as highly active communication residues the nodes with P > 0.75, using our previous studies showing that such residues typically correspond to important regulatory centers of signal transmission between protein domains. 96 The proposed methodology of network clustering was implemented as a Cytoscape plugin. 99 Dynamic-Based Modeling of Residue Interaction Network and Community Analysis. A graph-based representation of protein structures 100−102 is used to represent residues as network nodes and the inter-residue edges to describe residue interactions. The details of graph construction using residue interaction cutoff strength (I min ) were outlined in our previous studies. 103 −105 We constructed the residue interaction networks using both dynamic correlations 102 and coevolutionary residue couplings 105 that provide complementary descriptions of the allosteric interactions and yield robust network signatures of long-range couplings and communications. The ensemble of shortest paths is determined from the matrix of communication distances by the Floyd−Warshall algorithm. 106 Network graph calculations were performed using the python package NetworkX. 107 Using the constructed protein structure networks, we computed the residue-based betweenness parameter. The short path betweenness centrality of residue i is defined to be the sum of the fraction of shortest paths between all pairs of residues that pass through residue i: where g jk denotes the number of shortest geodesics paths connecting j and k, and g jk (i) is the number of shortest paths between residues j and k passing through the node n i . The betweenness centrality metric is also computed by evaluating the average shortest path length (ASPL) change by systematically removing individual nodes. 108,109 A Z-score is then computed for each node based on the change of ASPL compared to the initial one. The Girvan−Newman algorithm 110,111 is used to identify local communities. To characterize global bridges from a community structure, we employed the intercommunity bridgeness metric. 112−114 The algorithmic details were specified in our recent studies. 115, 116 The bridgeness parameter uses community detection as input: Here the sum is over communities J (different from the community of node i, denoted as I), δ iJ is equal 1 if there is a link between node i and community J and 0 otherwise. l ij corresponds to the distance between community I and community J. It is measured by the inverse of the number of links between them. Nodes that are only linked to nodes of their own community have G(i) = 0, while nodes involved in bridging two (or more) communities have a positive value of the index. All parameters were computed using the python package NetworkX. 107 The network parameters were also evaluated by Cytoscape package for network analysis. 117 The ModuLand program within the Cytoscape platform 118,119 was adapted to determine a hierarchical network structure and We employed multiple CABS-CG simulations followed by atomistic reconstruction and refinement to provide a detailed comparative analysis of dynamic landscapes that are characteristic of the major functional states of the SARS-CoV-2 S trimer. Atomistic structures of the SARS-CoV-2 S prefusion trimer mutants in multiple functional states were selected for this analysis. 36 These mutant variants of SARS-CoV-2 spike protein were originally chosen in structural studies to probe the dynamic equilibrium between functional states by targeting contact regions between the RBD and S2, the RBD and NTD, the SD1 and S2, and the SD2 and S2 (Figures 1, 2) . Simulations were performed for the closed form of a double cysteine mutant, S383C/D985C (RBD to S2 double mutant (rS 2d), pdb id 6X29), that introduced a covalent link between RBD and S2 regions at the C-terminus of HRR1 and was shown to stabilize the "all-down" conformation ( Figure 2A ). In addition we conducted CABS-CG simulations for the closed, partially open (1 RBD-up) and open forms (2 RBDs up) of the SARS-CoV-2 S SD1 to S2 quadruple mutant (u1S2q) (A570L T572I F855Y N856I) ( Figure 2B-D) . While all-atom MD simulations with the explicit inclusion of the glycosylation shield could provide a rigorous assessment of conformational landscape of the SARS-CoV-2 S proteins, such simulations remain technically challenging and computationally demanding due to the enormous size of a complete SARS-CoV-2 S system embedded onto the membrane. 56−59 The abundance and complexity of information in all-atom simulations of glycosylated S protein may also obscure the networking signatures and global modulation effect of the glycan shield 59 on the protein dynamics, energetics and particularly long-range allosteric couplings examined in our study. We opted to combine the efficiency and reproducibility of simplified CABS-CG simulations with more detailed energetic and network analysis of all-atom reconstructed trajectories that also incorporated the structurally resolved glycan environment in the key glycosite positions. The proposed approach was adopted to characterize the global dynamic and network signatures of the SARS-CoV-2 S structures that are generally robust to perturbations of the glycan microenvironment and driven mainly by the underlying topology and modularity of the S structures. We also used the proposed approach to identify hinge regions and key regulatory switch centers that can impose a primary control over conformational changes between distinct functional states of the SARS-CoV-2 S structures. These global dynamic and network features of the S protein landscapes are assumed to be largely determined by the underlying topology of the SARS-CoV-2 S structures and can be captured through a proposed combination of CABS-CG dynamics and subsequent atomistic refinement of trajectories with structurally resolved glycans. Using CABS-CG simulations, we examined how mutations in the interdomain contacts could affect the global dynamic profiles of the closed, partially open, and open states revealing the important regions of flexibility ( Figure 3) . We also explored connections between the intrinsic dynamics of the SARS-CoV-2 S protein and their propensity to access up and down RBD states. A comparative analysis of the dynamics profiles showed an appreciable and consistent difference between a considerably more flexible S1 subunit and generally structurally rigid S2 subunit in all functional states ( Figure 3 ). The S2 subunit exhibited very small thermal fluctuations with RMSF < 1.0 Å, which is consistent with a strong sequence conservation in these regions. 120−122 Indeed, the S2 domain displays considerably higher conservation as compared to the S1 domains, particularly in the FP region (residues 812−829), the HR1 (residues 910−985), and the CH regions (residues 986−1035). The CTD1 (residues 529− 591) ( Figure 2E ) is a structural relay between RBD and S2 regions that can arguably sense the functional movements in the S1 and S2 domains and communicate signal from and to the fusion peptide. 34 These segments of the S1 subunit are largely stable in all functional states of the S trimer ( Figure 3A−D) . The differences in mobility of the closed-down forms of the S trimer could be clearly seen from the RMSF profile of the S383C/D985C (RBD to S2 double mutant (rS 2d) ( Figure 3A ) and another "all-down" spike form of the quadruple mutant (u1S2q) (A570L T572I F855Y N856I) ( Figure 3B ). Strikingly, a significant stabilization of the S383C/D985C closed trimer can be contrasted with the mobility of the closed form of the quadruple mutant that may be triggered by rigid body shifts of SD1 domain and displacement of the A570L + T572I-containing loop from the native position. These findings are consistent with the experimental structural data suggesting that mutation-induced structural changes could trigger the increased flexibility in both the RBD and NTD regions. 36 Conformational dynamics of the quadruple mutant in the partially open form with 1 up RBD ( Figure 3C ) was similar to that in the closed mutant state ( Figure 3B) , showing an appreciable degree of mobility in the S1 subunit regions and a fairly sharp contrast with structurally stable S2 regions. In particular, S2 subdomains HR1, CH, and CD showed considerable stability during simulations. As expected, the mobility of the RBD-up protomer was greater than that of the other "down" RBDs ( Figure 3B ,C). Both closed and partially open states featured larger thermal fluctuations in the NTD and RBD regions, particularly of the exposed RBM segment involved in recognition of the host receptor. Interestingly, even in a fully closed conformation of the SARS-CoV-2 S quadruple mutant with all RBDs in the "closed-down" position, the flexibility of the RBDs can be maintained ( Figure 3B ). These results are consistent with our recent studies suggesting that closed and open prefusion forms of the SARS-CoV-2 spike trimer could become more dynamic than the "locked-closed" trimer conformations in which mobility of the NTD and RBD regions are largely suppressed. 64 Strikingly, conformational dynamics of the closed state for the S383C/D985C double mutant ( Figure 3A Hinge Regions in the SARS-CoV-2 Spike Structures Aligned with Mutational Sites That Drive Dynamic Shifts between Functional Forms. Here we examined the conjecture that SARS-CoV-2 S mutants can modulate the intrinsic equilibrium of the S protein by exploiting preexisting soft modes of collective motions defined by the protein topology. To characterize collective motions and determine the distribution of hinge regions in the SARS-CoV-2 S trimers, we performed PCA of trajectories derived from CABS-CG simulations and also determined the essential slow modes using ENM analysis. 89−92 The essential profiles for the slowest modes and reported the functional dynamics profiles were averaged over the first three major low frequency modes (Figure 4) . The conserved hinge regions that can regulate the interdomain movements between RBD and NTD are localized at F518 and V539 positions ( Figure 3A−C) . These hinges are responsible for the coordination of the functional movements of the NTD and RBD with respect to the S2 subunit and are shared by all states. The second hinge position V539 is localized in the CTD1 region (residues 529−591) and could coordinate the concerted dynamic changes between RBD and FP regions. Importantly, several additional hinge sites in the closed and partially open forms were precisely aligned with the sites of mutations. These hinge positions correspond to residues S383, K386, D428, T547, L570, I572, S591, S750, Y855, I856, T962, and D985 ( Figure 4A−C) . Notably, the slow mode shape for the rS 2d S mutant showed that RBD regions (residues 331−528) may be largely restricted in collective motions as these residues corresponded to minima along the essential mobility profile ( Figure 4A ). This reflected the enhanced stabilization of the RBD regions induced by the introduction of the disulfide bond S383C/D985C at the critical juncture between RBD and S2 domains. In contrast, the RBD regions are more prone to functional movements in the closed form of the u1s2q S mutant ( Figure 4B ) as the slow mode profile revealed that RBD residues are now aligned with local maxima and correspond to moving regions. The hinge centers in the closed state of the u1S2q S mutant corresponded to small regulatory clusters near F318, V539, 570−572, and 750−755 residues ( Figure 4B ). It is worth noting that L570, I572, Y855, I856, and S591 hinges are situated near the interdomain SD1-S2 interfaces that act as regulatory switch Figure 4D) . A direct comparison of ENM-based slow mode profiles for the quadruple mutant (u1S2q) (A570L T572I F855Y N856I) ( Figure 4D ) and S-2P (K986P/V987P) mutant ( Figure 4E ) showed that conformational plasticity of S protein and minor topological changes could lead to different collective dynamics signatures. In the S-2P mutant, we observed noticeable deviations in slow mode motions for the "standing-up" RBDs showing reduced mobility and similar amplitudes of displacements for both NTD and RBD regions. The S2 subunit in the S-2P mutant can feature moderate slow mode motions around local hinge sites ( Figure 4E ), while S2 is completely immobilized in the quadruple mutant ( Figure 4D ). Importantly, the collective dynamics in the open S quadruple mutant trimer with 2 up RBDs can be characterized by large movements of S1 subunit and RBDs around the largely immobilized S2 regions. These functional motions can promote access and binding with the host receptor and transition to the fully open state with all RBDs in the active position. On the basis of these findings, we argue that the quadruple mutant could also prime the S trimer for efficient binding with ACE2 and initiation of the membrane fusion and transition to the postfusion form. Structural mapping of the slow mode profiles illustrated that structurally immobilized region of functional motions in the quadruple trimer corresponded to the core of the S2 subunit, while the NTDs and RBDs in both forms are prone to large conformational changes ( Figure S1 By mapping major hinge centers onto the protomer conformations, we observed that regulatory centers of The Journal of Physical Chemistry B pubs.acs.org/JPCB Article functional motions in the closed form could coordinate rigid body rearrangements of RBDs and NTDs and global displacements of the S1 regions with respect to more rigid S2 subunit. The network of hinge sites in the closed form highlights the position of these residues at the border of SD1 and SD2 regions and between SD1 and RBD regions ( Figure 5A ). The functional movements of RBDs are dictated by main hinge centers located near F318, S591, F592, V539, L570, and I572 residues. Structural map of these hinge centers indicated that these sites can form a protective shield around one of the key salt bridges between D614 of one protomer and T859 of the adjacent protomer ( Figure 5B) . To probe a potential role of glycosylation sites in assisting functional dynamics and collective motions in the SARS-CoV-2 mutant structures, we highlighted the structural alignments of 16 glycosylation sites out of 22 most occupied N-glycans in each protomer (N122, N165, N234, N282, N331, N343, N603, N616, N657, N709, N717, N801, N1074, N1098, N1134, and N1158) projected onto conformational maps of the essential profiles in the closed and open forms ( Figure S1 , Supporting Information). Strikingly, this simple analysis showed a strong overlap of the glycosylation sites dispersed across both the S1 and S2 subunits with the regions of increased mobility in collective motions ( Figure S1 , Supporting Information). Among major conserved glycosites, eight are located in the NTD regions and two (N331 and N343) are located in the RBD but outside of the RBM region accessible to ACE2. Some of these glycosylation sites (N122, N149, N74, N165, N234, and N343) are aligned with peaks of the essential profile, especially in the open states ( Figure 4C,D) . In fact, the showed that N165 and N234 corresponded almost precisely to the local maxima in the NTD regions, while N331 of RBD is aligned with the local hinge site responsible for collective movements of RBD with respect to S1 subunit ( Figure 4C,D) . These findings may partly rationalize the experimental results showing that N165 and N234 glycans may not only shield the S protein from antibody-meditated immunity but also provide conformational control acting as stabilizers of the open state and inhibitors of transitions to the closed state. 49 According to our observations, by interacting with the key mobile elements of collective motions in the open state and restricting their movements, N165 and N234 glycans would likely impede conformational transitions to the closed states. Combined with the fact that the N331 glycosite is aligned with the local hinge positons, our results suggested a plausible scenario in which N331, N165, and N234 may be coupled during allosteric changes and act synchronously to modulate conformational preferences of the S protein. These findings are also consistent with microsecond all-atom MD simulations of a fully glycosylated SARS-CoV-2 S protein embedded in a viral membrane model showing a collective behavior of these glycans that could form topologically important glycan microdomains. 59 Recent experimental studies presented strong evidence for the presence of O-glycosylation at T323 position and indicated possible glycosylation at S325 of the S protein. 122, 123 A possible O-glycosylation at Thr323 of SARS-CoV-2 subunit S1 glycoprotein has been predicted by computational analysis. 124 Interestingly, our results showed a close overlap between glycosites N331, T323, and S325 and important hinge center in the RBD region anchored by F318 that is in the close proximity of these glycosites ( Figure S1 , Supporting Information). 125 Hence, glycans may provide an additional stabilizing contribution to immobilize the hinge site in the RBD region, thereby modulating collective conformational changes in the S protein. Ensemble-Based Alanine Scanning and Energetic Analysis of the Interprotomer Contact Sites in the SARS-CoV-2 Spike Mutants. We employed the equilibrium ensembles generated by CABS-CG simulations of the SARS-CoV-2 mutant to perform a systematic alanine scanning of the protein residues ( Figure 6 ). We also highlighted and examined stability and binding free energy changes for the interprotomer residues. The objective was also to identify whether mutational sites of the interdomain contacts and hinge positions in the SARS-CoV-2 trimers also corresponded to notable energetic hotspots. Although the general pattern of the residue-based alanine scanning profile in the SARS-CoV-2 S mutants is preserved, the structural differences and subtle packing alterations and rearrangements were revealed through this analysis. As expected, in the closed form of the rS 2d mutant ( Figure 6A ) quadruple mutant ( Figure 6B) , we observed the largest number of interprotomer sites that contribute significantly to protein stability. These hotspot cluster are broadly distributed, including the RBD residues and CTD1 domain (residues 529−591) but most of the hotspots are localized in the S2 subunit regions, including the HR1 (residues 910−985), CH (residues 986−1035), and UH regions (residues 736−781). These peaks included hydrophobic residues Y789, F759, F592, F562, Y756, S383, Y396, L570, and I572 ( Figure 6A,B) . Some of these hotspots included mutational sites and regulatory positions such as S383, L570, and I572. However, mutations in the key hinge positions F318 and V539 produced destabilization changes of 1.53 and 2.1 kcal/mol that are similar to a spectrum of hydrophobic tightly packed residues located away from regulatory regions. The experimental studies showed that the D614G mutation can shift the population of the SARS-CoV-2 S trimer toward open states. 42 According to the energetic analysis, the D614G mutation yielded a fairly moderate energy loss of 0.8 kcal/mol as this important residue is relatively exposed. These results suggested that the local energetic loss caused by the D614G mutation may not be sufficient to explain radical effects of this mutation on the dynamic equilibrium of S trimer that are likely linked with allosteric changes in the residue interaction networks and structural perturbations around important hinge centers. We found that the density of energetic hotspot clusters decreases in the partially open and open states ( Figure 6C,D) . In the open form of the S quadruple mutant, most of the interprotomer hotspots are localized in the S2 region, which is consistent with the conformational dynamics analysis showing the increased mobility of S1 regions. Interestingly, alanine scanning showed that the interprotomer contacts in S2 subunit of the open form can be weakened as compared to the closed state ( Figure 6D ). The characteristic peaks and energetic hotspots are mostly conserved, featuring L570 and I572 positions among the most important for protein stability of the interdomain and interprotomer interfaces in the quadruple S mutant. At the same time, mutations in other hinge positions F318, S383, V539, and S591 produced more moderate destabilization changes. Hence, alanine scanning analysis demonstrated that while the hinge centers and sites of engineered mutations (S383, K386, D428, T547, L570, I572, S591, S750, Y855, I856, T962, D985) can appreciably contribute to stability of the SARS-CoV-2 S trimer, only several regulatory points (L570, I572) corresponded to the energetic protein stability hotspots. These findings suggested that mutation-induced population shifts between distinct functional states observed for rS 2d and u1S2q S proteins arise primarily not from local interaction changes at the interprotomer interfaces but are driven by longrange allosteric changes and reorganizations in the allosteric interaction networks. Protein Contact Networks and Network Clustering in the SARS-CoV-2 Spike Mutants: Probing Interdomain Communicating Regions and Mediating Centers. We constructed protein contact network models for the SARS-CoV-2 mutant proteins using a model developed in our earlier studies. 96−98 The protein contact network is a network whose nodes are the protein residues and links are active contacts between residues in the protein structure. For the rS 2d (S383C/D985C) mutant, structural modularity is symmetrically arranged as evidenced by the distribution of two welldefined clusters ( Figure S3, Supporting Information) . The active region is in the central core of the structure and points to residues crucial for structural stability. In the closed state of the quadruple mutant we observed a similar pattern of clustering partition, but some relevant differences seemed to emerge as the clusters to be more distributed and less segregated ( Figure S3 Figure S3 , Supporting Information). We observed that the two clusters are well distributed and strongly intermingled as compared to the rS 2d down conformation. The participation coefficient P values describe the propensity of nodes to contribute to the intercluster communication. 96−98 In this model, we characterized the importance of functional regions in allosteric interactions based on the average participation coefficient values of the corresponding residues involved in signal transmission between different clusters (functional domains). Table 2 reports the general features of the four structures, including the number of detected intercluster communicating residues and their hydrophobicity according to the Kyte−Dolittle scale. 126 The number of highly communicating residues (P > 0. 75 (Table 2 ). These findings suggested that there may be a redistribution of local modules and partial rewiring of the protein contact network upon adopting a receptor-accessible open state in the S trimer mutants. The participation coefficient heat maps showed several interesting patterns of the intermodular communication in the SARS-CoV-2 S trimer mutants ( Figure S4 , Supporting Information). The actively participating residues (P > 0) in both closed forms are similarly distributed across the entire trimer ( Figure S4, Supporting Information) . However, this distribution is more biased toward the interdomain regions between RBDs and S2 unit in the S383C/D985C double mutant, while communicating sites form more densely populated clusters in the S2 regions of the quadruple mutant ( Figure S4, Supporting Information) . Structural mapping of communicating residues in the closed states showed a high density of these sites in the S2 regions The second, fourth and sixth columns report the mean hydrophobicity of the whole structures, of the communicating regions (P > 0), and of the highly active residues (P ≥ 0.75), respectively. The third and the fifth columns reports the number of communicating residues with P > 0 and of the number of the intermodular communicating hotspots (residues with P ≥ 0.75). The Journal of Physical Chemistry B pubs.acs.org/JPCB Article and the interdomain S1/S2 regions ( Figure 7A ,B). The potential communicating hotspots with high participation coefficient values (P > 0.75) in the closed S states included residues in the CTD1 region (N532, L533, N542, V539, T553), FP region (V826, T827), the HR1 (Q913, V915, Q926, A944, L948, C985) (Tables S1 and S2, Supporting Information). Notably, this list of communicating hotspots included C985 that is disulfide-bonded to S383 in the double mutant and residues in the FP and CTD1 regions mediating communication with the RBD regions (Table S1 , Supporting Information). In the partially open state with a single protomer in the "up" position, there is a noticeable narrowing of the intermodular communicating sites that become consolidated at the border of the two major clusters located near the hinge regions between RBD and S2 ( Figure 7C ) including positions F318, V539, Y855, I856, N960, V963 (Table S3 , Supporting Information). Some of these residues are involved in the interprotomer contacts near mutational sites (L570-V963, N960-L570). In the open form of the quadruple S mutant with 2 up RBDs, the clusters are rearranged, and clusters of communicating residues become less dispersed and more localized in the S2 regions and near the interdomain border of the S1 subinit (Figure 7 D) . The majority of strongly communicating residues in the open form are located in the CH and UH regions of the S2 subunit, including S555, VL585, D586, L752, and F759 residues (Table S4, Supporting Information) . Notably, the communicating sites near S1 and S2 borders (F318, V539, L570, and I572) that featured significant participation coefficients corresponded to the key hinge centers implicated in regulation of functional movements These findings suggested that different functional states of S mutant trimers can share the intercluster communicating sites that also correspond to major hinge centers of collective dynamics. At the same time, protein contact network analysis pointed to a partial redistribution of communication hubs in the open trimer that may reflect reorganization of the allosteric interaction networks and modulation of signal transmission mechanisms. Dynamic Network Modeling and Centrality Analysis Identify Regulatory Regions Mediating Allosteric Interaction Networks in the SARS-CoV-2 Spike Mutants. Understanding of the interplay between population-shift mechanism governing allosteric structural changes and "dynamics-driven" mechanism in which allosteric interactions can be mediated through alterations of functional motions and The Journal of Physical Chemistry B pubs.acs.org/JPCB Article rebalancing of rigid and flexible protein regions is central objective of mechanistic models underpinning allostery. 65−69 These models allow for quantitative analysis of allosteric molecular events in which conformational landscapes of protein systems can be remodeled by various perturbations such as mutations, ligand binding, or interactions with other proteins. In particular, conformational changes between open and closed forms of SARS-CoV-2 S protein and process of virus entry on the molecular level may exemplify a potential synergistic effect of population-shift and dynamically driven mechanisms, where allosteric structural changes between spike states are intimately coupled with an exchange of conformational mobility that accompany S protein binding with the host receptor. Structure-based network approaches have offered a simple and effective formalism for describing allosteric interactions, where the dynamic fluctuations are mapped onto a graph with nodes representing residues and edges representing weights of the measured dynamic properties. 65−69 Network-centric models of protein structure and dynamics used in the current study can also provide a complementary perspective to physics-based landscape models and allow for quantitative analysis of allosteric changes, identification of regulatory control centers and mapping of allosteric communication pathways. 108,127−129 Using this framework, the residue interaction networks in the SARS-CoV-2 spike trimer structures were built using a graph-based representation of protein structures 100−102 in which residue nodes are interconnected through both dynamic 102 and coevolutionary correlations. 105 By employing perturbation-based network modeling, we computed ensemble-averaged distributions of several residue-based topological network metrics: the short path residue centrality ( Figure 8A − C), the betweenness centrality based on average shortest path length (ASPL) change upon removing individual nodes ( Figure 8D−F) , and the intercommunity bridgeness ( Figure 8G −I). These network metrics were explored to establish consensus in identifying mediating centers of allosteric interactions in the SARS-CoV-2 mutant structures. We found that the high centrality residues are assembled in tight interaction clusters with the distribution of peaks that is generally similar for the closed and partially open forms of u1S2q S mutant ( Figure 8A ,B). In these structures, the network centrality peaks tend to align with the hinge center regions (residues 315−331, 569−572, 756−760), indicating that these key regulators of functional motions may also mediate communication in the residue interaction networks. Moreover, high centrality residues often precisely corresponded to the engineered mutational sites, suggesting that their functional role may be partly determined by global mediating propensities in the residue interaction networks ( Figure 8A−C) . In addition, there are a number of local peaks in the HR1 region (residues 910−985), CH regions (residues 986−1035), and UH (residues 736−781). Notably, the distribution of high centrality peaks was broad in the closed and partially open forms, indicating a large allosteric interaction network characteristic of these states. The number of local centrality peaks noticeably reduced in the open form with 2 up RBDs, even though distribution peaks were still seen in the CH regions of S2 and the interdomain borders of S1 and S2 subunits ( Figure 8C ). Characteristically, a noticeable density of high centrality sites in the CH region is present and shared in distinct conformational states of the SARS-CoV-2 spike protein ( Figure 8A−C) . A similar composition of high centrality peaks was also seen using a different betweenness centrality metric ( Figure 8D−F) . Using this network metric, we observed consolidation of local centrality peaks into three major clusters in the open form, corresponding to interdomain borders of RBD and S2, CTD1 regions, and CH residues ( Figure 8F ). The centrality profiles also highlighted a gradual narrowing of the betweenness centrality distribution as we proceed from the closed state to the open forms. These observations indicated that allosteric interaction network in the open state may be partially rewired and become more localized. Another important network signature observed in the open forms was the emergence of the largest and strongest peak aligned with residues 568−572 in the CTD1 region and strong centrality density in the NTD regions near glycosite positions ( Figure 8E ,F). These high centrality clusters mediate allosteric couplings between RBDs and NTDs as well as between RBDs and functional FP segment in the S2 subunit. Structural maps of high centrality sites highlighted a significant density of these residues in the interprotomer and interdomain regions, primarily in the HR1, CH, and UH regions of the S2 subunit ( Figure S5 , Supporting Information). Interestingly, the incorporation of glycan shield into atomistic models allowed for more detailed network analysis where the observed peaks and high centrality of glycosites in the NTD and RBD regions of the open states resulted from formation of tightly packed glycan microclusters ( Figure 8A − F). Indeed, in the open states glycans at N331 and N343 positions of "up" protomer B can interact with N122, N165, and N234 glycans from the NTD regions of another open protomer C. The increased betweenness centrality ( Figure 8B ,C) and Z-score ASPL network centrality ( Figure 8E ,F) for these regions suggested that the glycan environment at glycosites N165, N234, N331, and N343 enhances the mediating role of these positions in the intramolecular communications and may contribute to stabilization of the global allosteric interaction network. These results are consistent with recent experimental 49 and computational studies 56, 57 showing that N234 and N165 glycans can act as stabilizers of the open state and that alanine mutations at these positions destabilize the RBD up conformations. Furthermore, our findings indicated that N234 and N165 glycosites acting as mediating centers may be allosterically coupled with other regulatory elements and therefore exert a partial control on conformational transitions in the SARS-CoV-2 S mutant structures. While the density of high centrality sites is reduced in the partially open and open states, we observed mediating residues near glycosites in the NTD-RBD regions as well as high centrality clusters in the CTD1 region and at the borders between RBD and S2 regions ( Figure S5 , Supporting Information). Collectively, these findings provided several important insights. First, the observed narrowing in the distribution of centrality sites in the open states may reflect the mechanism of allosteric interactions in this form where signal transmission may proceed through fewer localized routes and functional centers connecting the S2 subunit with RBD regions. Second, several known regulatory sites that dictate dynamic switching between conformational states of the SARS-CoV-2 spike trimer overlap with the high centrality mediators of the interaction network. This may be important for facilitating therapeutic The Journal of Physical Chemistry B pubs.acs.org/JPCB Article intervention through targeting specific mediating centers involved in signaling pathways of the receptor-accessible states of SARS-CoV-2 protein. Community Analysis of Residue Interaction Networks and Intercommunity Bridgeness Analysis Suggest Principles of Signal Transmission in the SARS-CoV-2 Spike Mutants. Using community decomposition, the residue interaction networks were divided into local interaction modules in which residues are densely interconnected. A community-based model of allosteric communications is based on the notion that groups of residues that form local interacting communities are correlated and switch their conformational states cooperatively. In this model, allosteric communications can be transmitted through a chain of stable local modules connected via intercommunity bridges. By describing protein dynamics as a dynamically evolving network of interconnected modules, the topological regularities of the network structure can be identified. 65−69 We performed community analysis of the SARS-CoV-2 S trimers and computed the intermodular bridgeness profiles ( Figure 8G−I) . The bridgeness profiles reflect the differences in concentration of local communities in the S1 and S2 subunits as most of stable modules and corresponding connectors tend to be localized in the S2 regions, particularly in the HR1 (residues 910−985), CH regions (residues 986− 1035), and UH (residues 736−781). Interestingly, open forms featured a number of intercommunity bridges in the UH regions (residues 736−781), suggesting that UH residues could provide stepping stones for the global routes connecting the S2 and S1 subunits. Structural map of major communities illustrated a dense distribution of stable network modules that form hubs of an allosteric network that connects S2 and S1 subunits and mediated through the UH and CH regions of the protein core. The projection of largest and most stable communities on the protomer structures revealed the overall topological arrangement pointing to conservation of intra-and interprotomer stable modules in the S2 subunit, primarily localized in a core antiparallel β-sheet and UH and CH regions ( Figure 9A−C) . A conserved 3D organization of the core β-sheet and the UH and CH regions among different coronaviruses is a salient feature of the S glycoprotein structure. 21 The conserved communities in the S2 regions shared by all states included (Y1067-L1049-F906), (C1032-H1048-S1051-H1064), (M1029-V1033-P1053), and (F1042-E1031-R1039) that provide stabilization anchors for mediating allosteric interactions in the S2 subunit. Interestingly, the overall topology of local communities in the S2 subunit is mostly retained in all forms, suggesting that allosteric signaling mechanisms may be governed by functional regions in the core β-sheet and UH and CH regions. We also delineated conformation-specific community rearrangements that provide some insight into changes in allosteric interaction networks between different functional forms of the SARS-CoV-2 S trimer mutants. In the closed conformation, the intramolecular communities prominently ). In addition, we detected a number of local communities in the CH, FP, and CTD1 regions such as (F1042-E1031-R1039), (C1032-H1048-S1064), and (Q1054-S816-E819). A detailed structural map of local communities in the closed form is highlighted for a protomer ( Figure 9A ) and a schematic projection of largest communities is shown for the closed state of the u1S2q mutant trimer ( Figure 9D ). It can be seen that stabilizing network modules are broadly distributed across the structure, indicative of a broad and large allosteric network in the closed form. By comparing the community distributions in the closed ( Figure 9A ) and open forms ( Figure 9B ,C) we found a noticeable redistribution in the local network modularity that becomes particularly apparent in the open conformation with 2 up RBDs. When the S trimer adopts a partially open conformation, the RBD in the "up" position reduces the contacts with the S2 subunit of an adjacent protomer. Our community analysis reflected these structural changes as density of local communities in the RBD and NTD regions was appreciably reduced in the "up" protomer. At the same time, stable modules are preserved in the S2 subunit and are localized primarily in the UH and CH regions ( Figure 9B ,E). These communities include (T791-K795-L806), Figure 9F ) highlighted the reduced density and a partially rewired organization of stabilizing modules as compared to strongly interconnected closed form. The preservation of key stabilizing communities in the S2 subunit and CTD1 regions allows for the interprotomer connections and could arguably result in localized but still efficient signaling between S1 and S2 regions. Another emerging concept central to computational models of allostery is identification of regulatory control points that mediate long-range communications in the system and modeling of allosteric pathways between conformational switch centers. 65−69 The current understanding of communication paths in proteins is based on the ensemble-based statistical model that assumes existence of multiple paths where the validity and thermodynamic preferences of a particular route are often tested against experimental data using mutagenesis of key switch points that mediate specific communication. 65−69 In our study, we leveraged the results of community decomposition and used edge betweenness (or edge centrality) in the global interaction network as a proxy for modeling and analysis of allosteric communication pathways. This parameter is defined as the ratio of all the shortest paths passing through a particular edge to the total number of shortest paths in the network. Through structural mapping of the high centrality edges in most probable pathways of the ensembles, we obtained a global map of major communication routes in the interaction network of the SARS-CoV-2 trimer mutants ( Figure 10A −C). The topography of these maps showed significant differences in the distribution and density of major communication routes connecting S1 and S2 regions in different functional states of the S trimer. The communication highway system in the closed state featured a dense interconnected network with several alternative global routings that engage molecular switches connecting the protomers in key positions to enable allosteric signal transmission ( Figure 10A ,B). In the closed form, the major communication routes can proceed through key S2 segments, including primarily CH regions 986−1035 (S1021, Y1007, A1015) and UH regions 736−781 (S758, F759, Y759, N751, L754, L755, C749, E748). Notably, the UH residues provide an important part of the global route connecting S2 and S1 subunits that dominates allosteric paths. The high edge centrality links in the ensemble of communication paths featured by several key interprotomers are L570(chain A)-V963(chain B), S383(chain B)-D485(chain C), D571(chain B)-S967(chain C), P589(chain A)-Y855-(chain B), and P987 (chain A)-D427(chain C) ( Figure 10A ,B). By exploring links between the UH and CH regions, the communication paths within a single protomer can reach the regulatory switch points of P987 and D485 residues that enable the interprotomer connections ( Figure 10A,B) . We found that at this juncture communication paths can link P987 of the protomer with D427 of the one adjacent protomer, and D485 with S383 of the other adjacent protomer. Using these switch centers, allosteric paths can connect directly and efficiently to the RBD regions of the adjacent protomers ( Figure 10B ). These results highlighted the important role of S383, D485, and P987 residues that can function as critical regulators of the conformational equilibrium in S-2P and rS 2d spike mutants. 27,35−38 According to our findings, these strategically located residues can provide key mediating bridges for efficient allosteric communication paths between S2 and S1-RBD regions in the S trimer ( Figure 10A,B) . Moreover, these switch centers can act as hubs in the ensemble of allosteric communication paths in the closed form as a significant fraction of communication routes in the S2 subunit can utilize these hubs to reach RBD regions of the neighboring protomers ( Figure 10A,B) . Molecular mapping of allosteric pathways revealed a well-defined topography and significant density of communication routes with switch residues at the interprotomer contact interfaces. The central result of this analysis is that these experimentally confirmed regulatory hotspots (S383, D485 and K987P) 35−38 can also function as key network bridges of communication pathways and control switch points of signal transmission in the closed form of the SARS-CoV-2 spike trimer. It is worth noting that structural architecture of the closed form and strategic location of switch contacts can preserve the overall topology of communication routes in the closed states of the wild-type and mutant S proteins. However, not only do K986P and V987P mutations shift the equilibrium away from the "all-down" form 27 but also the efficiency of predicted pathways in the closed conformations can be compromised due to weakening of stabilizing interactions at the switch control points. The weakening of the key interprotomer switch contacts may significantly hamper the efficiency of allosteric communication in the S closed trimers. Structural mapping of hubs in the ensemble of most probable pathways indicated that a significant number of these mediating sites are located in the S2 regions and near the interprotomer and interdomain boundaries between S1 and S2 subunits ( Figure S6 , Supporting Information). The density of communication hubs is high in the closed form reflecting the broader allosteric network and diverse communication routes. However, the number of hubs can be markedly reduced in the open S trimer with the RBDs in the up conformation, making allosteric signaling more dependent on the fewer interprotomer bridges governing communications. We argue that the dense allosteric network of high centrality sites and communication hubs seen in the closed form could allow for diverse alternative routes between S2 and S1 regions, thereby protecting the robustness of allosteric signaling. At the same time, allosteric communications in the open form may become more sensitive and vulnerable to targeted mutations of the key mediating centers. In the partially open spike with one RBD in the up position (chain B), the high edge centrality links corresponded to key interprotomer contacts including L570(chain A)-V963(chain B), P589(chain A)-I856(chain B), L570(chain B)-N960(chain C), L570(chain B)-K964(chain C), and P987(chain B)-D427(chain C) (Figure 10 C,D). We found a partial redistribution of communication pathways in the partially open form, where allosteric routes in the "up" protomer proceeded via S2 subunit utilizing HR1 (residues 910−985) and passing through a hydrophobic stretch of 934−940, L948, V952, A956, V963 until connecting with the Y855/I856 switch to link to "down" protomer and reach the RBD regions ( Figure 10C ,D). An alternative path can pass through the S2 subunit of the "down" protomer by utilizing UH regions until reaching a key switch point Y855/I856 and connecting to the "up" protomer and its RBD regions. Hence, in the partially open state, the dominant pathways appeared to utilize Y855/I856 regulatory center to connect S2 regions with activated RBD in the "up" protomer. We also found that in the partially open state, a significant population of probable paths could proceed through HR1 region, including residues 934−940 that form a cluster targeted by naturally occurring mutations in SARS-CoV-2 protein. 130, 131 This region also undergoes significant conformational changes during postfusion transition 130 (Figure 10E ,F). Among high edge centrality links are a number of the interprotomer connections that could create a cascade of switches for propagating communication signal. These links included P589(chain A)-Y855(chain B), D1041 (chain B)-S1030 (chain C), G413 (chain A)-E988-(chain C), V963(chain A)-L570(chain C), I856(chain A)-F592(chain C), K964(chain A)-L570(chain C), and D427-(chain A)-P987(chain C). Importantly, some of these critical switch links dominate the ensemble of communication paths included sites L570, F592, Y855, I856, and K964. The ensemble of dominant pathways connects the FP residues I818, E819, V826 on one of the "up" protomers to the HR1 region (residues 910−985) via V952, A956, V963, and Q965, and then through the critical switch Y855/I856, the path links to the other "up" protomer by connecting to another potential switch center F592 of the CTD1 region ( Figure 10E ,F). In turn, the CTD1 region provides mediating bridges for communication routes to the RBD of the "up" monomer. Another optimal path along the "up" protomer can connect FP residues via I856 to L570/I572 in the CTD1 of the "down" protomer and subsequently to the RBD regions. Strikingly, a new network of efficient routes connecting the FP regions with RBD residues can emerge that leverages several conserved switches in the receptor-accessible "up" positons. A significant fraction of communication paths navigates through switch points Y855, I856, and F592 located in the close proximity of D614 at the interprotomer interface with T859 of the adjacent protomer ( Figure 5B ). On the basis of this analysis, we argue that D614G mutation not only abolishes the local interprotomer contacts and weakens stability of the trimer but also could adversely affect allosteric communication routes passing through switch Y855/I856 and F592 to reach the S1 subunit regions ( Figure 10E ,F). This is consistent with other studies showing that a close proximity of D614G to the hinge bending region could lead to the increased switch flexibility and promote transitions to the active "up" state, making the mutant more amenable for binding with host receptor ACE2. 132 It is worth noting that molecular mechanisms underlying the D614G effect are not well understood. 133 Several studies suggested that D614G may affect the intraprotomer energetics and favor a "one-up" open conformation by eliciting contract changes in the CTD1 and FP regions. 134, 135 In support of our results, recent cryo-EM structures of SARS-CoV-2 S mutants revealed an allosteric effect that seemingly small D614G changes in the SD2 domain have on RBD dynamics, suggesting that SD2 regions may be involved in modulating RBD transitions between up and down conformations. 136 Indeed, according to our analysis, the communication paths in the open form may often exploit FP and CTD1/CTD2 regions to connect S2 regions with activated RBDs in the up position ( Figure 10E,F) . Strikingly, even though some of the regulatory switch points are utilized by probable paths in all states, the composition of switching cascades and topography of signaling routes can considerably change between the closed and open forms. In the closed form the broad ensemble of pathways connecting S2 and S1 regions proceed through S383, D985, P987, DS427, L570, and V963 control points ( Figure 10A,B) . A different cascade of regulatory control points (L570, I572, F592, Y855, I856, and K964) may control signal transmission in the open form that featured paths navigating through FP and CTD1/ CTD2 regions to connect with RBDs. In addition, the density and width of the communication path ensemble can be progressively reduced in the partially open and open states. Although, the signaling routes become more localized in the open states, allosteric interactions and cross-talk of the S1 and S2 subunits can be maintained through a narrow ensemble of communication routes that are highly dependent on efficiency of switch points in the CTD1 region, including L570, I572, and F592. The greater diversity of allosteric communication paths in the closed-down form could ensure a proper balance of the network robustness and efficiency to maintain resilience against random mutational changes. In some contrast, efficiency of signal transmission in the open states could be more vulnerable to mutations of specific conformational switches, particularly in the CTD1 region. This study systematically examined functional mechanisms of the SARS-CoV-2 S mutants through the lens of allosteric regulation. We found that mutations of the interdomain residues can change the dynamic profiles and stability of the distinct functional states. Functional dynamics analysis revealed conserved hinge regions that appeared to be aligned with engineered mutational sites driving population shifts between closed and open states. This analysis supported the notion that the introduction of mutations at the border between SD1 and S2 domains engineered in the quadruple mutant may accelerate functional conformational changes in the SARS-CoV-2 S trimer, priming spike protein for binding with the host receptor and subsequent transitions to the postfusion state. The results of this study revealed key functional regions and regulatory centers that govern collective dynamics and allosteric interactions in the SARS-CoV-2 spike proteins. Dynamic modeling of residue interaction networks and community analysis identified stabilizing structural modules in the S2 subunit acting as anchoring regions for signal transmission in the SARS-CoV-2 S mutants. The central finding of this study is that the allosteric interactions and pathways in different functional states can be modulated through a versatile network of the interprotomer regulatory centers. We found that mutations at the key interdomain positions can differentially modulate distribution of conformational states and uniquely define topography of signal propagation pathways operating through state-specific cascades of control switch points. We argue that targeting the predicted mediating centers of allosteric interactions and communication hubs in the SARS-CoV-2 S proteins with chemical probes could be beneficial for probing molecular mechanisms underlying the virus infectivity profile and discovery of specific and broad-spectrum allosteric modulators of the SARS-CoV-2 protein functions. CoV-2 u1S2q quadruple mutants ( Figure S2 ), network clustering for the SARS-CoV-2 mutant structures ( Figure S3 ), structural mapping of the participation coefficient distributions ( Figure S4 , structural maps of high centrality sites in the SARS-CoV-2 mutant structures ( Figure S5 ), structural maps of communication pathway hubs for the SARS-CoV-2 mutant structures ( Figure S6 ), and the participation coefficients for highly communicating residues (P > 0.75) in the different functional forms of the SARS-C oV-2 double and quadruple mutants (Tables S1−S4) (PDF) Early transmission dynamics in Wuhan, China, of novel coronavirus-infected pneumonia A novel coronavirus outbreak of global health concern COVID-19: what has been learned and to be learned about the novel coronavirus disease Genome composition and divergence of the novel coronavirus (2019-nCoV) originating in China SARS and MERS: recent insights into emerging coronaviruses Clinical features of patients infected with 2019 novel coronavirus in Severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2): An overview of viral structure and host response Coronavirus envelope protein: current knowledge Characterization of the receptor-binding domain (RBD) of 2019 novel coronavirus: implication for development of RBD protein as a viral attachment inhibitor and vaccine SARS-CoV-2 cell entry depends on ACE2 and TMPRSS2 and is blocked by a clinically proven protease inhibitor A pneumonia outbreak associated with a new coronavirus of probable bat origin Genomic characterisation and epidemiology of 2019 novel coronavirus: implications for virus origins and receptor binding The SARS-CoV-2 spike glycoprotein biosynthesis, structure, function, and antigenicity: Implications for the design of spike-based vaccine immunogens Structural and functional basis of SARS-CoV-2 entry by using human ACE2 Receptor recognition by the novel coronavirus from Wuhan: An analysis based on decade-long structural studies of SARS coronavirus Cell entry mechanisms of SARS-CoV-2 Structure of SARS coronavirus spike receptor-binding domain complexed with receptor Receptorbinding domain of severe acute respiratory syndrome coronavirus spike protein contains multiple conformation-dependent epitopes that induce highly potent neutralizing antibodies Pre-fusion structure of a human coronavirus spike protein Cryo-electron microscopy structure of a coronavirus spike glycoprotein trimer Cryo-electron microscopy structures of the SARS-CoV spike The Journal of Physical Chemistry B pubs glycoprotein reveal a prerequisite conformational state for receptor binding Unexpected receptor functional mimicry elucidates activation of coronavirus fusion Cryo-EM structures of MERS-CoV and SARS-CoV spike glycoproteins reveal the dynamic receptor binding domains Cryo-EM structure of the SARS coronavirus spike glycoprotein in complex with its host cell receptor ACE2 Stabilized coronavirus spikes are resistant to conformational changes induced by receptor recognition or proteolysis Tectonic conformational changes of a coronavirus spike glycoprotein promote membrane fusion Cryo-EM analysis of the post-fusion structure of the SARS-CoV spike glycoprotein Neutralization of SARS-CoV-2 by destruction of the prefusion spike Veesler, D. Structure, Function, and Antigenicity of the SARS-CoV-2 Spike Glycoprotein Cryo-EM structure of the 2019-nCoV spike in the prefusion conformation Structural basis for the recognition of SARS-CoV-2 by full-length human ACE2 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 Tracking changes in SARS-CoV-2 Spike: Evidence that D614G increases infectivity of the COVID-19 virus Fontes-Garfias, C. R. Spike mutation D614G alters SARS-CoV-2 fitness SARS-CoV-2 D614G variant exhibits efficient replication ex vivo and transmission in vivo Structural and functional analysis of the D614G SARS-CoV-2 spike protein variant 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 Receptor binding and priming of the spike protein of SARS-CoV-2 for membrane fusion Cryo-EM analysis of a feline coronavirus spike protein reveals a unique structure and camouflaging glycans Vulnerabilities in coronavirus glycan shields despite extensive glycosylation Site-specific glycan analysis of the SARS-CoV-2 spike Glycans on the SARS-CoV-2 spike control the receptor binding domain conformation Preferred conformations of N-glycan core pentasaccharide in solution and in glycoproteins An atomistic perspective on ADCC quenching by core-fucosylation of IgG1 Fc N-glycans from enhanced sampling molecular dynamics Molecular-level simulation of pandemic influenza glycoproteins Analysis of the SARS-CoV-2 spike protein glycan shield reveals implications for immune recognition Virus-receptor interactions of glycosylated SARS-CoV-2 spike and human ACE2 receptor 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 potein A multiscale coarse-grained model of the SARS-CoV-2 virion AI-driven multiscale simulations illuminate mechanisms of SARS-CoV-2 spike dynamics Exploring dynamics and network analysis of spike glycoprotein of SARS-COV-2 Enhanced receptor binding of SARS-CoV-2 through networks of hydrogen-bonding and hydrophobic interactions Dynamics of the ACE2-SARS-CoV-2/SARS-CoV spike protein interface reveal unique mechanisms. Sci. Rep. 2020, 10, 14214. (62) Verkhivker, G. M. 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 The discovery of a putative allosteric site in the SARS-CoV-2 spike protein using an integrated structural/dynamic approach Molecular simulations and network modeling reveal an allosteric signaling in the SARS-CoV-2 spike proteins Protein allostery and conformational dynamics Controlling allosteric networks in proteins Allostery without a conformational change? Revisiting the paradigm A unified view of "how allostery works Allostery in its many disguises: From theory to applications Protein modeling and structure prediction with a reduced representation 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 Consistent view of protein fluctuations from all-atom molecular dynamics and coarse-grained dynamics with knowledge-based force-field Protocols for fast simulations of protein structure flexibility using CABS-Flex and SURPASS The protein data nank 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 ArchPRED: A template based loop structure prediction server The FALC-Loop web server for protein loop modeling Jr 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 3Drefine: an interactive web server for efficient protein structure refinement A fully automated taskoriented interface for the analysis of molecular dynamics trajectories Global dynamics of proteins: bridging between structure and function oGNM: online computation of structural dynamics using the gaussian network model DynOmics: dynamics of structural proteome and beyond 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 Modules identification in protein structures: The topological and geometrical solutions Shedding light on protein-ligand binding by graph theory: The topological nature of allostery Proteins as sponges: a statistical journey along protein structure organization principles GIANT: A Cytoscape plugin for modular tetworks A network representation of protein structures: Implications for protein stability Interaction energy based protein structure networks Dynamical networks in tRNA:protein complexes Atomistic simulations and network-based modeling of the Hsp90-Cdc37 chaperone binding with Cdk4 client protein: A mechanism of chaperoning kinase clients by exploiting weak spots of intrinsically dynamic kinase domains Dancing through life: Molecular dynamics simulations and network-centric modeling of allosteric mechanisms in Hsp70 and Hsp110 chaperone proteins Computational analysis of residue interaction networks and coevolutionary relationships in the Hsp70 chaperones: A community-hopping model of allosteric regulation and communication Algorithm 97: Shortest path Exploring network structure, dynamics, and function using NetworkX Residues crucial for maintaining short paths in network communication mediate signaling in proteins Comparing protein structures with RINspector automation in Cytoscape Community structure in social and biological networks Finding community structure in networks using the eigenvectors of matrices Detecting global bridges in networks Diversity and dissimilarity coefficients: A unified approach A general framework for analyzing diversity in science, technology and society Atomistic modeling of the ABL kinase regulation by allosteric modulators using structural perturbation analysis and community-based network reconstruction of allosteric communications Dynamic view of allosteric regulation in the Hsp70 chaperones by J-Domain cochaperone and post-translational modifications: Computational analysis of Hsp70 mechanisms by exploring conformational landscapes and residue interaction networks Cytoscape: A software environment for integrated models of biomolecular interaction networks Community landscapes: an integrative approach to determine overlapping network module hierarchy, identify key nodes and predict network dynamics ModuLand plug-in for Cytoscape: determination of hierarchical layers of overlapping network modules and community centrality Phylogenetic analysis and structural modeling of SARS-CoV-2 spike protein reveals an evolutionary distinct and proteolytically sensitive activation loop Sequence variation of SARS-CoV-2 spike protein may facilitate stronger interaction with ACE2 promoting high infectivity SARS-CoV-2, an evolutionary perspective of interaction with human ACE2 reveals undiscovered amino acids necessary for complex stability Identification of 22 N-glycosites on spike glycoprotein of SARS-CoV-2 and accessible surface glycopeptide motifs: implications for vaccination and antibody therapeutics Deducing the N-and O-glycosylation profile of the spike protein of novel coronavirus SARS-CoV-2 Glycoinformatics approach for identifying target positions to inhibit initial binding of SARS-CoV-2 S1 protein to the host cell A simple method for displaying the hydropathic character of a protein Small-world network approach to identify key residues in protein-protein interaction Network analysis of protein dynamics Integration of network models and evolutionary analysis into high-throughput modeling of protein dynamics and allosteric regulation: theory, tools and applications Spike mutation pipeline reveals the emergence of a more transmissible form of SARS-CoV-2 The impact of mutations in SARS-CoV-2 spike on viral infectivity and antigenicity Characterizations of SARS-CoV-2 mutational profile, spike protein stability and viral transmission Functional importance of the D614G mutation in the SARS-CoV-2 spike protein The SARS-CoV-2 spike variant D614G favors an open conformational state Conformational dynamics of SARS-CoV-2 trimeric spike glycoprotein in complex with receptor ACE2 revealed by cryo-EM D614G mutation alters SARS-CoV-2 spike conformational dynamics and protease cleavage susceptibility at the S1/S2 junction