The transmissible spongiform encephalopathies (TSEs), or prion diseases, are a group of fatal neurodegenerative disorders including scrapie in sheep,^[1,2] chronic wasting disease in deer,^[3] spongiform encephalopathy in bovine,^[4] Creutzfeldt-Jakob disease (CJD),^[5,6] fatal familial insomnia (FII),^[7,8] and kuru^[9,10] in human. The “prion hypothesis”, in other words, the TSEs are caused by the misfolding and aggregation of prion proteins (PrP), has dominated the field since 1982.^[11–16] This hypothesis has also been extended to cover other neurodegenerations such as Alzheimer’s disease, amyotrophic lateral sclerosis (ALS), and Parkinson’s disease, which are respectively associated with the aggregation of beta amyloids, TDP-43, and α-synucleins.^[17–23]

The key hallmark underlying the pathological misfolding and aggregation of PrPs is the conformational transition from a helix-rich conformer (PrP^C, in which C stands for cellular) to a sheet-rich pathogenic (or scrapie-like) conformation (PrP^Sc).^[24–26] The PrP^C exists as a glycoprotein composed of 209 residues including a highly disordered N-terminal segment and a globular folded region comprising residues 125–231. NMR^[27–29] and x-ray^[30] studies have shown that the folded domain of human PrP^C contains three α-helices (residues 144–154, 173–194, and 200–228, denoted respectively H1, H2, and H3) and two short antiparallel β-strands (residues 128–131 and 161–164, denoted respectively S1 and S2) (Fig. 1(a)). Although numerous works have suggested that PrP^C is related to a wide range of different cellular processes,^[31–36] its biological function remains largely unknown. On the other hand, the PrP^Sc isoform has characteristics of enriched β-sheet structures, proteases-resistance, and a high propensity to aggregate into amyloid fibrils.^[37–39] However, solving the atomic-level structure of PrP^Sc is very challenging, possibly due to its insolubility, short survival time, and fast aggregation into toxic fibrils.^[40–42]

About 10% to 15% of human prion diseases are caused by mutations of the PrP.^[43] These mutations lead to the spontaneous misfolding of the protein and the formation of neurotoxic protofibrils and fibrils. Examples of the more than twenty reported disease-related mutations include D178N,^[44,45] Q217R,^[46] and T188R,^[47] related respectively to FFI, Gerstmann–Sträussler–Scheinker syndrome (GSS), and CJD. On the contrary, a few mutations (M129V polymorphism,^[48] E219K polymorphism,^[49,50] E200K,^[51] and G127V) are reported to have protective effect on prion diseases. For example, Mead and coworkers reported a novel acquired prion disease resistance factor (G127V mutation) selected during the kuru epidemic in Papua New Guinea.^[52] Asante et al. further performed transgenic mice experiments and showed that the V127 variant completely prevents prion diseases and that its mechanism is distinct from the M129V polymorphism.^[53] In addition, Sabareesan and Udgaonkar showed that the G126V mutation in mouse PrP (equivalent to G127V in human PrP) can slow the initial fibril growth by extending the lag phase.^[54] In spite of these experimental advances, the molecular mechanisms underlying the complete resistance of the V127 variant to prion disease remain mostly unknown.

