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

We developed and applied a computational approach to simulate functional effects of the global circulating mutation D614G of the SARS-CoV-2 spike protein. All-atom molecular dynamics simulations are combined with deep mutational scanning and analysis of the residue interaction networks to investigate conformational landscapes and energetics of the SARS-CoV-2 spike proteins in different functional states of the D614G mutant. The results of conformational dynamics and analysis of collective motions demonstrated that the D614 site plays a key regulatory role in governing functional transitions between open and closed states. Using mutational scanning and sensitivity analysis of protein residues, we identified the stability hotspots in the SARS-CoV-2 spike structures of the mutant trimers. The results suggest that the D614G mutation can induce the increased stability of the open form acting as a driver of conformational changes, which may result in the increased exposure to the host receptor and promote infectivity of the virus. The network community analysis of the SARS-CoV-2 spike proteins showed that the D614G mutation can enhance long-range couplings between domains and strengthen the interdomain interactions in the open form, supporting the reduced shedding mechanism. This study provides the landscape-based perspective and atomistic view of the allosteric interactions and stability hotspots in the SARS-CoV-2 spike proteins, offering a useful insight into the molecular mechanisms underpinning functional effects of the global circulating mutations.

Introduction

Dissecting molecular basis and mechanisms of SARS-CoV-2 viral spike (S) glycoprotein is fundamental to understanding functional drivers of infectivity as S proteins interact with the host cell receptor.[1−8] The SARS-CoV-2 S protein consists of two subunits S1 and S2, where the subunit S1 includes an N-terminal domain (NTD) and the receptor-binding domain (RBD). The RBD undergoes stochastic transitions between closed and open conformations in which the acquisition of the “up” open form in individual protomers of the S functional trimer can promote binding with the host receptor ACE2.[9−12] Dissociation of S1 from S2 and binding with ACE2 can lead to dramatic structural rearrangements in the S2 subunit mediating the fusion of the viral and cellular membranes.[13−19] The cryo-EM structures of the SARS-CoV-2 S proteins unveiled unique conformational states of the S protein trimers in the prefusion form.[20−29] The range of functional S forms includes a structurally rigid closed form and more dynamic semiopen states that can be adopted as stable intermediates before transition to the fully open conformation.[30−34] The cryo-EM structures of the SARS-CoV-2 S trimers in situ on the virion surface and single-molecule fluorescence experiments revealed a spectrum of conformational states and transitions proceeding through functional intermediate in which all three closed RBD domains become oriented toward the viral membrane.[35,36] The cryo-EM structures of a closed SARS-CoV-2 S trimer with packed fusion peptide and an ACE2-bound S trimer showed concerted swing motions of the RBDs that accompany ACE2 recognition and determine the spectrum of S protein responses to the host receptor.[37] Conformational events associated with SARS-CoV-2 S binding with ACE2 were detailed by cryo-EM experiments showing transitions from a compact closed form of the SARS-CoV-2 S protein that becomes weakened after S1/S2 cleavage followed by acquisition of the partially open states and ACE2-bound open form priming the S protein for fusion.[38]

Th global circulating mutants of the SARS-CoV-2 S protein discovered through epidemiological surveillance prompted a significant effort in the biomedical community.[39−41] From studies of pseudotype viruses, it was initially conjectured that the D614G substitution can increase S glycoprotein incorporation into virions, reduce S1 loss, and promote enhanced infectivity. The competition assay studies revealed that the variant exhibits more efficient infection and competitive fitness but maintains similar neutralization as the ancestral wild-type virus.[40,41] D614G mutation is characterized by the enhanced transmissibility but this variant is not linked with the increased disease severity. Large-scale functional and clinical studies showed that the D614G mutation could be linked with the higher viral load and younger age of patients but does not result in the higher mortality.[42] The genomes of 5085 SARS-CoV-2 strains from two distinct COVID-19 disease waves were sequenced in the metropolitan Houston region, showing that strains with D614G mutation caused the vast majority of clinical cases in the massive second disease wave, likely due to the increased transmission.[43] Despite the dominance of circulating DG614 virus, mRNA-based vaccine candidates that encode D614 antigens have achieved 94–95% efficacies.[44]

The biochemical studies suggested a phenotypic advantage and the enhanced infectivity conferred by the D614G mutation. Contrary to the original hypothesis, a detailed structural and functional analysis found that the D614G mutation does not appreciably affect S protein synthesis and incorporation into SARS-CoV-2 particles but could reduce the binding affinity to ACE2 and lead to faster dissociation rates.[45] According to this investigation, the D614G mutation can shift the population of the SARS-CoV-2 S trimer from the closed form to an open topology of the “up” protomers with 39% of the population adopting two open protomers and 20% with all three protomers in the open conformation.[45] The increased population of the 1 RBD-up open form compared to that of the closed state was also observed in the S-GSAS/D614G structure.[46] Another study revealed the higher 84% percentage of the 1-up RBD conformation in the S-G614 protein.[47] The S-G614 mutant showed a reduced S1 shedding that was linked with the increased stabilization of the S-G614 mutant, ultimately resulting in a better incorporation of the S protein into the pseudovirion.[48] The increased stability of the S-G614 mutant was inferred from the cryo-EM structures that revealed three distinct functional states of the mutant representing a closed conformation, an intermediate closed form (or 1 RBD half-up) and a partially open 1 RBD-up conformation.[49] These structures confirmed that D614G eliminates a salt bridge between D614 of one subunit and K854 of the adjacent subunit but can promote ordering of the partly disordered loop (residues 620–640), which can strengthen the intra- and interdomain interactions and enhance the stability of the mutated S protein. Structure-based protein design and cryo-EM structural determination established that D614G and D614N mutations can result in the increased stability due to a decrease in the premature shedding of the S1 domain.[50] Other studies suggested that the D614G mutation increases entry efficiency with the enhanced ACE2-binding affinity by conferring the increased structural flexibility to S protein.[51,52] These studies advocated the mutation-induced shift to the open states favorable for the S-ACE2 interaction as primary mechanism associated with the increased infectivity of the D614G mutant. Another series of the cryo-EM structures of the S-G614 trimer showed that the S protein exists predominantly in the open, receptor-binding-competent conformation adopting a range of more open and flexible conformations than the wild-type S protein.[53] The greater openness and conformational flexibility of the open forms revealed by these S-G614 structures supported the “openness” mechanisms according to which preferences for the flexible open forms could promote ACE2 binding and lead to the higher infectivity.[53] Notably, in these structures, the important functional loops (residues 620–640 and 827–855) were completely disordered and missing in the structures, which may have enhanced the propensity for the S-G614 trimer to adopt highly flexible open states. While several alternative mechanisms were proposed to explain the experimental data, a consensus view on the mechanism underlying the functional effects and increased infectivity of the S-D614G spike mutant is yet to be established.[54]

Computer simulations and protein modeling played an important role in shaping up to our understanding of the dynamics and function of SARS-CoV-2 glycoproteins.[55−60] The initial all-atom molecular dynamics (MD) simulations of the SARS-CoV-2 S protein revealed that RBDs can undergo transitions through flexible intermediate states that reduce the free-energy barrier between closed and open forms.[55] The development of a fully glycosylated full-length SARS-CoV-2 S protein model[57] and MD simulations of the full-length SARS-CoV-2 S glycoprotein embedded in the viral membrane with a complete glycosylation profile[58] provided the foundation for rigorous atomistic modeling of SARS-CoV-2 S proteins. A “bottom-up” model of the SARS-CoV-2 virion integrated data from various experimental studies to build molecular models of SARS-CoV-2 proteins and assemble a complete virion model.[59] By providing valuable insights and establishing the blueprint for computational modeling, these studies paved the way for simulation-driven studies of SARS-CoV-2 spike proteins, also showing that conformational plasticity and the alterations of the SARS-CoV-2 spike glycosylation can synergistically modulate complex phenotypic responses to the host receptor and antibodies. Multi-microsecond MD simulations of a patch of viral membrane with four full-length, fully glycosylated S proteins allowed for mapping of generic antibody binding signatures and characterization of the antibody and vaccine epitopes.[60] A recent illumining account of computational simulation studies of SARS-CoV-2 S proteins highlighted the synergies between experiments and atomistic simulations advocating for launching community-based effort for mesoscale simulations of the complete virion.[61]

