Fluctuating Hydration Structure around Nanometer-Size Hydrophobic

J. Phys. Chem. C 2007, 111, 2861-2871

2861

ARTICLES Fluctuating Hydration Structure around Nanometer-Size Hydrophobic Solutes II. Caging and Drying around Single-Wall Carbon Nanotubes Takeshi Hotta† and Masaki Sasai*,‡,§ Department of Complex Systems Science and Department of Computational Science and Engineering, Nagoya UniVersity, Nagoya 464-8603, Japan, and CREST-JST, Nagoya 464-8603, Japan ReceiVed: May 16, 2006; In Final Form: October 10, 2006

The nanometer-scale hydration and hydrophobic interactions are studied by simulating aqueous solutions of single-wall carbon nanotubes (SWCNs) as example systems. Molecular dynamics simulation of the solutions reveals that hydration around the nanometer-size solute is strongly dependent on the topology of the solute configuration: When a single SWCN is solved in water, water molecules are attracted to the surface of a SWCN due to the dispersion interactions between carbons and water molecules. Water molecules in the first hydration layer around the SWCN move more slowly than in bulk. When a pair of SWCNs approach each other, the water density in the interstitial region between SWCNs fluctuates largely from trajectory to trajectory. The water-mediated force between SWCNs also fluctuates largely and becomes attractive when the water density in the interstitial region is small in fluctuation. When a triplet of SWCNs come close to each other, water molecules are expelled from the interstitial region to leave a vacant space, which drives SWCNs to aggregate. Correlation between the fluctuating caging and the fluctuating drying suggests that the topologydependent drying takes place when the water configuration becomes too much constrained to keep a consistent hydrogen bond network in the narrow region surrounded by the nanometer-size hydrophobic solutes.

1. Introduction How do hydrophobic solutes interact with each other in aqueous solution? The classical view of hydrophobic interaction has been based on the picture that water molecules form a clathrate-like configuration or “caging” around solutes and the gain of entropy due to the overlapping of cages may be an origin of the attractive interaction between them.1 Theoretical analyses showed that small apolar molecules such as rare gas atoms or butane can indeed fit into the hydrogen bond (HB) network and the picture of caging hydration has been supported.2-4 Recently, however, an alternative picture has emerged on the basis of the idea that water molecules are repelled from the large enough hydrophobic surfaces and attractive interaction between surfaces is induced to compensate the free energy cost arising from density depletion.5-9 This idea of drying hydration dates back to the insight given by Stillinger5 and was revived recently10 by the experimental observation of nanobubbles11-13 and density depletion on hydrophobic surfaces.14 Lum, Chandler, and Weeks15 developed a mean-field theory of hydration and suggested that the crossover from the caging hydration around small apolar solutes to the drying hydration around the large hydrophobic surfaces takes place for solutes of the nanometer size, which is just the size relevant to the biomolecular structural formation. Computer simulations of the nanometer-size systems have †

Department of Complex Systems Science, Nagoya University. Department of Computational Science and Engineering, Nagoya University. § CREST-JST. ‡

been performed16-20 to test the conjecture of Lum, Chandler, and Weeks. It has been shown that hydration sensitively depends on the interaction between solute atoms and water molecules. When the attractive interaction between solute atoms and oxygen in water molecules are artificially turned off in simulation, the clear drying phenomena have been confirmed around the nanometer-scale hydrophobic surfaces.16-22 When there are realistic strength of attractive dispersion interactions between solutes and water molecules, however, drying is much suppressed to be obscured.19-21 Huang et al. have examined the molecular dynamics (MD) simulation of paraffin plates and shown that drying does not take place when the paraffin and water have attractive interactions that are as strong as those between graphite and water molecules.20 On the other hand, simulations of biomolecules have been performed to show that drying can actually take place in the system of melittin tetramer.23 It is interesting to see that the hydrophobic surface of this molecule has some polar residues and the realistic dispersion interactions between hydrophobic residues and water molecules do not prevent the peptide surface from drying. Comparing with the absence of drying in a simulated twodomain protein system,21,22 drying around the simulated melittin tetramer suggested that hydration is very sensitive not only to interactions between atoms but also to topology of the solute configuration. In a previous paper we have studied hydrophobic hydration around the nanometer-size hydrophobic spheres, C60 and C60H60, with MD simulation by adopting realistic strength of dispersion interactions between carbons in solutes and oxygen in water

10.1021/jp062977h CCC: $37.00 © 2007 American Chemical Society Published on Web 01/30/2007

2862 J. Phys. Chem. C, Vol. 111, No. 7, 2007 molecules and shown that there is no definite drying transition around those spherical solutes.24 There are, however, both fluctuating drying and caging hydration around the solutes and the strength of drying or caging is largely different from one trajectory to the other in simulation, showing that the free energy barrier between different hydration structures is so large that each trajectory of 100 ps length or so does not surpass the barrier. It has been shown that when caging is more developed in the fluctuation, spheres in contact are more stabilized and when the hydration layer is more dried, it becomes easier for two spheres to come in contact. In this paper we examine the nanometer-size hydrophobic solute having a different topology; a single-walled carbon nanotube (SWCN), which has a cylindrical conformation. In previous works, hydration25,26 and the drying transition26 around a zigzag (16,0) SWCN have been studied by MD simulations. When we define the diameter of an SWCN as the diameter of the circle connecting the center of carbon atoms, the diameter of a zigzag (16,0) SWCN is 12.508 Å. Instead of using this thick SWCN, we here use an armchair (5,5) SWCN27 whose diameter is 6.932 Å. We focus on this thinner SWCN to compare the results with those obtained for the spherical solutes of similar diameter.24 A realistic strength of potential is used to represent the attractive dispersion interaction between carbons in solutes and oxygen in water molecules. We show that caging and drying coexist also in this system and fluctuating caging and drying decisively affect the solvent mediated interaction between solutes. We also show that the definite drying transition takes place when SWCNs take a certain configuration. Thus, we can clearly see that the balance between caging and drying in the nanometer length scale is strongly dependent on topology of the solute configuration. 2. Methods 2.1. Molecular Dynamics Simulation. MD simulations were performed by using the AMBER7 simulation package28 under an NVT ensemble with temperature set to be 298 K. Water was described by the TIP4P model29 and the Coulomb interaction was treated by the particle mesh Ewald method. Carbon atoms of the solute were modeled as uncharged particles interacting with the Lennard-Jones (LJ) potential with C-C ) 0.0556 kcal mol-1 (or 0.233 kJ mol-1) and σC-C ) 3.400 Å.30 The interaction between a carbon atom in a SWCN and a water molecule was modeled by the LJ potential by employing the Lorentz-Berthelot mixing rules;31 C-O ) (C-C·O-O)1/2 and σC-O ) (σC-C + σO-O)/2, where O-O and σO-O are the LJ parameters for oxygen atoms of TIP4P water, which yields C-O ) 0.0927 kcal mol-1 (or 0.388 kJ mol-1) and σC-O ) 3.277 Å. A variety of values have been adopted for the carbon-water LJ parameters.32-35 For example, C-O ) 0.114 kcal mol-1 and σC-O ) 3.275 Å,32 C-O ) 0.0930 kcal mol-1 and σC-O ) 3.280 Å,33 and C-O ) 0.0937 kcal mol-1 and σC-O ) 3.190 Å34 have been used previously. In ref 35 C-O ) 0.0873 kcal mol-1 (0.3651kJ mol-1) and σC-O ) 3.190 Å have been recommended to reproduce the appropriate contact angle of water droplets on graphite. Parameters used in this paper, C-O ) 0.0927 kcal mol-1 and σC-O ) 3.277 Å, are within the range of this variety. We performed simulations of three systems: system I was composed of a SWCN and 898 water molecules, system II was composed of a pair of SWCNs and 795 water molecules, and system III was composed of a triplet of SWCNs and 697 water molecules. Each of three systems was confined in a cubic box with a side length of Lbox ) 31.064 Å and the periodic boundary condition was applied. Densities around the edge of the box