Molecular dynamics (MD) simulations have been widely used to investigate the structural and dynamical properties of proteins.^[55–57] In 2016, Yao et al. reported the first computational study on PrP dimer and its V127 variant. Their 100-ns simulations showed that the constructed PrP dimer with G127V mutation is less stable than the WT dimer.^[58] Recently, Lin et al. solved the NMR structures of monomeric WT PrP and the G127V mutant and found G127V to be less stable than WT. They also performed 100-ns MD simulations starting from their NMR structures and showed that S1 and S2 in G127V are less populated than those in WT.^[59] However, due to the short time scales of their simulations, the atomic-level mechanism underlying the G127V-mutation-induced protective effect of PrP is not well understood. Herein, we have carried out six independent 2-μs all-atom MD simulations on WT PrP monomer and its V127 variant. Our simulations show that the G127V mutation stabilizes the PrP monomeric structure by greatly enhancing the structural rigidity of the S2–H2 loop and the stability of the C-terminal of H2. The increased rigidity of the S2–H2 loop and the enhanced H2 stability in V127 variant may inhibit the structural conversion from human PrP^C to PrP^Sc, thus prevent prion disease. Our hypothesis is supported by the fact that prions are poorly transmissible to animals with their PrP^C carrying a rigid loop and a recent finding showing that stabilization of H2 of prion protein prevents its misfolding and oligomerization. Interaction analysis and optimal allosteric path analysis provide the atomic-level mechanism of G127V in the enhancement of S2–H2 loop rigidity and the H2 C-terminal stability. Community network analysis further shows that G127V mutation also increases the global correlations of PrP.

We have performed three individual 2-μs-long simulations for both the globular domain of WT PrP (comprising residues 125–228) and its V127 variant. The initial structure of WT PrP is taken from a solution NMR structure,^[27] while the G127V mutation is obtained by replacing the sidechain of G127 with that of valine. For each system, we first calculate the time evolution of backbone root mean square deviations (RMSDs) with respect to the native WT NMR structure. The RMSD values of WT in the three simulations oscillate more largely than those of G127V. For example, RMSD values in WT-1 run experience a slow and continuous increase from 0.15 nm to 0.40 nm in the first 1.5 μs, followed by a sudden drop after which the RMSD values fluctuate around their equilibrium value of ∼ 0.35 nm (Fig. S1(a)). In contrast, the RMSD values in G127V-1 run fluctuate around ∼ 0.25 nm during the full period of the simulation (Fig. S1(d)). The probability density functions (PDF) of the WT and G127V RMSD values in four consecutive time windows are shown in Figs. 1(b) and 1(c). A shift from small RMSD values to larger ones can be clearly seen in the PDF of WT RMSD (Fig. 1(a)), while a similar G127V RMSD distribution is observed for all four time windows (a main peak centering 0.24 nm and a satellite peak centering 0.33 nm). These results show that the V127 variant has a relatively higher structural stability than WT.

We further investigate the influence of G127V mutation on the structural rigidity of PrP by calculating the root mean square fluctuation (RMSF) of both WT and G127V systems. As shown in Fig. 1(d), the RMSF values of almost all residues in G127V are lower than those of WT (except for a few residues on the S1–H1 loop). In the WT system, the residues in the S2–H2 loop (region 1, residues 164–174) and the residues in the C-terminal of H2 and those in H2–H3 loop (region 2, residues 192–197) have much higher RMSF values than other non-terminal residues, indicating their higher flexibilities than the other residues. The high flexibilities of these two regions are consistent with previous atomistic and coarse-grained simulation studies on WT PrP.^[60–64] It is noted that the RMSF values in our work are much higher than those reported in the previous all-atom simulation studies, probably due to the much longer simulation time in our work than previous studies (2 ns,^[60] 10 ns,^[61] 1.8 ns,^[65] 80 ns,^[62] 50 ns,^[63] 100 ns^[58,59]). To the best of our knowledge, our work is the first microsecond-scale all-atom MD simulation study on prion protein. However, the RMSF values of residues in these two regions of the G127V systems are remarkably lower than those of the WT systems. These findings demonstrate that the G127V mutation increases the rigidity of the S2–H2 loop, and enhances the stability of the C-terminal of H2.

The misfolding of PrP monomer and the consequent oligomerization/fibrilization are closely related to the prion disease pathology.^[26] The increased rigidity of S2–H2 loop as well as the enhanced stability of H2 C-terminal may stabilize the monomeric PrP^C structure, thus inhibiting the prion conversion of human PrP^C to PrP^Sc and finally preventing prion disease. Our interpretation is well supported by the following experimental results. Firstly, the relationship between the rigidity of the S2–H2 loop and the PrP’s resistance to prion diseases has been extensively discussed in previous studies.^[66–69] Accumulating evidence shows that a transmission barrier exists between prion proteins carrying rigid S2–H2 loop and those carrying flexible one. And mice with a flexible S2–H2 loop is much more easily infected by prions than rigid S2–H2 loop mice.^[70] Animals (such as bank vole, elk, horse, and rabbit) with a rigid S2–H2 loop PrP usually do not develop prion diseases.^[71–75] Secondly, the C-terminal region of H2 is identified as the first region to experience conformational conversion of PrP^C to PrP^Sc[76–78] and the stabilization of PrP H2 may prevent the misfolding and oligomerization of PrP.^[79]