Computational studies of the SARS-CoV-2 S binding with ACE2 provided important insights about molecular determinants of the binding affinity and selectivity.[62−65] Molecular mechanisms of the SARS-CoV-2 binding with ACE2 were analyzed in our recent study using coevolution and conformational dynamics.[64] Our recent studies combined coarse-grained and atomistic MD simulations with coevolutionary analysis and network modeling to present evidence that the SARS-CoV-2 spike protein can function as functionally adaptable allosterically regulated machine that exploits plasticity of allosteric hotspots to fine-tune response to antibody binding.[66−69] We found that examining allosteric behavior of the SARS-CoV-2 pike proteins may prove to be useful to uncover functional mechanisms and rationalize the growing body of diverse experimental data. This integrated approach is also highly complementary to more rigorous and detailed mesoscale simulations of SARS-CoV-2 S proteins by providing means for testing various models and allowing for systematic profiling of protein stability, binding interactions, and allosteric communications across distinct functional states of the S proteins.

Several computational studies examined the effects of global circulating mutations on dynamics and stability of the SARS-CoV-2 S proteins.[70−74] Microsecond all-atom simulations recently probed the effects of the D614G mutation showing a higher population of receptor-accessible open states in the S-G614 mutant forming through release of asymmetry present in the native interprotomer interactions.[70] The simplified physical model characterized the dynamics of the SARS-CoV-2 S protein and mutational variants, suggesting the increased population for the D614G mutation.[71] Other studies examined dynamics and energetics of the S-D614 and S-G614 proteins in the closed and partially open conformations, suggesting that the S-G614 mutant can improve the interprotomer interactions between S1 and S2 regions.[72−74] MD simulations of other circulating variants revealed that the N501Y mutation increases ACE2-binding affinity.[75,76] The growing body of computational studies investigating molecular mechanisms of circulating variants produced interesting but often inconsistent data that fit different mechanisms.

In this study, we examine molecular mechanisms underlying the functional effects of the D614G mutation by exploring conformational landscapes in different states of the SARS-CoV-2 S protein. We combined MD simulations with the ensemble-based mutational scanning and energetic sensitivity cartography analysis to simulate a number of cryo-EM structures of the SARS-CoV-2 D614 and G614 proteins in different conformational states. Using mutational scanning and sensitivity analysis of the spike residues, we identified the evolution of stability hotspots in the SARS-CoV-2 spike structures of the mutant trimers. The results suggest that the D614G mutation can induce the increased stability of the open form acting as a driver of conformational changes, which may result in the increased exposure to the host receptor and promote infectivity of the virus. The landscape-based network community analysis of the SARS-CoV-2 spike proteins showed that the D614G mutation can enhance long-range couplings between domains and strengthen the interdomain interactions in the open form. This study provides a novel insight into the molecular mechanisms underlying the effect of the D614G mutation by examining the SARS-CoV-2 S protein as an allosteric regulatory machine.

Results and Discussion

MD Simulations of the SARS-CoV-2 S-D614 and S-G614 Mutant Proteins

MD simulations were performed for the cryo-EM structure of SARS-CoV-2 S-D614 and SARS-CoV-2 S-G614 mutants in the closed and open states (Table and Figure ). The conformational dynamics profiles of the SARS-CoV-2 S-GSAS/D614 and SARS-CoV-2 S-GSAS/G614 showed a more flexible S1 subunit and very stable S2 subunit in both closed form (Figure A) and 1-up open state (Figure B). The thermal fluctuations of both S1 and S2 regions were smaller in the closed state with an RMSF <1.0 Å. In the closed state, the RBD (residues 331–528) and CTD1 (residues 528–591) are the most stable regions in the S1 subunit (Figure A). Only marginally larger fluctuations were seen in the CTD2 region (residues 592–686) that connects S1 and S2 subunits. Conformational dynamics profiles were similar for the closed forms, and no radical changes were detected in the dynamics profile of the D614G mutant (Figure A). A significant increase in the conformational mobility of the S-D614 and S-D614G open states was seen for the RBD-up protomer (Figure B). The “up” protomer can undergo fluctuations up to RMSF = 2.5–3.5 Å, also revealing signs of moderately smaller fluctuations and greater stability for the open form of the S-G614 mutant (Figure B). Of some notice were moderately greater residue displacements in the NTD regions of the S-D614 and S-G614 open structures, including the NTDs of the two closed-down protomers (Figure B,D). In general, the dynamic profiles of the S-D614 and S-G614 trimers remained similar in the native and mutant open forms. Although our results showed no substantial alteration of the dynamic signatures caused by the D614G mutation, simulations pointed to a marginally greater stabilization of the open form for the S-G614 mutant, which is consistent with previous assertions.

Table 1

Structures of SARS-CoV2 Spike Protein Structures Examined in This Study

PDB code	description	RBD position
6XR8	SARS-CoV-2 Spike Protein Trimer D614 locked closed	3-down
7KDG	SARS-CoV-2 Spike Protein Trimer (S-GSAS) D614	3-down
7KDH	SARS-CoV-2 Spike Protein Trimer (S-GSAS) D614	1-up
7KDK	SARS-CoV- 2 Spike Protein Trimer (S-GSAS) G614	3-down
7KDL	SARS-CoV- 2 Spike Protein Trimer (S-GSAS) G614	1-up
7KRQ	SARS-CoV- 2 Spike Protein Trimer G614 locked closed	3-down
7KRS	SARS-CoV- 2 Spike Protein Trimer G614 intermediate closed	1-half-up
7KRR	SARS-CoV- 2 Spike Protein Trimer G614 open	1-up

Figure 1

Cryo-EM structures of the SARS-CoV-2 S trimer structures used in this study. (A) Cryo-EM structure of SARS-CoV-2 S-GSAS/D614 in the closed state (pdb id 7KDG). The structure is in ribbons with protomers A, B, and C colored in green, red, and blue, respectively. A close-up of the key residues and interactions near the D614 position. D614, Q613, and N616 of protomer A (green sticks) and K854 and T859 of the adjacent protomer B are in red sticks. (B) S-GSAS/D614 in the 1 RBD-up open states (pdb id 7KDH). A close-up of the key residues and interactions near the D614 position is shown. (C) Cryo-EM structure of the SARS-CoV-2 S-GSAS/G614 mutant ectodomain in the closed form (pdb id 7KDK). A close-up of the key residues and interactions is shown. (D) Cryo-EM structure of the SARS-CoV-2 S-GSAS/G614 mutant in the 1 RBD-up open states (pdb id 7KDL). A close-up of the key residues and interactions near the D614 position is shown, and residues are annotated as in panel (A).

Figure 2