Hotta and Sasai averaged over trajectories were F ) 1.000 g/cm3 (system I), F ) 0.995 g/cm3 (system II), and F ) 0.998 g/cm3 (system III). We treated SWCN as a rigid body by fixing the position and orientation. As a reference system, we also calculated the pure water system composed of 1000 water molecules confined in the same periodic boundary box with the side length Lbox ) 31.064 Å, whose density was F ) 0.997 g/cm3. For system I, an armchair (5,5) SWCN was placed as a solute at the center of the box of the periodic boundary condition with the straight side of the SWCN perpendicular to a boundary wall of the box and 898 water molecules were placed around the SWCN. For system II, two paralleling armchair (5,5) SWCNs were placed symmetrically around the center of the box, where both SWCNs were oriented perpendicular to a boundary wall of the box, and 795 water molecules were placed around the SWCNs. The distance between the center line of one SWCN and the center line of another SWCN was fixed to be rc. For system III, three armchair (5,5) SWCNs were placed in isosceles triangle geometries. Three cylindrical solutes were aligned parallel to one another and perpendicular to a plane of the box. The distance between center lines of the first and the second SWCN was fixed to 10.20 Å as a dimer of SWCNs and both the center-line distance from the third SWCN to the first one and the center-line distance from the third to the second one were set to be rc.A total of 697 water molecules were placed around the SWCNs. After a short energy minimization for the initial configuration of each system, a high-temperature run at 498 K was performed for 1 ns duration for system I and for 2 ns duration for systems II and III to break the HB network in the initial configuration. We extracted a configuration at every 100 ps of these high-temperature runs. For system I or II, each of the thus obtained 10 or 20 configurations was equilibrated at 298 K for 100 ps duration, and for system III for 500 ps duration. The subsequent 100 ps trajectories were sampled at 298 K by saving snapshots at every 10 fs, which provides 10 000 snapshots for each trajectory. In all simulations a time step was 1 fs. To check whether the results depend on the choice of ensemble, we also performed simulations under an NPT ensemble. Starting from 10 structures obtained from the NVT ensemble, we calculated 10 equilibration runs for 500 ps under the NPT ensemble with P ) 1 atm and T ) 298 K. Subsequent to this equilibration, 100 ps length trajectories were calculated under the same NPT ensemble to sample the data. In the NPT ensemble Lbox fluctuates to keep the pressure constant. In systems II and III we allowed rc to fluctuate proportionally to Lbox around its average value of 〈rc〉. We further performed two additional types of simulations to check whether the results depend on the choice of potentials. One type of simulation was performed with a polarizable model of water molecules, POL3,36 under the NVT ensemble, and the carbon-water LJ parameters were adjusted by employing the Lorentz-Berthelot mixing rules for the POL3 parameters as C-O ) 0.0932 kcal mol-1 and σC-O ) 3.30 Å. The other type of simulation was performed by using the POL3 potential under the NVT ensemble, and parameters of ref 35, C-O ) 0.0873 kcal mol-1 and σC-O ) 3.190 Å, were used for the LJ interaction. For those two types of simulations the same equilibration and sampling procedures were carried out as described above for the case of the TIP4P model. In the following we mainly show the results of simulations calculated under the NVT ensemble with TIP4P potential, C-O ) 0.0927 kcal mol-1 and σC-O ) 3.277 Å unless specifically mentioned.

Fluctuating Hydration Structure in SWCNs

J. Phys. Chem. C, Vol. 111, No. 7, 2007 2863

Figure 1. Radial distribution function of oxygen atoms of water molecules around a single SWCN. The radius, r, is measured from the center of the SWCN in the radial direction. The distribution function is the average over 10 trajectories.

Figure 2. Radial distribution function of amplitude of site dipoles around a single SWCN. The line connects the values averaged over 10 trajectories. Error bars are drawn between maximum and minimum values in 10 trajectories.