The solution-state NMR structures of bank vole, elk, and horse PrPs have been reported (PDBID: 2K56, 1XYW, 2KU4).^[71,72,74] Figure 2(a) shows the cartoon representations of their S2–H2 loops (residues 165–174). Surprisingly, they have very similar conformations, comprising a helical N-terminal (residues 165–168), and a turn-like region (residues 169–174), as highlighted respectively in magenta and blue boxes. The helical nature of residues 165–168 can be further characterized by calculating the (Φ, Ψ) values of these residues. The (Φ, Ψ) values of D167 are given in Fig. 2(b), and they are distributed around (–80°, –10°), showing a propensity towards α-helical structure. Similar results are observed for residue M166 (Fig. S2).

To show whether the S2–H2 loop in V127 variant has a conformation similar to that in bank vole, elk, and horse PrPs, we perform cluster analysis on their conformations in the last 1 μs of the trajectories using a single-linkage algorithm^[80] with an RMSD cutoff of 0.05 nm. As shown in Fig. 2(c), the S2–H2 loop conformations in WT-1 run are clustered into 270 clusters, suggesting a wide variety of loop conformations. We show in Fig. 2(c) the representative configurations of the top four largest clusters. Interestingly, no helical N-terminal region is observed in the first and the second clusters. Although the third and the fourth clusters show a certain propensity towards forming helical-like structures, they only comprise respectively 7.6% and 7.4% of all conformations. Analyses of the other two WT systems give similar results (Fig. S3). In sharp contrast, for all of the three G127V MD runs, almost all S2–H2 loop conformations (containing respectively 97.5%, 99.7%, and 98.7% of conformations in G127V-1/2/3 system) are clustered into a single cluster, and all of them adopt a conformation similar to that of S2–H2 loop in bank vole, elk, and horse PrPs. We show in Figs. 2(d)–2(f) the representative conformations of the S2–H2 loop in the largest cluster with the helical-like N-terminal region highlighted in magenta boxes. This propensity towards helical structure is further characterized by calculating the potential mean force (PMF) of residue D167 in WT and in G127V PrPs as a function of its (Φ, Ψ) values. As shown in Fig. 2(g), the PMF of WT D167 has three energy minimum basins. Two of them are located at (Φ, Ψ) values of (–80°, –10°) and (–160°, –10°), while the third one is located at a (Φ, Ψ) value of (–160°, ± 180°). The deepest one at (–80°, –10°) corresponds to α-helical structures. In contrast, the D167 PMF of the G127V system presents only two basins. The α-helical basin of the G127V system is deeper and wider than that of the WT system, suggesting that the S2–H2 loop of the G127V system has a higher preference of forming helical N-terminal structures than that of WT. These results, together with the fact that prions are poorly transmissible to animals with their PrP^C carrying a rigid loop, provide explanation to the strong resistance of V127 variant against prion diseases.