Conformational dynamics of the SARS-CoV-2 spike (S) protein trimer mutants revealed by MD simulations. (A) Root-mean-square fluctuation (RMSF) profiles for the cryo-EM structure of SARS-CoV-2 S-GSAS/D614 in the closed-all-down state (pdb id 7KDG) and SARS-CoV-2 S-GSAS/G614 mutant ectodomain in the closed form (pdb id 7KDK). (B) RMSF profiles for S-GSAS/D614 in the 1 RBD-up open state (pdb id 7KDH) and SARS-CoV-2 S-GSAS/G614 mutant in the 1 RBD-up open states (pdb id 7KDL). The residue numbering is shown for three protomers A, B, and C in sequential order. (C) Structural map of the conformational mobility profiles for the S-GSAS/D614 trimer in the closed form. (D) Structural map of the conformational mobility profiles for the S-GSAS/D614 trimer in the partially open 1 RBD-up state. The colored-coded maps show the rigidity–flexibility changes in blue-to-red spectrum with the most stable regions in blue and most flexible residues in red.

A close inspection of the local dynamics profiles revealed only minor flexibility changes near the mutational site, where the loss of favorable interactions with T859 of the adjacent protomer may be compensated through the intraprotomer contacts, especially for the up protomer in the open state. These results revealed small and largely synchronous dynamic changes in the closed and open forms of the D614G mutant, showing no indication of a dramatic alteration of the dynamic signatures to suggest a clear mechanistic trigger for the dynamic preferences of the D614G mutant toward the open form, as was proposed in a computational study.[71]

We also performed atomistic MD simulations of the S-D614 timer in the locked closed form (pdb id 6XR8) and structures of the S-G614 mutant in the locked closed, intermediate, and open forms (pdb id 7KRQ, 7KRS, and 7KRR, respectively).[49] This analysis provided further evidence of the mutation-induced control of conformational dynamics in the S-G614 mutant (Supporting Information, Figure S1). In the closed form of the S-G614 mutant, the FPPR segments (residues 823–862) and loops (residues 620–640) are fully ordered, while in the S-D614 locked closed state, only the FPPR region is structured (Supporting Information, Figure S1A). In the intermediate conformation with a partly open single RBD, the loops are disordered in the RBD-shifted module and FPPR remains fully ordered (Supporting Information, Figure S1B).

The thermal fluctuations of both S1 and S2 subunits are very small in the locked closed state of the S-D614 trimer (Supporting Information, Figure S2A). The mobility of the RBD residues was largely suppressed, and only peripheral regions of the NTD showed some level of residual mobility, while thermal stabilization of the structurally ordered FPPR motif (residues 828–853) was evident. In the locked closed form and intermediate conformational state of the S-G614 mutant, there is an expansion of mobility regions in both the NTD and RBD regions (Supporting Information, Figure S2B,C). In the open conformation, the S-G614 mutant featured only moderately increased mobility in the RBD of the up protomer, while the two closed protomers exhibited enhanced rigidity with the essentially immobilized NTD and RBD regions (Supporting Information, Figure S2D).

The covariance matrixes of residue fluctuations highlighted a general similarity of correlated motions, but a more detailed comparison pointed to stronger couplings between S1 and S2 domains in the open S-G614 trimer (Supporting Information, Figure S2E–G). These observations suggested that long-range couplings between functional regions could become stronger in the open mutant form. We argue that by enhancing the interdomain communications in the open state, the D614G mutation could potentially improve stability of the S-G614 trimer and limit shedding of the S1 subunit.

Using principal component analysis (PCA) of MD trajectories, we examined functional dynamics of the S-G614 trimer states by computing the essential mobility profiles and identifying the hinge regions (Figure ). The regulatory regions that control collective motions are conserved as the key hinge positions are shared between different functional states and remained intact in the S-D614 trimer and G614 mutant trimer (Figure ).

Figure 3

Functional dynamics of the SARS-Cov-2 S-GSAS/D614 structures in the closed and open forms obtained using PCA of atomistic MD trajectories. The essential mobility profiles are averaged over the three lowest-frequency modes. (A) Essential mobility profiles for the SARS-Cov-2 S-GSAS/D614 in the closed form (pdb id 7KDG). (B) Slow mode profile for the SARS-Cov-2 S-GSAS/D614 structure in the 1-up open state (pdb id 7KDH). (C) Essential mobility profiles for the SARS-CoV-2 S-GSAS/G614 structure in the closed form (pdb id 7KDK). (D) Slow mode profile for the SARS-Cov-2 S-GSAS/G614 structure in the 1-up open state (pdb id 7KDL). The profiles for protomer chains A, B, and C are shown in green, red, and blue lines, respectively. The inset on each panel shows structural maps of the SARS-CoV-2 S trimers colored by the main-chain deformability where high deformability regions are colored in red.

The local minima regions in each protomer are localized around K310-F318 and S591-V618 regions in the S1 domain where the first cluster can modulate NTD–RBD motions, while the second hinge center is responsible for collective dynamics between RBD and S2 regions. These observations are in line with experimental studies, indicating that the D614G mutational site may be located in the immobilized structural region of the SD2 domain where the local environment of D614 combined with the β strand formed by residues 311–319 may correspond to a hinge center governing motions of NTD and RBD, as well as isolating the motions in S1 from the S2 subunit.[46] Importantly, the key hinge residues F318, S591, F592, F855, and N856 residues correspond to the local minima of the essential mobility profile for each of the protomer and are involved in the interprotomer interactions (Figure ). These sites are located in the close proximity of D614 and could form several interaction clusters in both the closed and open states of the SARS-CoV-2 S-D614 and SARS-CoV-2 S-G614 structures. Structural mapping of high deformability regions in collective motions showed the broader distribution of these sites in the closed form and a strong but localized density for the partially open state. The structural conservation of the interprotomer hinge centers is also seen in the slow mode profiles pointing to the immobilization of the key regulatory positions that control collective dynamics across distinct states of S trimers (Figure ).

We expanded this analysis by comparing the essential mobility profiles of the native S-D614 closed trimer with the S-G614 mutant in the locked/closed, intermediate, and open forms (Supporting Information, Figure S3). In the locked closed form of the S-D614 trimer, the movements of the RBD and NTD regions were completely suppressed (Supporting Information, Figure S3A). On the other hand, in the closed and intermediate states of the S-G614 mutant, we noticed moderate but appreciable displacements of the RBD regions along the slow modes for all protomers (Supporting Information, Figure S3B,C). These dynamic preferences of the RBDs may lead to more frequent shifts between the closed and open forms and arguably lower the kinetic barriers for the transitions between functional forms in the S-G614 mutant. Interestingly, in the open state of the S-G614 mutant trimer, only the up-positioned protomer featured large displacements of the NTD and RBD regions, while these regions in the closed-down protomers become mostly immobilized (Supporting Information, Figure S3D). The observed pattern of collective motions in the S-G614 mutant would imply that the RBD regions may experience functional displacements in the closed and intermediate forms that encourage movements toward the open form. At the same time, once the S-G614 trimer attains the 1-up open form, the closed protomers become rigidified with their RBD movements being largely curtailed. As a result, unless open RBD becomes fully engaged by the ACE2 receptor, this may arguably prompt the open protomer to return to the down conformation, as was proposed in structural studies.[49] These results may reconcile the experimental studies[45,49] showing that the D614G mutation confers the increased structural flexibility in the closed state, which may promote the exchange between the open and closed forms and the increased exposure of the RBDs for interactions with ACE2.

Mutational Sensitivity Cartography in the SARS-CoV-2 Spike Trimers Reveals D614G-Induced Stabilization of the Closed and Open States

