8 March 1996
ELSEVIER
OHISM|CAL PHYS|OS LEYTERS Chemical Physics Letters 250 (1996) 443-449
Effects of simple model solutes on water dynamics: residence time analysis F. Brug~
a.b,c,
E. Parisi
a,b,
S.L. Fornili a,b,c
a INFM, University of Palermo, Via A,chirafi 36. 90123 Palermo, Italy b Department of Physics, University of Palermo, Via Archirafi 36, 90123 Palermo. Italy c CNR-IAIF, Via Archirafi 36, 90123 Palermo. Italy
Received I November 1995
Abstract The water residence time around pairs of simple hydrophilic and hydrophobic solutes, kept immobilized to model groups on slowly diffusing macromolecules, has been investigated along 400 ps molecular dynamics simulations. Results show that the water mobility around one solute strongly depends on the polarity characteristics of the other solute of the pair.
1. Introduction Extensive molecular dynamicg (MD) simulations have been recently performed on aqueous solutions of the model solute pairs sketched in Fig. la [1,2] to study solute-solute solvent-induced forces (SIFs). These 'solutes' consist of immobilized 31'2 [3] water molecules (PHIL) o1" neon-like [3] Lennard-Jones (LJ) particles (PHOB). Results of these simulations suggest that SIFs can play a remarkable role in molecular mechanisms itwolving biomolecules, such as in the enzyme-substrate and receptor-ligand interactions [1,2]. Structural properties of the solvent correlated with characteristics of SIFs have already been investigated [1,2]. In this worg the water residence time within the first coordination shell of solutes is studied. Thanks to the systems simplicity, we obtain valuable information on the dynamic behavior of the solvent and interesting methodological suggestions for proper use of residence time analysis [6] even in
more complex cases, such as those concerning macromolecules.
2. Computational details The systems consisted of 'solutes' pairs (Fig. l a) in a bath of 727 ST2 water molecules confined into a cubic box of 28.52 ,g, side-length. A cut-off ~adius of 7.8 ,~, and reaction field were used [7]. The waterPHOB interaction was described through the LJ parameters of the ST2 model (e = 0.07575 kcal/mol, (r = 3.1 ~). Direct interactions between solutes were not considered. The simulations were performed in the microcanonical ensemble at an average temperature of 300 K. The Verlet integration scheme [8] was adopted with a time step of 1 fs using periodic boundary condaions under minimum image convention [9]. Intensive computing was carried out on a multitransputer system [10]. The concurrent MD program includes an efficieat distributed load balancer
0009-2614/96/$12.00 © 1996 Elsevier Science B.V. All rights leserved Pll S0009-2614(96)00038-3
444
F. Brug~ et al. / Chemical Physics Letters 250 (1996) 443-449
[11] and a parallel implementation of the "cell linked-list' algorithm [12]. The water residence time distributions, n(t) (called survival functions in what follows), were calculated using a serial implementation of an algorithm developed by Garcia and Stiller [6] following a concepttlal reasoning previously adopted to investigate the time permanence of water hydrogen bonds (FIB) [13], Nw
n(t) =n(jdt) = ~ . B , J ( N - j + 1), t=l
where is one if the molecule i resides during j contiguous time intervals dr within the solute first coordination shell, Nw is the water molecules number and N the total number of system configurations, which were stored every dt = 50 fs along the last 400 ps of each simulation. Otherwise B~s is zero.
3, Results and discussion The effective residence time, T, of water within the solute hydration shell is usually de~ermined by fitting an exponent:,al function to the survival func-
I1)
(2)
(1)
PHOB-PHOB
~
tion [14], which, however, is often a complicate function [6]. In Fig. 2a, a typical free water survival function is shown. Three different behaviors can be distinguished: a rapid initial decay, an approximately exponential decay and a final region, often very complicate. In what follows only the intermediate region of each survival functioe is analyzed, by considering time values larger than 5 ps and excluding the non-exponential tail, which is not statistically relevant. The initial decay is probably related to particles which, being c~ose to the shell border, rapidly cross it back and forth. This suggests that reliable residence time values cannot be obtained if the shell borders pass through regions where the water density is high. This is illustrated in Fig. 2b, where two survival functions are showr~. They have been calculated using spherical shell radii r --- 2.875 .A and r = 3.35 ~, which correspond to the first maximum and to the first minimum of the oxygen-oxygen pair correlation function, g(r), respectively. For the latter, the water mean residence time is 2.4 ps, in good agreement with previous findings for ST2 water, as shown in Table 1, where values are also reported for MCY
~2)
PHIL-PHIL
(1)
(2)
PHOB-PHIL
bulk
~ Z
Y
Fig. I. Upper. %u~ate' pairs referred to in the present work, and definition of the coordinate system. Tile PHIL mo~el solutes consist of ST2 [3] water molecules, while the PHOB solutes are neon-like [3] Lennard-Jones particles. Mutual distance: 4.2 ~. Lower left: cylindrical coordinates and volume used to evaluate tile number density of the water oxygen atoms. The cylinder has 28/~ length and 5 .~ radius. Lower right: spatial decomposition of the solvent region around each solute into 'interior' (dotted) and "exterior' (shaded) hemispheres.
F. Brugb et al. / Chemical Physics Letters 250 (/996)443-449
445
I02 4
I0 ~
I0' 1oo
100
104 104
= lO.Z lO-J
10-2 I04
r ffi2.875 ~. x=0.05 ps
lO-S
I
0
20
40
60
80
time (ps)
0
5
n
I
i0
I
10.3
I
15
2O
time (ps)
Fig. 2. (a) Typical 'survival' function obtained from a 400 ps molecular dynamics trajectory o f a system consisting of 729 ST2 water molecules all free In move. This curve is the average of 100 survival functions relative to randomly chosen water molecules. Spherical shell radius r ffi 5.875 ~ . (b) Water sun, ival functions (exponential parts) averaged over 100 free water molecules for two shells, whose radii correspond to the first maximum ( r = 2.875 ~ ) and to the first minimum ( r = 3.75/~) of the oxygen-oxygen pair correlation function g(r) (shown in the inset), respectively.
[15] and TIPS3 [16] models, and experimental data [17]. However, for the shell radius r--2.875 the residence time is abnormally small. For the pair systems shown in Fig. la, the number density of the water oxygen atoms has been evaluated along the last 400 ps of each simulation, considering the molecules within the cylinder shown in Fig. lb. It is calculated as follows: D( p, ~, y) - - ( N ( p, ~b, y,
t))lV,
where ( N ) indicates the time-averaged number of water oxygen atoms in the volume V ffi pApA,~Ay placed at the point (p, o~, y), with A p--0.05 /~, A ~ = 6 ° and Ayffi0.1A. In Fig. 3, contour maps are shown which repreTable I Mean residence time, f , of water within the first coordination shell of a water molecule
~-(ps)
Watermodel
Ref.
2.4 2.5 1.8 1.4 2.5-10.0
ST2 ST2 MCY TIPS3 exp.
present work 17 14 6 17
sent number density functions after integration with respect to the ~b coordinate. Arrows indicate the solute positions, and circles represent the spherical shells used to evaluate the water residence time. Each shell is subdivided into 'interior' (Fig. lc, dotted) and 'exterior' (Fig. lc, shaded) hemispherical regions, as previously done [4,18,19]. Radii are such that the shell borders lie prevalently where the number density function is minimal (lighter isodensity lines in Fig. 3). In Table 2, we report water mean residence times within the solutes first coordination shells. Approximately four and five water molecules occupy the interior and exterior regions, respectively, of each PHOB shell. The residence time for the former is 40% higher than for the latter, and 2.8 times larger than the free water value. In the PHIL-PHIL case, two water molecules reside within both interior and exterior regions around each solute. They correspond to molecules tetrahedrically coordinated to the PHIL solutes [l], which are represented in Fig. 4a by molecules 3 and 4 in the interior region, and by molecules 5 and 6 in the exterior region. The I" value for the former is 13 times larger than for free water. It is surprisingly
446
F. Brug~ et al./ Chemical Physics Letters 250 (1996) 443-449
8 6 4 2 0
' , ( b )
::
6 P 4 2 0
.......... i ............. i ...........
.............. ~.............. i ...........
i ...........
............
...........
.............
...............
.............
.............. ! ...........
: i p .............................................i............................................i.(c).. ..........i.............i............ i.............i.............~............ .L ................i...
6
4 2 0
'
....
.
........
..........
-14-12-10-8 -6 -4 -2 0 2 4 6 8 10 12 14 Y Fig. 3. Two-dimensional contour maps representing the number density distribution functions of the water oxygen atoms surrounding the PHOB-PHOB (a), PHIL-PHIL (b), and PHOB-PHIL (c) solute pairs, whose distance is 4.2 ,~. The functions have been integrated with respect to the ~ coordinate and plotted in the y-p plane. Arrows show the solutes positions. In (c) the PHOB and PHIL solutes are placed at - 2 . 1 ,~ and 2.1 /~, respectively. Circles represent first coordination shells, whose radii arc reported in Table 2.
close to the mean residence time (T = 36 ps) recently calculated tbr water surrounding polar atoms of the BPTI protein [20].
Two sections of Table 2 are reserved to the PHOB-PHIL system since, as expected [1], water behavior is different around PHIL and PHOB so-
Table 2 Mean residence time, I', and number of water molecules, no, around each solute pair shown in Fig. la. Radii of the spherical shells are indicated
Shell region
PHOB-PHOB
PHIL-PHIL
PHOB-PHIL a
PHOB-PHIL b
r = 4.75,~
r--4,~
r = 5.0,~
r = 4.0/~
no
• (ps)
no
• (ps)
no
• (ps)
no
• (ps)
interior
3.7
6.7
2.0
31.9
23.4
4.9
4.7
2.0
I 0.0
26.4 3.8 12.8 4.0
2.0
exterior
2.0 3.3 1.8 3.9
2.0
14.9
a PHOB site of the PHOB-PH1L pair. PHIL site of the PHOB-PH1L pair.
F. Brugb et aL / Chemical Physics l,otters 250 (1996) 443-449
447
b)
o) I
5
3
Fig. 4. Statistically relevant snapshots of water molecules surrounding the PHIL-PHIL (a) and the PHOB-PHIL (b) solute pairs.
Fig. 3c. The same peak occurs in the PHIL interior shell, for which the fitting yields n o = 2 and ~r= 23.4 ps. Therefore, in the space between the PHOB and PHIL two molecules are present with ~',-, 25. In the PHOB exterior shell two other molecules reside with ~-,,, 13 ps, and in the PHIL exterior shell two molecules with ~-,,, 15 ps. They correspond to the number density peak on the right of Fig. 3c. By comparing data for PHIL-PHIL and PHOBPHIL systems, one sees that the replacement of a PHIL with a PHOB, lute, which has already been shown to strongly affect S|Fs [2], causes a decrease
lutes. We also found that the water survival functions around the latter cannot be fitted to single exponontial functions, but pairs of exponential functions must be used. Moreover, the plane subdividing the PHOB shell into interior and exterior crosses the number density peak P (Fig. 3c). This causes an artificial decrease in the residence time value, as observed above, yielding 1---, 4 ps for n o = 3 molecules in the interior and for n o = 4 molecules in the exterior PHOB shell. In the PHOB interior shell two additional water molecules are statistically present with ~-= 26.4 ps, which correspond to the central peak of
P
8 6 4 2 0 6 4 2 0
,
,
()
I ~--~.~~ . ~ - - ~ . . ~ . - - . ~ ~,-~ _~ ~_,~.~--__-~
; . ~--~ ~-~ ~ - - - - ' - t " ~ ' ~ . , = .&-~..~...~ ~~,~,,~
;
i
!
:
i
-
i
( b )
........... ! ............. : .......... ! ............. ~.............. ~............ i ........ .............................................
:
~
_
~
~ ~. .e , ~7 ......~.-..... ~
~
-14-12-10 -8 -6 -4 -2
0
-I
~
'
~
:........
_
' ........
o-~
2 4 6 8 101214
Y Fig. 5. Two-dimensional contour maps representing the number density distribution functions of the water oxygen atoms surrounding a single PHOB (a), and a single PHIL (b) solutes which are placed at y = - 2 . I and y = 2.1, respectively.
448
F. Brug~ e, aL / Chemical Physics Letters 250 (1996) 443-449
in the water residence time to a value (¢'.-25 ps) practically coincident with the mean residence time around non-polar atoms of BPTI [20]. This might be fortuitous, since our model is oversimplified and the solute-solute distance of 4.2 ,g. has been arbitrarily chosen, though comparable with typical distances between protein surface groups exposed to the solvent [21]. However, our findings suggest that the water residence time around non-polar groups oa proteins could be remarkably affected by nearby polar groups, Rather surprisingly, the molecules residing within the PHIL exterior shells have a larger ¢ value in the PHOB-PHIL than in the PHIL-PHIL case. By inspection of the Fig. 4b, this result can be related with the enhancement of the HB network in the PHIL exterior side as a result of a corresponding weakening in the interior shell due to presence of the PHOB solute. This is also evident in the density maps (Fig. 3), where the central peak at y ~ 0 ~, which is symmetric in Fig: 3b, is shifted toward the PHIL position ( y -- 2.1 A) in Fig. 3c. To discriminate effects on the water residence time due to the immobilization of one solute from those induced by adding another solute, two further 500 ps MD simulations have been performed on a system of 728 water molecules and a single immobilized PHOB or PHIL solute placed in the same position and orientation as in the PHOB-PHIL pair (Fig. I a). The interior and exterior shells are conventionally defined in the same way. In Table 3 we report results for first coordination shells chosen considering the contour maps of Fig. 5.
Table 3 Mean residence time, ~'. and number of water molecules, n o, around a single PHOB and a single PHIL solute, which are kept immobile during 500 ps MD simulations in the same position and orientation as in the PHOB-PHIL pair sketched in Fig. la. Inxerior and exterior shells are defined as for the PHOB-PHIL pair, The r parameters indicale spherical shells radii Shell ragion
interior exterior
PHOB
PHIL
r -- 5,0 ~
r = 3.75 ~,
no
'r (ps)
no
~" (ps)
5.4 5, I
5.4
i .9 2,0
I 1,0 9.3
5,4
They show that, as expected, symmetry exists between interior and exterior regions of the PHOB shell, where -~ 5 molecules reside with T~ 5 ps. Around the immobilized water, instead, the HB donor side (interior region) affects more strongly the dynamics of the solvent water than the accepter side. Further, the water immobilization causes a four times increase in the water residence time as compared with the free water case. in qualitative agreement with previous tinddings [22]. Comparison of results reported in Tables 2 and 3 evidences that the addition of the PHOB solute at a distance of 4.2 ,~ from the PHIL solute doubles the residence time of the solvent molecules bound to the hydrogen atoms of the immobilized water, namely ~goes from 1 ! to 23.4 ps. Comparing Fig. 5b with Fig. 3c one sees that this is accompanied by a tighter localization of these water molecules, whose position is shifted towards the exterior shell, thus favoring their participation in the HB network, as observed above, and reducing the water mobility in that region. Instead, by adding a second PHIL solute to form the PHIL-PHIL solute pair, the ¢ value of the two moiccules (No. 3 and 4 in Fig. 4) present in the interior shells of both solutes increases almost three times, while in the exterior shells the water ~- changes by less than 10%. Work is in progress to investigate possible correlation between space-resolved properties of SIFs and solvent residence time.
4. Conclusions System configurations generated by MD simulations of aqueous solutions of simple solute models mimicking polar and non-polar groups attached to slowly diffusing macromolecules, have been analyzed to evaluate the water residence time within the first coordination shell of each solute. Results show that the mean residence time of the molecules nearest to the PHIL-PHIL solute pair is 32 ps, that is 13 times larger than the free water value (2.4 ps), and comparable with the water mean residence time (36 ps) around polar atoms of BPTI [20]. Although this agreement might be somewhat fortuitous, it seems significant that the same order of magnitude is obtained using such a simple model system.
F. Brugb et aL / Chemical Physics Letters 250 (1996~ 443-449
Analogously, the coincidence is remarkable of the residence time (25 ps) of the molecules closest to the PHOB-PHIL solute pair with the water mean residence time (25 ps) around non-polar atoms of BPTI [20]. Both these coincidences, if not accidental, would confirm that our modeling, though oversimplified, is able to capture essential features of the interaction between water and polar or non-polar groups belonging to macromolecules in solution, in analogy to the proven usefulness of using simple exact lattice models to explore the protein folding problem [23]. Further, the present results indicate that the solvent residence time depends remarkably upon the polarity characteristics not only of the single groups but also of the macromolecular environment to whom the groups belong. This could explain the rather large scattering of the residence time values reported for BPTI [20]. From a methodological point of view, the present findings suggest that the space region used to evaluate the solvent residence time should be carefully chosen in order to get reliable results. In particular, the border of this region should avoid any place where the solvent number density is large.
Acknowledgement We wish to thank Dr. A.E. Garcia for providing us with FORTRAN routines suitable to implement the residence time analysis, and G. Cottone, V. Martorana and R. Noto for useful discussions. This work has been partially supported by MURST (national and local funds) and by CRRN-SM.
449
References [1] F. Brugb, S.L. Fornili and M.B. Palma-Vittorelli, .L Chem. Phys. 101 (1994)2407. [2] F. Bmg~:, G. Cottone and S.L. Fornili, Chem. Phys. Letters 235 (1995) 402. [3] F.H. Stillinger and A. Rahman, J. Chem. Phys. 60 (1974) 1545. [4] C. Pangali, M. Rao and B.J. Berne, J. Chem. Phys. 71
0979) 2982. [5] [6] [7] [8] [9]
M. Mezei and A. Ben-Naim, J. Chem. Phys. 92 (1990) 1359. A.E. Garcia and L. Stiller, J. Comp. Chem. 14 (1993) 1396. O. Steinhauser, Mol. Phys. 45 (1982) 335. L. Verier, Phys. Rev. 159 (1967) 98. M.P. Allen and D.J. Tildesley, Computer simulation of molecular liquids (Clarendon Press, Oxford, 1987). [10] F. Brug~: and S.L. Fomili, Comput. Phys. Commun. 60 (1990) 31. [I I] F. Bmg~ and S.L. Fomili, Comput. Phys. Commun. 60 (1990) 39. [12] F. Brug~, J. Comput. Phys. 104 (1993) 263. [13] F. Sciortino and S.L. Fomili, J. Chem. Phys. 90 (1989) 2786. [14] R.W. Impey, P.A. Madden and I.R. McDonald, J. Phys. Chem. 87 (1983) 5071. [15] O. Matsuoka, E. Clemcnti and M. Yoshimine, J. Chem. Phys. 64 (1976) 1351. [16] W.L Jorgensen, J. Am. Chem. Soc. 103 (1981) 339. [17] P. Bopp, in: The physics and chemistry of aqueous ionic solutions, eds. M.C. Belissent-Funel and G.W. Neilson (Reidel, Dordrecht, 1987) p. 217. [18l A. Geiger, A. Rahman and F.H. S:illinger, J. Chem. Phys. 70 (1979) 263. [19] D.A. Zichi and P.J. Rossky, J. Chem. Phys. 83 (1985) 797. [20] R.M. Brunne, E. Liepinsh, G. Otting, K. W'fittich and W.F. van Gusteren, J. Mol. Biol. 231 (1993) 1040. [21] A. Ben-Naim, K.-L "ling and R.L. Jcmigan, Biopolymers 28 (1989) 1309. [22] R. Noto, M. Migliore, F. Sciortino and S.L. Fomili, Molec. Simul. I (198g)225. [23] K.A. Dill, S. Bromberg, K. Yue, K.M. Fiebig, D.P. Yee, P.D. Thomas and H.S. Chan, Protein Science 4 (1995) 561.