Although residues in the S2–H2 loop (residues 164–174) are distant from the mutation site (V127) in sequence, they are spatially close to each other (see Fig. 1(a)). We perform interaction analysis to understand how the G127V mutation improves the rigidity of the S2–H2 loop. Valine has a bulky hydrophobic side chain while glycine has no side chain, so we suspect that the hydrophobic interactions may play an important role. By calculating the contact probability between G/V127 and all residues in the S2–H2 loop, we find that the hydrophobic residue P165 has the highest contact probability with V127 (Fig. S5). The time evolution of P165–G/V127 contact number in WT-1 and G127V-1 runs shows that the P165–G127 contact number fluctuates around 2 throughout the full period of simulation, while the P165-V127 contact number increases from ∼ 2 to ∼ 6 within the first 0.2 μs and fluctuates around 4 in the remaining simulation time (Fig. 3(a)). Similar result is seen from the other four simulations (WT-2/3, and G127V-2,3), and the contact number PDF calculated using a combined trajectory of the three individual simulations (Figs. S6(b), S6(c), and 3(b)). These results suggest that the hydrophobic interaction between P165 and G/V127 is greatly enhanced as a result of the G127V mutation. Interestingly, a recent work demonstrated that substituting G127 into I127 which possesses higher hydrophobicity than both glycine and valine also significantly decreases the PrP cytotoxicity, indicating the important role of the V127-involved hydrophobic interaction in preventing prion diseases.^[81]

The enhanced V127–P165 hydrophobic interaction may facilitate other residues in PrP fragment ₁₂₅LGVYM₁₂₉ to interact with the S2–H2 loop (₁₆₂YYRPMDEY₁₆₉). We calculate the H-bond number between them. As shown in Fig. 3(c), the H-bond number in G127V-1 run increases with simulation time and reaches ∼ 5, while that in WT-1 run fluctuates around 3. Similar results are observed in WT-2/3 and G127V-2/3 systems (Figs. S6(e)–S6(f)), and a statistical analysis on all simulations for WT and G127V systems shows that G127V has a higher H-bond number than WT (Fig. 3(d)). A detailed analysis of the H-bonds shows that residue E168 plays an important role on the H-bond formation: it forms stable H-bonds with three residues (L125, G126, and G/V127), and these H-bonds are much more stable in G127V system than in WT system (Fig. 3(e)). In addition, the M129-Y163 H-bond is stable in both WT and G127V systems. The M129-Y163, and the E165-L125/G126/V127 H-bonds are shown in Figs. 3(h)–3(j). The enhancement of the E165-involved H-bonds may also contribute to the rigidity of the S2–H2 loop.

It has been reported that the R164–D178 salt bridge plays an important role on the structural stability of PrP.^[82,83] To explore the effect of G127V mutation on the stability of this salt bridge, we calculate the distance between the charged groups of sidechains in R164 and D178 in WT-1 and G127V-1 runs. Figure 3(f) shows that the R164-D178 distances in both systems are close to ∼ 0.4 nm in most of the simulation time, suggesting the existence of R164–D178 salt bridge in both systems. In addition, the distance in WT-1 run has a higher fluctuation than that in G127V-1. Analyses of the WT-2/3 and G127V-2/3 systems are given in Figs. S6(h)–S6(i), and a statistical analysis using all three simulations is given in Fig. 3(g). The PDF curve of R164–D178 distance in G127V systems has a single peak centering 0.38 nm, while that in WT systems has a dominant peak around 0.38 nm and a satellite peak centering at 0.45 nm. These data indicate that the R164–D178 salt bridge in G127V systems is slightly enhanced compared to that in WT.

Taken together, the G127V mutation enhances the V127–P165 hydrophobic interaction, the E165-L125/G126/V127 H-bonds, and the R164–D178 salt bridge. These interactions collectively rigidify the S2–H2 loop of PrP.

As mentioned above, G127V mutation also enhances the stability of H2 C-terminal region (residues 192–197). This region is far from the mutation site (G/V127) both in sequence and in space, suggesting that G127V improves the H2 C-terminal stability through an allosteric effect. The optimal dynamical path analysis method has been demonstrated as a useful tool for identifying allosteric paths in various proteins.^[84–88] Here we generate dynamic networks for WT and G127V systems by adding edges between residue pairs with a relatively high contact probability ( > 70%), and constructing the adjacency matrices using the absolute value of the inter-residue correlations (Fig. S4). We then calculate the optimal and suboptimal paths between residues G/V127 and G195 (the center residue in 192–197 region) to unravel how G127V mutation increases the stability of this distant region.