To provide further analysis of molecular drivers underlying functional effects of the D614G variant, we conducted the ensemble-based profiling and mutational scanning of protein stability in different functional forms of the S-G614 trimer. Through this analysis, we examine a hypothesis that the enhanced infectivity of the S-G614 mutant could be determined by the interplay between the increased thermodynamic stability of the S-G614 trimer in the open state and the enhanced kinetic accessibility of S-G614 RBD regions for interactions with the host receptor. We employed the equilibrium ensembles generated by MD simulations of the SARS-CoV-2 protein structures to perform deep mutational scanning and mutational sensitivity analysis of key functional sites in the native and mutant proteins (Figures and 5). We constructed mutational sensitivity heatmaps for the major hinge centers in the S-D614 closed trimer (Figure A) and partially open form (Figure B) by systematically substituting each residue to all possible variants and estimating the resulting effect on protein stability. Mutational cartography revealed several interesting trends and highlighted the energetic effect of the D614 site. First, we noticed that the energetic hotspots corresponded to hinge centers F592 and F855 showing the larger contributions in the closed state (Figure A). A similar heatmap was obtained for the S-D614 open form where mutations in the F855 hinge position caused large destabilization changes (Figure B). Importantly, the D614 position appeared to be tolerable to most modifications in both closed and open forms (Figure A,B).

Figure 4

Mutational sensitivity analysis for the SARS-CoV-2 S-D614 trimer in the closed and open states. (A) Mutational scanning map of the S-D614 trimer in the closed form. (B) Mutational scanning map of the S-D614 trimer in the open form. The rows of the heatmaps show the effect of all single mutations at the corresponding residue positions on the folding free-energy changes. The squares on the heatmap are colored by mutational effect in a 3-colored scale from red to light blue, with red indicating the largest destabilization effect. The data bars correspond to the computed binding free-energy changes. (C) Close-up of the interprotomer hinge cluster in the closed state. (D) Close-up of the interprotomer hinge cluster in the open state. D614, F318, and F592 belong to one of the protomers and are shown in green spheres. K854, F855, N856, and T859 belong to the other protomer and are shown in cyan spheres.

Figure 5

Mutational sensitivity analysis for the SARS-CoV-2 S-G614 trimer in the closed and open states. (A) Mutational scanning map of the S-G614 trimer in the closed form. (B) Mutational scanning map of the S-G614 trimer in the open form. The rows of the heatmaps show the effect of all single mutations at the corresponding residue positions on the folding free-energy changes. The squares on the heatmap are colored by mutational effect in a 3-colored scale from red to light blue, with red indicating the largest destabilization effect. The data bars correspond to the computed binding free-energy changes. (C) Close-up of the interprotomer hinge cluster in the closed state. (D) Close-up of the interprotomer hinge cluster in the open state. G614, F318, and F592 belong to one of the protomers and are shown in green spheres. K854, F855, N856, and T859 belong to the other protomer and are shown in cyan spheres.

Structural mapping showed that hinge residues F518 and F592 in one protomer can form the interprotomer cluster with K854, F855, N856, and T859 hinge sites of the adjacent protomer (Figure C,D). D614 packs tightly with T859, K854, F855, and N856 of the other protomer while also maintaining favorable contacts with F318 and F592 sites of the same protomer. Interestingly, D614 can anchor these clusters by connecting major hinge residues (Figure C,D). Notably, tight packing in this cluster is preserved in the closed and open forms, which is determined by the interactions between F592 and F855 (Figure C,D). Hence, structural and mutational changes in the D614 position may be well tolerated and accompanied by repacking of the interprotomer hinge cluster, thereby preserving the functional role of this region in regulation of collective motions.

Mutational sensitivity mapping of the hinge residues in the S-G614 mutant trimer revealed a redistribution of the energetic hotspots (Figure A,B). In the closed form of the G614 trimer, the stability hotspots corresponded to F592, F855, and G614 sites as mutations in these positions consistently resulted in large destabilization changes (Figure A). Of particular interest is the emergence of G614 as one of the stability hotspots that together with F592 and F855 strengthen the interprotomer hinge cluster (Figure C). Although D614G eliminates favorable interactions with T859 and K854, the tight packing is maintained between residues of the interprotomer hinge cluster (Figure C). In the open form of the S-G614 trimer, the stability hotspots F592 and G614 become dominant and anchor the interprotomer hinge cluster (Figure B). Importantly, the heatmap showed that G614 is more stable in the open state than in the closed state, which may drive conformational preference of the mutant trimer toward the 1-up form. Structural analysis of the interprotomer cluster reflected a partial rearrangement of packing interactions where contacts between F592, G614, and F855 residues determine the stability of the hinge center (Figure D). The central finding of the mutational sensitivity mapping is a significant energetic plasticity of D614 in the native trimer, while in the S-G614 mutant trimer, this position becomes the dominant energy hotspot of the interprotomer hinge center.

To provide a direct quantitative comparison of mutational scanning profiles for D614 and G614 positions in the SARS-CoV-2 S trimers, we plotted the ensemble-averaged folding free-energy changes that result from systematic modifications of these residues. In the closed form of the S-D614 trimer, most of D614 mutations can result in the improved stability, highlighting, in particular, an appreciable stabilization of the D614C, D614G, and D614N mutants (Figure A).

Figure 6

Mutational sensitivity analysis of the SARS-CoV-2 S-D614 and S-G614 trimers. Mutational sensitivity scanning of the D614 position in the closed form of the S-D614 protein (A) and open form of the S-D614 protein (B). Mutational sensitivity scanning of the G614 position in the closed form of the S-G614 protein (C) and open form of the S-G614 protein (D). The protein stability changes are shown in black-filled bars.

Notably, the D614G mutation caused the largest energetic change as compared to other substitutions. A similar pattern was seen in the mutational profile of the open state where D614G and D614N mutations displayed favorable and comparable stabilization changes (Figure B). These findings are consistent with the latest differential scanning fluorimetry studies showing that the D614G- and D614N-induced changes have a very similar effect, leading to a considerable improvement of thermal stability, which may be explained by a decrease in premature shedding of the S1 domain.[49] Consistent with this experimental study, we found that D614G/N mutations can moderate the repulsive charge interactions at the interface between S1 and S2 domains and allow for the tighter interdomain packing. The mutational sensitivity profile of the SARS-CoV-2 S-D614 in the locked closed form also showed moderate destabilization changes, whereas D614G, D614C, D614H, and D614N modifications resulted in the improved stability of the closed form (Supporting Information, Figure S4A). This was especially interesting given the rigidity of the locked closed S-D614 trimer, thus indicating that mutations in the D614 position can be accommodated even in this constrained state.

Strikingly, the mutational sensitivity analysis of the SARS-CoV-2 mutant structures revealed a far more significant destabilization effect of G614 substitutions in the closed and open forms (Figure C,D). The destabilization changes caused by mutations of G614 were greater in the open form, thereby indicating that the D614G mutation may have a stronger stabilization effect on the open state. In the locked closed form of S-G614 the largest destabilization changes were induced by the G614K, G614D, and G614E mutations (Supporting Information, Figure S4B), highlighting the fact that the reverse G614D mutation could significantly destabilize the S-G614 trimer. In the intermediate conformation, we observed larger destabilization changes for the moving protomer but also an appreciable loss of stability for the closed-down two protomers (Supporting Information, Figure S4C). An interesting pattern of protein stability changes was seen in the 1 RBD-up open state of the S-G614 trimer (Supporting Information, Figure S4D). The free-energy changes induced by mutations of G614 in the open protomer (protomer A) were somewhat smaller, reflecting a larger mobility and openness of this protomer near the site of mutation. At the same time, the losses of protein stability upon modifications of G614 in the closed protomers were more significant, indicating the increased stabilization of the S-G614 trimer (Supporting Information, Figure S4D). The results of the mutational scanning are consistent with the functional dynamics analysis, highlighting structural and energetic stability of the closed protomers in the open state of the S-G614 mutant. The preferential stabilization of the open mutant form may be driven by strengthening of the interdomain interactions and curtailing functional movements of the NTD and RBD regions in the closed-down two protomers. In this scenario, the reverse transition of the up RBD to the closed-down form can become energetically favorable unless the open protomer is engaged by recognition of the ACE2 receptor. This corroborates with the experimentally observed prevalence of 1 RBD-up conformation for the S-G614 mutant trimer.