2.2. Site Dipoles and Site Density. Individual water molecules move fast and lose their orientational or positional memories in several picoseconds at room temperature. These molecular motions, however, are not completely random but they fit the structure and movement of the surrounding HB network. Because such a surrounding network configuration moves more slowly than individual molecules, physical quantities averaged at each spatial position, or “site” in the simulation box, should have some long-time components, which are not expected in quantities obtained by tracing individual molecules.37 Such “site”-dependent quantities represent the tendency of water molecules passing by the site and are convenient to visualize the hydration structure.24,37 The periodic box was divided into 20 × 20 × 20 cubes, so that the side length of a cube was Lcube ) Lbox/20 ) 1.5532 Å. The position ri of the ith cube was represented by its body center. We defined the site dipole as dw(ri) ) Σtd(ri,t)/Ndetect(ri), where Ndetect(ri) was the number of snapshots when the oxygen atoms of water molecules were detected in the ith cube during the time duration w, and d(ri,t) was a unit vector pointing from the oxygen atom to the midpoint of two hydrogen atoms of the water molecule detected in the ith cube at time t. When no molecule was detected in the ith cube, d(ri,t) ) 0. In this paper, we sampled snapshots of the 10 fs interval in each trajectory during w ) 100 ps. We also introduced the site density, Fw(ri). The periodic box was divided into 60 × 60 × 60 cubes, so that the side length was L′cube ) Lbox/60 ) 0.5177 Å. The site density was defined by Fw(ri) ) Ndetect(ri)/(NsampleV′cubeF0) ) 0.0216Ndetect(ri), where Nsample ) 104 was the number of snapshots sampled during w ) 100 ps, and V′cube ) L′cube3. Ndetect(ri) was defined in the same way as explained above, and F0 ) 1000/Lbox3 was the averaged number density of water molecules in the reference pure water system. In the NPT ensemble Lbox fluctuates as Lbox(t). We defined Lcube(t) ) Lbox(t)/20 or L′cube(t) ) Lbox(t)/60 to calculate site dipole or site density.

The average density of water molecules in the cylindrical layer of 5.7 Å < r < 11.6 Å, Fhyd, is larger than the average water density of this system, Fav, with Fhyd/Fav ) 1.021. This tendency shows that water molecules are attracted to the solute in average and there is no evidence of drying around a single SWCN. We should note, however, that the water structure around two or three SWCNs is different from that around a single SWCN: As will be shown in section 3.2, the drying occurs as a trajectory-dependent fluctuation in a pinched region between a pair of SWCNs, and the dried cavity clearly appears in a region among a triplet of SWCN. The peaks and valleys found in Figure 1 suggest that hydrating water molecules form some structures around the solute, which can be more clearly seen with the site-dipole method. In Figure 2, the radial distribution of the amplitude of 〈Σi|dw(ri)|〉/(2πr∆rLbox) is shown together with the range of the trajectory-dependent fluctuation of Σi|dw(ri)|/(2πr∆rLbox), where summation is taken over sites in the layer between the cylinder with radius r and the one with radius r + ∆r, with ∆r ) 0.8Å, and 〈...〉 represents the average over 10 trajectories. Note that the length of each cube, Lcube, is 1.5532 Å, so that the point at 4.8 Å in Figure 2, for example, summarizes information of the region of approximately 4.8 Å - Lcube/2 < r < 5.6 Å + Lcube/2 ) 6.3766 Å and hence has a nonzero value though water molecules cannot penetrate into the region r < 5.7 Å. Figure 2 shows that water molecules passing by the solute have a larger tendency to be aligned than in the bulk, and the tendency is strong when molecules are in the region of 5.7 Å < r < 8.0 Å, which corresponds to the first hydration layer of SWCN, as shown in Figure 1. Figure 3 shows a typical pattern of site dipoles around a SWCN in an exemplified trajectory. We can see that water molecules tend to align parallel to the straight side of a SWCN in the first hydration layer. A distribution of the dipole orientation is shown in Figure 4, confirming the intuitive picture obtained from Figure 3: Dipoles of water molecules have a mild tendency to align parallel to the straight side of the SWCN. This tendency is consistent with the previously reported results that water molecules near SWCN tend to direct parallel to the cylindrical surface,25 but here we found that the dipole direction is more restricted, showing the tendency to be aligned to the straight side direction. This coherence in dipole orientation is shown in Figure 5 by drawing the site dipoles in the first hydration layer of SWCN in the “surface map” of the cylinder. Water molecules form a “vortex” like pattern on the surface with some bias to align parallel to the straight side of the cylinder. The large amplitude of site dipoles near the solute dose not imply that each water molecule is frozen on the hydrophobic

3. Results and Discussion 3.1. Fluctuating Caging around a SWCN. The solute configuration has a two-dimensional symmetry, so that the radial distribution function is calculated by sampling the number of atoms in the cylindrical layer of radius r from the center line of a SWCN. Figure 1 is such a radial distribution function of oxygen atoms of water molecules around a single SWCN in system I. Figure 1 shows that water molecules cannot penetrate into the cylinder of radius r < 5.7 Å from the center of SWCN in the radial direction. Peaks are at r ) 6.6 and 9.3 Å, which correspond to the first and second hydration layers, respectively, giving a density profile similar to the one reported in ref 25.

2864 J. Phys. Chem. C, Vol. 111, No. 7, 2007

Figure 3. Typical pattern of site dipoles, dw(ri), around a single SWCN in an example trajectory. Site dipoles with lengths larger than 55% of the largest value in the system are shown by bars. Each of site dipoles is oriented from red to white.

Figure 4. Distribution of the cosine of the angle between site dipoles and the vector parallel to the straight side of the cylindrical solute. Distribution is averaged over 10 trajectories. Values are normalized for their maximum to be 1.

Figure 5. Pattern of site dipoles, dw(ri), around a SWCN molecule in an example trajectory. Site dipoles in the cylinder of r < 8 Å from the center of SWCN are shown in the plane of the y coordinate and θ, where the coordinate along the straight side of the cylinder is denoted by “y coordinate” and the angle on the circular surface of the cylinder is denoted by θ. Site dipoles with lengths larger than 55% of the largest one in the system are shown.

surface with a fixed direction but implies that different molecules came into and went out of the site multiple times with a tendency to be aligned. Restless motions of water molecules can be confirmed by plotting the diffusive motions in Figure 6. Here, we should note that the estimated diffusion constant of water

Hotta and Sasai

Figure 6. Diffusive and rotational motions of water molecules found near the solute at the initial time. Averaged over water molecules which reside at time 0 ps within the cylinder of radius 8 Å from the center of SWCN in an example trajectory. (a) Square distance between the position of each water molecule at time 0 ps and that at time t in a trajectory. The solid line is fitted to the data of first 13 ps ,and the dotted line is fitted to the data later than 13 ps. (b) Relaxation of the cosine of the angle between the dipole of water molecule at time 0 ps and that at time t in a trajectory. Data are fitted by an exponential curve with the relaxation time τhyd ) 5.59 ps. To estimate the motion of water molecules near SWCN correctly, SWCN is not treated as a rigid body as in other calculations in this paper but is here treated as a network of carbon atoms connected by springs with the spring constant of 114.46 kcal/(Å2 mol).