In WT PrP system, we identify two shortest network paths from G127 to G195 (path length = 317) as well as 15 paths with slightly longer path lengths (between 317 and 327). As shown in Table S1, all of these paths are similar: they start from G127, propagate to residue M129 on S1, and to V161 on S2, and subsequently to C terminus of H3, then pass to its N-terminal through nodes in H3, and finally reach residue G195. We present in Fig. 4(a) a snapshot showing residues located in the shortest path, and in Fig. 4(b) the correlation values of all edges alongside the path. As expected, all edges have high inter-residue correlation values ( > 0.6), and the edges within H3 have the highest correlation values, probably due to the existence of the stable H3 helix. Most of the edges (especially those located on the helices) connect two residues which are close in sequence, expect for the S1–S2 (M129–V161) and the S2–H3 (V161–V210) edges. These two edges are stabilized respectively by H-bonds between S1 and S2, and the V161-V210 hydrophobic interactions. To identify the important residues in this allosteric path, we remove the node corresponding to each residue and examine its effect on the optimal path length by calculating the path lengths after the residue removal. As shown in Fig. 4(c), removal of residue M129 results in the largest increase of optimal path (14%), revealing its highest importance in the G127–G195 allosteric path. In contrast, the optimal allosteric path from V127 to G195 in G127V system is much shorter than that in WT system (path length = 257). Similar to the WT system, the V127–G195 allosteric paths start from V127, propagate to residue M129 on S1, and to Y162 on S2. However, after that, instead of propagating to C-terminal of H3, the allosteric paths propagate to the N-terminal of H2, and pass to its C-terminal through nodes in H2, and finally reach residue G195 (Table S2). We show in Fig. 4(d) all the residues on the shortest path, and in Fig. 4(e) the correlation values of all edges alongside the path in both WT and G127V systems. Except for the V127–M129 edge, most edges in G127V system have a higher correlation value than those in WT system, suggesting that the correlations in G127V system are much stronger than those in WT system. The main difference between the G127V allosteric path and the WT path is that S2 connects to H2 in G127V, while it connects to H3 in WT. As residues R164 and D178 are located respectively in S2 and H2, the S2–H2 correlation in G127V system may result from the enhanced R164–D178 salt bridge. The importance of the S2–H2 correlation can be revealed by the alternation of optimal path length upon node removal. As shown in Fig. 4(f), the removal of the node corresponding to residue Y162 results in the largest increase in the optimal path length (4%). More interestingly, removing the S2–H2 interactions by cutting all edges connecting Y162 and H2 (Y162-C192/I182/T183) leads to a suboptimal V127–G195 path (Fig. 4(d) and Table S3) that is similar to the G127–G195 optimal path in WT system, but with a shorter path length (282). These results provide insights into the allosteric effect by which G127V mutation enhances the stability of the remote region (H2 C-terminal). The optimal allosteric path from G127 to H2 C-terminal in WT PrP passes edges connecting S2 and H3. After the G127V mutation, the interaction between S2 and H2 N-terminal is enhanced (through stabilized R164–D178 salt bridge), thus generates a new optimal allosteric path with a much shorter path length. This path passes through edges connecting S2 and H2. In addition, a path same to the WT optimal path also exists, but becomes the suboptimal path in G127V. The new-generated optimal path and the shortened suboptimal path from V127 to G195 explain the increased long-range correlation between the mutation site and the C-terminal of H2, thus stabilize the structure of H2 C-terminal.

We further investigate the effect of G127V mutation on the global correlation of PrP by calculating the correlation values between each two residues. As shown in Fig. 5(a), the WT PrP correlation map contains three diagonal positively-correlated regions corresponding to the inter-residue interactions in the three helices (red circles), an off-diagonal region corresponding to the anti-parallel H2–H3 interaction (blue circle), and two regions corresponding to H2/H3–S2 interactions (green circles). The G127V correlation map has a similar pattern to that of WT, but the correlation values of the six positively-correlated regions are mostly higher than those of WT (Fig. 5(b)). In addition, more positively-correlated regions are observed in the correlation map of G127V system. These results demonstrate that the inter-residue correlations in PrP significantly increase after G127V mutation.