We also employed the FoldX approach with the all-atom representation of protein structure[77−82] to evaluate residue-based folding free-energy contributions in the D614 and G614 S trimers (Figure ). The protein stability ΔΔG changes were computed by averaging the results of computations over 1000 samples obtained from MD simulation trajectories.[83,84] In this analysis, positive folding free-energy contributions indicate stability weakness sites and highly negative folding free-energy contribution characterized sites of local stability in the S protein.

Figure 7

Residue-based folding stability analysis of the SARS-CoV-2 S-D614 and S-G614 trimers. The folding free-energy profile for the S-D614 protein in the closed form (A), the S-D614 protein in the 1-up open form (B), the S-G614 protein in the closed form (C), and the S-G614 protein in the 1-up open form (D). The folding stability residue profiles are shown in color-coded bars (chain A is in gray bars, chain B is in maroon bars, and chain C is in orange bars). The positions of the RBD sites K417, E484, N501, and D614/G614 that are subjected to circulating mutational variants are shown in yellow-colored filled squares.

In the locked closed form of the D614 S trimer, the stability hotspots are distributed across multiple regions including the NTD residues, the RBD core, and the S2 domain regions (Figure A). It is worth noting that the most prominent stability centers were found in the RBD core and the region near the hinge hotspot cluster (residues 591–600) that featured stability peaks for hinge sites S591, F592, and G593 (Figure A). In both closed and open forms of the D614 S trimer, the receptor-binding motif (RBM) of the RBD region (residues 470–491) showed only marginal stability as the energetic plasticity of this region is required for functional interactions with ACE2.

In the open form of the D614 S trimer, the stability of the RBD core and S2 regions was partially reduced while maintaining the overall stabilizing signature in these positions (Figure B). By highlighting positions of circulating mutations K417, E484, N501, and D614, we noticed that D614 corresponds to a moderately stable position, while other sites are marginally destabilizing (Figure A,B). Hence, mutations in these sites can be easily tolerated without impairing protein stability of the S trimer. The overall stability pattern in the G614 S trimer remained preserved, highlighting similar regions of high structural stability (Figure C,D). However, the D614G mutation can incur a more significant stabilization in the site of mutation and in the neighboring regions, particularly enhancing stability of the hinge site cluster formed by F318, F592, G593, G614, K854, F855, N856, and T859 residues (Figure C,D). Notably, the stability of the G614 position is stronger in the open state of the S-G614 trimer (Figure D). Significantly, the open form of the G614 mutant showed a stronger stabilization of the RBD regions and S1 and S2 interdomain regions (Figure D) as compared to the open form of the D614 S protein (Figure B).

Our findings support the notion that the D614G mutation in the SD2 domain could strengthen stability of the distal NTD and RBD regions in the open state[46,47] and potentially promote exposure to the host receptor. These results offer support to the experimental observations that the enhanced stability of the S-D614G mutant may be linked with the mechanism of the reduced S1 shedding.[50]

Network Modeling and Community Analysis Suggest D614G-Induced Reorganization of the Residue Interaction Networks and Improved Allosteric Signaling in the Open States

The residue interaction networks in the SARS-CoV-2 spike trimer structures were built using a graph-based representation of protein structures[85] in which residue nodes are interconnected through both dynamic[86] and coevolutionary correlations.[87] Using community decomposition, the residue interaction networks were divided into local interaction modules. 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. We explored community analysis as a network proxy for stability assessment. The number of communities for the S-D614 protein was greater in the closed form (Figure A) as compared to that in the open form (Figure B). In some contrast, the network analysis of the S-G614 mutant showed a moderate redistribution of communities as the open state of the mutant harbored more communities (Figure C,D). Interestingly, the total number of local communities is moderately increased in both closed and open states of the S-G614 mutant. These results are consistent with the latest experimental data that demonstrated the improved stability of the D614G mutant as compared to that of the S-D614 protein, allowing for reduction in a premature shedding of the S1 domain.[49]

Figure 8

Community analysis and structural community maps in the SARS-CoV-2 S-D614 and S-D614G mutant structures. (A) Structural mapping of local communities is projected onto a single protomer for the S-GSAS/D614 in the closed-all-down state (pdb id 7KDG). (B) Structural mapping of local communities for the S-GSAS/D614 in the open state (pdb id 7KDH). (C) Structural mapping of local communities in the closed form of S-GSAS/G614 (pdb id 7KDK). (D) Structural mapping of local communities for the S-GSAS/G614 in the open state (pdb id 7KDL). The cryo-EM structures of the SARS-CoV-2 S-D614 and S-G614 trimers are shown in ribbons (chain A in green, chain B in red, and chain C in blue). The local communities are highlighted in spheres. The detailed close-ups of the local communities are shown for a single protomer, with the protomer in green ribbons and communities in red spheres.

A more detailed inspection of local communities in the closed and partially open conformational states of the native S-D614 protein provided some important insight into the proposed mechanisms (Tables and 3). We found that there are important changes in the communities situated in the RBD–CTD1, RBD–CTD2, and CTD1–CTD2 regions. The local modules in these regions that are unique to the closed form included Q314–S596–Q613, I693–V656–Y660, F543–L598–V576, and I666–L650–I670–T645 clusters (Table ). In the partially open form, the state-specific communities in these regions included P579–P330–N544, F541–L552–I587, and R328–F543–P579 (Table ). Characteristically, the key community in the closed form of the S-D614 protein is anchored by Q613, which is the immediate neighbor of D614, forming a tight stable cluster with Q314 in the NTD and S596 in CTD2. Structural mapping of local communities in the S-D614 states illustrated subtle differences in the distribution and density of stable modules (Figure A,B). A general comparison of structural maps indicated the better connectivity of local communities in the closed form of S-D614 forming a broad network linking the S2 regions with the NTD and RBD regions (Figure A,B). In addition, a number of unique communities are localized in the CTD2 region (I693–V656–Y660 and I666–L650–I670–T645), suggesting the stronger S1–S2 interfacial interactions and tighter packing between S1 and S2 domains in the closed form of the S-D614 protein. We also detected the larger number of communities in the NTD and RBD S1 regions in the closed form, indicating that the allosteric interaction network in the S1 and S1–S2 regions is stronger in the closed form, protecting all-down closed form in the S-D614 protein (Figure A). A weaker connectivity of local communities near S1–S2 interfaces and in the NTD/RBD regions was observed in the open state of S-D614 (Figure B). These findings suggested that efficiency of allosteric communication in the residue interaction network could favor the closed-down form of the S-D614 spike protein.

Table 2

Local Interacting Communities of the Protomer A in the Structure of SARS-CoV-2 S-GSAS/D614 in the Closed State (pdb id 7KDG)