molecules fluctuates from one trajectory to the other: In an example trajectory of Figure 6a the diffusion constant of watermolecules which were within 8.0 Å from the center of a SWCN at time 0 ps is Dhyd ) 1.60 Å2/ps during the first 13 ps and Desc ) 2.17 Å2/ps for the period later than 13 ps. The latter value of Desc is close to that of water molecules found in the separated region of further than 8.0 Å from the center of a SWCN at time t ) 0 ps in the same system, Dsep ) 1.95 Å2/ps or that in pure water system, Dpure ) 2.11 Å2/ps. After averaging over 10 trajectories, we obtain Dhyd ) 1.70 Å2/ps and Desc ) 2.05 Å2/ps and Dsep ) 1.96 Å2/ps. Diffusion of water molecules in the first hydration layer is slower than that in pure water as Dhyd < Dpure as was noted in refs 18, 22, and 38. Average hydrated water molecules diffuse about 14 Å, over 9 cubes, during the period to evaluate a site dipole in a single cube, w ) 100 ps. A similar alignment of orientation of rapidly moving water molecules was found around the simulated C60 and C60H60 by calculating the site-dipole distribution.24 Alignment of water molecules shown in Figures 3-5 is not the rigid fixed one but represents the statistical tendency of moving molecules and hence we should refer it to as the “fluctuating caging” around SWCN. The averaged relaxation of cosine of the angle between the dipole of a molecule within 8.0 Å in the radial direction from the center of a SWCN at time 0 ps and that of the same molecule at time t is plotted in Figure 6b. By fitting an exponential curve to the data and averaging it over 10 trajectories, we find that the relaxation time is τhyd ) 4.96 ps, which is larger than that for water molecule separated from the solute, τsep ) 3.54 ps, in



Figure 7. Potential of mean force acting between a pair of SWCNs or among a triplet of SWCNs: V(rc) calculated by eq 1 in the text (solid line); VLJ(rc) derived by summing the Lennard-Jones interactions between atoms in solutes (dotted line); Vhyd(rc) that is the solvent contribution to the potential of mean force (dashed line). (a) Potentials between a pair of SWCNs and (b) potentials among a triplet of SWCNs, separated by the distance rc.

the same system, or τpure ) 3.47 ps in the reference system of pure water. This slower rotational relaxation should be an origin of the larger amplitude of site dipoles near the solute, but it alone does not explain the coherent pattern of site dipoles: If randomly oriented water molecules would visit the site with a slower rotational relaxation than bulk water molecules, the site should have a larger site dipole than in the bulk, but the resultant pattern of site dipoles should not be spatially correlated, as shown in Figure 5. Coherent patterns of site dipoles lasting for w . τhyd imply that multiple water molecules visiting the observed site near the surface tend to accommodate their directions to the surrounding structure of the HB network, which brings about the spatial correlation among site dipoles. In this way the fluctuating cage, which was expressed as the hydration structure around C60 or C60H60 in our previous study,24 also exists around a SWCN. In summary, diffusion of water molecules in the first hydration layer of SWCN is slower than in pure water. Water molecules near a SWCN tend to align to form patterns of dipoles around the hydrophobic cylinder. This is not the rigidly frozen iceberg-like cage but is a fluctuating cage, which can be seen via statistical sampling of directions of water molecules during the period much longer than the rotational relaxation time. The drying phenomenon, however, is not evident around a single SWCN. 3.2. Fluctuating Drying and Caging around a Pair of SWCNs and a Triplet of SWCNs. Figure 7 is the potential of mean force (PMF), V(rc), between a pair of solutes or among a triplet of solutes, which was calculated by integrating the force acting on solutes whose positions and orientations were fixed during the simulation, n

V(rc) ) -

〈F(ri)‚u〉δr ∑ i)0

(1)

Here, ri is distance between the center of one solute and that of the other solute; ri ) ro - iδr, with δr ) 0.2 Å and r0 ) 15 Å, and rc ) r0 - (n + 1/2)δr. For a pair of SWCNs F(ri) was obtained from F(ri) ) (F1(ri) - F2(ri))/2, where F1(ri) was the total force acting on one SWCN and F2(ri) was that acting on the other SWCN. u is a unit vector pointing from the center of one solute to that of the other solute, and 〈...〉 is average over 20 trajectories. V(rc) has a minimum at distance rcontact ≈ 10.2 Å at which two SWCNs are in contact and the additional minimum at rseparated ≈ 13.5 Å at which two SWCNs are separated by a layer of one interstitial water molecule. For a triplet of SWCNs the pair of the first and second SWCNs were treated as a dimer pair with the fixed distance between their centers to be 10.2 Å. The third SWCN was placed on the top of the isosceles triangle. rc was defined by the distance between the center of the third SWCN and the center of the first SWCN or equivalently between the center of the third SWCN and the center of the second SWCN. Then, F(ri) ) (F1(ri) - F2(ri))/2 was calculated from the force acting on the dimer pair, F1(ri), and the force acting on the third SWCN, F2(ri). V(rc) has a minimum at distance rcontact ≈ 10.2 Å at which three SWCNs are in contact and a very shallow minimum at rseparated ≈ 13.3 Å. The mean force, F(rc)‚u, consists of the contribution from the direct Lennard-Jones interactions between atoms in solute molecules, FLJ(rc), and the contribution from forces acting between water molecules and solutes, Fhyd(rc); F(rc)‚u ) FLJn (rc) + Fhyd(rc). In Figure 7, VLJ(rc) ) -Σi)0 FLJ (ri)δr and Vhyd n (rc) ) -Σi)0 Fhyd (ri)δr are plotted together with V(rc). We can find that Vhyd(rc) is small at rc ≈ rcontact and V(rcontact) ≈ VLJ(rcontact), so that the interaction between solutes in contact is dominated by the dispersion interaction both for a pair of SWCNs and for a triplet of SWCNs. For a pair of SWCNs, Vhyd(rc) and Fhyd(rc) are significant in the range rcontact < rc < rseparated, yielding a free-energy barrier separating rcontact and rseparated, as shown in Figure 7. For a triplet of SWCNs, on the other hand, absolute values of Vhyd(rc) and Fhyd(rc) are small in the rage of rcontact < rc < rseparated, leading to an extremely low barrier separating rseparated and rcontact; The height of the barrier is less than the thermal energy, which can be easily surmounted by the thermal fluctuation. In Figure 8, F(rc)‚u and FLJ(rc) are plotted, showing that for a pair of SWCNs, Fhyd(rc) at around the free-energy barrier has a large trajectory-to-trajectory fluctuation, which implies that the chance for the separated pair to be in contact or the probability for the pair in contact to be separated is strongly affected by the trajectory-dependent fluctuation. For a triplet of SWCNs, Fhyd(rc) also has a large trajectory-to-trajectory fluctuation, implying that the force among SWCNs at rc ≈ rseparated can be either attractive or repulsive. Insight into these fluctuations should be obtained by analyzing the correlation among Fhyd(rc), site densities, and site dipoles around the solutes. A rectangular box with the size rrec ) 8.0 Å is defined as in Figure 9 to examine the interstitial region pinched by SWCNs. As shown in Figure 10, Fint, the site density averaged over the region within this box, has a large trajectory-dependent fluctuation. We can note that the trajectory-dependent drying, or “fluctuating drying”, takes place in the region pinched by SWCNs. The large fluctuation of Fint implies that the free-energy barrier between different water configurations is rather large and is not often surmounted by the system during a single trajectory of 100 ps both for a pair of SWCNs and for a triplet of SWCNs. Fluctuation of Fint is strongly anticorrelated to the force between solutes as shown in Figure 10. The correlation