Community network method has become a useful strategy for analyzing global correlations for biomolecules.^[89–93] Based on the dynamical networks of WT and G127V, we calculated the optimal community distributions using the Girvan–Newman algorithm.^[94] The residues belonging to the same community are more strongly and densely interconnected to each other, and have weaker connections to other residues in the protein. As shown in Fig. 5(c), the WT community network comprises ten communities, including a single community for H1, two separate communities for N and C terminals of H2 (H2N and H2C), three separate communities for N-terminal, middle part, and C-terminal of H3 (H3N, H3M, and H3C), and one single community for the S2–H2 loop. In sharp contrast, the G127V community network comprises only six communities (Fig. 5(d)). The three WT communities corresponding to S1/S2, H2N, and H3M have combined into one large community, and the two WT communities corresponding to H3N and H2C have also combined into one large community. In addition, the inter-community connectivity in G127V is stronger than that in WT (the width of edges connecting communities corresponds to the inter-community connectivity strength). The decrease in cluster number, the combination of several clusters into single large cluster, as well as the increase of inter-community connectivity strength indicate that the G127V mutation enhances the global correlations and interactions in PrP, thus stabilizes the PrP^C structure and inhibits the prion conversion to PrP^Sc.

We have performed six 2-μs-long MD simulations on the globular domain of WT human PrP and its V127 variant to investigate the effect of G127V mutation on the structural and dynamical properties of monomeric human PrP, with the aim of understanding the mechanism of G127V’s disease-prevention. RMSD and RMSF analyses show that G127V mutation increases the structural stability and rigidity of PrP. Especially, (1) the rigidity of the S2–H2 loop is greatly increased, and (2) the stability of the C-terminal of H2 is enhanced. The increased rigidity and stability of these two regions may inhibit the prion conversion of PrP^C into toxic PrP^Sc, thus prevents prion diseases. Our speculation is supported by the fact that prions are poorly transmissible to animals with rigid S2–H2 loop and recent findings showing that stabilization of H2 of prion protein prevents its misfolding and oligomerization. Interaction analysis shows that the increased rigidity of the S2–H2 loop results from the enhanced V127–P165 hydrophobic interaction, E165-L125/G126/V127 H-bonds, and R164–D178 salt bridge. The mutation-induced stabilization of the R164–D178 salt-bridge enables the formation of a new allosteric path from the mutation site to the H2 C-terminal with a much shorter path length than the optimal path in WT. The WT optimal path still exists in G127V and becomes the suboptimal path. The new optimal path and the shortened suboptimal path lead to an enhanced stability of the H2 C-terminal. At last, community network analysis shows that the global correlations and interactions in PrP are strengthened by the G127V mutation. Our findings provide structural and dynamical basis for understanding the role of human V127 variant in preventing prion conversion and propagation, which may be helpful for the rational design of potent anti-prion therapies.

The initial structure of WT PrP globular domain was taken from a solution NMR structure comprising residues 125–228 (PDBID: 1HJM).^[27] We built the G127V structure by replacing the sidechain of glycine with that of valine using tools implemented in the Pymol package.^[95] Three individual 2-μs-long MD simulations were performed for each system using the AMBER99SB-ILDN force field^[96] in combination with the TIP-3P water model. All simulations were preformed using the GROMACS-5.1.2 package^[97] in the NPT ensemble. Electrostatic interactions were treated using the Particle Mesh Ewald (PME) method^[98] with a real space cut-off of 1.2 nm. The vdW interactions were calculated using a cutoff of 1.2 nm. The solute and solvent were separately coupled to an external temperature bath using a velocity rescaling method^[99] and a pressure bath using the Parrinello–Rahman method.^[100] The temperature and pressure were maintained at 310 K and 1 bar with a coupling constant of 0.1 ps, respectively. The neighbor-list was updated every 10 steps with a cut-off distance of 1.2 nm using a Verlet buffer.^[101]

