PHYSICAL
PRL 114, 187802 (2015)
REVIEW
week ending 8 MAY 2015
LETTERS
Curvature Dependence of Hydrophobic Hydration Dynamics R. Gregor WeiB,1’2 Matthias Heyden,2 and Joachim Dzubiella1,2' ^ Institut fiir Physik, Humboldt-Universitdt zu Berlin, Newtonstrasse 15, D-12489 Berlin, Germany "Soft Matter and Functional Materials, Helmholtz-Zentrum Berlin, Hahn-Meitner Platz 1, D-14109 Berlin, Germany Max-Planck-Insitut fiir Kohlenforschung, Kaiser-Wilhelm-Platz 1, D-45470 Miilheim an der Ruhr, Germany (Received 25 September 2014; revised manuscript received 25 February 2015; published 5 May 2015) We investigate the solute curvature dependence of water dynamics in the vicinity of hydrophobic spherical solutes using molecular dynamics simulations. For both the lateral and perpendicular diffusivity, as well as for H-bond kinetics of water in the first hydration shell, we find a nonmonotonic solute-size dependence, exhibiting extrema close to the well-known structural crossover length scale for hydrophobic hydration. Additionally, we find an apparent anomalous diffusion for water moving parallel to the surface of small solutes, which, however, can be explained by topology effects. Our findings regarding the intimate connection between solute curvature and water dynamics has implications for our understanding of hydration dynamics at heterogeneous biomolecular surfaces. DOI: 10.1103/PhysRevLett.l 14.187802
PACS numbers: 61.20.Ja, 05.20.Jj, 61.25.Em, 82.30.Rs
One of the greatest advances in our understanding of the hydrophobic effect is the recognition that the hydration structure and thermodynamics of apolar solutes is qualita tively length scale dependent [1-11 ]. The microscopic reason is that water structures very differently at small (convex) solutes, where the bulk H-bond network is only moderately deformed, as compared to large solutes, which significantly distort the tetrahedral bulk structure. The structural crossover happens at subnanometer length scales and has important implications for the interpretation of the structure and thermodynamics of hydrophobically driven assembly proc esses [1,12], such as protein folding and association [13-15], The dynamics of the hydration layer that surrounds molecular self-assemblies and proteins in solution has attracted plentiful interest in the past decade [16], Solute fluctuations and hydration dynamics are understood to be highly coupled with important consequences to biological function, such as enzyme catalysis and molecular recog nition in binding [16-26]. Despite the obvious importance of the solute chemical composition, apparently the intrinsic topological and geometric features play an important role as well [19-21,27], possibly even leading to anomalous dif fusion behavior [20]. Therefore, and due to the established fact that water considerably restructures at radii of curvature close to the subnanometer scale, a natural question to ask is, how does the water structural crossover affect the dynamics of water in the hydration layers of solutes? One of the first simulation studies of curvature effects on hydration dynamics was performed by Chau, Forester, and Smith for three solute radii between 0.35 and 0.8 nm [28], They found a slowdown of the diffusion of water in the first hydration shell relative to the bulk with an apparent minimum for the intermediate solute size. This interesting finding was not commented on, probably due to the little amount of data and statistical uncertainty of the results.
0 0 3 1 -9 0 0 7 /1 5 /1 14(18)/187802(5)
Furthermore, a slowdown of water reorientational dynam ics was found compared to bulk, an effect that decreased with solute size [28], The reorientational slowing-down has been recently explained by excluded-volume effects on the H-bond exchange dynamics [29,30], In disagreement with the simulations with Chau, Forester, and Smith, however, the excluded-volume concept predicts a mono tonic increase of the reorientation times [29], an apparent controversy which has not been addressed in literature, yet. In this Letter, we report on systematic molecular dynamics (MD) simulations of water around hydrophobic spherical model solutes of varying radius between 0.3 and 2.1 nm and investigate intrinsic curvature effects on water hydration dynamics. We find that the perpendicular and lateral diffu sivity of hydration water exhibit a nonmonotonic curvature dependence, with temperature-dependent minima located close to the structural crossover length scale [7,31]. Both curvature and temperature dependence are strikingly similar to that of entropy during hydrophobic solvation [31] and thus allow strong conjecture that the hydrophobic effect extends beyond thermodynamic into dynamic anomalies. Furthermore, we find that the intriguing curvature depend ence of the diffusivity is found to be related to a nonmonotonic dependence of H-bond lifetimes on curvature, exhibiting a maxima at the crossover length scale. This finding reconciles previous, apparently contradicting results from simulations [28] and theory [29,30] on reorientation times, as we find both predicted behaviors but in different length scale regimes. The intimate connection between solute curvature and water structure and dynamics should be of fundamental importance for biomolecular function mediated by heterogeneous bio molecular surfaces [16,17,19-21,27]. Our MD simulations are performed each containing a single and fixed hydrophobic model solute solvated in the simple point charge extended (SPC/E) model of water [32]
187802-1
©2 0 1 5 American Physical Society
PRL 114, 187802 (2015)
PHYSICAL
REVIEW
by using the Gromacs simulation package [33]. The solutewater interaction is mediated by a shifted Lennard-Jones potential Uu(r') — 4e,[(a /r/) 12 - (a/V )6], whereas r7 = r — r0 describes a coordinate shift by r0. With such a shifted potential, it is guaranteed that by changing the solute size only intrinsic curvature effects are probed, not, e.g., those of an additionally varying dispersion attraction. Values of a — 3.768 A and e = 1.197 kJ/ mol are chosen from a model for methane to mimic a reasonable dispersion attraction with a well-defined first solvation shell [6], We choose shift radii r0 = 0 ,1 ,2 ,4 ,6 , 8,10, 12.5,15, and 17.5 A in ten separate simulations. As a limiting case of a hydrophobic surface with zero curvature, the box delimiting walls in the x-y plane in a pseudo-2D simulation are chosen to interact with water by a 12-6 potential in the z direction Ux2- 6(z) — 4e[(o/ z)xl - (o/ z)6]. After 100 ps of Gibbs ensemble (NP T ) equilibration, canonical (N V T ) production runs are performed up to 200 ns at a temperature of T = 300 K and with N = 6000 to 12 000 water molecules. Further details on the simulation technique can be found in Supplemental Material [34], A simulation snapshot and water oxygen radial distribution functions (RDFs) around all the solutes are presented in Fig. 1. We first characterize the water diffusion parallel to the solute surface. For this, we calculate the lateral mean square displacement (MSD) of the arc length a water oxygen has traveled in time f, via
(W H J& gW + O -W )* -
week ending 8 MAY 2015
LETTERS
The lateral MSDs for water around all solutes are shown in Fig. 2(a) along with those of water at the planar hydrophobic surface and in bulk. Comparison to the bulk MSD indicates a crossover from below to above bulk water self-diffusion with increasing slopes with growing solute size. Additionally, the hydrating water of the smallest solutes exhibits nonlinear, anomalous behavior in t. This can be better recognized by the negative slope for small solutes in Fig. 2(b), where we plot the MSD divided by time: (S(t))/4t. Linearity of the MSDs is restored for larger solutes, converging towards the limiting MSD r0 -*■ oo of the planar hydrophobic surface. The apparent anomalous diffusion behavior for small solutes can be explained by intrinsic curvature effects, which modify the standard 2D diffusion law [37]. Here, the probability distribution function (PDF) of diffusion on spherical surfaces reads P{6, t) =
T
v /0sin(0)e_(02/2T)>
(2)
where r = 2D ||t/R avg and N(z) is a normalization constant. By expanding Eq. (2) up to second order for small 0, we find that the second moment of the PDF can be written for small solutes as
4-Dh2
(02)R2vg« 4 D ||t - - - L ,
(3)
3 ^avg
establishing a first-order curvature correction to the 2D diffusion law. On surfaces with high curvature
(i)
Here Ravg = f g ' dr rgso(r) denotes the average distance of the water molecules inside the first hydration layer, which is delimited by the location R , of the first minimum in the RDF (cf. Fig. 1), to the solute center. The time average (...)R over starting times t’ is taken only for water molecules that reside in the first hydration shell. The variable 6(t) is the azimuthal angle of the water oxygensolute distance vector at time t. For water molecules diffusing at a wall, Eq. (1) reduces to the usual twodimensional MSD ([A.r(t) + Ay(f)]2)X] in the x-y plane.
r [nm] FIG. 1 (color online). Left: MD simulation snapshot of a hydrophobic sphere with shift radius r0 = 8 A in explicit water. Right: Radial distribution function (RDF) gso(r) of water around differently sized hydrophobic model solutes.
FIG. 2 (color online), (a) MSD of water molecules parallel to differently curved hydrophobic surfaces and MSD of water in bulk (thick black dashed line), (b) The data from the top panel (a) divided by At. The inset shows the curves for water at the five smallest solutes along with the corresponding fits from Eq. (3).
187802-2
PRL 114, 187802 (2015)
PHYSICAL
REVIEW
(/?aVg < ZD||r), the second term on the right-hand side of Eq. (3) slows down the MSD as found in our simulations. In the limit Rmg -* oo, this term vanishes and the MSD on planar surfaces, oc.4D^t, is restored, consistent with the vanishing nonlinearity as observed in Fig. 2. By fitting Eq. (3) to our MSDs, we obtain the lateral diffusion constant Dy; see the inset in Fig. 2(b) for the fits. The results are plotted versus the Gibbs radius of the solutes [we defined where gso (r = RG) = 0.5] in Fig. 3 at T — 300 K normalized by the bulk diffusion constant Dbulk of water. For small solutes, we find that the parallel diffusion D\\ slows down with increasing solute sizes, which leads to parallel mobilities smaller than in bulk. This trend changes at a minimum at a length scale of RG « 0.45 nm, beyond which the lateral diffusivity monotonically rises with decreasing surface curvature to become bulklike at about R*G k 0.75 nm. Finally, it saturates to a water mobility higher than in bulk, as is well known for smooth and planar hydrophobic surfaces [38]. The significant length scales observed here match well the structural crossover length scale of —0.5 nm found for SPC/E water estimated from solvation free energies of spherical model solutes [7], To calculate the water diffusivity perpendicular to the hydrophobic surfaces, we employ the mean first-passage time (MFPT) analysis introduced by Hinczewski et al. [39,40], The MFPT Tfp(r,rt) describes the mean time required by a molecule to travel from distance r to the solute to a target distance r, > r and has an exact solution in a Smoluchowski description in terms of the free energy F(r) and diffusivity profile D± (r) [41], Its inversion leads to an expression for the perpendicular diffusivity profile [39,40]
FIG. 3 (color online). Parallel and perpendicular diffusivities D || (green dashed lines) and D (red solid lines) scaled by bulk water diffusion D bulk against solute size Rc for T = 300 K. For the four smallest cavities, the temperature effect on D ± and D is shown for T = 260, 280, and 320 K. The gray dashed horizontal line represents bulk water diffusivity. The values corresponding to diffusivity at the limiting case of zero curvature (Rc —> oo, only T = 300 K) are drawn as horizontal arrows with respective color and line coding.
week ending 8 MAY 2015
LETTERS
ePF(r) D±(r)
dTf p (r, rt) / d r j
/>
dde kF(r)
(4)
with / L 1 — kBT and rmin being a reflective boundary close to the solute where F(rmm) = l0kBT. From our MD simulations, we extract MFPTs T(r, r, = r0 + 14 A) as shown in Fig. 4(b) along with free energy profiles F(r) in Fig. 4(a). The latter result from a simple Boltzmann inversion of the solute-water RDFs plus the entropic con tribution F(r) = - k BT\ngS0(r) - 2kBT\nr [40]. Figure 4(c) plots the resulting diffusivity profiles F)L(r)/ Dhu]k normal ized by bulk diffusivity for all solute sizes. The general shape of the profiles describes a maximum inside the first solvation shell followed by a minimum in diffusivity which is reached at the outside margin of the solvation layer, whereafter diffusivity converges towards bulk diffusion in an oscillatory fashion. (The final rise of the profiles for radii close to the target distance is an artifact of non-Markovian contributions [39,40]; see also [34] for technical details.) In our calculated profiles, an interesting nonmonotonic curvature dependence of the perpendicular water mobility near hydrophobic surfaces is visible in the extrema of the profiles. With growing solute size, the maximum of diffusivity in the first hydration shell continues to increase, which leads to high
r [nm] FIG. 4 (color online), (a) Free energy profile F(r) for water around differently sized hydrophobic solutes obtained by the Boltzmann inversion of the water RDFs in Fig 1. (b) MFPT curves TfP(r,r,) for water molecules to reach the target distance at r, =
14 A given they started at a distance r < r,. (c) Perpendicular diffusivity profiles D±(r) obtained from Eq. (4).
187802-3
PHYSICAL
PRL 114, 187802 (2015)
REVIEW
perpendicular diffusivity in the vicinity of weakly curved surfaces, whereas the minimum becomes less pronounced and almost vanishes in the profile at the planar limit (see [34]). We now define the perpendicular diffusion coefficient D i (Rmax) m the first hydration shell by the value of the diffusivity curve at Rmax, the position of the first peak of the RDF. The scaled value Oj_(Nmax) /D buIk is presented next to the scaled values of Dy in Fig. 3 versus solute size Rc for T = 300 K. Like D p the perpendicular diffusivity also shows the crossover from below-to-above bulk mobility but at a smaller distance of about 0.4 nm, indicating a preceding nonmonotonic solute-size dependence for small solutes as it is interpolated towards the limiting (bulk) case of a solute with vanishing radius. Since the magnitude of the perpendicular diffusivity slightly depends on the definition of what radius defines the first hydration shell, or if over all water molecules in this shell should be averaged, we compare those different definitions in Supplemental Material [34]. We find similar values; in particular, all tested definitions rigorously reproduce the crossover from belowto-above bulk mobility behavior of water versus solute size close to the structural crossover length scale. Additionally, Fig. 3 plots diffusivity changes upon three different temperatures T — 260,280, and 320 K for solutes of sizes at which the dynamic anomaly occurs. Decreasing temperature shifts the curves and the dynamic anomaly towards larger radii of curvature. Hence, the observed temperature dependence of diffusive dynamics obeys the same trend as the thermodynamic and structural crossover length scale [31]. Taking into account entropy scaling laws for diffusion [42,43], which scale exponentially with excess entropy, suggests entropy to be a constitutive measure for the dynamic anomaly, thus firmly corroborating our conclusions. Ultimately, water hydrogen bond (HB) kinetics and mobility are intimately linked because translational diffu sion is accompanied by processes breaking, forming, and reforming HBs. The cooperation of elemental dynamic processes in water is successfully described with reactiondiffusion models [44-46] and was probed by simulations [44-47], Here, we probe HB lifetimes z which are estimated by the correlation function c(t) = (h(t)h(0 ))/ {h2) of the HB operator h(t), which is 1 for a specific pair of water molecules while they are bonded and 0 otherwise [44-46]. The negative derivative k(t) — - d c ( t ) / d t is the reactive flux hydrogen bond correlation function. We estimate “intermittent” HB lifetimes by the zero frequency part [45], the integral of c(t);
/
co r co tk(t)dt — / c(t)dt.
o
(5)
Jo
It is significantly influenced by diffusion leading to the separation of initially bonded water molecule pairs after
LETTERS
week ending 8 MAY 2015
HB breaking. Without this separation, breaking a HB between two water molecules is reversible, and the kinetics associated with this process are much faster than the irreversible rearrangements of the water hydrogen bond network analyzed here [45]. The correlation functions and details on the H-bond definition and numerical evaluation can be found in Supplemental Material [34]. The results for z (RG) scaled by their bulk value Tbuik = 3.4 ps are plotted in Fig. 5. They exhibit a nonmonotonic size dependence with a maximum retardation of about 25% versus bulk at roughly RG — 0.45 nm, close to the structural cross over length scale. After the maximum, the retardation decreases for solute sizes R G > 0.5 nm down to a remaining «5% near the planar surface. Concurrently, the average number of HBs per water molecule n / n bu{k, also shown in Fig. 5, decreases monotonically. Hence, the nonmonotonic translational mobility behavior of water is solely reflected in the kinetic behavior of the HBs. The observed behavior reconciles previous, apparently contradicting results from simulations [28] and theory [29,30] on water reorientation times, as we find both predicted trends (increasing and decreasing lifetimes), but in different solute size regimes, roughly separated by the important structural crossover scale. Figure 5 also plots the T dependence on the HB correlation times where a significant slowing-down with decreasing temperature can be observed with not much change of the position of the maximum on the solute-size axis. We note that the latter behavior cannot be related one to one to our observed diffusivity behavior (in Fig. 3) due to the multiple time scales entering z [34,45],
FIG. 5 (color online). Mean HB lifetime of first-solvation shell water t (Rc ) versus solute size scaled by the bulk value rbu|k (orange symbols). The average number of HBs per water molecule n(Rc ) scaled by its bulk value nbulk (blue errorbars). The values corresponding to the respective measures at the limiting case of zero curvature (RG -> oo, only T = 300 K) are drawn as horizontal arrows with respective color coding. The error bars for n/ nbu]k were estimated from block averages as described in Ref. [34],
187802-4
PRL 114, 187802 (2015)
PHYSICAL
REVIEW
In summary, we have established a firm link between the structure and dynamics of hydration water around hydrophobic solutes with a novel dynamic anomaly happening at the well-known crossover length scale. Because of the fundamental importance of surface water to biomolecular processes and function [16,26], in particular, at topologi cally heterogeneous protein surfaces [19-21,27], our find ings imply that nature has the means to employ local surface topology to mediate biological function. Hence, our results will help in the interpretation of experimentally found dynamic heterogeneities on biomolecular surfaces [16-21]. Locally slowed-down water, for instance, could thus fine-tune the folding kinetics of hydrophobic polymers and peptides [16] or may mediate the appropriate time scales for the association of ligands to catalytically active sites or binding pockets [16,22,48], The authors thank Roland R. Netz for helpful discus sions and the Deutsche Forschungsgemeinschaft (DFG) for financial support. M. H. is grateful for financial support by the Cluster of Excellence RESOLV (EXC-1069) funded by the DFG.
[1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17]
To whom correspondence should be addressed. [email protected] D. Chandler, Nature (London) 437, 640 (2005). C. Y. Lee, J. A. McCammon, and R J. Rossky, J. Chem. Phys. 80, 4448 (1984). A. Nicholls, K. A. Sharp, and B. Honig, Proteins 11, 281 (1991). K. Lum, D. Chandler, and J. D. Weeks, J. Phys. Chem. B 103, 4570 (1999). D. M. Huang, P. L. Geissler, and D. Chandler, J. Phys. Chem. B 105, 6704 (2001). D. M. Huang and D. Chandler, J. Phys. Chem. B 106, 2047 ( 2002 ). S. Rajamani, T. M. Truskett, and S. Garde, Proc. Natl. Acad. Sci. U.S.A. 102, 9475 (2005). H. S. Ashbaugh and L. R. Pratt, Rev. Mod. Phys. 78, 159 (2006). L. R. Pratt, Annu. Rev. Phys. Chem. 53, 409 (2002). Y. S. Djikaev and E. Ruckenstein. J. Chem. Phys. 139, 184709 (2013). J. G. Davis, K. P. Gierszal, P. Wang, and D. Ben-Amotz, Nature (London) 491, 582 (2012). B.J. Beme, J. D. Weeks, and R. Zhou, Annu. Rev. Phys. Chem. 60, 85 (2009). D. M. Huang and D. Chandler, Proc. Natl. Acad. Sci. U.S.A. 97, 8324 (2000). M. V. Athawale, G. Goel, T. Ghosh, T. M. Truskett, and S. Garde, Proc. Natl. Acad. Sci. U.S.A. 104, 733 (2007). I. T. S. Li and G. C. Walker, Proc. Natl. Acad. Sci. U.S.A. 108, 16527 (2011). B. Bagchi, Chem. Rev. 105, 3197 (2005). G. Niehus, M. Heyden, D. A. Schmidt, and M. Havenith, Faraday Discuss. 150, 193 (2011).
LETTERS
week ending 8 MAY 2015
[18] O. H. Kwon, T. H. Yoo, C. M. Othon, J. A. V. Deventer, D. A. Tirrell, and A. H. Zewail, Proc. Natl. Acad. Sci. U.S.A. 107, 17101 (2010). [19] L. Zhang, L. Wang, Y.-T. Kao, W. Qiu, Y. Yang, O. Okobiah, and D. Zhong, Proc. Natl. Acad. Sci. U.S.A. 104. 18461 (2007). [20] F. Pizzitutti, M. Marchi, F. Sterpone, and P. J. Rossky, J. Phys. Chem. B 111, 7584 (2007). [21] A. C. Fogarty and D. Laage, J. Phys. Chem. B 118, 7715 (2014). [22] P. Setny, R. Baron, P. Kekenes-Huskey, J. A. McCammon, and J. Dzubiella, Proc. Natl. Acad. Sci. U.S.A. 110, 1197 (2013). [23] J. Mittal and G. Hummer, Proc. Natl. Acad. Sci. U.S.A. 105, 20130 (2008). [24] S. Vaikuntanathan and P. L. Geissler, Phys. Rev. Lett. 112, 020603 (2014). [25] A. J. Patel, P. Varilly, S. N. Jamadagni, H. Acharya, S. Garde, and D. Chandler, Proc. Natl. Acad. Sci. U.S.A. 108, 17678 (2011). [26] S.N. Jamadagni, R. Godawat, and S. Garde, Annu. Rev. Chem. Biomol. Eng. 2, 147 (2011). [27] Y. K. Cheng and P. J. Rossky, Nature (London) 392, 696 (1998). [28] P.L. Chau, Mol. Phys. 89, 1033 (1996). [29] D. Laage, G. Stirnemann, and J. T. Hynes, J. Phys. Chem. B 113, 2428 (2009). [30] N. Galamba, J. Phys. Chem. B 118, 4169 (2014). [31] H. S. Ashbaugh, Chem. Phys. Lett. 477, 109 (2009). [32] H. J. C. Berendsen, J. R. Grigera, and T. P. Straatsma, J. Phys. Chem. 91, 6269 (1987). [33] B. Hess, C. Kutzner, D. van der Spoel, and E. Lindahl, J. Chem. Theory Comput. 4, 435 (2008). [34] See Supplemental Material at http://link.aps.org/ supplemental/10.1103/PhysRevLett.l 14.187802 for details on the MD simulation setup as well as more technical details on the MFPT method and the HB analysis, which includes Refs. [35,36], [35] S. Nose, J. Chem. Phys. 81, 511 (1984). [36] W. G. Hoover, Phys. Rev. A 31, 1695 (1985). [37] A. Gosh, J. Samuel, and S. Sinha, Europhys. Lett. 98, 30003 ( 2012). [38] P. Kumar, S. V. Buldyrev, F. W. Starr, N. Giovambattista, and H.E. Stanley, Phys. Rev. E 72, 051503 (2005). [39] M. Hinczewski, Y. von Hansen, J. Dzubiella, and R. R. Netz, J. Chem. Phys. 132, 245103 (2010). [40] F. Sedlmeier, Y. von Hansen, L. Mengyu, D. Horinek, and R.R. Netz, J. Stat. Phys. 145, 240 (2011). [41] G. H. Weiss, Adv. Chem. Phys. 13, 1 (1966). [42] Y. Rosenfeld, Phys. Rev. A 15, 6 (1977). [43] M. Dzugutov, Nature (London) 381, 137 (1996). [44] A. Luzar and D. Chandler, Nature (London) 379, 55 (1996). [45] A. Luzar, J. Chem. Phys. 113, 10663 (2000). [46] O. Markovitch and N. Agmon, J. Chem. Phys. 129, 084505 (2008). [47] R. G. Pereyra, A. J. B. di Lorenzo, D. C. Malaspina, and M. A. Carignano, Chem. Phys. Lett. 538, 35 (2012). [48] M. Grossman, B. Bom, M. Heyden, D. Tworowski, G. B. Fields, I. Sagi, and M. Havenith, Nat. Struct. Biol. 18, 1102 ( 2011 ).
187802-5
Copyright of Physical Review Letters is the property of American Physical Society and its content may not be copied or emailed to multiple sites or posted to a listserv without the copyright holder's express written permission. However, users may print, download, or email articles for individual use.
PRL 114, 187802 (2015)
REVIEW
week ending 8 MAY 2015
LETTERS
Curvature Dependence of Hydrophobic Hydration Dynamics R. Gregor WeiB,1’2 Matthias Heyden,2 and Joachim Dzubiella1,2' ^ Institut fiir Physik, Humboldt-Universitdt zu Berlin, Newtonstrasse 15, D-12489 Berlin, Germany "Soft Matter and Functional Materials, Helmholtz-Zentrum Berlin, Hahn-Meitner Platz 1, D-14109 Berlin, Germany Max-Planck-Insitut fiir Kohlenforschung, Kaiser-Wilhelm-Platz 1, D-45470 Miilheim an der Ruhr, Germany (Received 25 September 2014; revised manuscript received 25 February 2015; published 5 May 2015) We investigate the solute curvature dependence of water dynamics in the vicinity of hydrophobic spherical solutes using molecular dynamics simulations. For both the lateral and perpendicular diffusivity, as well as for H-bond kinetics of water in the first hydration shell, we find a nonmonotonic solute-size dependence, exhibiting extrema close to the well-known structural crossover length scale for hydrophobic hydration. Additionally, we find an apparent anomalous diffusion for water moving parallel to the surface of small solutes, which, however, can be explained by topology effects. Our findings regarding the intimate connection between solute curvature and water dynamics has implications for our understanding of hydration dynamics at heterogeneous biomolecular surfaces. DOI: 10.1103/PhysRevLett.l 14.187802
PACS numbers: 61.20.Ja, 05.20.Jj, 61.25.Em, 82.30.Rs
One of the greatest advances in our understanding of the hydrophobic effect is the recognition that the hydration structure and thermodynamics of apolar solutes is qualita tively length scale dependent [1-11 ]. The microscopic reason is that water structures very differently at small (convex) solutes, where the bulk H-bond network is only moderately deformed, as compared to large solutes, which significantly distort the tetrahedral bulk structure. The structural crossover happens at subnanometer length scales and has important implications for the interpretation of the structure and thermodynamics of hydrophobically driven assembly proc esses [1,12], such as protein folding and association [13-15], The dynamics of the hydration layer that surrounds molecular self-assemblies and proteins in solution has attracted plentiful interest in the past decade [16], Solute fluctuations and hydration dynamics are understood to be highly coupled with important consequences to biological function, such as enzyme catalysis and molecular recog nition in binding [16-26]. Despite the obvious importance of the solute chemical composition, apparently the intrinsic topological and geometric features play an important role as well [19-21,27], possibly even leading to anomalous dif fusion behavior [20]. Therefore, and due to the established fact that water considerably restructures at radii of curvature close to the subnanometer scale, a natural question to ask is, how does the water structural crossover affect the dynamics of water in the hydration layers of solutes? One of the first simulation studies of curvature effects on hydration dynamics was performed by Chau, Forester, and Smith for three solute radii between 0.35 and 0.8 nm [28], They found a slowdown of the diffusion of water in the first hydration shell relative to the bulk with an apparent minimum for the intermediate solute size. This interesting finding was not commented on, probably due to the little amount of data and statistical uncertainty of the results.
0 0 3 1 -9 0 0 7 /1 5 /1 14(18)/187802(5)
Furthermore, a slowdown of water reorientational dynam ics was found compared to bulk, an effect that decreased with solute size [28], The reorientational slowing-down has been recently explained by excluded-volume effects on the H-bond exchange dynamics [29,30], In disagreement with the simulations with Chau, Forester, and Smith, however, the excluded-volume concept predicts a mono tonic increase of the reorientation times [29], an apparent controversy which has not been addressed in literature, yet. In this Letter, we report on systematic molecular dynamics (MD) simulations of water around hydrophobic spherical model solutes of varying radius between 0.3 and 2.1 nm and investigate intrinsic curvature effects on water hydration dynamics. We find that the perpendicular and lateral diffu sivity of hydration water exhibit a nonmonotonic curvature dependence, with temperature-dependent minima located close to the structural crossover length scale [7,31]. Both curvature and temperature dependence are strikingly similar to that of entropy during hydrophobic solvation [31] and thus allow strong conjecture that the hydrophobic effect extends beyond thermodynamic into dynamic anomalies. Furthermore, we find that the intriguing curvature depend ence of the diffusivity is found to be related to a nonmonotonic dependence of H-bond lifetimes on curvature, exhibiting a maxima at the crossover length scale. This finding reconciles previous, apparently contradicting results from simulations [28] and theory [29,30] on reorientation times, as we find both predicted behaviors but in different length scale regimes. The intimate connection between solute curvature and water structure and dynamics should be of fundamental importance for biomolecular function mediated by heterogeneous bio molecular surfaces [16,17,19-21,27]. Our MD simulations are performed each containing a single and fixed hydrophobic model solute solvated in the simple point charge extended (SPC/E) model of water [32]
187802-1
©2 0 1 5 American Physical Society
PRL 114, 187802 (2015)
PHYSICAL
REVIEW
by using the Gromacs simulation package [33]. The solutewater interaction is mediated by a shifted Lennard-Jones potential Uu(r') — 4e,[(a /r/) 12 - (a/V )6], whereas r7 = r — r0 describes a coordinate shift by r0. With such a shifted potential, it is guaranteed that by changing the solute size only intrinsic curvature effects are probed, not, e.g., those of an additionally varying dispersion attraction. Values of a — 3.768 A and e = 1.197 kJ/ mol are chosen from a model for methane to mimic a reasonable dispersion attraction with a well-defined first solvation shell [6], We choose shift radii r0 = 0 ,1 ,2 ,4 ,6 , 8,10, 12.5,15, and 17.5 A in ten separate simulations. As a limiting case of a hydrophobic surface with zero curvature, the box delimiting walls in the x-y plane in a pseudo-2D simulation are chosen to interact with water by a 12-6 potential in the z direction Ux2- 6(z) — 4e[(o/ z)xl - (o/ z)6]. After 100 ps of Gibbs ensemble (NP T ) equilibration, canonical (N V T ) production runs are performed up to 200 ns at a temperature of T = 300 K and with N = 6000 to 12 000 water molecules. Further details on the simulation technique can be found in Supplemental Material [34], A simulation snapshot and water oxygen radial distribution functions (RDFs) around all the solutes are presented in Fig. 1. We first characterize the water diffusion parallel to the solute surface. For this, we calculate the lateral mean square displacement (MSD) of the arc length a water oxygen has traveled in time f, via
(W H J& gW + O -W )* -
week ending 8 MAY 2015
LETTERS
The lateral MSDs for water around all solutes are shown in Fig. 2(a) along with those of water at the planar hydrophobic surface and in bulk. Comparison to the bulk MSD indicates a crossover from below to above bulk water self-diffusion with increasing slopes with growing solute size. Additionally, the hydrating water of the smallest solutes exhibits nonlinear, anomalous behavior in t. This can be better recognized by the negative slope for small solutes in Fig. 2(b), where we plot the MSD divided by time: (S(t))/4t. Linearity of the MSDs is restored for larger solutes, converging towards the limiting MSD r0 -*■ oo of the planar hydrophobic surface. The apparent anomalous diffusion behavior for small solutes can be explained by intrinsic curvature effects, which modify the standard 2D diffusion law [37]. Here, the probability distribution function (PDF) of diffusion on spherical surfaces reads P{6, t) =
T
v /0sin(0)e_(02/2T)>
(2)
where r = 2D ||t/R avg and N(z) is a normalization constant. By expanding Eq. (2) up to second order for small 0, we find that the second moment of the PDF can be written for small solutes as
4-Dh2
(02)R2vg« 4 D ||t - - - L ,
(3)
3 ^avg
establishing a first-order curvature correction to the 2D diffusion law. On surfaces with high curvature
(i)
Here Ravg = f g ' dr rgso(r) denotes the average distance of the water molecules inside the first hydration layer, which is delimited by the location R , of the first minimum in the RDF (cf. Fig. 1), to the solute center. The time average (...)R over starting times t’ is taken only for water molecules that reside in the first hydration shell. The variable 6(t) is the azimuthal angle of the water oxygensolute distance vector at time t. For water molecules diffusing at a wall, Eq. (1) reduces to the usual twodimensional MSD ([A.r(t) + Ay(f)]2)X] in the x-y plane.
r [nm] FIG. 1 (color online). Left: MD simulation snapshot of a hydrophobic sphere with shift radius r0 = 8 A in explicit water. Right: Radial distribution function (RDF) gso(r) of water around differently sized hydrophobic model solutes.
FIG. 2 (color online), (a) MSD of water molecules parallel to differently curved hydrophobic surfaces and MSD of water in bulk (thick black dashed line), (b) The data from the top panel (a) divided by At. The inset shows the curves for water at the five smallest solutes along with the corresponding fits from Eq. (3).
187802-2
PRL 114, 187802 (2015)
PHYSICAL
REVIEW
(/?aVg < ZD||r), the second term on the right-hand side of Eq. (3) slows down the MSD as found in our simulations. In the limit Rmg -* oo, this term vanishes and the MSD on planar surfaces, oc.4D^t, is restored, consistent with the vanishing nonlinearity as observed in Fig. 2. By fitting Eq. (3) to our MSDs, we obtain the lateral diffusion constant Dy; see the inset in Fig. 2(b) for the fits. The results are plotted versus the Gibbs radius of the solutes [we defined where gso (r = RG) = 0.5] in Fig. 3 at T — 300 K normalized by the bulk diffusion constant Dbulk of water. For small solutes, we find that the parallel diffusion D\\ slows down with increasing solute sizes, which leads to parallel mobilities smaller than in bulk. This trend changes at a minimum at a length scale of RG « 0.45 nm, beyond which the lateral diffusivity monotonically rises with decreasing surface curvature to become bulklike at about R*G k 0.75 nm. Finally, it saturates to a water mobility higher than in bulk, as is well known for smooth and planar hydrophobic surfaces [38]. The significant length scales observed here match well the structural crossover length scale of —0.5 nm found for SPC/E water estimated from solvation free energies of spherical model solutes [7], To calculate the water diffusivity perpendicular to the hydrophobic surfaces, we employ the mean first-passage time (MFPT) analysis introduced by Hinczewski et al. [39,40], The MFPT Tfp(r,rt) describes the mean time required by a molecule to travel from distance r to the solute to a target distance r, > r and has an exact solution in a Smoluchowski description in terms of the free energy F(r) and diffusivity profile D± (r) [41], Its inversion leads to an expression for the perpendicular diffusivity profile [39,40]
FIG. 3 (color online). Parallel and perpendicular diffusivities D || (green dashed lines) and D (red solid lines) scaled by bulk water diffusion D bulk against solute size Rc for T = 300 K. For the four smallest cavities, the temperature effect on D ± and D is shown for T = 260, 280, and 320 K. The gray dashed horizontal line represents bulk water diffusivity. The values corresponding to diffusivity at the limiting case of zero curvature (Rc —> oo, only T = 300 K) are drawn as horizontal arrows with respective color and line coding.
week ending 8 MAY 2015
LETTERS
ePF(r) D±(r)
dTf p (r, rt) / d r j
/>
dde kF(r)
(4)
with / L 1 — kBT and rmin being a reflective boundary close to the solute where F(rmm) = l0kBT. From our MD simulations, we extract MFPTs T(r, r, = r0 + 14 A) as shown in Fig. 4(b) along with free energy profiles F(r) in Fig. 4(a). The latter result from a simple Boltzmann inversion of the solute-water RDFs plus the entropic con tribution F(r) = - k BT\ngS0(r) - 2kBT\nr [40]. Figure 4(c) plots the resulting diffusivity profiles F)L(r)/ Dhu]k normal ized by bulk diffusivity for all solute sizes. The general shape of the profiles describes a maximum inside the first solvation shell followed by a minimum in diffusivity which is reached at the outside margin of the solvation layer, whereafter diffusivity converges towards bulk diffusion in an oscillatory fashion. (The final rise of the profiles for radii close to the target distance is an artifact of non-Markovian contributions [39,40]; see also [34] for technical details.) In our calculated profiles, an interesting nonmonotonic curvature dependence of the perpendicular water mobility near hydrophobic surfaces is visible in the extrema of the profiles. With growing solute size, the maximum of diffusivity in the first hydration shell continues to increase, which leads to high
r [nm] FIG. 4 (color online), (a) Free energy profile F(r) for water around differently sized hydrophobic solutes obtained by the Boltzmann inversion of the water RDFs in Fig 1. (b) MFPT curves TfP(r,r,) for water molecules to reach the target distance at r, =
14 A given they started at a distance r < r,. (c) Perpendicular diffusivity profiles D±(r) obtained from Eq. (4).
187802-3
PHYSICAL
PRL 114, 187802 (2015)
REVIEW
perpendicular diffusivity in the vicinity of weakly curved surfaces, whereas the minimum becomes less pronounced and almost vanishes in the profile at the planar limit (see [34]). We now define the perpendicular diffusion coefficient D i (Rmax) m the first hydration shell by the value of the diffusivity curve at Rmax, the position of the first peak of the RDF. The scaled value Oj_(Nmax) /D buIk is presented next to the scaled values of Dy in Fig. 3 versus solute size Rc for T = 300 K. Like D p the perpendicular diffusivity also shows the crossover from below-to-above bulk mobility but at a smaller distance of about 0.4 nm, indicating a preceding nonmonotonic solute-size dependence for small solutes as it is interpolated towards the limiting (bulk) case of a solute with vanishing radius. Since the magnitude of the perpendicular diffusivity slightly depends on the definition of what radius defines the first hydration shell, or if over all water molecules in this shell should be averaged, we compare those different definitions in Supplemental Material [34]. We find similar values; in particular, all tested definitions rigorously reproduce the crossover from belowto-above bulk mobility behavior of water versus solute size close to the structural crossover length scale. Additionally, Fig. 3 plots diffusivity changes upon three different temperatures T — 260,280, and 320 K for solutes of sizes at which the dynamic anomaly occurs. Decreasing temperature shifts the curves and the dynamic anomaly towards larger radii of curvature. Hence, the observed temperature dependence of diffusive dynamics obeys the same trend as the thermodynamic and structural crossover length scale [31]. Taking into account entropy scaling laws for diffusion [42,43], which scale exponentially with excess entropy, suggests entropy to be a constitutive measure for the dynamic anomaly, thus firmly corroborating our conclusions. Ultimately, water hydrogen bond (HB) kinetics and mobility are intimately linked because translational diffu sion is accompanied by processes breaking, forming, and reforming HBs. The cooperation of elemental dynamic processes in water is successfully described with reactiondiffusion models [44-46] and was probed by simulations [44-47], Here, we probe HB lifetimes z which are estimated by the correlation function c(t) = (h(t)h(0 ))/ {h2) of the HB operator h(t), which is 1 for a specific pair of water molecules while they are bonded and 0 otherwise [44-46]. The negative derivative k(t) — - d c ( t ) / d t is the reactive flux hydrogen bond correlation function. We estimate “intermittent” HB lifetimes by the zero frequency part [45], the integral of c(t);
/
co r co tk(t)dt — / c(t)dt.
o
(5)
Jo
It is significantly influenced by diffusion leading to the separation of initially bonded water molecule pairs after
LETTERS
week ending 8 MAY 2015
HB breaking. Without this separation, breaking a HB between two water molecules is reversible, and the kinetics associated with this process are much faster than the irreversible rearrangements of the water hydrogen bond network analyzed here [45]. The correlation functions and details on the H-bond definition and numerical evaluation can be found in Supplemental Material [34]. The results for z (RG) scaled by their bulk value Tbuik = 3.4 ps are plotted in Fig. 5. They exhibit a nonmonotonic size dependence with a maximum retardation of about 25% versus bulk at roughly RG — 0.45 nm, close to the structural cross over length scale. After the maximum, the retardation decreases for solute sizes R G > 0.5 nm down to a remaining «5% near the planar surface. Concurrently, the average number of HBs per water molecule n / n bu{k, also shown in Fig. 5, decreases monotonically. Hence, the nonmonotonic translational mobility behavior of water is solely reflected in the kinetic behavior of the HBs. The observed behavior reconciles previous, apparently contradicting results from simulations [28] and theory [29,30] on water reorientation times, as we find both predicted trends (increasing and decreasing lifetimes), but in different solute size regimes, roughly separated by the important structural crossover scale. Figure 5 also plots the T dependence on the HB correlation times where a significant slowing-down with decreasing temperature can be observed with not much change of the position of the maximum on the solute-size axis. We note that the latter behavior cannot be related one to one to our observed diffusivity behavior (in Fig. 3) due to the multiple time scales entering z [34,45],
FIG. 5 (color online). Mean HB lifetime of first-solvation shell water t (Rc ) versus solute size scaled by the bulk value rbu|k (orange symbols). The average number of HBs per water molecule n(Rc ) scaled by its bulk value nbulk (blue errorbars). The values corresponding to the respective measures at the limiting case of zero curvature (RG -> oo, only T = 300 K) are drawn as horizontal arrows with respective color coding. The error bars for n/ nbu]k were estimated from block averages as described in Ref. [34],
187802-4
PRL 114, 187802 (2015)
PHYSICAL
REVIEW
In summary, we have established a firm link between the structure and dynamics of hydration water around hydrophobic solutes with a novel dynamic anomaly happening at the well-known crossover length scale. Because of the fundamental importance of surface water to biomolecular processes and function [16,26], in particular, at topologi cally heterogeneous protein surfaces [19-21,27], our find ings imply that nature has the means to employ local surface topology to mediate biological function. Hence, our results will help in the interpretation of experimentally found dynamic heterogeneities on biomolecular surfaces [16-21]. Locally slowed-down water, for instance, could thus fine-tune the folding kinetics of hydrophobic polymers and peptides [16] or may mediate the appropriate time scales for the association of ligands to catalytically active sites or binding pockets [16,22,48], The authors thank Roland R. Netz for helpful discus sions and the Deutsche Forschungsgemeinschaft (DFG) for financial support. M. H. is grateful for financial support by the Cluster of Excellence RESOLV (EXC-1069) funded by the DFG.
[1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17]
To whom correspondence should be addressed. [email protected] D. Chandler, Nature (London) 437, 640 (2005). C. Y. Lee, J. A. McCammon, and R J. Rossky, J. Chem. Phys. 80, 4448 (1984). A. Nicholls, K. A. Sharp, and B. Honig, Proteins 11, 281 (1991). K. Lum, D. Chandler, and J. D. Weeks, J. Phys. Chem. B 103, 4570 (1999). D. M. Huang, P. L. Geissler, and D. Chandler, J. Phys. Chem. B 105, 6704 (2001). D. M. Huang and D. Chandler, J. Phys. Chem. B 106, 2047 ( 2002 ). S. Rajamani, T. M. Truskett, and S. Garde, Proc. Natl. Acad. Sci. U.S.A. 102, 9475 (2005). H. S. Ashbaugh and L. R. Pratt, Rev. Mod. Phys. 78, 159 (2006). L. R. Pratt, Annu. Rev. Phys. Chem. 53, 409 (2002). Y. S. Djikaev and E. Ruckenstein. J. Chem. Phys. 139, 184709 (2013). J. G. Davis, K. P. Gierszal, P. Wang, and D. Ben-Amotz, Nature (London) 491, 582 (2012). B.J. Beme, J. D. Weeks, and R. Zhou, Annu. Rev. Phys. Chem. 60, 85 (2009). D. M. Huang and D. Chandler, Proc. Natl. Acad. Sci. U.S.A. 97, 8324 (2000). M. V. Athawale, G. Goel, T. Ghosh, T. M. Truskett, and S. Garde, Proc. Natl. Acad. Sci. U.S.A. 104, 733 (2007). I. T. S. Li and G. C. Walker, Proc. Natl. Acad. Sci. U.S.A. 108, 16527 (2011). B. Bagchi, Chem. Rev. 105, 3197 (2005). G. Niehus, M. Heyden, D. A. Schmidt, and M. Havenith, Faraday Discuss. 150, 193 (2011).
LETTERS
week ending 8 MAY 2015
[18] O. H. Kwon, T. H. Yoo, C. M. Othon, J. A. V. Deventer, D. A. Tirrell, and A. H. Zewail, Proc. Natl. Acad. Sci. U.S.A. 107, 17101 (2010). [19] L. Zhang, L. Wang, Y.-T. Kao, W. Qiu, Y. Yang, O. Okobiah, and D. Zhong, Proc. Natl. Acad. Sci. U.S.A. 104. 18461 (2007). [20] F. Pizzitutti, M. Marchi, F. Sterpone, and P. J. Rossky, J. Phys. Chem. B 111, 7584 (2007). [21] A. C. Fogarty and D. Laage, J. Phys. Chem. B 118, 7715 (2014). [22] P. Setny, R. Baron, P. Kekenes-Huskey, J. A. McCammon, and J. Dzubiella, Proc. Natl. Acad. Sci. U.S.A. 110, 1197 (2013). [23] J. Mittal and G. Hummer, Proc. Natl. Acad. Sci. U.S.A. 105, 20130 (2008). [24] S. Vaikuntanathan and P. L. Geissler, Phys. Rev. Lett. 112, 020603 (2014). [25] A. J. Patel, P. Varilly, S. N. Jamadagni, H. Acharya, S. Garde, and D. Chandler, Proc. Natl. Acad. Sci. U.S.A. 108, 17678 (2011). [26] S.N. Jamadagni, R. Godawat, and S. Garde, Annu. Rev. Chem. Biomol. Eng. 2, 147 (2011). [27] Y. K. Cheng and P. J. Rossky, Nature (London) 392, 696 (1998). [28] P.L. Chau, Mol. Phys. 89, 1033 (1996). [29] D. Laage, G. Stirnemann, and J. T. Hynes, J. Phys. Chem. B 113, 2428 (2009). [30] N. Galamba, J. Phys. Chem. B 118, 4169 (2014). [31] H. S. Ashbaugh, Chem. Phys. Lett. 477, 109 (2009). [32] H. J. C. Berendsen, J. R. Grigera, and T. P. Straatsma, J. Phys. Chem. 91, 6269 (1987). [33] B. Hess, C. Kutzner, D. van der Spoel, and E. Lindahl, J. Chem. Theory Comput. 4, 435 (2008). [34] See Supplemental Material at http://link.aps.org/ supplemental/10.1103/PhysRevLett.l 14.187802 for details on the MD simulation setup as well as more technical details on the MFPT method and the HB analysis, which includes Refs. [35,36], [35] S. Nose, J. Chem. Phys. 81, 511 (1984). [36] W. G. Hoover, Phys. Rev. A 31, 1695 (1985). [37] A. Gosh, J. Samuel, and S. Sinha, Europhys. Lett. 98, 30003 ( 2012). [38] P. Kumar, S. V. Buldyrev, F. W. Starr, N. Giovambattista, and H.E. Stanley, Phys. Rev. E 72, 051503 (2005). [39] M. Hinczewski, Y. von Hansen, J. Dzubiella, and R. R. Netz, J. Chem. Phys. 132, 245103 (2010). [40] F. Sedlmeier, Y. von Hansen, L. Mengyu, D. Horinek, and R.R. Netz, J. Stat. Phys. 145, 240 (2011). [41] G. H. Weiss, Adv. Chem. Phys. 13, 1 (1966). [42] Y. Rosenfeld, Phys. Rev. A 15, 6 (1977). [43] M. Dzugutov, Nature (London) 381, 137 (1996). [44] A. Luzar and D. Chandler, Nature (London) 379, 55 (1996). [45] A. Luzar, J. Chem. Phys. 113, 10663 (2000). [46] O. Markovitch and N. Agmon, J. Chem. Phys. 129, 084505 (2008). [47] R. G. Pereyra, A. J. B. di Lorenzo, D. C. Malaspina, and M. A. Carignano, Chem. Phys. Lett. 538, 35 (2012). [48] M. Grossman, B. Bom, M. Heyden, D. Tworowski, G. B. Fields, I. Sagi, and M. Havenith, Nat. Struct. Biol. 18, 1102 ( 2011 ).
187802-5
Copyright of Physical Review Letters is the property of American Physical Society and its content may not be copied or emailed to multiple sites or posted to a listserv without the copyright holder's express written permission. However, users may print, download, or email articles for individual use.