Hotta and Sasai

Figure 8. Trajectory-dependent fluctuations of force, F(rc)‚u, acting between a pair of SWCNs or among a triplet of SWCNs. Error bars are drawn between the minimum and maximum forces of 10 trajectories. FLJ(rc) derived from sum of the Lennard-Jones interactions between atoms in solutes is shown with a dotted line. The sign of forces is defined to be positive when attractive and negative when repulsive. (a) Forces between a pair of SWCNs and (b) forces among a triplet of SWCNs.

Figure 9. Illustration of configurations of a pair of SWCNs and a triplet of SWCNs separated by the distance rc. The interstitial region between the pair or among the triplet is analyzed by defining a rectangular box of size rrec.

coefficient is -0.986 in Figure 10a and is -0.986 in Figure 10b. The same anticorrelation was observed when we performed the simulation under the NPT ensemble: The correlation coefficient is -0.818 when 〈rc〉 ) 12.52 Å in system II, and -0.997 when 〈rc〉 ) 12.94 Å in system III. The results do not depend much on the choice of the water model either: In the POL3 model with C-O ) 0.0932 kcal mol-1 and σC-O ) 3.30 Å, the correlation function is -0.895 when rc ) 12.60 Å in system II and -0.999 when rc ) 13.04 Å in system III. Values of force and site density, however, depend on the choice of parameters in the LJ potential. In the POL3 model with C-O ) 0.0873 kcal mol-1 and σC-O ) 3.190 Å, the force

Figure 10. Correlation between site densities averaged within the rectangular box defined in Figure 9 and the mean force F(rc)‚u. Both site densities and F(rc)‚u show large trajectory-dependent fluctuation. (a) Site densities and the mean force between a pair of SWCNs separated by the distance rc. Calculated with the TIP4P model under the NVT ensemble at rc ) 12.60 Å (red), the TIP4P model under the NPT ensemble at 〈rc〉 ) 12.52 Å (green), the POL3 model under the NVT ensemble with C-O ) 0.0932 kcal mol-1 and σC-O ) 3.30 Å at rc ) 12.60 Å (blue), and the POL3 model under the NVT ensemble with C-O ) 0.0873 kcal mol-1 and σC-O ) 3.190 Å at rc ) 12.60 Å (pink). (b) Site densities and the mean force among a triplet of SWCNs separated by the distance rc. Calculated with the TIP4P model under the NVT ensemble at rc ) 13.04 Å (red), the TIP4P model under the NPT ensemble at 〈rc〉 ) 12.94 Å (green), the POL3 model under the NVT ensemble with C-O ) 0.0932 kcal mol-1 and σC-O ) 3.30 Å at rc ) 13.04 Å (blue), and the POL3 model under the NVT ensemble with C-O ) 0.0873 kcal mol-1 and σC-O ) 3.190 Å at rc ) 13.04 Å (pink) and at rc ) 12.85 Å (light blue). F(rc)‚u is defined to be positive when attractive. The negative correlation shown here implies that the repulsive force is weakened or turned into attractive when the interstitial region is dried.

among a triplet of SWCNs at rc ) 13.04 Å and 12.85 Å is more attractive than in other calculations with the larger σC-O as shown in Figure 10b. With this small σC-O parametrization the correlation coefficient is -0.977 when rc ) 12.60 Å in system II, -0.941 when rc ) 13.04 Å in system III, and -0.980 when rc ) 12.85 Å in system III. Thus, although the precise values of density and force depend on the parameters of the LJ potential, the features of the large trajectory-dependent fluctuation in Fint and the strong anticorrelation between the force and Fint can be also found in simulations with different LJ parameters. Analyses of Figure 10 are further extended by varying rc by using the TIP4P model under the NVT ensemble. In Figure 11a the correlation coefficient between Fint and F(rc)‚u is plotted


Figure 11. Correlation coefficients between site densities and the mean force and between site dipoles and the mean force. (a) Correlation coefficients between site densities averaged within the rectangular box and the mean force. (b) Correlation coefficients between the relative strength of site dipoles in the cylinder defined in Figure 12 and the mean force. Coefficients for the system of a pair of SWCNs (dashed line) and coefficients for the system of a triplet of SWCNs (solid line).