Trajectory analysis was performed using our in-house-developed codes and the facilities implemented in the GROMACS-5.1.2 software package.^[97] We discarded the first 1 μs data for each MD trajectory to remove the bias of the initial states. Therefore, the structural and dynamical properties of each system were based on the simulation data generated in the last 1 μs. The cluster analysis for the S2–H2 loop was performed using a single linkage method^[80] with an RMSD cutoff of 0.05 nm. An inter-residue contact was defined when the aliphatic carbon atoms of two nonsequential side chains (or main chains) come within 0.54 nm or any other atoms of two nonsequential side chains (or main chains) lie within 0.46 nm. One H-bond is taken as formed if the N … O distance is less than 0.35 nm and the N–H … O angle is greater than 150°. Graphical analysis and visualizations were performed using the Pymol package.^[95]

Dynamical networks of WT and G127V were generated using our in-house-developed codes. Each C_α atom was selected as a representative node for the corresponding residue. Edges were added between two nodes if the corresponding residues are in contact during a majority of the simulation time (> 70%). Nearest neighbors in sequence are not considered to be in contact as they lead to a number of trivial allosteric paths in the weighted network. The weight of each edge was defined as –log | C_ij|, where C_ij stands for the dynamic cross correlation of two nodes (i and j), so that the network distance between two nodes connected by an edge decreases as the correlation of the two corresponding residues increases.

On the basis of the dynamical networks, the allosteric signal transmission from the mutation site (G/V127) to the C-terminal of H2 (residues 192–197) was analyzed by calculating the allosteric paths between residue G/V127 and residue G195. The length of a path D_ij between distant residues i and j was defined as the sum of the edge weights between the consecutive nodes k, l along the path: D_ij = Σ_k,lw_kl. The path between residues i and j with the shortest distance Dij(0) was found by using the Floyd–Warshall algorithm.^[102] Because the high correlation values between two nodes may not necessarily guarantee strong physical interactions between the corresponding residues (it could also result from strong interaction between nearby residue pairs), the shortest path as well as the paths whose lengths lie within a certain limit of the shortest distance were equally important.^[88] We defined the “optimal allosteric path” between residues i and j as the shortest paths between residues i and j (path length =Dij(0)) as well as the slightly longer paths (path length <Dij(0)+10). The importance of each residue in the allosteric path was evaluated by calculating the change in optimal path length upon removing the node corresponding to this residue while keeping all other contacts in the network intact. Although the optimal allosteric path is the most dominant mode of communication between the two strongly-correlated remote residues, other paths with longer path length (corresponding to weak correlations/interactions) may also play a role. However, the optimal path analysis method^[88] cannot identify these longer paths. Thus, we identified the most important residue in the allosteric effect by the aforementioned residue-removal technique, removed all edges connecting to this residue (these edges correspond to the strongest interaction along the optimal path), and calculated the path between residues i and j with the shortest distance Dij(1) in this modified network. By removing these edges, the network path from the mutation site to its strongly-correlated remote residue has to pass through edges corresponding to weak interactions. The “suboptimal allosteric path” between residues i and j was thus defined as the path with path length Dij(1) as well as the slightly longer paths (path length <Dij(1)+10). Identification of both optimal and suboptimal allosteric paths enables us to reveal the mechanism of the mutation-induced enhancement of global structural stability.

We calculated the shortest network path between each residue pair (i, j, for i ≠ j). The betweenness of each edge was then defined as the number of shortest paths that cross that edge. The optimal community distribution was calculated using the Girvan–Newman algorithm,^[94] which iteratively removed the edge with the highest betweenness and recalculated the betweenness of all remaining edges until the modularity of the community network was maximized. The modularity is a measure of the quality of a particular division of a network, and the bigger the modularity, the better the division quality.^[103] The optimal paths and dynamical networks were obtained using codes developed by Eargle and Sethi^[88] with our modifications. Visualization of the community networks was performed using the VMD package^[104] (snapshots) and the Omni Graffle package (schematic diagrams).