7584
J. Phys. Chem. B 2007, 111, 7584-7590
How Protein Surfaces Induce Anomalous Dynamics of Hydration Water Francesco Pizzitutti† and Massimo Marchi Commissariat a l’Energie Atomique, DSV-DBJC-SBFM, Centre d’EÄ tudes, Saclay, 91191 Gif-sur-YVette Cedex, France
Fabio Sterpone*,† and Peter J. Rossky Department of Chemistry and Biochemistry, UniVersity of Texas at Austin, 1 UniVersity Station, CM A 5300, Austin, Texas 78712 ReceiVed: March 2, 2007; In Final Form: April 26, 2007
Water around biomolecules slows down with respect to pure water, and both rotation and translation exhibit anomalous time dependence in the hydration shell. The origin of such behavior remains elusive. We use molecular dynamics simulations of water dynamics around several designed protein models to establish the connection between the appearance of the anomalous dynamics and water-protein interactions. For the first time we quantify the separate effect of protein topological and energetic disorder on the hydration water dynamics. When a static protein structure is simulated, we show that both types of disorder contribute to slow down water diffusion, and that allowing for protein motion, increasing the spatial dimentionality of the interface, reduces the anomalous character of hydration water. The rotation of water is, instead, altered by the energetic disorder only; indeed, when electrostatic interactions between the protein and water are switched off, water reorients even faster than in the bulk. The dynamics of water is also related to the collective structuresa` Voir the hydrogen bond (H-bond) networksformed by the solvent enclosing the protein surface. We show that, as expected for a full hydrated protein, when the protein surface offers pinning sites (charged or polar sites), the superficial water-water H-bond network percolates throughout the whole surface, hindering the water diffusion, whereas it does not when the protein surface lacks electrostatic interactions with water and the water diffusion is enhanced.
I. Introduction Almost all biomolecules function in aqueous environments and their activity is critically modulated by the structure and dynamics of water. A biomolecule and its solvatation water molecules, or biological water,1 can be treated as a coupled system whose dynamics occurs on several time scales:2,3 from picoseconds to nanoseconds and longer. To elucidate the extent of this coupling, experiments, theoretical models, and computer simulations have recently focused on water dynamics around biomolecule surfaces, including small globular proteins,2,4-8 multidomain protein,9 DNA strands,10 direct11,12 and reverse micelles,13-15 and bilayers.16 It has been appreciated for some time that the water dynamics around the surface of biomolecules is retarded with respect to bulk water.17 For example, close to the surface of a globular protein or of a micelle with small hydrophilic groups, the water dynamics slows down from 4 to 7 times.2,5,6,18-20 Around non-ionic micelles with long hydrophilic tails the retardation could become larger, up to 1-2 orders of magnitude,5,21 although on this issue there is not complete agreement.22 In the interior of reverse micelles, water structure and dynamics increasingly differs from the bulk as the size of the water pool decreases, and depends on the distance from the surfactant hydrophilic part;13,15,23 in the proximity of the interface, water rotation has been estimated to be more than 5 times slower than in the interior of the pool.15 In the case of proteins, theoretical * Corresponding author. Current address: F. S. Caspur, via dei Tizii 6B, 00185 Rome, Italy. † These authors contributed equally to this work.
efforts also have been made to assess a relationship between surface composition (hydrophilic/hydrophobic amino acids),7,24 secondary structure,25 3D conformation,26 and water dynamics. Of more recent interest, the dynamical slowdown is accompained by an anomalous behavior, reminescent of a glassy state. Molecular dynamics (MD) simulations have shown that the water mean square displacement (MSD) in the proximity of a biomolecule is sublinear in time, and that its dipolar relaxation, as well as its survival probability on the protein surfaces, is described through stretched exponentials.6,18,27 Moreover, neutron scattering measurements28 and MD simulations on protein crystals29 have shown the presence of the Boson peak in the single particle dynamic spectra and of a stretched decay of the intermediate scattering function. All of these findings have led the scientific community to depict the solvation water as a glassy system characterized by energetic and topological disorder. It is worth noting that the origin of the sublinear diffusion is as intriguing as it is elusive. It is a common assumption to trace the subdiffusive dynamics back to either the fractal geometry of the biomolecular surface (geometrical disorder) or the dynamics on a multiple minima energetic landscape generated by the biomolecule surface (temporal disorder).27,30 In the former case, the water dynamics on the surface can be viewed as described by a random walk in a geometrically disordered environment, characterized by a fractal dimension df and a spectra dimensionality ds which are related to the MSD through the following relationship: 〈∆r2〉 ≈ tdf/ds. In the latter case, the water motion occurs in a disordered energetic landscape that generates a long tail in the distribution of the residence times
10.1021/jp0717185 CCC: $37.00 © 2007 American Chemical Society Published on Web 06/12/2007
Anomalous Dynamics of Hydration Water on the surface: ψ(t) ≈ t-µ. This tail produces a MSD proportional to tµ-1. At the present, however, there is no clear evidence as to which process dominates the hydration water dynamics. Recently, Bizzarri, Rocchi, and Cannistraro18,27,31 have carefully tackled the problem, studying the water diffusion around a hydrated plastocyanin (PC) at ambient temperature and around a low hydrated PC in a wide range of temperatures, above and below the dynamical temperature transition (Td = 150 K). They found that formally both the models can describe the results obtained by the simulated trajectories. Others have suggested that water diffusion at the protein surface is determined by geometrical disorder.6 Nanosecond simulations of a fully hydrated globular protein (lysozyme) were performed by using two different models of water (SPC/E and TIP3P) ,and for both models, the water dynamics at the protein surface is slower with respect to the bulk solution and the MSD can be fitted by a power law tR with R ) 0.6. The two models should create a different protein-solvent energetic landscape, and possibly (although this is not demonstrated) a different spreading of the water-protein residence times. For this reason, the geometrical disorder scenario seemed favored. Simulations of other complex environments such as reverse23 and direct micelles21 have, however, suggested that both disorders contribute in slowing down water dynamics. Recently, Mukherjee and Bagchi32 have also proposed a simple atomistic model in which the sublinear diffusion follows from the exchange dynamics between free and bound water environments around a biomolecule. The model was created on the basis of NOE measurements1 and simulation data (see references in ref 17). Despite the fact that a large number of experimental and simulation results are now available, the origin of the so-called anomalous behavior of biological water is still puzzling. In particular, the role of the biomolecule’s multiple time scale dynamics must be considered carefully. Indeed, it is well-known that the nanosecond exchange dynamics of water is coupled to the slow rearrangment of the protein matrix.33-35 Exchange dynamics is shown also in the intermediate temporal scale of hundred of picoseconds (7% of hydration water) and is coupled to faster protein relaxations.33 Such a multiple time scale coupling must affect the solvation dynamics in protein folding, aggregation, and function.3 Moreover, as discussed first by Nilsson and Halle4 and recently by Qiu et al.8 and Golosov and Karplus,36 the dynamics of the protein per se must be carefully considered when spectroscopic data collected from a chemical probe attached to a protein are used to investigate solvation dynamics. Therefore, concerning the water diffusion, the theoretical explanations invoking either the geometrical or the energetic disorder of a biomolecule surface have to take into account the time evolution of such disorder. In this paper, we focus on the basic ingredients that underlie the slowdown of the hydration water dynamics around a globular protein, namely, (i) the role of protein motion and (ii) the role of protein topological and energetic disorder. For this aim, we have performed ad hoc computer simulations of different protein-water systems specifically targeted to tackle these issues. First, we have analyzed how the protein motion affects the anomalous dynamics of the solvatation water molecules in a fully hydrated lysozyme. We performed a simulation of a frozen protein matrix surrounded by a water solution thermalized at T ) 300 K and have compared the results with those obtained for a system where the protein and water are in thermal equilibrium at T ) 300 K. The approach we adopted follows
J. Phys. Chem. B, Vol. 111, No. 26, 2007 7585 that of other studies where ad hoc computer simulations based on dual bath dynamics, which thermalizes the protein and the solvent at different temperatures, have been used to study the coupling of protein-water motion in the context of the protein/ water glass transition.29,37,38 In a second step, we turned our attention to the role of protein topological and energetic disorder. Hence, we switched off the electrostatic interactions between the protein and water to analyze the effect of the static geometrical disorder without interference of the energetic disorder generated by the polar and charged groups of the protein surface. Such strategy was successfully used by others in the study of hydrophobic collapse in multidomain protein folding.39 The paper is organized as follows: In section II, we present the systems, the models employed in the simulations, and the tools used for analyzing the MD trajectories. In section III, the results are presented and discussed. The discussion is focused on the water residency dynamics and on the rotational and translational water dynamics at the protein surface. In the last section the main results are summarized and final conclusions are drawn. II. Methods We present results from simulations of several systems constituted by a lysozyme protein surrounded by 3954 water molecules. The first of these systems (LYS-F) is built by freezing the protein matrix and thermalizing the water solution at T ) 300 K and P ) 1 atm. The initial configuration of lysozyme was obtained from a previous simulation described in ref 6. In this case, protein and water interact with their full potential. In contrast, in the second system considered here, LYS-VDW, we switched off the water-protein electrostatic interaction but kept the protein frozen as in LYS-F. In what follows, all the results for LYS-F and LYS-VDW are compared to results of a previous simulation of a lysozyme/water solution at T ) 300 K and P ) 1 atm (LYS-T).6 We used the Force Field developed and parametrized by MacKerell et al.40 and the SPE/E model for water.41 The simulations of LYS-F and LYS-VDW have been carried out in the NVE ensemble. A multiple time step scheme was used to integrate the equations of motions, see ref 33 for details. The SPME algorithm with the same parameters as in ref 33 was used to handle the electrostatic interactions. The 9 ns trajectories were sampled every 240 fs for statistical analysis. The analysis of the hydration water dynamics is performed by computing the survivial probability function, Nw(t), the restricted mean square displacement (MSD), and dipolar relaxation. Nw(t) gives the probability that a water molecule is in contact with the protein for a time t. To obtain Nw(t), we compute the conditional probability Pj(tn,t) for each jth water molecule of the system and then average it over the length of the run:
Nw(t) )
1
Nt
∑ ∑j Pj (tn,t)
Nt n)1
(1)
Here, Pj(tn,t) takes the value of 1 if the jth water is adjacent to any of the residues of the protein between times tn and tn + t, and it has the value of zero otherwise. Nt is the number of the simulation time-frames of length t. The adjacent waters, or hydration waters, are selected by using a cutoff radius; thus a water molecule is considered in the hydration shell of the protein if the internuclear distance to the water oxygen is less than Rc:
Rc ) f(rw + rp)
(2)
7586 J. Phys. Chem. B, Vol. 111, No. 26, 2007
Pizzitutti et al. TABLE 1: Residence Times Computed by Fitting the Survival Probability Nw(t) with One Stretched Exponential (〈τs〉 and ns) and Two Exponentials (τi and ni with i ) 2, 3)a survival probability system
〈τs〉 (ps)
LYS-F LYS-VDW LYS-T
15.6 23.6 14.6
ns
γ
τ2 (ps)
n2
490.6 0.473 250.5 30.8 342.8 0.382 4.5 181 516.2 0.474 301.1 26
τ3 (ps) 1930 4487 2110
n3
np
23.5 14 9.6 8 12.1 3
a The last column reports the number of permanently attached water molecules np. See eq 5.
Figure 1. Survival probability for water in the hydration shell of lysozyme: black squares, the frozen system LYS-F; maroon diamonds, LYS-VDW; and red circles, LYS-T. The solid lines show the fit functions (see Table 1). In the inset, the interval for t < 100 ps is magnified.
where rw and rp are the van der Waals radii of the water oxygen and of the protein atoms, respectively. The coefficient f is set to 1.1 as in ref 6. We remark that Nw(t) takes important meanings for t ) 0 and for the length of the simulation, t ) tsim. Indeed, Nw(0) is the average number of hydration waters, whereas Nw(tsim) corresponds to the number of waters “permanently” attached to the solute, on the 9 ns time scale. The dynamics of water at the protein surface can be characterized by studying the translational and rotational diffusion of water molecules residing on the protein surface. Hence, the mean-square displacement (MSD) 〈|r(t) - r(0)|2〉, as a function of time, can be restricted to the water molecules in contact with the protein surface. In the same way, one can compute the first- and second-rank dipole-dipole correlations:
P1(t) ) 〈cos θ(t)〉 P2(t) )
〈23cos θ(t) - 21〉 2
(3) (4)
where θ(t) is the angle between the water dipole of the same molecule at time τ and τ + t, and 〈...〉 stands for statistical average. The comparison of P1(t) and P2(t) computed for the hydration waters with a simple kinetic relaxation model makes it possible to estimate their characteristic times, τ1R and τ2R. III. Results and Discussion A. Hydration State and Exchange Dynamics Time Scales. To extract the time scales of the water/protein dynamics, we fitted the survival probability curves Nw(t), shown in Figure 1, with the following function:
Nw(t) ) nse-(t/τs)γ + n2e-(t/τ2) + n 3e-(t/τ3) + np
(5)
where np indicates the number of water molecules that were attached to the protein for the whole simulation time, np ) Nw(tsim). The Nw(t) curves were computed for the water molecules in the hydration shell of the protein, and the obtained fit parameters are reported in Table 1 for the systems LYS-F, LYSVDW, and LYS-T described in the last section. For the stretched exponential part, the averaged characteristic time is computed as:
〈τs〉 )
τ
∫0∞ e-(t/τ )γ dt ) γs Γ(γ1)
with Γ as the gamma function.
s
(6)
The frozen system LYS-F exhibits the same time scales already discerned for the water dynamics around the protein in the LYS-T system. The fast decay, characterized by a time scale in the range of picoseconds, is well reproduced by a stretched exponential and describes the majority of the hydration water, ns = 490 of the total Nw(t)0) ) 531. The exponent γ of the stretched exponential decay measures the manifold dynamics of hydration water and thus provides a quantitative indication of the complex landscape underlying the escape dynamics of water from the protein surface. It is important to note that the computed γ values in LYS-F and LYS-T have the same value, γ = 0.47; thus, at short time, the protein motion does not affect the escape process. This is shown graphically in the inset of Figure 1 wherein the survival probility is plotted for t′ < 100 ps. As seen from the parameters in Table 1, the two slower time scales, labeled τ2 and τ3, have characteristic times in the 100-1000 ps range and are experienced by only =11% of the hydration water. However, the freezing of the protein motion induces an increase in the amount of “τ3 waters”, compared to the LYS-T system. Furthermore, the number of water molecules remaining in contact with the protein (np) increases from 3 to 14. The data in Table 1 highlight the rather different dynamic behavior of the LYS-VDW system, charaterized by the disappearence of the time scale around 300 ps: Thus, at faster times we find two contributions to the decay of τ2 ) 4.5 ps and 〈τs〉 ) 23.6 ps which have an exponential and a stretched exponential nature, respectively. The exponent γ ) 0.382 of LYS-VDW is smaller than that in the LYS-F and LYS-T cases, indicating that the topological disorder has a strong effect on the escape dynamics at shorter time. It should be noted that the weighted average of the two shorter times yields to a characteristic time of 16.9 ps, which is smaller than 29.2 and 28.7 ps obtained with the same averaging for LYS-F and LYS-T, respectively. This finding is consistent with the basic idea that electrostatic interactions strengthen the water-protein coupling. The fast decay is indeed quite enhanced compared to the other two systems, as is shown in the inset graph of Figure 1, whereas in the nanosecond scale the survival probability seems to be more similar to that of the LYS-F system. Finally we note that for the nonpolar frozen system, the number of permanent waters is np ) 8. B. The Permantly Attached Water. The 14 permanently attached water molecules found in LYS-F are found to be located in the superficial clefts or the internal pores of the protein, and are hydrogen bonded to either the backbone oxygens or the polar lateral groups. Among this class of molecules, we recognized the conserved group of water molecules buried in the hydrophilic pore of the protein, found in previous simulations33 and in X-ray experiments.42 The other locations for the fixed water molecules are defined by the initial conformation of the protein. Indeed, the water molecules remain buried in superficial clefts because no local fluctuations of the
Anomalous Dynamics of Hydration Water
J. Phys. Chem. B, Vol. 111, No. 26, 2007 7587
Figure 2. The first- and second-rank dipolar relaxations for the simulated systems: LYS-T in black, LYS-F in red, and LYS-VDW in blue. For LYS-F we also computed the relaxation of the water molecules residing in the hydration layer for less than 1 ns, LYS-FNoLong in green.
protein allow the exchange with the bulk. Therefore, a partially open pocket thatsthanks to the flexibility of the protein matrixs would fluctuate between an open and a closed conformation, in the frozen system is blocked in its initial configuration. Switching off the protein-water electrostatic interactions, all the H-bonds existing in LYS-F and LYS-T are suppressed. On the other hand, the water molecules buried in closed spaces cannot escape from the protein surface due to geometrical confinement. This explains the high number of permanent water molecules detected in the LYS-F and LYS-VDW system and the long time decay of the associated survival probability functions. C. Residence Time, Dipolar Relaxation, and Diffusion. The survial probability Nw(t) estimates the several time scales of the water exchange dynamics, but does not provide insight into the dynamics occurring at the protein surface. To explore this, we computed the dipolar relaxation and the translational diffusion of water residing on the protein surface and then estimated the residence time of the hydration water. Here, we limit ourselves to two distinct definitions of the residence time. The first consists of defining the residence time as the time required for a water molecule to lose its orientational memory. Such a definition can be exploited in Nuclear Magnetic Resonance Dispersion (NMRD) experiments:2,43 by assuming all water rotation as isotropic the residence time is given by 2 τR ) 3FτR-Bulk
(7)
where the factor F is the second-rank rotational retardation, F 2 2 ) (τ2R)/τR-Bulk , with τ2R and τR-Bulk denoting the second-rank relaxation of surface and bulk water, respectively; finally 3 is the proportionality factor relating the characteristic times of the first- and second-rank dipole-dipole correlation function for the isotropic model. The second definition assumes that the residence time is the time that a water molecule takes to move away from its residence site (τw). τw can be extracted from the MSD of the hydration water molecules, as the time needed by a diffusing particle to cover one water diameter, i.e., ∼3 Å. The dipolar relaxation (P1(t) and P2(t)) and the translational diffusion are shown in Figures 2 and 3, respectively. The P1(t) and P2(t) are fitted for t < 25 ps with a stretched exponential, Pn(t) ) exp(-(t/τnR)γ), where n is the rank of the dipole-dipole correlation function. By using eq 6 we computed the average time 〈τ1R〉 and 〈τ2R〉. The parameters extracted from the fits of the first- and second-rank dipole-dipole correlation functions are reported in Table 2. Table 3 compares the residence times, τR and τw. For LYS-F we find 〈τ2R〉 ) 15.5 ps, which is about 1.6 time larger than that computed for LYS-T, 〈τ2R〉 ) 9.8 ps. Thus, we
Figure 3. The mean square displacement (MSD) of hydration water computed for the systems LYS-T (black circles), LYS-F (red squares), and LYS-VDW (blue diamonds) is compared to the MSD for bulk water (black dashed line). The MSD for the water molecules residing less than 1 ns in the hydration shell of LYS-F is also reported in green triangles (LYS-FNoLong).
TABLE 2: Parameters for the Fit of the First-Rank, P1(t), and Second-Rank, P2(t), Dipole-Dipole Correlation Functiona P1(t)
P2(t)
system
τ1R (ps)
γ
〈τ1R〉
τ2R (ps)
γ
〈τ2R〉
LYS-F LYS-VDW LYS-T
28.00 1.67 19.50
0.42 0.69 0.44
81.9 (4.9) 2.1 48.5
3.50 0.45 2.55
0.36 0.60 0.38
15.5 (1.9) 0.7 9.8
a The fit is performed only for the initial decay (t < 25 ps) with a stretched exponential. The average time 〈τiR〉 is computed by using eq 6. The value in parentheses is the computed average relaxation time for bulk SPC/E water.
TABLE 3: Residence Time Computed for a Rotationally Isotropic Model, τR, and from the MSD, τwa system
τR (ps)
τw (ps)
R
LYS-F LYS-VDW LYS-T
46.5 2.1 29.4
62.8 29.0 26.8
0.46 0.6 0.6
a In the last column we report the parameters of the fit of the computed MSD with a power law function, i.e., 〈|r(t) - r(0)|2〉 ∝ tR.
find a retardation of F = 8 with respect to bulk SPC/E water where 〈τ2R〉 ) 1.9 ps, and the rotational residence time τR of 46.5 ps. On the contrary, hydration water in LYS-VDW shows a rotational relaxation faster than that occurring in bulk water, specifically 〈τ2R〉 ) 0.7 ps, and an estimated residence time τR ) 2.1 ps. This faster relaxation is a direct consequence of the lack of H-bonds between the water molecules and the protein surface. For the frozen system LYS-F, we find τw ) 62.8 ps, which is about 12 larger than τw = 5.2 ps for bulk. When the electrostatic interactions are switched off, the water diffusion is still slower than that in the bulk, but the retardation is considerably smaller than that detected in LYS-F, namely, τw ) 29 ps is about 5.5 larger than τw for bulk. We point out that in the LYS-T, the residence time defined through the dipolar relaxation, τR, and that extracted from translational diffusion, τw, are very similar,6 whereas they significantly differ (25%) in the LYS-F system and more so in the case of LYS-VDW, see Table 3. Thus the energetic landscapes that govern the rotational and the translational diffusions in our frozen systems are distinct and are poorly described by a simple isotropic model.
7588 J. Phys. Chem. B, Vol. 111, No. 26, 2007 The difference between the translational and rotational mechanism, also detected in a low hydrated protein simulated at low temperatures38 and in complex micellar interfaces,21 must be related to the confinement effect created by the frozen protein state and to the energetic competition among the solventsolvent and solvent-protein H-bonds. The LYS-VDW represents an extreme situation. In this case, the water molecules can rotate in their residence sites without interference from the water-protein H-bonds. At the same time, they have a restrained translational displacement due to the static disordered protein surface. To further characterize the nature of the water rotation on the protein surface, for LYS-F and LYS-T we computed the ratio between the first-rank (〈τ1R〉) and the second-rank (〈τ2R〉) dipolar relaxation. We find values considerably larger than 3, β ) 〈τ1R〉/〈τ2R〉 )5-5.3. A β close to 3, as those obtained for a bulk water solution and for water on the surface of LYS-VDW systems, is the signature of a quasi-isotropic landscape for the dipolar relaxation. Generally, a β different from 3 measures the anisotropy of the dipolar rotation. In particular a value larger than 3 reflects a rotation constrained onto a conic surface.44 Anisotropic rotations of water have been previously found in biosolutions. Smith and co-worker44 computed a β close to 1 in a high concentration PEO/water solution for the vector perpendicular to the water plane. Also, a ratio close to 1 was computed for the dipolar relaxation of water at the interface of a non-ionic C12E6 micelle.21 In those cases, the rotation is spatially constrained onto a plane because of the structural H-bonds formed by the water molecules and the ether oxygen groups. These findings show that the electrostatic interactions between water and the protein surface alter locally the water rotations with respect to those in the bulk. Furthermore, the freezing of the protein motion increases the activation energy of the boundto-free transition of the H-bonds formed between the water molecules and the protein surface, causing a slower rotational relaxation. Last but not least, we highlight the fact that the γ exponent of the stretched exponential fit of the dipolar correlation functions takes very similar values for LYS-F and LYS-T: 0.40 and 0.44 for the first rank P1(t), and 0.36 and 0.38 for the second rank P2(t), respectively. On the contrary for LYS-VDW, γ is higher, 0.6-0.7, a value very close to that obtained for water rotation on the large hydrophobic surface of a LDAO micelle.21 To summarize, our results on the dipole-dipole dynamics show that the protein-water electrostatic contribution has an important effect on both the energetics, changes in 〈τ1R〉 and 〈τ2R〉, and the manifold dynamics, changes in γ values, of water rotation at the protein surface. D. The Origins of the Anomalous Diffusion. In Figure 3, we report the MSD of the hydration water molecules in proximity of the protein surface. The curve associated with LYS-F can be fitted by a power low, tR with R = 0.46, a value smaller than that previously computed for LYS-T (R = 0.6). Such a value demonstrates that the protein matrix fluctuations act to reduce the anomalous character of hydration water translational diffusion. As observed for water motion in a glassy environment,45 when the protein is thermally activated the translational diffusion is guided by two mechanisms: (i) water jumps in sites previously occupied by other water molecules and (ii) water jumps in sites previously occupied by protein groups. For a frozen system, the second channel is blocked. The effective dimensionality of the space where the hydration
Pizzitutti et al. water translation occurs is then reduced with respect to a system wherein all protein atoms move. To discriminate between the geometrical and the energetic disorder, it is, however, crucial to also consider the system LYSVDW. The profile of the MSD computed in this case is fitted by a power law with R ) 0.6 and the curve overlaps with that of LYS-T, see Figure 3 and Table 3. The diffusion is then faster compared to LYS-F but still very far from the brownian regime. Thus the intrisic fractal dimensionality of the protein surface (an estimate of the dimensionality gives df ) 2.1746) contributes significantly to the anomalous behavior of water. This finding quantifies for the first time the influences of the geometric and energetic disorder on the water diffusion around a static protein surface. Both of these two properties of the protein matrix acting on water molecules contribute to the behavior of the local hydration dynamics. Nevertheless, we also expect that for LYSVDW, if the atoms of the protein were allowed to move the water would reduce its slowdown, eventually leading to a bulkwater-like regime. This is supported by our finding on a micelle, LDAO, with strongly hydrophobic surfaces, where we showed that the hydration water dynamics is almost unperturbed compared to its bulk behavior.21 Finally, we turn our attention to the specific contribution of water exchanging with the bulk in the nanosecond time scale. Specifically, we are interested in the question whether the sublinearity depends on the set of molecules tightly bound to the protein and on their interconversion from or to a free state in the hydration shell, as postulated in some models.32 With this in mind, we have computed the MSD and the dipolar relaxation for hydration water of LYS-F, excluding all those molecules in contact with the protein for more than 1 ns. This allows us to precisely and completely exclude the conversions between free and bound environments that occur in the protein hydration shell on the nanosecond time scale, and also exclude the direct exchange with bulk occurring in the same time scale. We find that these waters affect only weakly the translational diffusion, which remains sublinear in nature, as is evident in Figure 3. Similarly, the dipolar relaxation for times smaller than 500 ps is almost unaffected by the exclusion of long residing molecules. This shows clearly that the anomalous behavior of water is fully determined by water molecules that are in contact with the protein for less than a nanosecond. E. Hydrogen Bond Networks. The different behavior of water around the two frozen systems, LYS-F and LYS-VDW, can be rationalized if one considers the water-water hydrogen bond network connectivity at the protein surface. Using a geometrical criterion53 we first counted the number of H-bonds formed by water residing in the hydration shell. The analysys of the H-bond connectivity has been carried out over a part of each trajectory, 2 ns. The number of H-bonds per molecule in the hydration shell is nhb = 2 and nhb = 2.7 for LYS-F and LYS-VDW, respectively. Thus when the electrostatic interactions are switched off the water molecules tend to compensate for the lack of connectivity with the protein surface increasing the water-water H-bonding. However, since the number of hydration waters is almost identical in all the systems, 〈Nw〉 = 530, the enhanced water-water connectivity observed in the LYS-VDW system suggests that water molecules at the surface tend to increase their H-bonding with waters that are outside the first shell.47 To describe the H-bond connectivity around the protein surface, it is essential to consider the behavior of clusters of water-water H-bonded molecules.54 The probability distribution for the maximum cluster size, P(Smax), computed over 2 ns of
Anomalous Dynamics of Hydration Water
Figure 4. Distribution of the maximal cluster size, P(Smax), computed for the systems LYS-F, LYS-VDW, and LYS-T.
each trajectory is reported in Figure 4. The plot shows that there is a drastic difference between the two systems LYS-F and LYSVDW. The former shows the characteristic profile of systems above a percolation threshold. In this case, the average maximum cluster size is 〈Smax〉 = 860 molecules. The number of water molecules found in the hydration shell being about 530, the largest cluster on the surface includes also molecules outside the first hydration shell H-bonded to the hydration water molecules. This finding shows, as is expected for a full hydrated protein,48 that when the water interacts electrostatically with the protein surface, the solvent is able to create a network of H-bonds spanning the whole protein. On the contrary, we find that the peak in the P(Smax) for LYS-VDW is well below the percolation threshold. Hence, despite the higher connectivity of the single molecule, hydration waters, lacking the support of pinning residues on the protein surface, cannot encage the protein in a single hydration bonded network. Indeed, the solvent forms a set of small clusters with a far smaller average maximum size, on the order of Smax = 60 molecules. In such hydrogen-bonded clusters, each particle displacement implies the reorganization of essentially all bonds connecting the particles. Moreover, in the case of LYS-F the water-water H-bond network is also anchored to the protein surface and the single particle diffusion slows down because of the size, 〈Smax〉 = 860, and the structure of the percolating cluster. It should also be noted that vibrational spectroscopy experiments showed that the water percolating cluster around a protein is formed by very strong water-water H-bonds48 and that simulations have appreciated that water-water H-bonds nearby a protein surface or a small polypeptide as well as the pinning H-bonds between water and protein are more stable than the H-bonds in bulk water.49,50 This helps in understanding the strong slowdown observed in LYS-F. In the case of LYS-VDW, the retardation of water diffusion is less severe because of the smaller size of the clusters, 〈Smax〉 = 60, and the lack of anchoring with the protein surface, hence the water motion has fewer reorganizational requirements, on average. In particular at the interface between disconnected clusters the particle displacement is weakly constrained by the H-bond connectivity. The thermally activated protein, LYS-T, is also found to be surrounded by a percolating network as in the LYS-F system, but the translational diffusion of hydration water is faster than in the frozen protein case. A similar percolating structure has also been observed in the hydration shell of staphilococcal nuclease (SNase) at T ) 300 K.51 These results suggest that (i) the H-bond network around the protein surface is flexible and rearranges with protein motion without breaking the overall
J. Phys. Chem. B, Vol. 111, No. 26, 2007 7589 connectivity and (ii) the rearrangement is associated with an increase in the water translational diffusion. Moreover, our findings are complementary to the recent work of Oleinikova et al.52 in which the single molecule connectivity, or nhb, at the protein surface has been demostrated to be a percolating parameter. Also in agreement with results on LYS-F and LYS-T (for both nhb = 2) discussed in this study, for low hydrated Lysozyme the critical value for the existence of a spanning H-bond network has been estimated in the range nhb = 2.0 to 2.3, depending on the system and temperature. Here, we show also that at high local connectivity, nhb = 2.7, the hydrobobic nature of the protein (LYS-VDW) induces a nonpercolating structure of the hydration layer and this is correlated with the lack of H-bond connectivity between the protein and the proximal hydration water molecules. IV. Conclusion In this paper we studied the dynamics of water molecules found adjacent to the surface of protein models. To explain the origin of the water anomalous dynamical behavior at the solutesolvent interface, two different protein-water systems have been considered. The first consists of a frozen protein surrounded by water molecules thermalized at 300 K. In the second system, the protein matrix is still kept frozen and, additionally, the electrostatic interactions between the protein and water are eliminated. Results from these systems are compared to a solution where protein-water interactions are included in full, dynamics is unconstrained, and the entire system is thermalized at 300 K. The obtained results allow us to consider separately the contributions from geometrical and energetic disorder, and from protein motion. We have demonstrated for the first time that around a static protein both types of disorder act on the water translational diffusion, contributing both to the average retardation and also to the glassy-like anomalous translational diffusion of water at the protein surface. We have also shown that when the protein surface offers electrostatic H-bond pinning sites to water molecules, the hydration water forms a percolating H-bond cluster that surrounds the whole protein and hinders the water dynamics whereas when the electrostatic interactions between the protein and the solvent are eliminated, water only forms relatively small H-bonded clusters and superficial water diffusion is enhanced. We observed that the rotational dynamics on the protein surface is basically shaped by electrostatic interactions alone and that the H-bonds formed by water with the protein surface break the quasi-isotropic nature of the dipolar rotation that is found in the bulk. Also, for the fully thermalized protein, LYS-T, we find a ratio between the characteristic times of the first and the second dipole-dipole correlation function, β ) (τ1R)/(τ2R), of about 5 that is at variance with the isotropic assumption, β ) 3, used in NMRD estimates of the translational residence time. In addition, we find that the protein motion reduces the retardation of the water dynamics, because the dimension of the water translational space is increased and at the same time the orientational decorrelation is accelerated. In spite of this accelerated dynamics, hydration water diffusion remains anomalous for a thermalized protein such as LYS-T. Finally, we have demonstrated unequivocally that the dynamical anomalies are not related to the set of water molecules tightly bound to the protein (τ > 1 ns) and to their behavior in the nanosecond time scale, whether due to internal shell conversion or direct exchange with the bulk.
7590 J. Phys. Chem. B, Vol. 111, No. 26, 2007 Acknowledgment. P.J.R. gratefully acknowledges support by the U.S. National Science Foundation (Collaborative Research in Chemistry, CHE0404695) and by the R. A. Welch Foundation (F0019). References and Notes (1) Nandi, N.; Bagchi, B. J. Phys. Chem A 1997, 101, 10954. (2) Halle, B. In Hydration Process in Biology; Bellissent-Funel, M. C., Ed.; IOS Press, Ohmsha: Amsterdam, The Netherlands, 1999. (3) Fenimore, P. W.; Fraunfelder, H.; McMahon, B. H.; Young, R. D. Proc. Natl. Acad. Sci. U.S.A. 2004, 101, 14408-14431. (4) Nilsson, L.; Halle, B. Proc. Natl. Acad. Sci. U.S.A. 2005, 102, 13867-13872. (5) Pal, S.; Peon, J.; Zewail, A. H. Proc. Natl. Acad. Sci. U.S.A. 2002, 99 (4), 1763. (6) Marchi, M.; Sterpone, F.; Ceccarelli, M. J. Am. Chem. Soc. 2002, 124, 6787. (7) Russo, D.; Muraka, R. K.; Copley, J. R. D.; Head-Gordon, T. J. Phys. Chem. B 2005, 109, 12966-12975. (8) Qiu, W.; Kao, Y.-T.; Zhang, L.; Yang, Y.; Wang, L.; Stites, W. E.; Zhong, D.; Zewail, A. H. Proc. Natl. Acad. Sci. U.S.A. 2006, 103, 13979-13984. (9) Hua, L.; Huang, X.; Zhou, R.; Berne, B. J. J. Phys. Chem. B 2006, 110, 3704-3711. (10) Pal, S.; Zhao, L.; Zewail, A. H. Proc. Natl. Acad. Sci. U.S.A. 2003, 100, 8113-8118. (11) Kumbhar, M.; Goel, T.; Mukherjee, T.; Pal, H. J. Phys. Chem. B 2004, 108 (4), 19246-19254. (12) Vass, S.; Grimm, H.; Ba´nyai, I.; Meier, G.; Gila´nyi, T. J. Phys. Chem. B 2005, 109, 11870. (13) Abel, S.; Sterpone, F.; Badyopadhyay, S.; Marchi, M. J. Phys. Chem. B 2004, 108, 19458-19466. (14) Faeder, J.; Ladanyi, B. M. J. Phys. Chem. B 2005, 109, 67326740. (15) Dokter, A.; Wourtersen, S.; Bakker, H. J. Proc. Natl. Acad. Sci. U.S.A. 2006, 103, 15355-15358. (16) Lopez, C. F.; Nielsen, S. O.; Klein, M. L.; Moore, P. B. J. Phys. Chem. B 2004, 108, 6603-6610. (17) Bagchi, B. Chem. ReV. 2005, 105, 3197. (18) Bizzarri, A. R.; Cannistraro, S. J. Phys. Chem. B 2002, 106, 66176633. (19) Balasubramanian, S.; Bagchi, B. J. Phys. Chem. B 2001, 105, 12529. (20) Bruce, C. D.; Senapati, S.; Berkowitz, M. L.; Perera, L.; Forbes, M. D. E. J. Phys. Chem. B 2002, 106, 10902-10907. (21) Sterpone, F.; Marchetti, G.; Pierleoni, C.; Marchi, M. J. Phys. Chem. B 2006, 110, 11504-11510. (22) Carlstro¨m, G.; Halle, B. J. Chem. Soc., Faraday Trans. 1 1989, 85, 1049. (23) Faeder, J.; Ladanyi, B. M. J. Phys. Chem. B 2000, 104, 10331046. (24) Makarov, V. A.; Andrews, B. K.; Smith, P. E.; Pettitt, B. M. Biophys. J. 2000, 79, 2966. (25) Bandyopadhyay, S.; Chakraborty, S.; Balasubramanian, S.; Bagchi, B. J. Am. Chem. Soc. 2005, 127, 4071. (26) Bandyopadhyay, S.; Chakraborty, S.; Bagchi, B. J. Phys. Chem. B 2006, 110, 20629.
Pizzitutti et al. (27) Rocchi, C.; Bizzarri, A. R.; Cannistraro, S. Phys. ReV. E 1998, 57, 3315. (28) Settles, M.; Doster, W. Faraday Discuss. 1996, 103, 269. (29) Tarek, M.; Tobias, D. Phys. ReV. Lett. 2002, 89, 275501. (30) Klafter, J.; Zumofen, G.; Blumen, A. Chem. Phys. 1993, 177, 821. (31) Bizzarri, A. R.; Paciaroni, A.; Cannistraro, S. Phys. ReV. E 2000, 62, 3991. (32) Mukherjee, A.; Bagchi, B. Phys. Chem. Commun. 2003, 6, 2831. (33) Sterpone, F.; Ceccarelli, M.; Marchi, M. J. Mol. Biol. 2001, 311, 409-419. (34) Garcia, A. E.; Hummer, G. Proteins 2000, 38, 261-272. (35) Damjanovic, A.; Garcia-Moreno, B.; Lattman, E. E.; Garcia, A. E. Proteins 2005, 60, 433-449. (36) Golosov, A. A.; Karplus, M. J. Phys. Chem. B 2007, 111, 14821490. (37) Vitkup, D.; Ringe, D.; Petsko, G. A.; Karplus, M. Nature 2000, 7, 34-38. (38) Tournier, A. L.; Xu, J.; Smith, J. C. Biophys. J. 2003, 85, 18711975. (39) Zhou, R.; Huang, X.; Margulis, C. J.; Berne, B. J. Science 2004, 305, 1605-1608. (40) MacKerell, A., D. Jr.; Bashford, D.; Bellot, R. L.; Dunbrack, R. L., Jr.; Evanseck, J. D.; Field, M. J.; Fischer, S.; Gao, J.; Guo, H.; Ha, S.; Joseph-McCarthy, D.; Kuchnir, L.; Kuczera, K.; Lau, F. T. K.; Mattos, C.; Michnick, S.; Ngo, T.; Nguyen, D. T.; Prodhom, B.; Reiher, W. E., III; Roux, B.; Schlenkrich, M.; Smith, J. C.; Stote, R.; Straub, J.; Watanabe, M.; Wiorkiewicz-Kuczera, J.; Yin, D.; Karplus, M. J. Phys. Chem. B 1998, 102, 3586-3616. (41) Berendsen, H.; Grigera, J.; Straatsma, T. J. Phys. Chem. 1987, 91, 6269-6271. (42) Kundrot, C. E.; Richards, F. M. J. Mol. Biol. 1987, 193, 157-170. (43) Halle, B. In Hydration Process in Biology; Bellissent-Funel, M. C., Ed.; IOS Press, Ohmsha: Amsterdam, The Netherlands, 1999. (44) Borodin, O.; Bedrov, D.; Smith, G. J. Phys. Chem. B 2002, 106, 5194-5199. (45) Molinero, V.; Goddard, W. A., III Phys. ReV. Lett. 2005, 95, 457011. (46) Pfeifer, P.; Welz, U.; Wippermann, H. Chem. Phys. Lett. 1985, 113, 535-540. (47) Lee, C. Y.; McCammon, J. A.; Rossky, P. J. J. Chem. Phys. 1984, 80 (4448), 2299. (48) Khoshtariya, D. E.; Hansen, E.; Leecharoen, R.; Walker, G. J. Mol. Liq. 2003, 105, 13. (49) Xu, H.; Berne, B. J. J. Phys. Chem. 2001, 105, 11929. (50) Bandyopadhyay, S.; Chakraborty, S.; Bagchi, B. J. Am. Chem. Soc. 2005, 127, 16660. (51) Brovchenko, I.; Krukau, A.; Smolin, N.; Oleinikova, A.; Geiger, A.; Winter, R. J. Chem. Phys. 2005, 123, 2249051. (52) Oleinikova, A.; Brovchenko, I.; Smolin, N.; Krukau, A.; Geiger, A.; Winter, R. Phys. ReV. Lett. 2005, 95, 247802. (53) Two water molecules are considere H-bonded if the distance between the two oxygens is less than 3.5 Å and the angle formed between the OH bond and the OO axis is less than 30°. We then count the number of H-bonds formed between two water malecules of which at least one is adjacent to the protein surface. (54) A cluster of water molecules on the protein surface is defined as a set of water molecules all connected by H-bonds as defined previously.