for various values of rc. We can find the strong anticorrelation for the range 12.6 Å < rc < 12.8 Å in the case of a pair of SWCNs, and for the range 12.8 Å < rc < 13.4 Å in the case of a triplet of SWCNs. For a pair of SWCNs the maximum of PMF strongly diminishes when the interstitial region becomes dried, which allows the two solutes to approach each other across the PMF maximum. For a triplet of SWCNs the PMF barrier is low in average and is easily surmounted. When the interstitial region is dried in fluctuation, this tendency to overcome the barrier is further enhanced. Thus, the fluctuating drying in the interstitial region decisively determines the stability and lifetime of the separated solutes in both systems. To see the development of the site-dipole ordering in the interstitial region, we define a cylindrical region as illustrated in Figure 12 and define the relative weight of strength of site dipoles in the cylinder, dcyl ) ∑i∈cylinder |dw(ri)|/∑i∈outside |dw(ri)|, where ∑i∈cylinder is sum over sites in the interstitial cylinder of radius rcyl ) 1.5532 Å and ∑i∈outside is sum over sites which are further than 8.0 Å in the radial direction from any center of SWCN. Figure 13 shows that the relative strength of site dipole, dcyl, is enhanced in the range of rc for which Fint begins to decrease as solutes approach each other. The correlation coefficient between F(rc)‚u and dcyl is shown in Figure 11b. The attractive force between solutes and the strength of site dipole are positively correlated in the range of rc for which dcyl is large. Distributions of site dipoles are visualized in Figure 14. In Figure 14a,b we can see that site dipoles are enhanced in the interstitial region between a pair of SWCNs or among a triplet of SWCNs, as was expected from Figure 13. The repulsive force between or among SWCNs is reduced or becomes attractive when the site dipole is enhanced and the site density is reduced in these interstitial regions. It is interesting to see in Figure 14c that site dipoles are diminished in the interstitial region of system III when SWCNs come close to rc ) 12.85 Å. This is due to the drying in system


Figure 12. Illustration of configurations of a pair of SWCNs and a triplet of SWCNs separated by the distance rc. The interstitial region between the pair or among the triplet is analyzed by defining a cylinder of radius rcyl.

Figure 13. Dependence of the averaged site densities and the relative strength of site dipoles on the distance between SWCNs. Averaged over 10 trajectories. (a) Site densities averaged within the rectangular box of size rrec = 8.0 Å defined in Figure 9. (b) Relative strength of site dipoles in the cylinder of radius rcyl = 1.5532 Å defined in Figure 12. A pair of SWCNs (dashed line) and a triplet of SWCNs (solid line).

III: Site dipoles are not detected because water molecules are expelled from the interstitial region. Drying in system III is clearly seen in Figure 15 of the solvent configuration. In system II (Figure 15a) there are water molecules filling the interstitial region with rc ) 12.80 Å and drying is observed only infrequently in the trajectory-dependent fluctuation but in system III (Figure 15b) we can find water molecules expelled from the interstitial region in all the trajectories examined. Such a difference between system II and system III was also confirmed


Hotta and Sasai

Figure 15. Snapshots of MD simulation: (a) pair of SWCNs with the distance rc ) 12.80 Å; (b) triplet of SWCNs with the distance, rc ) 12.85 Å.

Figure 14. Typical patterns of site dipoles, dw(ri), in examples of trajectories. Patterns calculated (a) around a pair of SWCNs separated by the distance rc ) 12.60 Å, (b) around a triplet of SWCNs with rc ) 13.04 Å, and (c) around a triplet of SWCNs with rc ) 12.85 Å. Site dipoles with the length larger than 55% of the largest value in the system are shown by bars. Each of site dipoles is oriented from red to white.

under the NPT ensemble: A clear drying takes place in system III with 〈rc〉 ) 12.81 Å whereas water molecules reside between a pair of SWCNs in system II with 〈rc〉 ) 12.71 Å. The results do not depend on the water model either: For simulations in the POL3 model with C-O ) 0.0932 kcal mol-1 and σC-O ) 3.30 Å, a clear dying transition takes place in system III with rc ) 12.85 Å, whereas water molecules reside between a pair of SWCNs in system II with rc ) 12.80 Å. In the POL3 model with C-O ) 0.0873 kcal mol-1 and σC-O ) 3.190 Å, on the other hand, drying was not evident at rc)12.85 Å in system III. With this small σC-O the tendency of drying appears at smaller rc than in other calculations. With this parametrization a clear drying transition takes place in system III with rc ) 12.67 Å, whereas a clear drying transition is absent in system II and water molecules reside between a pair of SWCNs with rc ) 12.60 Å in system II. Thus, the probability that the narrow interstitial region is dried is strongly dependent on the topology of the solute configuration. This tendency of drying can be found in the rc dependence of Fint shown in Figure 13a: In system II Fint gradually becomes smaller as rc decreases and Fint ≈ 0 only when rc ≈ 11 Å, for which it is geometrically impossible to insert a water molecule in the interstitial region, whereas in system III Fint ≈ 0 at rc ≈ 12.5 Å, for which a water molecule can enter into the interstitial



Figure 16. Schematic phase diagrams of water pinched by solutes. In the phase of fluctuating caging, the relative strength of the site dipole is large, as shown in Figure 13b, and in the dried phase, the average site density is small, as shown in Figure 13a. In the fluctuating drying phase, the site density and F(rc)‚u are strongly anticorrelated, as shown in Figure 11a. rseparated is the distance at the PMF minimum at which SWCNs are separated by a layer of water molecules, and rpeak is the distance at the top of the PMF barrier.

region geometrically but the entrance was prohibited in the simulated results due to its free energy cost. 3.3. Drying Transition around a Triplet of SWCNs. Differences in hydration in system II and system III are summarized in Figure 16. When rc is large enough, drying is not evident in either system II or system III and the hydration is dominated by the fluctuating caging as in system I. In system II a pair of SWCNs are stabilized by the fluctuating caging at rc ≈ rseparate. When a pair of SWCNs approach each other, the fluctuating drying begins to play important roles to determine the interaction. The free energy barrier that stabilizes the separated pair of SWCNs becomes lower when water molecules are expelled from the interstitial region in the trajectorydependent fluctuation. In system III, on the other hand, the fluctuating drying takes place at around rc ≈ rseparate and the separated SWCNs are not very stable. This predominance of drying in system III is consistent with the low free energy barrier in PMF. When SWCNs further approach in system III, the definite drying transition takes place. Analyses of the site-dipole distribution in Figures 12-14 showed that the fluctuating caging is enhanced in the interstitial region for which the fluctuating drying takes place. A consistent explanation for the simultaneous enhancement of caging and drying in the interstitial region is that the low-density water molecules in the interstitial region form a stable HB network with a relatively large tendency of the orientational alignment, which should enhance site dipoles. This constrained water configuration should be energetically favorable but need certain entropic cost. This cost should become larger as rc decreases, so that the low-density constrained configuration is replaced by the state with much lower density. In this way the fluctuating caging and drying coexistent at around rc ≈ rseparate should turn into the definite drying for smaller rc. If this explanation of the change in hydration is correct, the density fluctuation in the interstitial region should be observed as a cooperative transition between the state with Fint < 1.0 and the state with Fint , 1.0. Distributions of Fint of 20 trajectories are shown in Figure 17. For a pair of SWCNs, the distribution of Fint has a single peak. On the other hand, for the system of a triplet of SWCNs, the distribution has two peaks that are separated by a narrow well. Existence of multiple peaks for a triplet of SWCNs suggests that the density fluctuation is cooperative and there is large fluctuation between different hydration states.