number	contributing structural domains	local community interacting residues
1	HR1–HR2	A_1042_PHE A_1031_GLU A_1039_ARG
2	S2–HR1–HR2	A_1067_TYR A_1049_LEU A_906_PHE A_909_ILE A_911_VAL A_915_VAL
3	S2–HR1–HR2	A_1054_GLN A_816_SER A_819_GLU
4	NTD	A_215_ASP A_266_TYR A_64_TRP
5	NTD	A_275_PHE A_290_ASP A_58_PHE
6	NTD	A_279_TYR A_44_ARG A_47_VAL A_49_HIS
7	CTD2–NTD	A_664_ILE A_312_ILE A_598_ILE
8	RBD	A_464_PHE A_355_ARG A_514_SER A_429_PHE A_425_LEU A_426_PRO
9	RBD	A_509_ARG A_442_ASP A_438_SER
10	CTD2	A_693_ILE A_656_VAL A_660_TYR
11	S2	A_888_PHE A_789_TYR A_880_GLY
12	S2–HR1	A_1001_LEU A_1005_GLN A_759_PHE
13	HR1–HR2	A_1024_LEU A_1028_LYS A_1042_PHE A_1032_CYS
14	S2–HR1–HR2	A_1062_PHE A_1029_MET A_1033_VAL A_1053_PRO A_877_LEU
15	HR1–HR2	A_1032_CYS A_1043_CYS A_1048_HIS A_1051_SER A_1064_HIS
16	S2–HR1–HR2	A_819_GLU A_1055_SER A_874_THR
17	HR1–HR2	A_1075_PHE A_1096_VAL A_1110_TYR A_714_ILE
18	HR1–HR2	A_1106_GLN A_1109_PHE A_915_VAL
19	NTD	A_231_ILE A_130_VAL A_168_PHE
20	NTD	A_193_VAL A_204_TYR A_37_TYR
21	NTD	A_285_ILE A_279_TYR A_38_TYR
22	NTD–CTD1–CTD2	A_299_THR A_315_THR A_597_VAL
23	RBD	A_338_PHE A_342_PHE A_368_LEU
24	RBD	A_451_TYR A_401_VAL A_442_ASP
25	RBD	A_353_TRP A_398_ASP A_464_PHE
26	RBD–CTD1	A_543_PHE A_585_LEU A_576_VAL
27	CTD2	A_666_ILE A_650_LEU A_670_ILE A_645_THR
28	HR1	A_977_LEU A_749_CYS A_993_ILE A_997_ILE
29	S2	A_878_LEU A_806_LEU A_882_ILE
30	S2	A_906_PHE A_923_ILE A_916_LEU
31	HR1	A_1000_ARG A_977_LEU A_996_LEU
32	S2–HR1–HR2	A_741_TYR A_1004_LEU A_962_LEU A_858_LEU
33	S2–HR1–HR2	A_902_MET A_1050_MET A_898_PHE
34	RBD	A_342_PHE A_511_VAL A_374_PHE A_436_TRP
35	RBD	A_353_TRP A_400_PHE A_423_TYR A_512_VAL A_410_ILE
36	RBD	A_392_PHE A_515_PHE A_395_VAL
37	RBD–CTD2	A_314_GLN A_596_SER A_613_GLN

Table 3

Local Interacting Communities of the Protomer A in the Structure of SARS-CoV-2 S-GSAS/D614 in the Open State (pdb id 7KDH)

number	contributing structural domains	local community interacting residues
1	S2–HR1	A_731_MET A_1014_ARG A_955_ASN
2	NTD	A_101_ILE A_240_THR A_265_TYR A_92_PHE
3	S2–HR1–HR2	A_1028_LYS A_1062_PHE A_727_LEU
4	HR1–HR2	A_1042_PHE A_1031_GLU A_1039_ARG
5		A_1032_CYS A_1048_HIS A_1051_SER A_1064_HIS
6	S2–HR1–HR2	A_1067_TYR A_1049_LEU A_909_ILE
7	S2–HR1–HR2	A_1054_GLN A_816_SER A_819_GLU
8	HR1–HR2	A_1102_TRP A_1081_ILE A_1135_ASN
9	NTD	A_238_PHE A_92_PHE A_267_VAL
10	NTD	A_275_PHE A_290_ASP A_58_PHE
11	NTKD	A_279_TYR A_44_ARG A_47_VAL A_49_HIS
12	NTD–CTD2	A_299_THR A_315_THR A_597_VAL
13	CTD2	A_664_ILE A_312_ILE A_598_ILE
14	RBD	A_418_ILE A_495_TYR A_453_TYR
15	RBD	A_464_PHE A_429_PHE A_425_LEU A_426_PRO A_514_SER
16	RBD	A_454_ARG A_457_ARG A_467_ASP
17	HR1–HR2	A_1062_PHE A_1029_MET A_1033_VAL
18	S2–HR1–HR2	A_819_GLU A_1055_SER A_874_THR
19	HR1–HR2	A_1106_GLN A_1109_PHE A_915_VAL
20	NTD	A_193_VAL A_204_TYR A_37_TYR
21	NTD	A_285_ILE A_279_TYR A_38_TYR
22	RBD–CTD1	A_579_PRO A_330_PRO A_544_ASN
23	RBD	A_351_TYR A_454_ARG A_492_LEU
24	RBD–CTD1	A_365_TYR A_387_LEU A_515_PHE A_432_CYS
25	CTD1	A_541_PHE A_552_LEU A_587_ILE
26	S2–HR1	A_997_ILE A_749_CYS A_993_ILE
27	S2	A_898_PHE A_802_PHE A_797_PHE
28	S2	A_878_LEU A_806_LEU A_882_ILE
29	S2	A_818_ILE A_804_GLN A_935_GLN
30	S2–HR1	A_923_ILE A_916_LEU A_906_PHE
31	S2–HR1	A_905_ARG A_1050_MET A_898_PHE A_902_MET
32	S2–HR1–HR2	A_1065_VAL A_1052_PHE A_802_PHE A_805_ILE A_878_LEU A_927_PHE
33	RBD–CTD1	A_328_ARG A_543_PHE A_579_PRO
34	RBD	A_342_PHE A_511_VAL A_374_PHE A_436_TRP
35	S2–HR1–HR2	A_805_ILE A_1054_GLN A_818_ILE

The analysis revealed the increased number of stable communities for the S-G614 mutant in both the closed form (Table ) and partially open state (Table ). Both forms of the S-G614 mutant featured the appreciable number of communities in the CTD1 and CTD2 regions involved in stabilization of the S1–S2 interfaces (Figure C,D). These interdomain local communities bridging S1 and S2 regions included T315–T299–V597, R328–S530–Q580, Q314–S596–Q613, I664–I312–I598, L611–P6645–L650, and F898–F802–F797. The stabilizing communities in the S1–S2 interfaces appeared to strengthen the stability of both closed and open S-G614 forms, which is consistent with the experimentally observed increase in protein stability of the D614G mutant.[49] However, the balance in the number of communities may be shifted toward the open state of the S-G614 mutant, showing a more significant increase of local stable modules and promoting the preferential stabilization of the open mutant form. The key community near mutational site Q314–S596–Q613 is uniquely present only in the open form of the S-G614 mutant, likely pointing to state-specific rearrangements of stable interactions induced by the D614G mutation (Table ). The important subtle rearrangements were found in the local communities localized in the N2R linker region (residues 306–334) that connects the NTD and RBD regions stacking against the CTD1 and CTD2 domains. Interestingly, the open form of the S-G614 mutant featured some of these modules (T299–T315–V597 and Q314–S596–Q613) and several additional N2R communities unique to this state (R328–F543–P579, P579–P330–N543, F338–F342–L368) that strengthen connections of the N2R linker with the NTD and RBD regions (Figure D and Table ). This local subnetwork may be instrumental in bridging the individual domains in the S1 subunit and ensure a more efficient interconnectivity between S2 regions, N2R, and S1-RBD regions in the S-G614 mutant.

Table 4

Local Interacting Communities of the Protomer A in the Structure of SARS-CoV-2 S-GSAS/G614 in the Closed State (pdb id 7KDK)

