Reconciling Mediating and Slaving Roles of Water in Protein Conformational Dynamics

Proteins accomplish their physiological functions with remarkably organized dynamic transitions among a hierarchical network of conformational substates. Despite the essential contribution of water molecules in shaping functionally important protein dynamics, their exact role is still controversial. Water molecules were reported either as mediators that facilitate or as masters that slave protein dynamics. Since dynamic behaviour of a given protein is ultimately determined by the underlying energy landscape, we systematically analysed protein self energies and protein-water interaction energies obtained from extensive molecular dynamics simulation trajectories of barstar. We found that protein-water interaction energy plays the dominant role when compared with protein self energy, and these two energy terms on average have negative correlation that increases with increasingly longer time scales ranging from 10 femtoseconds to 100 nanoseconds. Water molecules effectively roughen potential energy surface of proteins in the majority part of observed conformational space and smooth in the remaining part. These findings support a scenario wherein water on average slave protein conformational dynamics but facilitate a fraction of transitions among different conformational substates, and reconcile the controversy on the facilitating and slaving roles of water molecules in protein conformational dynamics.


Introduction
Protein dynamics is critical for their functions [1][2][3][4] and evolvability [5], and is to a great extent determined by the roughness of their potential energy surface (PES). Solvents play an indispensable role in shaping dynamic behaviour of proteins through molecular interactions that contribute to protein PES. Energy transferred from the first hydration shell to surface residues of cyclophilin A was computationally demonstrated to influence catalysis through network fluctuations [2]. It is well known that, due to the hierarchical nature of PES [6], conformational dynamics of native proteins is hierarchical and occurs on many different time scales corresponding to different types of molecular motions, including bond stretch and bending motions on femtoseconds, rotations of small groups (e.g. methyl) on picoseconds, side chain and backbone dihedral rotations on subnanoseconds to microseconds, and major domain motions up to multiple milliseconds. It is likely that water molecules play different roles in the above mentioned various type of dynamic processes. Coupling between the function and internal motions of proteins and their water environment has been intensively studied [7][8][9][10][11][12][13]. Two lines of evidences that have been presented by many experimental [14][15][16][17][18][19][20] and computational [21][22][23][24][25][26][27] reports supporting either mediating or slaving roles of water molecules are briefly summarized below.
It was demonstrated by both experimental [14] and computational studies [22] that below glass transition temperature, protein dynamics is slaved(or caged) by surrounding frozen water molecules. Protein dynamics on 10 ps to 100 ps time scales was found to be closely correlated with dynamics of surrounding water hydrogen bond network and thus slaved by water molecules [23]. Water molecules' relaxation was observed to correlate well with conformational transitions of myoglobin among statistical substates [16], which occur on time scales ranging from sub-nanoseconds to microseconds at the physiological temperature, suggesting the slaving role of water molecules on corresponding time scales.
At physiological or room temperature, certain hydration level is essential for functions of many proteins [13,28]. Based on the analysis of crystallographic water molecules, it was proposed that water molecule ''lubricate'' folding of proteins through three bond centre hydrogen bonds [21]. Theoretical protein structure prediction studies [24] revealed that addition of water mediated potential in protein design facilitated search of the free energy minimum (i.e. native state), water molecules were also found to mediate native state dynamics of eglin C [26,27]. Raman optical activity studies [15] support the role of water molecules as ''lubricant of life''. From an energy landscape perspective, observations along this line were explained with the belief that water molecules facilitates protein dynamics through effectively smoothing their PES. However, direct evidence supporting this concept is lacking.
In this study, we generated collectively 5 microsecond molecular dynamics (MD) trajectories for a small globular protein barstar [29], which is synthesized by the bacterium Bacillus amyloilyquefaciens as an inhibitor of the ribonuclease protein barnase. By systematically analysing the time series (or evolution in conformational space) of relevant energy terms (protein self energy (E p ), protein-water interaction energy (E p{w ) and their sum (E tot )) obtained from MD trajectories, we found that while the negative correlation between E p and E p{w in most parts of conformational space provides possibility for PES smoothing, the dominance of s Ep{w over s Ep (s stands for standard deviation) resulting in a rougher PES on average, especially for picoseconds and longer time scales. These two aspects contribute to the end effects of water molecules on the PES of proteins, that is roughening the majority and smoothing the remaining part of the potential energy landscape. Thus, the conflicting roles of water in the protein conformational dynamics are reconciled with an energy landscape perspective. It is noted here that due to the constraint of computational resource, our analysis were limited to submicrosecond time scale, correspond to transitions covering a few hierarchies of statistical substates. The impact of water molecules on major domain motions that occur on milli-seconds and longer time scales and folding dynamics is beyond the scope of this study.

Results
The PES that underlies the dynamics of a protein molecule can be decomposed into two components, protein self energy (E p ) and protein-solvents interaction energy (E p{s ). Since deciphering roles of water molecules in protein dynamics is the goal of this study, we focus our attention on E p and protein-water interaction energy (E p{w ). For a single protein molecule travelling in the conformational space as in a typical MD simulation or a single-molecule experiment, the PES can be represented as a time series of potential energies with its roughness represented by corresponding standard deviations.
where s stands for standard deviation, Var stands for variance and Cov stands for covariance. A brief explanation of using standard deviation to represent PES roughness is given as follows. Unlike folding/unfolding and large scale conformational change between major conformations that occur on milliseconds and longer time scales, where one (or a few ) free energy barriers dominate, conformational transitions among a large number of hierarchical statistical substates on time scales up to microseconds involve many barriers on each specific time scale (e.g. nanoseconds) that have similar heights and are distributed over many degrees of freedoms. Therefore, we think standard deviations (s) of relevant potential energy terms are a reasonable representation of PES roughness for a given time scale dt. Expressing variances in eq.1 with standard deviations and the Pearson correlation coefficient r between E p and E p{w , we have: where E(dt) p stands for the average of the consecutive n potential energy values that have dt intervals in a time series E p . From eq.2, it is apparent that if on a given time scale dt, water molecules indeed smooth PES of a protein (i.e. s(dt) E tot vs(dt) Ep ), it is essential that the sum of the last two terms s(dt) 2 Ep{w z2s(dt) Ep s(dt) Ep{w r being non-positive, as standard deviations are always a non-negative number, this amounts to one of necessary conditions rv0. Therefore, E p{w need to negatively correlate with E p (e.g. when a protein molecule change into a higher energy configuration, water molecules compensate energetically by exerting a lower E p{w ). MD simulation provides great convenience in dissecting different potential energy terms on any time scales that is accessible by available computational power. To this end, we performed collectively 5 microsecond MD simulations of barstar and analysed resulting trajectories to obtain the Pearson correlation coefficients r between E p and E p{w . As shown in Fig. 1a, for time scales varying from 10 fs to 100 ns (see methods for specific procedures used to calculate energy correlation on a give time scale), E p{w on average negatively correlate with E p , with the correlation being the minimal on the shortest time scale analysed in our study (10 fs), becoming stronger on longer time scales up to 100 ps and levelling off beyond that point. However, the distribution of calculated correlation coefficient (r) in Fig. 1b exhibit both positive and negative correlations between the two energy terms (E p and E p{w ), only that the negative correlation has a larger probability of occurrence.
The immediate question one would ask is that does E p{w , that negatively correlate with E p , indeed smooths the PES of barstar. As mentioned above, negative correlation between E p and E p{w is only one of the necessary conditions. The sufficient condition is s 2 Ep{w z2s Ep s Ep{w rv0 (i.e. s Etot vs Ep ). To verify this condition, we calculated the averages for the standard deviations s Ep , s Ep{w and s Etot , and plotted them as a function of time scales (Fig. 2a). It is apparent that for all the time scales studied, addition of proteinwater interaction energies increased the roughness of the protein PES, and the roughening effects increases with increasing time scales. One interesting observation is that on short time scales (10 fs and 100 fs), due to small negative correlation between E p{w and E p , s Etot is greater than both s Ep{w and s Ep . On longer time scales (10 ps and longer), with increased negative correlation, s Etot becomes smaller than s Ep{w but remains greater than s Ep . Distributions of relevant energy standard deviations for two time scales (10 fs and 100 ps, representing short and long time scales) are shown in Fig. 2c and Fig. 2d respectively. When dt = 10 fs, s Etot exhibits the largest spread while s Ep{w has the largest spread for dt = 100 ps. For all time scales, s Ep has the smallest spread. Distributions of these three energy terms for all time scales analysed are available in Fig. S1. As ensemble averaged observables from molecular simulations has established correspondence with ensemble experimental measurements, our average PES roughness data demonstrated net roughening effects on all time scales studied, therefore unequivocally support slaving theory. Distributions of correlation coefficient r between E p and E tot (Fig. 1b) demonstrate the complex relationship between these two different PES compo-nents. Distributions of standard deviations for the three energy terms s Ep , s Ep{w and s Etot (Fig. 2c and 2d) indicate that in some region of the conformational space it is possible for s Etot to be smaller than s Ep . This observation suggests the possibility that water molecules do smooth corresponding part of protein PES. To validate this speculation, we compared the standard deviations calculated from each set of potential energy data (representing a specific region of PES) for E p and E tot , and plotted the probability that D~s Etot {s Ep being greater or equal to (P z , indicating roughening of PES) and smaller than 0 (P   demonstrate that in most parts of observed conformational space, protein-water interactions effectively roughen protein PES and play a smoothing role in the remaining part. The relative importance of the smoothing and roughening roles of water molecules shows weak dependence on time scales. The smoothing probability are larger on shorter time scales (10 fs to 1 ps) than on longer time scales (10 ps to 100 ns).
As mentioned above (see also eq. 2 and 3), the net effect of water molecules on the roughness of protein PES (s Etot {s Ep ) dependent on both correlation coefficient r between E p and E p{w and relative magnitude of s Ep and s Ep{w . To reveal the relationship between r and the net effects of water molecules on barstar PES, r vs. (s Etot {s Ep ) plots were generated for all the time scales that we studied and data for dt~10fs,10ps and 10 ns are displayed in Fig. 3a, b and c (data for other time scales are shown in Fig. S6). Each point in these plots represents a local region on a given time scale dt in the configurational space of barstar. In the four quadrants (noted I, II, III and IV in Fig. 3), quadrant II is always empty as it is a mathematically impossible region (from Eq. 2, if rw0, s Etot {s Ep w0), points in quadrant III corresponds to the smoothing role of water molecules, while points in quadrant I and IV correspond to roughening effects. Is it interesting to see that when dt change from 10 fs to 10 ps, the relative weight (noted as percentage in quadrants I, III and IV) of quadrant IV increased at the expense of quadrant I, and when dt change from 10 ps to 10 ns, the relative weight of quadrant IV increased at the expense of both quadrant I and III, but mainly quadrant III. This observation demonstrates that on short time scale (dt~10 fs), when the amplitude of s Ep{w is comparable with that of s Ep r becomes the major factor of smoothing/roughening effects.

Discussion
On very short time scales (dt~10{100 fs), dielectric relaxation experiments revealed that protein dynamics is more or less independent of water behaviour, while slaving is mainly observed for longer time scales. This is in qualitative agreement with our data (Fig. 2a) that the difference between s Etot and s Ep is the smallest on these time scales and monotonically increase with increasingly longer time scales. Additionally, the probability of occurrence for smoothing by water is larger on these short time scales (Fig. 2b and Fig. 3). However, the difference between s Etot and s Ep is very significant even for dt~10fs and should not be negligible in experimental dynamic measurements. This puzzle may be explained by the following fundamental physical causes that are not embodied in Eq. 2. All the bonding and bending degrees of freedom (DOFs) simultaneously contributing to the PES on femtoseconds time scales, when the limited increase of PES roughness (20 to 30 kcal=mol) are distributed among so many DOFs (1447 bonds and 2622 angles for barstar), the net effect (0:01{0:02 kcal=mol per DOF) is negligible considering large force constants of these motions (,50 kcal=mol=rad 2 for bending and 200-500 kcal=mol= Å 2 for bonding). However, although the total number of rotatable dihedral angles are not so small for a protein (3845 for barstar), on time scales longer than subnanosecond, transitions among various statistical substates mainly involve rotation of side chains (x 1 ) and backbone dihedral angles (w and y) in flexible regions of protein. The total number of these dihedral angles are about three times the number residues (,267 for barstar) and for any given region on PES, most of them are not rotatable on large scales (e.g. trans to gauche). When relatively large increase of PES roughness (,120 kcal=mol) are distributed over a small number of DOFs, the net effects (,1 kcal=mol per DOF) are significant considering the small force constants (0:1{5 kcal=mol) of most dihedral rotations.
Our data support both slaving and mediating roles of water molecules when different regions of PES were investigated separately. However, due to the fact that the slaving role has a larger probability of occurrence, and the fact that ensemble based experimental characterizations (e.g. photolysis analysis and dielectric relaxation measurements) are sensitive to the ensemble average of observables, thus only majority events (slaving) were seen. Conformational analysis performed on MD trajectories and PDB structures, should in principle be able to overcome such issues. However, the conformational states correspond to mediating roles of water (e.g. bridging residues of the same charge) are much easier to identify due to their relatively higher stability and structural simplicity, while a much larger number of conformational states correspond to slaving role of water molecules have lower stability and higher structural complexity. Therefore, in structural analysis of MD trajectories, the later tend to be neglected as random and featureless (which they are despite their significance in contribution to overall PES) events.
The role of water molecules is apparently even more important for the folding process of proteins as hydrophobic interactions are the most important driving force for protein folding. The critical role of individual water molecules is observed in both folding of protein molecules [30] and in triggering folding of model polymers with hydrophobic and hydrophilic monomers [31]. Both mediating [21,24] and slaving [17] roles of water in protein folding dynamics are reported. However, due to the constraint of computational resource, we were not able to perform similar analysis for the protein folding processes.
Sampling and accuracy of force fields are the two fundamental limitations that compromise the prediction power of molecular simulation methods. In this study, we have used collectively 5 microseconds of MD simulation trajectories. Based on the success of many simulation studies we believe that our trajectories should explore a significant and representative part of the concerned protein phase space for time scales varying from multiple femtoseconds to multiple nanoseconds. The time scales investigated here cover a few hierarchies of statistical substates transitions ranging from rapid bending motion to rotation of many side chains and backbone dihedral angles. As dynamic behaviour of proteins on longer time scales (micro-to milliseconds) have been shown to be coupled with shorter time scale dynamics [32], negative correlation and dominance of protein-water interactions on these time scales likely to impact long time dynamics in some way. We hope to address this issue in the immediate future. Additionally, our analysis is mainly based upon the fluctuation of relevant energetic terms. Therefore, the absolute value of these energetic terms, which absorb a significant part of inaccuracy of the adopted molecular mechanical force fields, is not a major concern.
Two different approaches have been utilized to analyse the roles of water molecules in protein conformational dynamics. One is to monitor the behaviour of fully solvated proteins via experimental [14][15][16][17][18][19][20] or computational methodologies [21][22][23][24][25][26][27], and this is what we adopted in the current study. In this scenario, the PES (of solvated proteins) consists of two components, one is the intramolecular contribution (E p ) and the other is the intermolecular contribution (E p{w ). Our analyses indicate that E p{w contribute more to the PES roughness than E p . In contrast, the other approach has compared dynamical behaviour of proteins with various extents of hydration [33][34][35][36], where the observed difference is indisputably resulted from different extents of hydration/solvation. In studies adopting the later approach, it was found that above 250 K (below that temperature, water molecules effectively freeze up most interesting protein motions), hydration significantly enhances protein dynamics. From a PES point of view, that unequivocally leads to the conclusion that water molecules smooth protein PES. However, different conclusions from these two distinct approaches do not necessarily constitute a direct contradict. In the former approach, the relative contribution of two components of the solvated protein PES is considered, one structural ensemble (the solvated native ensemble, this is a very approximate term as different solvation conditions correspond to different ensembles) is the focus of investigation. In the later case, totally different protein PESs (that of dry proteins, fully hydrated proteins or something in between) are compared, the two extreme structural ensembles (dry and fully hydrated ensembles) are different with the extent of differences (shared and distinct conformations) being unknown. Future investigations that compare dry and solvated proteins using experimental and/or computational techniques will provide more insights.
In conclusion, by decomposing the PES E tot of barstar into E p{w and E p , and analysing their correlations and roughness on time scales varying from femto-seconds to sub-microseconds, we found energetic evidence for both the slaving and mediating roles of water molecules. These analysis revealed that on the above mentioned time scales, in most part of the configurational space, water molecules slave protein dynamics through effectively roughening local PES and in the remaining part, water molecules may facilitate conformational dynamics by smoothing local PES.
Here we carefully studied the impact of protein-water interactions on the PES of barstar from an energy landscape perspective. It is possible that other proteins with different folds may have qualitatively different behaviour. Based on the experimental reports of similar slaving behavior of many proteins [16], our conclusion is likely to be qualitatively applicable for many other globular proteins as well. It is important to note that our study focuses on the PES of fully solvated proteins, the change of dynamical behavior from dry proteins to solvated ones is not covered. As exhaustive studies of all protein folds is not achievable due to prohibitive computational cost, we hope that this study may stimulate the community's interest to utilize available trajectories of different proteins to quantitatively answer this question, and accurately gauge the role of water molecules on protein conformational dynamics in a general sense.

Molecular Dynamics Simulations
All MD simulations were performed with NAMD software package [37],version 2.7 using CHARMM27 force fields. barstar (pdb code:1bta) was solvated with TIP3P water model. 100 mM N z a and Cl { were added to neutralize net charges of our simulation systems. Bond-lengths involving hydrogen atoms were constrained using the SHAKE algorithm, and the integration time step is set to 2 fs. Periodic boundary conditions were used, a switch distance of 10Å and a cutoff distance of 12 Å were used for nonbonded interactions. Particle Mesh Ewald (PME) were used to calculate the long-range electronic interactions. All systems were minimized and then heated to 310 K with heavy atoms restrained, water molecules were equilibrated with 200-ps runs in NVT ensemble. After that, restraints for protein heavy atoms were released, and the whole system was equilibrated in the NPT ensemble for another 4 ns. A frame with the volume value that is closest to the average volume obtained from NPT equilibration run was selected for the next production runs which were performed in the NVT ensemble at 310 K. 10 500-ns trajectories were generated. Coordinates were saved every ps for analysis. To generate potential energy statistics on short time scales (10 and 100 fs). 10 100-ps trajectories originating from snapshots taken every 100 ns from arbitrarily selected long trajectories were generated and coordinates were saved every 10 fs.

Energy Correlation Analysis
To measure the correlation between E p and E p{w of barstar, mean value and distributions of Pearson correlation coefficients r were calculated on time scales ranging from 10 fs to 100 ns. For a given time scale dt, consecutive n data points (t s , t s zdt, …, t s z(n{1)dt) were picked from each potential energy time series(E p and E p{w ) of available trajectories to calculate one Pearson correlation coefficient, 50000 coefficients were obtained with all the t s uniformly distributed in our trajectories, and were used to generate the mean and distributions. The presented data were obtained with n~5, a larger n(e.g. 7,8) would mix neighbouring time scales and are therefore not used. The data for n~3,4 were also calculated and presented in Fig. S2. As expected, different n generate similar results.