Figure 17. Distribution of frequency of trajectories to have the site density, Fint, averaged within the rectangular box defined in Figure 9: (a) distribution of frequency of trajectories to have the site density, Fint, between a pair of SWCNs with the distance, rc ) 12.80 Å (gray), rc ) 12.60 (solid line), and rc )12.40 (light gray); (b) distribution of frequency of trajectories to have the site density, Fint, among a triplet of SWCNs with the distance, rc ) 13.22 Å (gray), rc ) 13.04 (solid line), and rc )12.85 (light gray).

4. Topology and Hydration In this paper we investigated hydration around a single SWCN, a pair of SWCNs, and a triplet of SWCNs. For a single SWCN, drying is suppressed due to the attractive dispersion interactions between water molecules and carbons in SWCN. Instead, the fluctuating caging is found in the first hydration layer of SWCN. Water molecules are less mobile in the first hydration layer than in bulk. When a pair of SWCNs are separated enough, their hydration structure is similar to that around a single SWCN: drying is not evident and the first hydration layer is structured with the fluctuating caging. When two SWCNs approach, however, the interstitial region between SWCNs tend to dry. The density fluctuates largely from trajectory to trajectory and the repulsive force between two SWCNs are weakened when the interstitial region has lower density of water molecules. This is the fluctuating drying found in the interstitial region. In ref 26 a rather clear drying was observed between a pair of SWCNs. The SWCN observed in this paper has a smaller diameter than that in ref 26, and the definite drying does not take place as in ref 26 for a pair of SWCNs. Instead, we observed the strong trajectory-dependent fluctuation in water density in the interstitial region. This difference showed that the drying fluctuation is sensitive to the size and the curvature of the solute surface. The pair of SWCNs separated by a layer of one water molecule thickness should be stable until the free energy barrier in PMF is diminished by this fluctuating drying or the excess thermal energy is added to overcome the barrier. When a triplet of SWCNs approach each other, on the other hand, the fluctuating drying appears even when the SWCNs are separated enough and by further reducing the distance among SWCNs, water molecules are expelled from the interstitial region to leave a vacant pipe pinched by three SWCNs. Corresponding to the

2870 J. Phys. Chem. C, Vol. 111, No. 7, 2007 predominance of this drying, the free energy barrier in PMF is small and there is strong tendency of aggregation for a triplet of SWCNs. In this way we showed that the hydration and hydrophobic interaction are largely dependent on the topology of the configuration of SWCNs. This is in contrast to the hydration of small hydrophobic molecules such as methane. Raschke et al. conducted MD simulations of methane molecules solved in aqueous solution.39 They showed that the free energy of adding a single methane molecule to an existing cluster of n molecules of methane is a decreasing function of n and suggested that aggregation of methane is cooperative. Shimizu and Chan40 and Chan et al.,41 on the other hand, discussed that aggregation may look cooperative even when interactions among solute molecules are pairwise additive. In ref 42 Shimizu et al. compared the interaction between a pair of methane molecules and interactions among a triplet of them and showed that the PMF is almost additive without definite evidence of cooperativity: PMF(rpeak) - PMF(rseparated) ) 0.61 kcal/mol in the system of a triplet of methane molecules is about two times as large as PMF(rpeak) - PMF(rseparated) ) 0.31 kcal/mol in the system of a pair of methane molecules, where rpeak is the distance at the top of the barrier in PMF.42 The temperature dependence of PMF is different for different solute configurations, showing that the degree of this additivity varies depending on the geometrical configuration of small solutes.43 For the nanometer scale solutes of our case, the interactions are far from additive: PMF(rpeak) - PMF(rseparated) ) 0.547 kcal/mol in system III is 0.186 times as large as PMF(rpeak) - PMF(rseparated) ) 2.93 kcal/mol in system II. The large dependence of hydration on topology of the solute configuration can be also found for spherical solutes. When a pair of C60 spheres are solved in bulk water, there is no definite drying transition but the water density between two solutes fluctuates largely from trajectory to trajectory, showing the coexistence of the fluctuating drying and the fluctuating caging. In Figure 18 two of C60 molecules are placed in between two plates of graphite. Figure 18 shows that water molecules are expelled from the interstitial region between a pair of C60s when the center-of-mass distance between C60 spheres are fixed to be rc ) 12.6 Å. This drying phenomenon found between the sandwiched C60 molecules was not seen in a pair of C60 molecules solved in a box.24 These examples of SWCNs and C60s showed that the definite drying takes place when a narrow region is enclosed by some nanometer scale hydrophobic objects. A plausible explanation for this topology-dependent drying found in the present paper is that the structure restriction on water molecules determines the free energy cost of hydration. The constrained water configuration in such a narrow region should largely reduce entropy to form the consistent HB network. When the free energy cost arising from the entropy reduction in such configuration is too large, the configuration should not be realized and water molecules are expelled from the region to leave the vacancy. In this work we have found the topology-dependent drying phenomena in simulated nanometer scale systems and proposed the hypothesis that the drying is induced by the free energy cost to form the HB network in a narrow region. It should be interesting to further study the topology-dependent drying in the folding process of proteins and the binding process of proteins to form protein complexes: The process that water molecules are squeezed out from the inner region of a protein in the last stage of protein folding and the process of expelling

Hotta and Sasai

Figure 18. Snapshot of MD simulation of the system of 207 water molecules and a pair of C60 molecules bound by two plates of graphite. (a) Side view of the snapshot, where water molecules and carbons in graphite that locate in front of C60 molecules are not shown. (b) Top view of the snapshot, where a top layer of graphite is not shown. Graphite and C60s are treated as rigid bodies whose orientation and position are fixed during the simulation. Two plates of graphite are separated by the distance 12.9 Å. The center-of-mass distance between a pair of C60s is 12.6 Å. Simulation was performed under an NVT ensemble with 298 K in a box of the periodic boundary condition with the side length Lbox ) 31.064 Å. The averaged water density in the simulated result was 0.995 g/cm3. Coulomb interaction was treated with the particle mesh Ewald methd.