number	contributing structural domains	local community interacting residues
1	NTD	A_101_ILE A_190_ARG A_94_SER
2	HR1–HR2	A_1042_PHE A_1031_GLU A_1037_SER A_1039_ARG A_1032_CYS
3	S2–HR1–HR2	A_1067_TYR A_1049_LEU A_906_PHE A_909_ILE A_911_VAL A_915_VAL
4	S2–HR1–HR2	A_1054_GLN A_816_SER A_819_GLU
5	S2–HR1–HR2	A_1063_LEU A_724_THR A_934_ILE
6	S2–HR1–HR2	A_1075_PHE A_1110_TYR A_714_ILE
7	NTD	A_265_TYR A_92_PHE A_240_THR
8	NTD	A_275_PHE A_290_ASP A_58_PHE
9	NTD	A_279_TYR A_44_ARG A_47_VAL A_49_HIS
10	NTD–CTD2	A_315_THR A_299_THR A_597_VAL
11	NTD–RBD–CTD1	A_328_ARG A_530_SER A_580_GLN
12	CTD1	A_579_PRO A_544_ASN A_564_GLN
13	S2	A_888_PHE A_789_TYR A_880_GLY
14	S2–HR1	A_1001_LEU A_1005_GLN A_759_PHE
15	HR1–HR2	A_1024_LEU A_1028_LYS A_1042_PHE
16	S2–HR1–HR2	A_1062_PHE A_1029_MET A_1033_VAL A_1053_PRO A_877_LEU
17	S2–HR1–HR2	A_905_ARG A_1050_MET A_898_PHE A_901_GLN
18	S2–HR1–HR2	A_819_GLU A_1055_SER A_874_THR
19	NTD	A_194_PHE A_238_PHE A_106_PHE A_117_LEU A_201_PHE A_235_ILE A_86_PHE A_231_ILE A_90_VAL
20	NTD	A_193_VAL A_204_TYR A_37_TYR
21	NTD–RBD–CTD1	A_328_ARG A_533_LEU A_578_ASP
22	RBD	A_338_PHE A_342_PHE A_368_LEU
23	RBD	A_353_TRP A_398_ASP A_464_PHE
24	RBD	A_365_TYR A_387_LEU A_515_PHE
25	RBD	A_406_GLU A_403_ARG A_495_TYR
26	RBD	A_509_ARG A_401_VAL A_451_TYR A_442_ASP A_438_SER A_507_PRO
27	CTD1	A_543_PHE A_585_LEU A_576_VAL
28	S2–HR1	A_977_LEU A_749_CYS A_993_ILE A_997_ILE
29	S2	A_878_LEU A_806_LEU A_882_ILE
30	HR1	A_973_ILE A_984_LEU A_992_GLN
31	S2–HR1	A_1000_ARG A_977_LEU A_996_LEU
32	S2–HR1	A_741_TYR A_1004_LEU A_962_LEU A_858_LEU
33	HR1–HR2	A_1065_VAL A_1052_PHE A_802_PHE A_927_PHE
34	NTD–RBD–CTD1	A_328_ARG A_543_PHE A_579_PRO
35	RBD	A_342_PHE A_511_VAL A_374_PHE A_436_TRP A_347_PHE
36	RBD	A_353_TRP A_400_PHE A_423_TYR A_512_VAL A_410_ILE
37	NTD	A_104_TRP A_194_PHE A_238_PHE A_92_PHE
38	NTD–CTD2	A_664_ILE A_312_ILE A_598_ILE

Table 5

Local Interacting Communities of the Protomer A in the Structure of SARS-CoV-2 S-GSAS/G614 in the Open State (pdb id 7KDL)

number	contributing structural domains	local community interacting residues
1	S2–HR1	A_1001_LEU A_1005_GLN A_759_PHE
2	HR1–HR2	A_1042_PHE A_1031_GLU A_1037_SER A_1039_ARG
3	S2–HR1–HR2	A_1067_TYR A_1049_LEU A_906_PHE A_909_ILE A_911_VAL
4	S2–HR1–HR2	A_1054_GLN A_816_SER A_819_GLU
5	S2–HR1–HR2	A_1075_PHE A_1110_TYR A_714_ILE
6	HR1–HR2	A_1102_TRP A_1081_ILE A_1135_ASN
7	NTD	A_186_PHE A_264_ALA A_66_HIS
8	NTD	A_189_LEU A_210_ILE A_217_PRO
9	NTD	A_220_PHE A_288_ALA A_36_VAL
10	NTD	A_265_TYR A_92_PHE A_240_THR
11	NTD	A_265_TYR A_65_PHE A_82_PRO
12	NTD	A_275_PHE A_290_ASP A_58_PHE
13	NTD	A_279_TYR A_44_ARG A_49_HIS
14	NTD	A_299_THR A_315_THR A_597_VAL
15	CTD2	A_664_ILE A_312_ILE A_598_ILE
16	RBD	A_342_PHE A_511_VAL A_374_PHE A_436_TRP A_347_PHE A_399_SER A_509_ARG
17	RBD	A_351_TYR A_452_LEU A_454_ARG A_492_LEU
18	NTD	A_91_TYR A_35_GLY A_56_LEU
19	RBD	A_464_PHE A_429_PHE A_425_LEU A_426_PRO A_514_SER
20	RBD	A_451_TYR A_401_VAL A_497_PHE A_448_ASN
21	RBD	A_454_ARG A_457_ARG A_467_ASP
22	CTD2	A_611_LEU A_666_ILE A_650_LEU
23	S2	A_888_PHE A_789_TYR A_880_GLY
24	HR1–HR2	A_1000_ARG A_977_LEU A_996_LEU
25	HR1–HR2	A_1024_LEU A_1028_LYS A_1042_PHE A_1032_CYS
26	S2–HR1–HR2	A_1062_PHE A_1029_MET A_1033_VAL A_1053_PRO A_877_LEU
27	HR1–HR2	A_1032_CYS A_1043_CYS A_1048_HIS A_1051_SER A_1064_HIS
28	S2–HR1–HR2	A_819_GLU A_1055_SER A_874_THR
29	HR1–HR2	A_1115_ILE A_1104_VAL A_1119_ASN
30	S2–HR1–HR2	A_1106_GLN A_1109_PHE A_915_VAL
31	NTD	A_193_VAL A_204_TYR A_37_TYR
32	NTD	A_104_TRP A_194_PHE A_238_PHE A_92_PHE A_267_VAL A_84_LEU
33	NTD	A_238_PHE A_106_PHE A_117_LEU A_201_PHE A_235_ILE A_86_PHE A_90_VAL
34	RBD–CTD1	A_579_PRO A_330_PRO A_544_ASN
35	RBD	A_338_PHE A_342_PHE A_368_LEU
36	RBD	A_353_TRP A_398_ASP A_423_TYR A_464_PHE
37	RBD	A_365_TYR A_387_LEU A_515_PHE
38	CTD1	A_543_PHE A_585_LEU A_576_VAL
39	S2–HR1	A_997_ILE A_749_CYS A_993_ILE
40		A_878_LEU A_806_LEU A_882_ILE
41	S2	A_906_PHE A_902_MET A_923_ILE A_916_LEU
42	S2–HR1	A_741_TYR A_1004_LEU A_962_LEU A_858_LEU
43	S2–HR1–HR2	A_905_ARG A_1050_MET A_898_PHE A_902_MET
44	HR1–HR2	A_1065_VAL A_1052_PHE A_802_PHE A_927_PHE
45	HR1–HR2	A_1081_ILE A_1115_ILE A_1137_VAL
46	NTD	A_285_ILE A_38_TYR A_279_TYR
47	NTD–RBD–CTD1	A_328_ARG A_543_PHE A_579_PRO
48	S2	A_898_PHE A_802_PHE A_797_PHE
49	RBD–CTD2	A_314_GLN A_596_SER A_613_GLN

