Proteins are important bio-macromolecules because they perform a vast variety of functions and involve virtually every process within cells. However, proteins are not static entities but undergo complex internal motions, which are essential for their functions [1
]. Therefore, a deeper understanding of protein functions needs to investigate their dynamic behaviors.
In a protein–solvent system or under physiological conditions, protein molecules are immerged in an aqueous solvent and thus the dynamic behaviors of the water molecules and the protein–solvent interactions are fundamental to the structural, dynamic, and functional properties of the proteins [3
], including facilitating protein folding [8
], maintaining structural stability/flexibility [10
], mediating protein–ligand binding [13
], accelerating enzymatic catalysis [15
], and inducing protein glass transition [17
The FEL, which is defined as a function describing the conformational states/substates of a protein with the corresponding Gibbs free energy of the protein–solvent system, is a seminal concept in characterizing and explaining many aspects of the protein dynamics and functions [2
]. Due to an absolute advantage of the water molecules in terms of both quantity and mass compared to the solute molecules, the changes in entropy and enthalpy of the solvent will make a substantial contribution to the change in Gibbs free energy of the protein–solvent system [22
], and the protein–solvent interactions also play a crucial role in shaping the FEL [26
]. In the first fast stage of protein folding, i.e.
, the hydrophobic collapse that results in the molten globule intermediate, the solvent entropy maximization contributes substantially to the lowering of the system free energy [27
]. In the following bottleneck stage, the formation of water network on the protein surface and the formation of HBs between the protein and water also contribute to lowering the system free energy via an enthalpy decrease [22
]. The protein native state residing at the bottom (the global free energy minimum region) of the funnel-like FEL (or folding funnel) is not a single, rigid conformer but exists as an ensemble of statistically populated conformational states/substates [31
]. The hierarchically organized free energy wells (i.e.
, an organization that the low-tier wells reside within the high-tier wells, and the lower tier the wells are, the lower are barrier heights between these wells [33
]) determines the hierarchy of protein dynamics—different structural components feature differentiated amplitudes and timescales of the fluctuations [2
]. Of particular note is that the solvent is able to influence protein dynamics at various hierarchical levels. It has been shown that the fast residue side chain rotations occurring on picosecond (ps) timescale are coupled to the hydration-shell fluctuations, and that the relatively slow loop motions occurring on nanosecond (ns) timescale are slaved to bulk-solvent motions [4
]. Furthermore, the complete exchange of protein-bound water molecules arising from water translational motions have been proposed to be necessary for protein structural relaxation that allows anharmonic and collective motions in proteins [36
]. Because the hierarchical arrangement of free energy wells determines that the higher-tier protein motions result from the lower-tier fluctuations [34
], it is reasonable to postulate that the protein collective motions, which involve the entire protein structure and occur on the timescales of microsecond (μs) to millisecond (ms), are a consequence of the accumulation of loop motions and side chain rotations, in which the solvent motions play an important role [22
The mobility of the atoms, and further, of the atomic groups and the molecules, is temperature dependent because the temperature is a fundamental determinant of the atomic thermal energy. A generally observed phenomenon is that the protein conformational flexibility increases/decreases with the elevated/decreased temperatures of the protein–solvent system. The fact that solvent fluctuations affect protein dynamics at various hierarchical levels raises a question of whether the temperature of the solvent, or that of the protein itself, is a primary factor in determining the internal motions of the protein. The second law of thermodynamics determines that the thermal energy always distributes as evenly as possible over the entire protein–solvent system, thus making it difficult to distinguish experimentally the effect of the solvent temperatures from that of the protein temperatures on protein dynamics [41
]. The method of MD simulation can be a valuable alternative, since it can create a thermodynamic system with the protein at one temperature and the solvent at a different temperature [41
]. Using such a simulation procedure, Vitkup et al.
] have shown that, in the temperature range 180–300 K, the elevated solvent temperature, and hence the increased solvent mobility, plays an essential role in enhancing the amplitude of atomic fluctuations in the protein, but that the temperatures of the protein itself only weakly affect the amplitude of protein atomic fluctuations. Despite this pioneering work, the detailed information about effects of the solvent temperatures on the molecular motions and conformational space sampling of the protein, as well as on the FEL of the protein–solvent system, remains unclear. Moreover, the question of how the solvent mobility dictates the hierarchical dynamics, from the low-tier group fluctuations to the high-tier collective motions of the protein, remains unanswered.
In order to resolve the above issues, in the present work the proteinase K (which is a subtilisin-like serine protease from Tritirachium album
limber with a broad-spectrum degradation capability to degrade proteins [43
]) in explicit solvent was subject to a series of MD simulations with temperatures of the protein and the solvent coupled to either 180 or 300 K. Proteinase K was chosen because it is a medium-sized globular α/β protein (279 amino acids) with an available high-resolution X-ray crystal structure [45
] and, in particular, its dynamic properties, molecular motions, and structure–dynamics–function relationship have been investigated by us in a series of simulation/modeling studies [46
]. The long, multiple-replica MD simulations with different initial atomic velocities [53
] were performed on each system of different protein/solvent temperature combinations to guarantee a more effective sampling of the conformational space accessible to proteinase K. Apart from comparative analyses of several structural properties, the differences in the sampled protein conformational spaces, collective motions, and constructed FELs between different temperature combinations were compared. The reasons for the enhanced mobility/flexibility of proteinase K caused by the elevated solvent temperature were discussed in depth based on the competitive HB interactions between the protein and solvent and between atoms within the protein. Finally, a refined FEL model was proposed to explain the role of the solvent in determining the hierarchical dynamics of proteins.
2.1. Evaluation of the Conformational Sampling
The time evolution of backbone RMSD relative to the starting structure was computed for each replica of the four simulation systems. The results (Supplementary Figure S1
) show that for the systems P180/S180 and P300/S180, proteinase K requires only a few ps to reach stable RMSD values, and for the systems P180/S300 and P300/P300, the protein requires about 600 ps and about 1 ns to reach an approximate equilibrium, respectively. In order to guarantee the consistency of our data over time among the four simulation systems, the first 1-ns trajectories were discarded and only the equilibrated portions (1–15 ns) of each replica were concatenated, thus yielding four 84-ns single joined trajectories.
The cosine contents for the first three eigenvectors derived from ED analyses of the replicas and of the joined trajectories were calculated to assess the conformational sampling. The results (Supplementary Table S1
) show that, for the independent replicas of each simulation system, most of their eigenvectors have significantly higher cosine contents (close to 1) than the corresponding eigenvectors of the single joined trajectory (~0.006 to ~0.2). This suggests that the regions of conformational space sampled by these independent replicas are displaced in different directions from the starting structure, resulting in different but partially overlapping samplings. Our strategy of performing the multiple, long simulations therefore improves the ergodicity of the system and achieves a relatively higher degree of sampling convergence in the joined trajectory than in the individual replicas. Consequently, the single joined trajectories of the four simulation systems were used for further analyses except where otherwise specified.
2.2. Structural Properties
lists the average values and standard deviations (in parentheses) for several structural properties of proteinase K calculated based on the joined trajectories at different combined temperatures. Under the combination of P180/S180, proteinase K has the largest average values for NNC and SSE while the smallest average values for RMSD and SASA, whereas under the P300/S300 combination, proteinase K has the smallest average values for NNC and SSE but the largest values for RMSD and SASA. This implies that the simulations at the combined temperature P180/S180 lead to the most stable structure and the most compact packing of proteinase K, whereas the opposite situations occur when simulations were performed at P300/S300. Interestingly, although the differences in structural properties between P300/S180 and P180/S300 are minor, proteinase K has a greater number of inter-atomic contacts, more compact packing, and smaller structural deviation from the starting structure when the solvent temperature is at 180 K than when at 300 K. In conjunction with larger standard deviations of all structural properties at P180/S300 than at P300/S180, these results indicate that the effect of the solvent temperature overwhelms that of the protein temperature, resulting in a higher rigidity/flexibility at the low/high solvent temperature.
When the solvent temperatures increase from 180 to 300 K while the protein temperatures remain unchanged, e.g., at 300 K (or 180 K), the NNC and number of residues in SSE decrease by 1.5% (0.9%) and 6.6% (3.3%), respectively, and the RMSD and SASA increase by 82.3% (45.2%) and 2.7% (1.8%), respectively; when the solvent temperatures remain at 300 K (or 180 K) while the protein temperatures increase from 180 to 300 K, the NNC and number of residues in SSE decrease by 1.2% (0.6%) and 5.4% (2.1%), respectively, and the RMSD and SASA increase by 37.8% (9.7%) and 1.9% (1.0%), respectively. These results reveal that, although the changes in protein temperatures have effects on structural properties of proteinase K, their effects are weaker and quantitatively smaller than those of the solvent temperature changes.
The Rg of proteinase K is 16.6 and 16.7 Å when the solvent is at 180 and 300 K, respectively, independent of the protein temperature. In addition, for all of the structural properties, their standard deviations follow the order P180/S180 < P300/S180 < P180/S300 < P300/S300, implying that proteinase K experienced larger structural variations and hence, a stronger internal mobility, during simulations at the high solvent temperatures than at the low solvent temperatures.
2.3. Structural Flexibility
root mean square fluctuation (RMSF) was computed from the joined trajectories to evaluate and compare the structural flexibility of proteinase K under the four different temperature combinations (Figure 1
A). As shown in Figure 1
A, although proteinase K exhibits the highest and lowest overall structural flexibility at the combined temperatures P300/S300 and P180/S180, respectively, the low (P300/S180) and high (P180/S300) solvent temperatures significantly reduce and elevate the overall protein flexibility, respectively, leading to very similar flexibility profiles between P300/S180 and P180/S180, and between P180/S300 and P300/S300. Therefore, the intensity of atomic fluctuations, and hence the protein internal flexibility, is strongly related to the solvent temperatures. Of particular note is that the solvent temperatures have a greater influence on the flexibility of loop regions than on that of the SSEs. For instance, the intensity of loop fluctuations is drastically increased and decreased when the temperature combination goes from P180/S180 to P180/S300 and from P300/S300 to P300/S180, respectively, whereas the corresponding changes in the intensity of SSE fluctuations are smaller, implying that the solvent-temperature-dependent changes in protein flexibility originate from the loop motions.
In order to further determine the effect of the solvent temperatures on flexibility of the entire protein structure, we computed the Cα
RMSF values as a function of distance to the protein surface (Figure 1
B). As shown in Figure 1
B, whether at the protein surface or in the protein core, the RMSF values are in an order P300/S300 > P180/S300 > P300/S180 > P180/S180 and, for these four combinations, a consistent trend is observed: RMSF values increase with increased distance from the protein core to the surface. However, the degree of this increase is different, with the most dramatic increase observed for P300/S300, large increase for P180/S300, and small increase for P300/S180 and P180/S180, ultimately leading to the amplitudes of atomic fluctuations, which are distributed over the entire protein structure, in the order P300/S300 > P180/S300 > P300/S180 > P180/S180.
For all four temperature combinations, the largest amplitudes of atomic fluctuations occur invariably on the protein surface. This is mainly due to: (i) the lack of structural packing environment for the surface-exposed atoms; and (ii) the direct interactions/contacts between the protein surface atoms and the water molecules. Nevertheless, the ordinal increase in RMSF values of the protein core from the low to the high solvent temperatures indicates that the mobility of the solvent can also determine the amplitude of atomic fluctuations in the protein core. The small increases in RMSF values from P180/S180 to P300/S180 reveal that elevating the temperature of the protein does make a contribution, to some extent, to its own mobility. Moreover, the degrees of these RMSF increases are relatively consistent over the entire protein structure (from the core to the surface), which may reflect a pure thermal relaxation of the protein arising from the elevated protein temperature. When compared to P300/S180, under the P180/S300 combination the protein shows not only significantly higher surface and near-surface RMSF values but also slightly higher core RMSF value. This result reveals that the high surface mobility caused by the high solvent temperature is not only transmitted to the protein core, but also improves the mobility of protein core.
2.4. ED and Collective Motions
In order to investigate the effect of the solvent temperatures on the collective motions (or global large-scale concerted motions) of proteinase K, ED analyses were performed on protein Cα atoms using the joined trajectories of the four simulation systems.
The total mean square fluctuation (TMSF) values obtained by diagonalization of the covariance matrices are 2.97, 1.33, 0.64, and 0.47 nm2 for the combinations of P300/S300, P180/S300, P300/S180, and P180/S180, respectively, indicating that proteinase K experiences significantly larger atomic fluctuations when the solvent temperatures are at 300 K than when at 180 K, in agreement with the above structural flexibility analysis. When the temperature combination goes from P300/S180 to P300/S300 and from P180/S180 to P180/S300, the TMSF values increase by 3.64- and 1.83-fold, respectively; the TMSF values increase by 1.23- and 0.36-fold when the temperature combination goes from P180/S300 to P300/S300 and from P180/S180 to P300/S180, respectively. These results reveal that elevating the temperature of the solvent makes a substantial contribution to increasing protein fluctuations. Although elevating protein temperature also contributes to protein fluctuations, such a contribution is quantitatively smaller, especially when the solvent temperature is at 180 K.
shows the eigenvalues as a function of eigenvector index. The four simulation systems share a common trend: the largest eigenvalue is for the first eigenvector and then the eigenvalues decrease with increased eigenvector index until they reach a minimum plateau. However, the degree of the steepness for these four eigenvalue curves is different and follows an order P300/S300 > P180/S300 > P300/S180 > P180/S180. Apparently, the first few eigenvectors of P300/S300 have the highest eigenvalues among those of the four combinations. The second highest, low, and the lowest eigenvalues of the first few eigenvectors are observed for the P180/S300, P300/S180 and P180/S180, respectively. These results indicate that the high solvent temperatures enhance greatly the anharmonic atomic motions and thus promote the global collective motions of the protein along the first few eigenvectors. In contrast, at the low solvent temperatures (i.e.
, P300/S180 and P180/S180), the first few eigenvectors have eigenvalues only slightly higher than those of eigenvectors with higher index, implying that the insignificant anharmonic atomic motions predominate along the first few eigenvectors, and as thus, the collective motions of the protein are greatly suppressed.
shows the projection extremes of the first eigenvector, which represent the most significant collective motions of proteinase K at each combined temperature. It should be noted that the linear interpolations between these two extremes are not the conformational transition pathway between the two extremes but emphasize merely the conformational differences between them. Under the P300/S300 combination (Figure 3
D), proteinase K exhibits the most significant collective motions that involve conformational displacements over the entire structure, with the large and small displacements observed in the surface-exposed loops and SSEs, respectively. The second significant collective motions of proteinase K are observed for the P180/S300 combination (Figure 3
C), which involve conformational displacements in most of the surface loops and some of the SSEs. Under the P300/S180 (Figure 3
B) and P180/S180 (Figure 3
A) combinations, most of the structural regions in the protein do not exhibit apparent displacements and only a very limited number of surface-exposed loops exhibit the large displacements. These observations are consistent with the eigenvalue order of the first eigenvector described above, both suggesting that the high solvent temperatures promote the global collective motions of the protein.
Combined ED analyses were performed to investigate the differences in ED properties. There are six possible pairwise combinations among the four joined trajectories; for purpose of clarity, we choose only the pair of P300/S180–P180/S300 to demonstrate the effect of the solvent temperatures on ED properties. Figure 4
shows the projections of the merged trajectories onto the combined eigenvectors and the properties of these projections. Only in the case of the first combined eigenvector can the projection be well separated (top panel, Figure 4
A), which is reflected clearly by the distinctly different projection distributions (top panel, Figure 4
B) and average values (Figure 4
C) between P180/S300 and P300/S180. This reveals that the changes in the protein and solvent temperatures lead to different equilibrium fluctuations and average structures of proteinase K along the first combined eigenvector. Starting from the combined eigenvector 2, the two equal halves of the projection show gradually increasing overlap between their distributions, as well as increasing similarity between the average values, e.g., when reaching the eigenvector 30, the projection distributions overlap into a near-Gaussian distribution and the average values are almost identical. These results indicate a gradually increasing similarity between the motion modes and average structures of the protein under these two temperature combinations, implying that the degree of harmonicity in protein motions increases with increased eigenvector index.
Of particular note is that for each of the first few combined eigenvectors, the right half of the projection (corresponding to P180/S300) exhibits a larger amplitude of the fluctuations than the left half of the projection (corresponding to P300/S180) (Figure 4
A). This reflects that at the high solvent temperature, proteinase K experiences more dramatic conformational changes along these combined eigenvectors than at the low solvent temperature. Figure 4
D shows the MSD values for the two halves of the projection as a function of eigenvector index, which can be used to quantitatively compare the difference in degree of the protein conformational shifts at different combined temperatures. Apparently, most of the eigenvectors (with the exception of eigenvectors 5 and 13) have significantly larger MSD values when the combined temperature is at P180/S300 than when at P300/S180, revealing that the high solvent temperature greatly facilitates the conformational shifts/changes of the protein in a subspace spanned by these eigenvectors.
2.5. Sampled Conformational Subspaces
In order to examine the influence of the solvent temperatures on the extent of conformational space explored by proteinase K, the four joined trajectories were projected onto a two-dimensional conformational subspace spanned by the first two eigenvectors (Figure 5
). As shown in Figure 5
, proteinase K explores similar shapes of the sampled regions under the temperature combinations P180/S180, P300/S180, and P300/S300, with the smallest and largest regions observed for P180/S180 and P300/S300, respectively. The second largest region is observed for the P180/S300; although this region exhibits a different shape with respect to the other three ones, it not only covers the entire regions of P180/S180 and P300/S180, but also has the largest overlap with that of P300/S300, revealing that the simulations at the high solvent temperatures enlarge the sampled conformation space of the protein.
As described above, elevating solvent temperature facilitates not only greatly the protein surface atomic fluctuations but also to a certain extent the fluctuations of atoms in the protein core (Figure 1
B). In order to further determine the influence of the solvent temperatures on the sampling degrees of the protein core, core periphery, and surface loops, these three structural parts were projected onto the subspace spanned by the first two eigenvectors (Supplementary Figure S2
). Visualization of the sampled conformational spaces reveals that: (i) for all four temperature combinations, the protein surface loops, core periphery, and protein core sample the largest, moderate, and smallest regions, respectively; (ii) the surface loops sample significantly more extended regions at the high solvent temperatures (P180/S300 and P300/S300) than at the low solvent temperatures (P300/S180 and P180/S180); and (iii) the core periphery and protein core explore relatively larger regions at the high solvent temperatures than at the low solvent temperatures. The above results reveal that the increased solvent temperatures facilitate not only the conformational sampling of the surface loops but also those of the structural parts that have no direct contacts/interactions with the solvent, in agreement with the comparative analysis of RMSF values as a function of distance to the protein surface.
2.6. Free Energy Landscapes
shows the FELs of the protein–solvent systems at the four combined temperatures as a function of two-dimensional conformational subspace (Figure 6
A–D) and of one-dimensional eigenvector projection (Figure 6
E,F). It should be pointed out that: (i) the present ns-timescale MD simulations sample only a limited portion of the protein accessible conformational space (e.g., the conformational space close to the native state but containing no the unfolded states); (ii) despite the reduction of dimensionality to the two- or one-dimensional space, protein conformational space is highly multi-dimensional; and (iii) the eigenvector projection biases towards second-order correlations of protein motions that result in conformational populations with mixed properties [55
]. Therefore, the eigenvector-projection-based FELs constructed from the metadynamics simulations are incomplete, representing only a major portion of the landscape but containing no detailed characteristics of the free energy surface (e.g., multiple minima or hierarchical arrangement of the free energy wells). Nevertheless, comparison between these FELs is still useful in probing characteristic differences resulting from alternation of the solvent/protein temperatures because the differences in statistical trends are still preserved.
As shown in Figure 6
A–D, all four FELs present a funnel-like shape. However, it is apparent that the FELs of the low solvent temperatures (P180/S180 and P300/S180) exhibit a perfectly smooth surface, whereas the surfaces of P180/S300 and P300/S300 feature a weak and strong rough character, respectively. This implies that the high solvent temperatures are able to produce local free energy minima either on the funnel wall or at the funnel bottom, thus increasing the protein conformational diversity or the number of conformational states/substates near the native state and hence the flexibility of the protein. The width of the FELs, from the bottom to the top outer surface, is found to be narrowest for P180/S180 and P300/S180, whereas for P180/S300 and P300/S300, the funnel width is obviously wide and the widest, respectively (Figure 6
E,F). This reveals that the high solvent temperatures greatly increase the conformational entropy of the protein, whereas the protein temperatures have only minor effect on the protein conformational entropy. Another interesting result is that the FEL corresponding to the P300/S300 has a pronounced higher minimum free energy level than the other three FELs. Although the other three FELs have similar minimum free energy values, the value of the P180/S300 is slightly higher than those of P180/S180 and P300/S180. This suggests that proteinase K has lower thermostability (or structural stability) at the high solvent temperatures than at the low solvent temperatures.
2.7. HB Interactions
HBs are the most important class of non-covalent interactions/forces that can contribute not only to the stability [56
] but also to the flexibility [58
] of the protein. In order to ascertain the origin of the characteristic differences in FELs under the four temperature combinations, we calculated the HB properties for both static and dynamic HBs within the protein and between the protein and solvent based on the joined MD trajectories.
lists the properties of the static and dynamic HBs formed inside the protein. The static HB number, which refers to that averaged over all frames in a trajectory, is 239.2 ± 5.3, 232.9 ± 6.1, 229.0 ± 6.8, and 215.6 ± 7.8 for P180/S180, P300/S180, P180/S300, and P300/S300, respectively. A clear trend is that at the low solvent temperatures the protein has higher numbers of internal HBs than at the high solvent temperatures, although lowering the protein temperature (when the solvent temperatures remain unchanged) can still increase the HB number, e.g., 6.3 from P300/S180 to P180/S180 and 13.4 from P300/S300 to P180/S300. The static HBs play a crucial role in stabilizing protein structure because the HB formation is an exothermic electrostatic interaction that contributes to lowering the system free energy via a negative enthalpy change. In this context, the static HB number can be considered as indicative of the structural stability/rigidity. Therefore, the above results show that the low solvent temperatures increase the protein structural stability/rigidity through increasing the static HB number. Interestingly, the number of the dynamic HBs, which refers to the total number of HBs that have ever appeared in a MD trajectory, exhibits an opposite trend with respect to that of the static HBs. As shown in Table 2
, significantly higher dynamic HB numbers occur at the high solvent temperatures (P180/S300 and P300/S300) than at the low solvent temperatures (P180/S180 and P300/S180), indicating that elevating the solvent temperature increases the number of the dynamic HBs inside the protein. The dynamic HB number is inversely correlated with the average HB persistency: the higher the number of the dynamic HBs, the lower is the average value of the bond persistency, and vice versa
); this implies that it is the increase of the low-persistence HBs that leads to an increase in the dynamic HB number. A low-persistence HB forms transiently and is easy to break during a simulation so that its donor and/or acceptor can form new HBs with the neighboring donor and/or acceptor atoms, thus leading to competitive HB interactions and, ultimately, the increase in the number of the dynamic HBs. Because the competitive interaction is able to cause conformational transition, the dynamic HB number can be considered as indicative of the structural flexibility or even as an approximate measure of the protein conformational entropy. To this end, the high solvent temperatures increase the protein flexibility/mobility through increasing the dynamic HB number.
Protein internal HBs involve either the main chains or the side chains of amino acid residues and therefore can be divided into three classes: the main chain–main chain, main chain–side chain, and side chain–side chain HBs. Although elevating the solvent temperature increases the numbers of all these three classes of HBs, the bonds involving the side chains (including main chain–side chain and side chain–side chain HBs) make a greater contribution than do the main chain–main chain HBs to the increase in the dynamic HB number (Table 2
). This is not surprising, because the side chains have a larger degree of freedom than the main chains and as thus are more sensitive to changes of the solvent temperature.
lists the properties of the static and dynamic HBs formed between the protein and the solvent, which directly reflect the changes in protein–solvent interactions among the four temperature combinations. The high solvent temperatures (P180/S300 and P300/S300) result in the lower/higher numbers of the static/dynamic protein–solvent HBs than do the low solvent temperatures (P180/S180 and P300/S180), resembling the situations of the protein internal HBs. The increased static protein–solvent HB interactions at the low solvent temperatures suggest the role of the low-temperature water molecules in rigidifying/freezing the protein structure. Because of the significantly enhanced mobility of the water molecules when the solvent temperatures go from 180 to 300 K, the numbers of the dynamic HBs at the high solvent temperatures are at least two orders of magnitude higher than those at the low solvent temperatures (Table 3
). This is also reflected by the two orders of magnitude reduced HB persistency from the low- to the high-solvent temperatures, suggesting that it is the high-frequency competitive interactions between different water molecules and the protein atoms that give rise to the increased dynamic HB number and thus, the increased mobility of the protein. For all four temperature combinations, the HBs formed between the amino acid side chains and the water molecules (side chain-solvent HBs) account for ~90% of the corresponding dynamic HBs and exhibit a slightly lower persistency than the main chain–solvent HBs, revealing that the side chains are the primary participant of the protein–solvent HBs or, in other words, water molecules affect the protein dynamics mainly through their interactions with the amino acid side chains.
“Glass transition” (or dynamical transition) of proteins, like that of glass-forming liquids [60
], is a well-established phenomenon for all proteins manifesting as a dramatic change in their dynamical properties at approximately 200 K [17
, at temperatures above 200 K, the protein dynamic behavior is dominated by the large-scale collective motions involving the entire protein structure, while at temperatures below 200 K, the simple harmonic atomic vibrations predominate, thus resulting in the cease of the protein activity. Because the dried proteins do not show any dynamical transition [36
] and the glass transitions in the protein and its hydration water occur at roughly the same temperature [17
], the protein glass transition is believed to be induced by the glass transition of the solvent. In order to provide detailed information about the effect of the solvent temperatures on the changes in dynamics of the protein and, further, to probe the mechanism by which the solvent mobility dictates the hierarchical dynamics of the protein, we take advantage of the molecular simulation technique in which the protein and solvent temperatures are coupled to either 300 K (above the protein glass transition temperature) or 180 K (below the transition temperature). In contrast to the previous MD study performed on a relatively small protein carboxy-myoglobin (153 amino acids) with a short production simulation time (only 100 ps) [42
], our simulations were performed on a relatively large globular protein (279 amino acids), the serine protease proteinase K, with application of the long, multi-replica simulation strategy to ensure a more effective and extensive conformational sampling of the protein. Indeed, the evaluation of cosine contents for the first few eigenvectors shows that the single joined trajectories from the multiple simulations achieve a higher degree of sampling convergence than the individual replicas, and therefore our comparative analyses based on these joined trajectories can reflect differences in intrinsic dynamic properties of proteinase K caused by different protein/solvent temperature combinations.
As a necessary complement to this work, we have also evaluated the influence of protein temperatures on the solvent mobility by calculating the mean square displacement (MSD) of the solvent in the P300/S180 and P180/S300 trajectories. The results (Supplementary Figure S3
) show that under the P180/S300 combination, the water molecules demonstrate a free diffusion behavior, and thus, are not “frozen” by the low protein temperature, whereas under the P300/S180 combination the MSDs of the hydration water and bulk water are characterized by an initial rise followed by a plateau formation, reflecting that the water molecules are in a structurally arrested, glass-like state. Of interest is that in the P300/S180 system the hydration water molecules are found to have slightly higher MSD values than the bulk water during both the rising and plateau stages. This indicates that although the hydration water molecules remain still localized or arrested around the protein surface (i.e.
, they are not “melted” by the high temperature of the protein), the high temperature of the protein can improve, to a certain extent, their diffusivity due to the physical contacts between the protein and the hydrated water molecules. Taken together, it can be concluded that the protein temperatures have only a very limited effect on the diffusivity (mobility) of the solvent.
The comparative analyses of several static structural properties and of conformational flexibility reveal that the solvent temperatures have an overwhelming effect, when compared to the protein temperatures, on these properties and atomic fluctuations of proteinase K. For instance, the changes of the solvent temperatures can lead to significant changes in average structural properties investigated (Table 1
RMSF profiles (Figure 1
A), and TMSF values, whereas the changes of the protein temperatures only result in minor changes in the structural properties and flexibility of proteinase K. Our results are consistent with previous MD simulation on myoglobin [42
], both revealing the dominant role of the solvent temperature in determining the internal mobility and flexibility/rigidity of the protein.
The ED analyses reveal that at the high solvent temperatures, the motion mode described by the first eigenvector is the global collective motions involving displacements of the entire protein structure. On the contrary, at the low solvent temperatures, the large displacements were observed only in very limited structural regions, indicating that the low solvent temperatures greatly suppress the collective motions of proteinase K. Further combined ED analyses reveal that the dramatic conformational shifts occur only at the high solvent temperatures, while at the low solvent temperatures, the simple harmonic fluctuations predominate even along the first few eigenvectors. In our previous MD studies [46
], the collective motions of proteinase K have been shown to not only lead to the breathing/twisting of the substrate-binding channel and pockets, but also play a role in modulating dynamic behaviors of the structural regions not directly involved in the catalysis but indirectly influencing the enzymatic function. Thus, the global collective motions of proteinase K are related to the substrate binding, orientation, catalysis, product release, and the regulation of these processes. The disappearance of the global collective motions below glass transition temperature explains why the proteins will cease to be active upon the glass transition.
The comparison between RMSF values as a function of distance to the protein surface, in conjunction with the comparison between sampling degrees of the conformational space explored by the protein core, core periphery, and surface loops at the four combined temperatures, reveals that: (i) the surface mobility can transmit to the core periphery and further to the core; and (ii) the mobility and sampling degrees of the non-surface-exposed regions are proportional to those of the protein surface. In light of the above findings, we can propose the following rough scenario: the high mobility of the water molecules, which arises from the high solvent temperatures, can lead to more frequent collisions/interactions with protein surface when compared to water molecules at the low temperatures, and therefore results in larger atomic fluctuations and structural displacements of the protein surface, especially of the solvent-exposed loop regions that lack the structural constraints. Such a high mobility can further transmit via specific mechanic mechanism to the protein core periphery and finally to the protein core, thus causing the global collective motions of the protein.
The number of dynamic HBs can be considered as a representative of competitive interactions within the protein and between the protein and the solvent because of the universal distribution of HBs in the protein–solvent system. Comparison between protein–solvent HB numbers reveals that the high solvent temperatures result in more intense competitive interactions between the water molecules and the protein surface atoms than do the low solvent temperatures. In addition, elevating solvent temperature also facilitates greatly the competitive interactions between atoms within the protein. In a previous theoretical study, Agarwal has investigated the transfer of energy from the solvent to the protein by adding kinetic energy to solvent molecules in the hydration shell, and the results show that within a short time, the solvent kinetic energy is transferred to the protein residues even more than 8 Å away from the protein surface, resulting in increased mobility of the protein [35
]. Taken together, it is reasonable to conclude that it is the competitive interactions that transmit solvent kinetic energy/mobility to the protein surface and further to the protein interior, ultimately leading to an overall increased conformational flexibility and the collective motions of the protein.
Comparison between our constructed FELs reveals that elevating the solvent temperature increases the conformational substates and conformational entropy while reducing the thermostabilility of the protein. These results, in conjunction with the data presented in this work and available from the literature [2
], allow us to propose a more refined FEL model to explain the effect of the solvent temperatures on the protein dynamics (Figure 7
). It should be point out that: (i) the FEL model plotted in Figure 7
is not original, but is mainly based on the model documented in references [2
] and [66
]; (ii) the information about the conformational transitions occurring on timescales of μs to ms comes from the experimental data rather than the MD data presented in this work; and (iii) the model presented here is to highlight the role of the solvent mobility in determining the hierarchical dynamics of proteins as described follows. First, the elevated solvent temperature increases the kinetic energy of water molecules and hence their mobility. This will result in the high-frequency collisions/competitive interactions of the water molecules with the protein surface (especially with the surface-exposed amino acid side chains) and thus cause the rotational motions of the side chains. As shown in Figure 7
, side chain rotations correspond to the tier-2 dynamics, which occur on the timescale of ps and involve crossing among the tier-2 free energy wells. Second, the accumulation of side chain rotations will break the local non-covalent interactions involved in the loop regions that lack the structural constraints, which in turn exacerbates the competitive interactions between residues and water molecules and as thus leads to the loop motions on the protein surface (i.e.
, the ns-timescale tier-1 dynamics that involve crossing among the tier-1 free energy wells). Finally, by perturbing the water network (hydration shell) around the protein surface [6
] and increasing interactions of the loops with the other protein structural parts [35
], the large loop fluctuations will facilitate the competitive interactions between atoms within the protein, thus transmitting the solvent kinetic energy over the entire protein structure and leading to the global collective motions of the protein (i.e.
, the tier-0 dynamics that occur on the timescales of μs to ms and involve crossing among the tier-0 free energy wells). Therefore, the high mobility of the solvent, which arises from its high temperature, plays a crucial role in facilitating the cascade amplification of microscopic motions of atoms and atomic groups into the protein global collective motions. However, the low solvent mobility, which results from its low temperature, will inhibit such cascade amplification through suppressing the competitive interactions between atoms within the protein and between the protein and solvent, ultimately leading to insignificant anharmonic motions of the protein as observed in this study.