water from the region sandwiched by multiple proteins should depend on the topology of those protein conformations. The probability whether the region is dried or wetted should be regulated by the probability of the HB network formation through the flexible closing or opening of the narrow region surrounded by protein chains. To check this hypothesis, it is necessary to further investigate the relation between the HB network formation and the fluctuating drying or the drying transition in systems of nanometer scale, and we believe that such investigation should open a new way to analyze water mediated interactions and to clarify the mechanisms of biomolecular structure formation. Acknowledgment. This work was supported by grants from the Ministry of Education, Culture, Sports, Science, and

Fluctuating Hydration Structure in SWCNs Technology, Japan, and from Japan Society for Promotion of Science and by grants for the 21st century COE program for Frontiers of Computational Science. References and Notes (1) Frank, H. S.; Evans, M. W. J. Chem. Phys. 1945, 13, 507. (2) Lazaridis, T. Acc. Chem. Res. 2001, 34, 931. (3) Southall, N. T.; Dill, K. A.; Haymet, J. Phys. Chem. B 2002, 106, 521. (4) Shimizu, S.; Chan, H. S. J. Chem. Phys. 2000, 113, 4683. (5) Stillinger, F. H. J. Solution Chem. 1973, 2, 141. (6) Be´rard, D. R.; Attard, P.; Patey, G. N. J. Chem. Phys. 1993, 98, 7236. (7) Wallqvist, A.; Berne, B. J. J. Phys. Chem. 1995, 99, 2885. (8) Wallqvist, A.; Berne, B. J., J. Phys. Chem. 1995, 99, 2893. (9) Hummer, G.; Garde, S. Phys. ReV. Lett. 1998, 80, 4193. (10) Ball, P. Nature 2003, 423, 25. (11) Parker, J. L.; Claesson, P. M.; Attard, P. J. Phys. Chem. 1994, 98, 8468. (12) Tyrrell, J. W. G.; Attard P. Phys. ReV. Lett. 2001, 87, 176104. (13) Steitz, R.; et. al. Langmuir 2003, 19, 2409. (14) Jensen, T. R.; et al. Phys. ReV. Lett. 2003, 90, 086101. (15) Lum, K.; Chandler, D.; Weeks, J. D. J. Phys. Chem. B 1999, 103, 4570. (16) Huang, X.; Margulis, C. J.; Berne, B. J. Proc. Natl. Acad. Sci. U.S.A. 2003, 100, 11953. (17) Wallqvist, A.; Gallicchio, E.; Levy, R. M. J. Phys. Chem. B 2001, 105, 6745. (18) Li, J.; Liu, T.; Li, X.; Ye, L.; Chen, H.; Fang, H.; Wu, Z.; Zhou, R. J. Phys. Chem. B 2005, 109, 13639. (19) Choudhury, N.; Pettitt, B. M. J. Am. Chem. Soc. 2005, 127, 3556. (20) Huang, X.; Zhou, R.; Berne, B. J. J. Phys. Chem. B 2005, 109, 3546. (21) Zhou, R.; Huang, X.; Margulis, C. J.; Berne, B. J. Science 2004, 305, 1605.

J. Phys. Chem. C, Vol. 111, No. 7, 2007 2871 (22) Hua, L.; Huang, X.; Zhou, R.; Berne, B. J. J. Phys. Chem. B 2006, 110, 3704. (23) Liu, P.; Huang, X.; Zhou, R.; Berne, B. J. Nature 2005, 437, 159. (24) Hotta, T.; Kimura, A.; Sasai. M. J. Phys. Chem. B 2005, 109, 18600. (25) Walther, J. H.; Jaffe, R.;Halicioglu, T.; Koumoutsakos, P. J. Phys. Chem. B 2001, 105, 9980. (26) Walther, J. H.; Jaffe, R. L.; Kotsalis, E. M.; Werder, T.; Halicioglu, T.; Koumoutsakos, P. Carbon 2004, 42, 1185. (27) Hamada, N.; Sawada, S.; Oshiyama, A. Phys. ReV. Lett. 1992, 68, 1579. (28) Pearlman, D. A.; et al. Comput. Phys. Commun. 1995, 91, 1. (29) Jorgensen, W.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L. J. Chem. Phys. 1983, 79, 926. (30) Steele, W. A. Surf. Sci. 1973, 36, 317. (31) Ree, F. H. J. Chem. Phys. 1984, 81, 1251. (32) Hummer, G.; Rasaiah, J. C.; Noworyta, J. P. Nature 2001, 414, 188. (33) Gordillo, M. C.; Marti, J. Chem. Phys. Lett. 2000, 329, 341. (34) Werde, T.; Walther, J. H.; Jaffe, R. L.; Halicioglu, T.; Koumoutsakos, P. J. Phys. Chem. B 2003, 107, 1345. (35) Jaffe, R. L.; Gonnet, P.; Werder, T.; Walther, J. H.; Koumoutsakos, P. Mol. Simul. 2004, 30, 205. (36) Caldwell, J. W.; Kollman, P. A. J. Phys. Chem. 1995, 99, 6208. (37) Higo, J.; Sasai, M.; Shirai, H.; Nakamura, H.; Kugimiya, T. Proc. Natl. Acad. Sci. U.S.A. 2001, 96, 5961. (38) Choudhury, N.; Pettitt, B. M. J. Phys. Chem. B 2005, 109, 6422. (39) Raschke, T. M.; Tsai, J.; Levitt, M. Proc. Natl. Acad. Sci. U.S.A. 2001, 98, 5965. (40) Shimizu, S.; Chan, H. S. J. Chem. Phys. 2001, 115, 3424. (41) Chan, H. S.; Shimizu, S.; Kaya, H. Mehods Enzymol. 2004, 380, 350. (42) Shimizu, S.; Chan, H. S. Proteins 2002, 48, 15. (43) Moghaddam, M. S.; Shimizu, S.; Chan, H. S. J. Am. Chem. Soc. 2005, 127, 303.

Fluctuating Hydration Structure around Nanometer-Size Hydrophobic

Recommend Documents