The observed difference in the total number of local modules favoring the open mutant state is largely due to emerging communities in the NTD and RBD regions. Indeed, while the open form of the S-D614 spike protein is characterized by mobile NTDs that could provide flexible access to the ACE2 receptor, the increased stabilization of the NTD and RBD regions in the S-G614 mutant may accompany strengthening of the S1–S2 interdomain interactions and could contribute to a decrease in premature shedding of the S1 domain. Hence, the D614G mutation may also exert allosteric effect by partly immobilizing the distal NTD and RBD regions. These findings support the previously suggested notion that the D614G mutation in the SD2 domain could allosterically strengthen the stability of the distal NTD and RBD regions in the open state and therefore potentially promote exposure to the host receptor and greater infectivity. Our findings of the network analysis support the latest experimental data[48] by showing that the D614G mutation may strengthen the stability of the local community Q314–S596–Q613 and promote hydrogen-bonding interactions between Q613 and T859 of adjacent protomer afforded by the local backbone flexibility at the mutational site. The presented results also partly reconciled several scenarios offered to explain functional effects of the D614G mutation. Indeed, the community network analysis suggested that preferential stabilization of the open mutant form can be determined by two main factors: (a) through strengthening of the S1–S2 interactions and (b) by reducing functional movements of the NTD and RBD regions exposed to binding with the host receptor. Accordingly, the D614G mutation may exert its effect through allosteric stabilization of the S1–S2 interfaces limiting shedding of the S1 domain and by reducing flexibility and enhancing thermodynamic preferences of the open state. Hence, dynamic network modeling and community analysis of the S-D614 and S-G614 proteins revealed that the D614G mutation can induce a partial rearrangement of the residue interaction networks and promote the larger number of stable communities in both the closed and open forms by enhancing the S1–S2 interdomain interactions. Furthermore, the network analysis suggested a differential stabilization of the S-G614 mutant, favoring the open form in which strengthened allosteric couplings between mutational sites, S1–S2 regions, and NTD/RBD regions in S1 could contribute to a decrease in premature shedding of the S1 domain. Although recent studies indicated that the D614G variant did not itself drive to escape from antibody binding, it was found that D614G can remarkably potentiate escape mutations at some positions in certain patients, supporting an allosteric mechanism of action triggered by this mutation on dynamics and function in remote regions exposed to interactions with antibodies. In light of these newly emerging experimental data, our results argue that the D614G mutation may exert its global impact on other sites by acting as an important mediating center governing regulation of the SARS-CoV-2 machine.

Conclusions

We combined simulation-based approaches with dynamic network modeling and community analysis to quantify the effect of the D614G mutation on dynamics, stability, and network organization of the SARS-CoV-2 S proteins. The results of this study provide an insight into the molecular mechanisms underlying the effect of the D614G mutation by examining the SARS-CoV-2 S protein as an allosteric regulatory machine. The functional dynamics analysis suggested that the S-G614 trimer may promote an exchange between the open and closed forms, which could increase the exposure of RBDs to ACE2 binding but also force the return to the down conformation unless RBD is locked by the interactions with the host receptor. Using mutational sensitivity analysis of the SARS-CoV-2 S-D614 and S-G614 proteins, we demonstrated the improved stability of the D614G mutant as compared to the S-D614 protein. The differential stabilization of the S-G614 closed and open forms may help to reconcile mechanisms of mutation-induced preferences toward the open state and the reduced S1 shedding scenario underlying functional effects of the D614G circulating mutation. By examining the dynamic and network properties of the SARS-CoV-2 S trimer proteins, we characterized the distribution of allosteric hotspots in the S-D614 and S-G614 mutant structures revealing consolidation of the communication hotspots in the CTD1–CTD2 regions of the S-G614 mutant that may promote the greater stability and efficient allosteric couplings between S1 and S2 regions. Our results suggest that functional effects of the globally circulating D614G mutation may be determined by confluence of both local and global stability changes that collectively exert control over conformational transformations and differential thermodynamic stabilization of closed and open states. In the proposed mechanistic view, the enhanced infectivity of the S-G614 mutant could be determined by the interplay between the increased thermodynamic stability of the S-G614 trimer in the open state and the enhanced kinetic accessibility of S-G614 RBD regions for interactions with the host receptor. Examining functions of the SARS-CoV-2 pike proteins through the prism of an allosterically regulated system[88] may prove to be useful to uncover functional mechanisms and rationalize the growing body of diverse experimental data.

Methods

Structure Preparation and MD Simulations

All structures were obtained from the Protein Data Bank.[89] Hydrogen atoms and missing residues were added using the WHATIF program web interface.[90] The missing loops in the studied cryo-EM structures of the SARS-CoV-2 S protein were reconstructed and optimized using the template-based loop prediction approach ArchPRED server.[91] The side-chain rotamers were reconstructed and optimized with the aid of the SCWRL4 tool.[92] The protein structures were also subjected to an additional refinement and global optimization using the 3Drefine method.[93] Structurally resolved glycan residues in the structures of the SARS-CoV-2 S proteins were checked and added. In addition, the glycosylated microenvironment for simulation trajectories was enriched by considering 16 most occupied N-glycans in each protomer as determined in the cryo-EM structures of the SARS-CoV-2 spike S trimer in the closed state (pdb id 6VXX) and open state (pdb id 6VYB).

The structures of the SARS-CoV-2 S proteins were simulated in a large box size of 85 Å with a buffering distance of 12 Å. Sodium (Na+) and chloride (Cl−) counterions were added to achieve charge neutrality and placed at least 8 Å away from any protein atoms and from each other. All-atom MD simulations were performed for an N, P, T ensemble in explicit solvent using the NAMD 2.13 package[94] with the CHARMM36 force field.[95] The van der Waals interactions were computed using an atom-based cutoff of 12 Å. Long-range electrostatic interactions were calculated using the particle mesh Ewald method.[96] MD simulations were executed using a 2 fs integration time step. Energy minimization after addition of solvent and ions was carried out using the steepest descent method for 100,000 steps. An NPT production simulation was run on the equilibrated structures for 500 ns, keeping the temperature at 310 K and constant pressure (1 atm). In simulations, the Nose–Hoover thermostat[97] and isotropic Martyna–Tobias–Klein barostat[98] were used to maintain the temperature at 310 K and pressure at 1 atm, respectively. PCA of simulation trajectories was based on the set of backbone heavy atoms (N, Cα, Cβ, C, O) to determine the essential dynamics of the chaperone structures. PCA computations were conducted using the CARMA package.[99]

Mutational Scanning and Sensitivity Mapping

Mutational scanning of protein residues was performed using the BeAtMuSiC approach.[100−102] If a free-energy change between a mutant and the wild-type (WT) proteins ΔΔG = ΔG (MT – ΔG (WT)) > 0, the mutation is destabilizing, while when ΔΔG < 0, the respective mutation is stabilizing. The protein stability changes were averaged using conformational ensembles obtained from MD trajectories. A total of 10 000 samples from equilibrium simulation trajectories were used for rapid evaluation of protein stability changes.

Modeling of Residue Interaction Networks

The residue interaction networks are constructed using dynamic cross-correlations[86] and coevolutionary residue couplings[87] that can adequately describe the allosteric interactions, as detailed in our previous studies.[67] The Floyd–Warshall algorithm is employed for computation of shortest paths.[103] Network graph calculations were performed using NetworkX[104] and the Girvan–Newman algorithm;[105−109] the topological parameters were computed using NetworkX[104] and Cytoscape software.[110,111] The network community analysis is used to evaluate protein stability and allosteric communications in the SARS-CoV-2 S-D614 and S-G614 mutant structures.[112]