5462
J. Phys. Chem. B 2000, 104, 5462-5470
Molecular Dynamics Simulations of Octyl Glucoside Micelles: Structural Properties Stephen Bogusz, Richard M. Venable, and Richard W. Pastor* Biophysics Laboratory, Center for Biologics EValuation & Research, Food & Drug Administration, RockVille, Maryland 20852-1448 ReceiVed: January 12, 2000; In Final Form: March 22, 2000
Molecular dynamics simulations were used to explore the effect of aggregate size on structural properties of octyl glucoside (OG) micelles. Systems included micelles of 1, 5, 10, 20, 27, 34, 50, and 75 surfactant molecules in water, as well as an OG bilayer, and neat octane. Micelles consisting of 5 lipids were unstable, while those of 10 or more remained intact (except for rare single lipid escapes) during the 4 ns simulations. Aggregate shape and internal properties (tail length, dihedral angle distributions, and isomerization rates) change little with size. In contrast, surface properties (hydrophobic accessible surface area and headgroup cluster size) vary with size, a consequence of the decreasing surface area-to-volume ratio in larger aggregates. Micelles in the studied size range are nonspherical with rough, uneven surfaces. Their tail lengths are comparable to those in the bilayer and neat octane, indicating that packing into a compact highly curved shape does not significantly distort the tails. Instead, the preference of the tails for an extended conformation may drive the micelle to exhibit a rough surface with large hydrophobic patches.
Introduction Amphiphilic molecules consist of a hydrophilic head (either polar or charged) and a hydrophobic tail (usually a hydrocarbon chain). They have long been of interest to biophysicists because of their ability to dissolve hydrophobic molecules in aqueous environments and to self-assemble into a variety of structures such as bilayers, vesicles, and micelles. Amphiphiles that form micelles usually have a large head and a small tail. Above a certain concentration, defined as the critical micelle concentration (cmc), these amphiphiles will self-assemble into micelles with their headgroups on the exterior where they shield the hydrophobic tails from interactions with water. 1-O-n-octyl β-D-glucopyranoside (commonly referred to as octyl glucoside or OG) is a nonionic detergent consisting of an octane attached to glucose (Figure 1). It has been used to solubilize, reconstitute, purify, and crystallize membrane proteins and membrane associated-protein complexes without denaturation,1-3 and has been of interest as an emulsifier, a cleaning agent, and a drug carrier due in part to its nontoxic and biodegradable nature.4 Along with β-linked OG, an R-linkage of the hydrocarbon tail to the glucose headgroup has been synthesized (Figure 1). Both form micelles, though R-OG has a lower solubility, forms larger aggregates, and has a lower cmc than β-OG.4,5 Many experimental measurements of β-OG micelles appear to conflict with one another. For example, aggregation numbers of 21,1 27,6 65,7 68-84,8 70,9 and 27-9710 have been measured near the cmc. These differences can be partly explained by different experimental temperatures and concentrations.1,10-12 The micelle sizes measured near the cmc and at the same temperature as the simulations described here (25 °C) range from 211 to 65.7 Commercial preparations of OG have been found to contain impurities which significantly affect measured properties.10 Micelle structure has been explored by a number of lattice simulations13-31 in which simple model lipids are placed on a grid and equilibrated using Monte Carlo techniques. These 10.1021/jp000159y
Figure 1. R-octyl glucoside (left) and β-octyl glucoside (right) with heavy atom numbering scheme.
methods benefit from rapid calculation of system energy, allowing a much larger conformational space to be explored. This speed comes at the cost of possible artifacts induced by system simplification such as constraining headgroups to the outermost spherical shell.14-16 Molecular dynamics (MD) simulations have also been used to explore micelle properties. MD simulations of micelles by other groups include simple chainlike model systems,32-34 and atomic level models of real surfactants such as sodium dodecyl sulfate (SDS),35,36 n-decyltrimethylammonium chloride,37 lyso-
This article not subject to U.S. Copyright. Published 2000 by the American Chemical Society Published on Web 05/19/2000
Structural Properties of OG Micelles by MD TABLE 1: Potential Parameters for Region Linking Glucose Head and Octane Taila bond
b0 (Å)
kb (kcal Å-2 mol-1)
C7-O1 C1-O1
1.4450 1.4066
296.7000 384.0792
angle O1-C7-C8 O1-C7-H7 O1-C1-C2 O5-C1-O1 H1-C1-O1 C1-O1-C7
θ0 (°) 105.0000 107.2400 107.6019 115.7322 109.3850 107.5000
kθ (kcal mol-1 radian-2) 110.0000 45.2000 112.2085 74.2586 52.5070 110.0000
dihedral O1-C1-O5-C5 C3-C2-C1-O1 O1-C1-C2-H2 O1-C1-C2-O2 C7-O1-C1-O5 C2-C1-O1-C7 C8-C7-O1-C1 X-C1-O1-X X-C7-O1-X
periodicity 1 2 3 1 2 3 1 2 3 1 2 3 1 3 1 3 1 3 3 3
kφ (kcal mol-1 radian-2) 1.9193 1.0102 0.7294 -1.914 -0.3739 -0.0340 0.0000 0.0000 0.1472 -4.9362 0.2907 0.4638 -1.3500 0.2000 0.4500 0.1500 0.4500 0.1500 0.1000 0.1000
a Equilibrium bond length b , bond force constant k , equilibrium 0 b angle value θ0, angle force constant kθ, and dihedral force constant kφ. The O1 oxygen is an ether oxygen type from Reiling et al. (ref 46) with charge -0.300/e. Carbon C7 is a methylene carbon identical to carbons 8 to 13 except for a net charge of -0.110/e. All other charges are as in glucose and octane residues.
phosphatidylethanolamine,38 and sodium octanoate.39-42 Of the simulations on nonmodel systems, all included a charged headgroup except for lysophosphatidylethanolamine,38 which is zwitterionic. MD of OG micelles, which have nonionic glucose headgroups, significantly extends the breadth of micelle simulations. To investigate the structural and dynamic properties of OG micelles we have run MD simulations on a variety of OG aggregate sizes. The sizes chosen span much of the experimental range of 20-100 lipids/micelle, and include micelles of 5 and 10 lipids to explore stability of small aggregates. In addition, simulations of neat octane, a single β-OG, an R-OG micelle, and a β-OG bilayer were run for comparison. This work considers the stability of various aggregate sizes, micelle shape, tail structure, clustering of glucose headgroups on the micelle surface, and a comparison with the results of lattice simulations. Methods Parameters. Simulations were carried out with the CHARMM (Chemistry at HARvard Molecular Mechanics)43 simulation package version 26. Parameters used were from CHARMM version 26, which includes glucose parameters developed by Ha et al.44 and alkane parameters developed by MacKerell et al.45 Parameters for linking the glucose head to the octane tail were based either on those of Reiling et al.46 or on similar CHARMM parameters (see Table 1). Water molecules were modeled by a modified rigid TIP3P potential.47 Electrostatic interactions were evaluated with the smooth particle-mesh ewald
J. Phys. Chem. B, Vol. 104, No. 23, 2000 5463 (PME) method,48 with a real space cutoff of 10 Å, a Gaussian width of 0.32 Å-1, and a sixth order cubic spline for interpolation between grid points (described in detail in Feller et al.49). Grid sizes were between 1.0 and 1.5 Å/grid. A 2 fs time step was used, with covalent bonds to hydrogen atoms held rigid with the SHAKE algorithm.50 The Lennard-Jones interactions were truncated between 8 and 10 Å with a smooth switching of the potential. Accessible surface area (ASA) was calculated using the Lee and Richards method51 with a probe radius of 1.6 Å. Error Estimation. Estimates listed for error of measured properties are the standard errors of the mean. The standard error (equal to the standard deviation divided by the square root of the number of independent samples) was calculated in one of two ways. For properties of individual lipids (tail length, angle between tail and the radial vector, number of dihedral transitions), the number of independent samples equaled the number of lipids in the micelle. For properties of the entire micelle (radius of gyration, ASA, moments of inertia, eccentricity) the number of independent samples were determined from the statistical inefficiency.52-54 Initial Aggregate Structure. To generate a series of random configurations, a single β-OG lipid was run in a vacuum at 293 K using Langevin dynamics with a collision frequency of 5 ps-1. After initial equilibration for 500 ps, coordinates were saved every 20 ps until 500 conformations were obtained. A series of random micelle conformations were then built for 5, 10, 20, 27, 34, 50, and 75 lipids (referred to here as a 5mer, 10mer, 20mer, etc.). To build these micelles, first a conformation from the vacuum set of 500 monomer configurations was chosen at random and placed with the O1 atom a set distance from the origin along a random direction. This distance (14 Å for the 34mer) was scaled with size. The lipid was rotated randomly until the tail was pointing roughly inward (the O1 atom at least 2.0 Å farther from the origin than the center of the tail). If the lipid had excessive overlap energy (> 50,000 kcal/mol) with the other placed lipids it was given successive small translations and rotations (keeping the O1 distance to the origin fixed) until the overlap was removed. The tails were then harmonically restrained within a sphere around the origin (12 Å for the 34mer) and the system energy minimized with 50 steps of steepest descents (SD) followed by 2000 steps of the adopted basis Newton-Raphson (ABNR) algorithm. The spherical restraints were removed and the system reminimized with 50 SD and 2000 ABNR. The micelle was then immersed in a previously equilibrated box of water molecules. All water molecules within 2.5 Å of the micelle were deleted and the box size set to be approximately 15 Å larger than the micelle diameter (25 Å for the monomer, 35 Å for the 5mer, 45 Å for the 10mer, 20mer and 27mer, 50 Å for the 34mer, 55 Å for the 50mer, and 60 Å for the 75mer). The system was minimized again with 50 steps of SD and 2000 of ABNR. This process was repeated to generate 100 micelles of each size. The conformation with the most favorable water/lipid interaction energy served as the initial condition for each of the MD simulations, except for those of the 10mer and 20mer where the conformation with the lowest tail accessible surface area was used. Simulation Technique. Simulations consisted of a heating phase in which the system temperature was increased from 98 to 298 K in 10 K increments over 20 ps, a 100 ps equilibration, followed by 4 ns of production (where averages were accumulated). All simulations (other than the bilayer) were carried out in the constant pressure and temperature (NPT) ensemble using the Langevin Piston method55 with a piston mass of 250
5464 J. Phys. Chem. B, Vol. 104, No. 23, 2000
Bogusz et al.
amu and a collision frequency of 20 ps-1. The temperature was maintained at 298 K by coupling to a Hoover thermostat56 with a mass of 1000 kcal-ps2. The R-OG 27mer was built and simulated using the same method as the β-OG micelles. For all micelle simulations other than the 75mer, nearly all C4-C5-C6-O6 dihedral values (dihedral angles will henceforth be referred to only by their 4 atoms) were close to 180° at the beginning of the simulation and 60° at the end. This situation resulted from the vacuum monomer trajectory used to produce the initial lipid conformations, where C4-C5-C6O6 was almost exclusively 180°. Initial C4-C5-C6-O6 values were set to 60° for the 75mer simulation, which was generated after this population shift was discovered. This did not have a significant effect on the resulting properties, as expected for a peripheral hydroxyl group with limited structural impact. The 1 ns bilayer simulation included 200 lipids and 584 water molecules in a 60.09 × 60.09 × 28.8 Å unit cell. These dimensions and water content are consistent with experimental measurements on the β-OG lamellar phase4 which yielded a repeat distance of 28.8 Å and an area per surfactant of 38.1 Å2 for 84.7 wt % β-OG. The cell lengths were fixed along the x and y directions while the cell length along z was allowed to fluctuate (the NPAT ensemble). The collision frequency was set to 25 ps-1, and the piston mass for the z direction to 215 amu. The mass of the Hoover thermostat equaled 2000 kcal-ps2. The octane simulation consisted of 100 molecules and was run for 8 ns using the same parameters as the micelle simulations. The average cell size during the simulation corresponded to a density of 0.66 g/cm3, compared to the experimental density of 0.70 g/cm3. Ellipsoidal Fit. The micelles were approximated as ellipsoids in order to characterize their shape. The equation for an ellipsoid with its longest axis aligned along the x direction, the second longest along y, and the shortest along z is
x2 y2 z2 + + )1 a2 b2 c2
(1)
The lengths of the semi-axes along the x, y, and z directions are a, b, and c respectively. For a prolate ellipsoid b ) c, and for an oblate, b ) a. The eccentricity, e, is defined as
e)
x
1-
c2 a2
(2)
For a sphere e ) 0, and e f 1 for flat or needlelike shapes. The principal moments of inertia for an ellipsoid (assuming constant density) are
1 1 1 I1 ) M(a2 + b2), I2 ) M(a2 + c2), I3 ) M(b2 + c2) (3) 5 5 5 where M is the total mass and I1 g I2 g I3. For point masses
I1 )
∑mi(xi2 + yi2),
I2 )
∑mi(xi2 + zi2), I3 ) ∑mi(yi2 + zi2)
(4)
After the MD run was completed, the coordinates at the end of every ps were aligned as specified for eq 1. The principal moments of inertia were calculated from eq 4, the equivalent axes lengths for a uniform ellipsoid from eq 3, and the eccentricity from eq 2. Random Model Micelles. Ellipsoidal fits were also used to create random model micelles for comparison with MD
structures. First, a best fit ellipsoid was generated for each set of O1 coordinates used in the cluster calculations (described in the Results Section) with a Levenberg-Marquardt algorithm.57 N random points (where N is the number of lipids) were chosen on the surface of this ellipsoid. To simulate the irregular nature of the surface each point’s radial distance was then subjected to a Gaussian blur with σ ) 0.3, with values less than 0.4 and greater than 2.3 truncated. Each lipid was translated to such a point, giving a random model micelle. The interaction energy between each surfactant and its neighbors was calculated and all with greater than 5 kcal/mol were given new random coordinates to remove unrealistic contacts. The energy cutoff and the Gaussian parameters were chosen to best reproduce the interaction energy and distribution about a best fit ellipsoid found for the MD frames. To ensure that conclusions were not too strongly influenced by small changes in the preceding construction procedure, the hydroxyl hydrogens and oxygens were energy minimized after removing water molecules and fixing the positions of other lipid atoms. Results are presented for the random and the dynamics structures, minimized and unminimized. A Comment on Effective Concentration. The simulations presented here model micelles at sizes found near the cmc, which is approximately 0.025 M.10 The concentrations of OG range from 0.069 M for the 5mer to 0.28 M for the 75mer; i.e., 3 to 11 times the cmc, and a concentration range where large aggregates are observed.11,58 Nearly all effects of high concentration (such as creation of large aggregates) result from direct interaction of the surfactants. Periodic boundary conditions restrict many of these interactions (e.g., a micelle cannot collide with itself and thereby grow in size unless it is longer than the cell length). Hence, these simulations do not model OG micelles at their nominal concentrations. A system at the cmc contains 1700 water molecules per lipid, or 17 times the ratio in the present systems. Although it would be interesting to simulate a system explicitly at the cmc (at the cost of an approximately 17-fold increase in computer time), the same resources could instead be used to run multiple trajectories at the present concentrations in order to create an ensemble, explore different systems, or to carry out a single trajectory in order to observe longer time scale properties. To obtain 4 ns simulation lengths, the systems described here were limited to a maximum cell length of 60 Å. An average simulation required 6 weeks on a cluster of 16 IBM SP2 nodes. Since the micelles cannot aggregate, the high concentrations of OG should have little effect on the physical properties studied here. Simulations that attempt to explore the effects of concentration on micelle size and formation will, however, require mixtures of micelles and free monomers. Results Stability. Although numerous lipid near escapes were observed, micelles with 10 or more lipids were stable over the 4 ns simulation time. Once or twice per nanosecond a lipid protruded almost completely from the micelle, only to reinsert itself within a few hundred picoseconds. In two cases, once in the 27mer and once in the 75mer, a lipid fully escaped and diffused away, suggesting that lipid exchange with the bulk solution may be a fairly rapid process. These escaped lipids were not included in the subsequent calculations of micelle properties. The micelle in the initial 5mer simulation, however, disassociated after 2.2 ns. (A 5mer was considered disassociated when the minimum distance between any single lipid or pair of lipids and the remainder of the micelle became 5.0 Å or
Structural Properties of OG Micelles by MD
J. Phys. Chem. B, Vol. 104, No. 23, 2000 5465
Figure 2. Trajectory frames from a 5mer simulation after 0.2 ns (top left), 0.8 ns (top right), 1.4 ns (bottom left), and 1.8 ns (bottom right). Tail carbons - dark gray; head carbons, polar hydrogens, and oxygens - light gray. Water molecules are omitted for clarity.
greater and remained so for more than 100 ps.) The stability of the 5mer was further explored with nineteen additional 1.2 ns simulations. Of the twenty 5mer simulations, five lost a single lipid and three lost two or more within 1.2 ns. Figure 2 shows an example of a nearly complete disassociation. These results suggest that octyl glucoside aggregates require a minimum size larger than 5 but probably less than 10 for stability on the nanosecond time scale. The 5mer simulation’s nominal lipid concentration is three times the cmc. As discussed in the final subsection of the Methods Section, the micelles cannot interact with their images, restricting the system from modeling a true solution of that concentration. Shape. Representative dynamics frames for the 20mer and 50mer simulations are shown in Figure 3. Both underwent significant shape changes on the hundreds of ps time scale. The average radius of gyration, , as a function of micelle size is shown in Figure 4. We define an effective micellular radius, RS, based on the relationship between the radius of a solid sphere of uniform density and its radius of gyration,
RS )
x53 G
(5)
The linear fit of RS vs N1/3 yields a slope of 4.71 ( 0.07 Å/N1/3 and a y-intercept of 2.73 ( 0.22 Å. These parameters yield, in turn, RS ) 16.8 Å for a 27mer, which is slightly larger than the experimental result of 15 ( 1 Å for purported 27mers.6,10 Given the approximate nature of the spherical assumption (see Figure 3 and the discussion below) and the underlying variability in aggregation number, the agreement of simulation and experiment is satisfactory. The radial density functions for the 20mer and 50mer are plotted in Figure 5. Density values at small radii are poorly sampled due to the low volume of these shells, and thus fluctuations in this region are relatively large. The 50mer curves are shifted toward larger radii compared to the 20mer, as expected for a larger micelle. Mixing of water and headgroups decreases as the micelle size increases. This is illustrated in Figure 5, where there is greater overlap of water and headgroup densities for the 20mer than for the 50mer. Note that the
nonspherical character of the micelles causes some broadening and overlap of these density functions. Compared to the atomic level micelle MD simulations noted in the Introduction, the overlaps of the tail and water regions are larger for the OG micelles (10 to 15 Å for OG vs 6 to 12 Å for the other systems35-37,39-42), as are the spreads of the carbon positions for most systems (11 to 18 Å for OG compared with 10 to 15 Å,35,37-41 or 20 Å36). As shown in Figure 6, the OG micelles deviate significantly from the spherical (e ≈ 0.6 for N > 10), but do not become less spherical as the number of lipids increase. (triangles in Figure 6) varies between 1.3 and 1.7 for different aggregation numbers. The micelle shapes fluctuate between prolate and oblate-like over the 4 ns simulations, as illustrated by the histogram of for the 34mer (Figure 7). The 1.3-1.7 range of calculated here is consistent with most results from MD of other micelles: 1.2 to 2.0 for model systems;32-34 1.02 to 1.13 for sodium dodecyl sulfate;35,36 1.3 for n-decyltrimethylammonium chloride;37 and 1.3 to 2.3 for sodium octanoate.39-42 The averages of the instantaneous axial ratios calculated from eq 3 range from 1.3 to 1.8. Experimental values for for OG micelles range from 1.5 to 20.4,7,58,59 However, most of these experiments were carried out at concentrations 5 to 20 times the cmc, and therefore measured micelles much larger than those simulated here. La Mesa et al. found that higher concentrations lead to rod-shaped micelles, with an axial ratio of 1.5 for concentrations 5 times the cmc (consistent with the simulations described here) and increasing to 5.7 at 20 times the cmc.11 In contrast, D’Aprano et al.7 deduced an axial ratio of 5 for concentrations near the cmc (a cylinder of radius of 12 Å and a length of 60 Å, consisting of approximately 65 lipids). Their fitting procedure, however, was based on a static rod model, and it is probable that a model incorporating the ensemble of shapes found here would lead to different results. The inconsistencies between results presented here and experiment, and between the experiments themselves may also be due to conditions such as ionic strength, or even impurities in the lipid samples, which have been found to significantly alter micelle size and properties for OG.10 Hydration. As shown in Figure 8 (squares) the total ASA varies linearly with N2/3 (as expected for spheres and ellipsoids), with a slope of 620.3 ( 9.7 Å2/N2/3. The area of the micelle surface covered by the headgroups (circles in Figure 8) increases linearly at nearly the same rate as the total ASA, while the tail ASA (triangles) increases only slightly over the range of micelle sizes. As is evident from Figure 9 (and consistent with Figure 8), both the ASA/lipid ratio and the relative contribution of tails to the ASA vary inversely with micelle size. The surfaces of the 20mer and 50mer include large hydrophobic patches, the relative size of which decrease as size increases (Figure 3). Water molecules (not shown in Figure 3) are found on these hydrophobic patches amongst the headgroups, leading to part of the overlap in the water and headgroup radial density functions. This overlap and the relative size of the hydrophobic patches vary inversely with micelle size. An estimate of 44 Å2/headgroup obtained experimentally for OG micelles4 with a measured radius of gyration of 40 Å is much less than the ASA/lipid found for any of the micelles simulated here (e.g., 130 Å2 for the 75mer, where ) 17.6 Å). The experimental value is, however, close to the 38.1 Å2/lipid measured for bilayers.4 Hence, we may conclude that the large micelle probably has bilayerlike packing. As is shown in Figure 9, the ASA/lipid for the bilayer MD is 59 Å. This
5466 J. Phys. Chem. B, Vol. 104, No. 23, 2000
Bogusz et al.
Figure 3. The 20mer (top row) and 50mer (bottom row) micelles after 1, 2, and 3 ns (left to right). Shading as for Figure 2. Water molecules are omitted for clarity.
Figure 4. Average radius of gyration vs number of lipids (N) to the 1/3 power for β-OG (solid symbols) and R-OG (open symbols). Linear fit is for β-OG micelles only. Block lengths of 250 ps were used to estimate standard errors. Error bars (all less than 0.09 Å) are smaller than symbol size.
Figure 6. Eccentricity (top, diamonds), and ratios for moments of inertia (bottom) vs number of lipids calculated as explained in text. - squares, - triangles, - circles. Results for R-OG and for β-OG are shown with open and solid symbols, respectively. Block lengths of 100 ps were used to estimate standard errors. Most error bars are smaller than symbol size.
Figure 5. Radial density distributions from micelle center of mass of the 20mer (top) and 50mer (bottom) for all micelle atoms (solid line), tail (dotted), head (dot-dash), and water (dashed).
55% increase over the fixed surface area of the system is associated with surface roughness that is included in the ASA calculation.
To better characterize the hydration of the micelle, water molecules were divided into “bulk” (those further than 4 Å from any lipid atom) and “contact” (not bulk). Contact waters were then subdivided into “nonisolated” and three types of “isolated”: Type I: Those not within 3.5 Å of a bulk water, and with an ASA of less than 35.64 Å2 (1/4 the average water molecule’s ASA). Type II: Those not within 3.5 Å of a headgroup oxygen atom. Type III: Those not within 3.5 Å of any oxygen atom. Many of the Type I waters are found among the headgroups or between the headgroups and the lipid tails, others are imbedded in pockets in the hydrophobic regions, and a few are found in the interior of the hydrophobic region itself. Type II are found in hydrophobic pockets on the surface, or part of a “finger” extending into the interior. Type III waters are completely buried in the micelle’s hydrophobic interior.
Structural Properties of OG Micelles by MD
J. Phys. Chem. B, Vol. 104, No. 23, 2000 5467 TABLE 2: Average Number (with Standard Errora) of Isolated Water Molecules (Defined in Text) average water molecules per frame type I type II
micelle size
Figure 7. Histogram of the instantaneous ratios (calculated each ps) between the moments of inertia along the minor and major axes for the 34mer simulation.
10 20 27 27alpha 34 50 75
10.1 ( 0.3 26.8 ( 0.5 39.3 ( 0.7 30.1 ( 0.6 53.1 ( 1.1 76.0 ( 1.0 125 ( 1
average water molecules per frame per lipid type I type II
1.98 ( 0.11 4.99 ( 0.18 6.64 ( 0.23 5.04 ( 0.18 8.60 ( 0.39 11.0 ( 0.4 16.8 ( 0.4
1.01 ( 0.03 1.34 ( 0.02 1.51 ( 0.03 1.12 ( 0.02 1.56 ( 0.03 1.52 ( 0.02 1.69 ( 0.02
0.198 ( 0.011 0.250 ( 0.009 0.256 ( 0.009 0.187 ( 0.007 0.253 ( 0.011 0.220 ( 0.009 0.228 ( 0.005
a The number of isolated waters was counted every 10 ps and averaged over the entire 4 ns. Block lengths of 100 ps were used for standard error calculations.
TABLE 3: Lipid Length, Tail Length, and Angle Between Tail and Micelle Center of Mass with Standard Error
Figure 8. Accessible surface area (ASA) vs number of lipids (N) to the 2/3 power for the tail (triangles), head (circles), and total micelle (squares). Linear fit is for β-OG micelles only (solid symbols). R-OG micelle results are shown with open symbols. Block lengths of 250 ps were used to estimate standard errors. Error bars (all less than 55 Å2) are smaller than symbol size.
Figure 9. Accessible surface area per lipid vs number of lipids for the tail (triangles), head (circles), and total micelle (squares). R-OG results are shown in open symbols, β-OG in solid symbols. Bilayer results (gray symbols) are arbitrarily placed at N ) 0. Block lengths of 250 ps were used to estimate standard errors. Error bars (all less than 2.5 Å2) are smaller than symbol size.
Table 2 lists statistics on Types I and II. Included are the average number of isolated waters on a per micelle and per lipid basis. The average numbers of Type I waters per lipid are in the range 1.0 to 1.7, and Type II waters are 0.19 to 0.26. Thus 80% of these waters are between headgroups or between a headgroup and the hydrophobic core. The number of Type I waters per lipid increases as size increases, while the number of Type II waters per lipid are relatively constant. Type III water molecules, which are in the hydrophobic core and have no nearby oxygens, are extremely rare; only between one and three were found for any of the trajectories, and none remained in the interior for more than a few ps. Tail Structure. As listed in Table 3, the average tail length (distance between C7 and C14) and lipid length (end-to-end distance along the long axis of the moment of inertia tensor) show little dependence on micelle size. The tails are close in
system
lipid length (Å)
tail length (Å)
tail angle (°)
monomer 10mer 20mer 27mer 27alpha 34mer 50mer 75mer bilayer
14.98 ( 0.11 15.05 ( 0.10 15.14 ( 0.04 15.14 ( 0.03 14.03 ( 0.07 15.15 ( 0.04 15.23 ( 0.03 15.26 ( 0.03 15.18 ( 0.05
8.09 ( 0.05 8.21 ( 0.02 8.21 ( 0.02 8.20 ( 0.01 8.19 ( 0.02 8.20 ( 0.01 8.25 ( 0.01 8.24 ( 0.01 8.24 ( 0.02
s 48.0 ( 1.8 46.6 ( 1.3 44.7 ( 1.5 44.5 ( 1.4 43.6 ( 1.3 41.5 ( 1.2 38.0 ( 0.9 41.7 ( 1.4
length to those in the bilayer, suggesting similar environments. They are significantly longer than that of the monomer, due to unfavorable hydrophobic interaction between the monomer tail and water. All the micelle tails are slightly longer than neat octane (8.166 ( 0.003 Å). Lipid length is indistinguishable between systems except for the R-OG micelle, as expected from its less extended shape (shown in Figure 1). In contrast to the tail and lipid length, the average angle between the tail and radial vector from the micelle center of mass decreases with size (last column of Table 3). The average angle for the 75mer is indistinguishable from the average angle between the lipid tails and the bilayer normal. The nonradial nature of the micelle tails is illustrated in Figure 10, which provides a cutaway view of the micelle interior. Some of the tails in each micelle point toward the center of mass, but many do not, resulting in the large average angle presented in Table 3. Figure 10 also illustrates the phenomenon of alignment of small clusters of adjacent lipid tails. It is striking how different these structures are compared to the classical micelle diagram, with lipids arranged like the spokes of a wheel, or the gridbased micelle models. Table 4 lists the trans to gauche rate constant ktg60 and the fraction of tail dihedral angles in the trans conformation, ft. The C11-C12-C13-C14 dihedral angle has the lowest ft for most sizes, due to the lower steric conflict found at the terminus of a hydrocarbon chain. The monomer chain has, in almost all cases, a smaller ft than that of the corresponding micelle dihedral, consistent with its shorter average tail length. ft is roughly constant for each dihedral across the micelle size range. The one exception is O1-C7-C8-C9 which has a much larger gauche content in the 10mer than in the other micelle sizes or the monomer. This may be due to the conformational constraints of packing the headgroups around such a small micelle. Values of ft for both the bilayer and octane are indistinguishable from the corresponding micelle dihedrals. The trans to gauche transition rates follow (in most cases) the pattern one would expect; ktg is smallest for dihedrals nearest to the glucose head (due to steric constraints inhibiting both
5468 J. Phys. Chem. B, Vol. 104, No. 23, 2000
Bogusz et al. cluster. At this point a new seed was chosen from the remaining lipids not yet assigned to a cluster, and the process repeated until all lipids had been placed in a cluster (even if the cluster only contained a single lipid). The cluster sizes were computed every 100 ps and averaged over the length of the simulation. Both number averages (Mn) and weight averages (Mw) were calculated according to eq 6:
NiMi NiMi2 ∑ ∑ Mn ) , Mw ) ∑Ni ∑NiMi
Figure 10. Micelle interiors at 4 ns for the 10mer (left), 34mer (center), and 75mer (right). Bottom views are created from 90° rotations about the horizontal. Only lipids whose centers of mass lie within 7 Å of the micelle midplane parallel to the plane of the page are included. For clarity, the radii of the tail carbons are reduced, and headgroup atoms for each lipid are replaced by a single sphere whose diameter is approximately that of the headgroup. The micelle centers are, in fact, densely packed; i.e., the apparent voids are filled by the deleted lipids.
head rotation and rotation of the whole hydrocarbon tail) and increases toward the hydrocarbon chain terminus. The monomer has larger ktg values, especially for O1-C7-C8-C9, due to lack of steric conflicts with other lipids. Isomerization rates do not vary with micelle size and are very similar to those in the bilayer, suggesting that the interior is similar for all of the OG systems simulated here. The ktg values for the octane dihedrals are larger than those in the micelles or bilayer, suggesting the OG tails have decreased mobility caused by tethering to the glucose head. To determine if the simulations have an order gradient in the interior as a function of radius as proposed by Dill and Flory,14 dihedral transitions for the tails of the 34mer were evaluated separately for the three (overlapping) regions from the micelle center of mass: 0 to 10 Å, 5 to 15 Å, and greater than 10 Å. The trajectory was divided into 100 ps segments, and a dihedral was defined as part of a region if, at the start of the 100 ps, all four atoms were within the region’s cutoff values. The average numbers of transitions per dihedral per 100 ps for the 3 regions are 1.72 ( 0.19, 1.57 ( 0.06 and 1.50 ( 0.21 respectively; i.e., the values are statistically indistinguishable. Values of ft were evaluated over 1 Å intervals along the radial vector. They are constant at 0.80 for all micelles except for the 75mer, where ft is 0.85 for dihedrals less than 7 Å from the center of mass and 0.80 otherwise. Glucose Headgroup Clusters. As illustrated in Figure 3, the headgroups do not appear to be arranged randomly across the surface; instead, they form linear clusters surrounding large regions of hydrophobic tail segments exposed to the solution. Headgroup clusters were quantified in the following way. First a lipid was arbitrarily chosen as a cluster seed. All lipids whose headgroups interacted with the headgroup of this seed more favorably than a given cutoff energy were assigned to the cluster. The cluster was then extended by including all lipids whose interaction energy with the lipids of the cluster were less than the cutoff, and so on until no new lipids were added to the
(6)
where Mi is the cluster size in number of lipids, and Ni the number of clusters of that size. For a cutoff value of -0.5 kcal/ mol the 50mer has a Mn of 12.1 and a Mw of 33.8 lipids/cluster. Upon minimization of the hydroxyl groups (as described in the Methods Section), these averages increased to 13.3 and 34.9 lipids/cluster. A certain level of clustering would be expected for a random placement of headgroups on the micelle surface. A series of random model micelles were generated (as described in the Methods Section) to estimate this effect. The cluster sizes and standard deviations for a random model 50mer are: 2.8 ( 0.5 (Mn, unminimized); 7.1 ( 2.8 (Mw, unminimized); 4.9 ( 2.9 (Mn, minimized); 13.4 ( 8.9 (Mw, minimized). Hence, the minimization affected the randomly generated structures more than those from the dynamics simulation. Nevertheless, the cluster sizes for the MD are larger by multiple standard deviations, and thus statistically significant. Similar trends and conclusions were obtained for other energy based and distance based cutoffs. The length of time clusters remain intact (defined as the time a cluster contains more than half of its original members and new lipids joining the cluster do not outnumber original cluster members) was calculated for the 50mer. If a cluster failed to meet these criteria temporarily (25 ps or less) and then reformed, it was considered intact. A new set of clusters was defined at each 500 ps time point. Each cluster of 2 or more lipids was analyzed every 5 ps to determine its persistence. The time the clusters remained intact varied from less than 5 ps to greater than 2 ns, with most less than 50 ps. Persistence is a function of cluster size. Clusters with 10 or fewer lipids lasted from less than 5 ps up to 100 ps, while those of more than 10 lipids lasted from 10 ps to greater than 2 ns. Analysis of the hydrogen bonds between headgroups with interaction energies more favorable than -1.0 kcal/mol reveals a large variety of hydrogen bond donor/acceptor combinations. O2 and O6 provide 70% of the donors and 50% of the acceptors. The number of hydrogen bonds per lipid increases approximately linearly with the size, with only an average of 0.25 hydrogen bonds per lipid existing for the 10mer and 0.87 for the 75mer. Octyl r-D Glucoside Micelle. Results from the simulation of the 27 lipid R-OG micelle were indistinguishable from that of the β-isomer for many properties: the moments of inertia ratios and eccentricity (Figure 6), tail length and angle (Table 3), dihedral distributions, and rate constants. The R-isomer, however, has a slightly smaller ASA (Figures 8 and 9) and 25 to 40% fewer isolated water molecules than the β-isomer. As evident from Figure 1, the β-isomer molecule is relatively linear, while the R-isomer contains a virtual right angle bend between head and tail. Hence, the β-isomer’s head may extend into the solution while the head of the R-isomer is more constrained to the surface of the micelle. The radial distribution function for the R-isomer has an almost identical distribution of tail carbons,
Structural Properties of OG Micelles by MD
J. Phys. Chem. B, Vol. 104, No. 23, 2000 5469
TABLE 4: Trans to Gauche Rate Constants ktga and Fraction of Dihedrals in the Trans Conformation ftb for OG Tail Dihedrals (or Corresponding Octane Dihedrals) 10mer
75mer
Bilayer
Octane
dihedral
ktg (ns-1)
ft
ktg (ns-1)
ft
ktg (ns-1)
ft
ktg (ns-1)
ft
O1 C7 C8 C9 C7 C8 C9 C10 C8 C9 C10 C11 C9 C10 C11 C12 C10 C11 C12 C13 C11 C12 C13 C14
1.24 ( 0.21 2.53 ( 0.28 2.43 ( 0.27 4.13 ( 0.36 2.94 ( 0.30 6.48 ( 0.47
0.691 ( 0.027 0.790 ( 0.016 0.833 ( 0.017 0.809 ( 0.013 0.840 ( 0.008 0.748 ( 0.012
0.67 ( 0.05 1.60 ( 0.08 1.71 ( 0.08 2.87 ( 0.11 3.74 ( 0.12 5.71 ( 0.16
0.813 ( 0.013 0.793 ( 0.009 0.846 ( 0.007 0.813 ( 0.006 0.837 ( 0.004 0.765 ( 0.005
1.12 ( 0.04 2.42 ( 0.06 1.89 ( 0.05 3.16 ( 0.07 4.27 ( 0.08 5.85 ( 0.10
0.759 ( 0.019 0.794 ( 0.013 0.838 ( 0.009 0.798 ( 0.009 0.836 ( 0.006 0.776 ( 0.007
6.22 ( 0.10 6.42 ( 0.07 7.66 ( 0.08
0.802 ( 0.003 0.811 ( 0.002 0.756 ( 0.002
a Calculated using correlation functions as described by Zhang and Pastor (ref 60). Dihedral transitions were modeled as a Poisson process. Consequently, the error in ktg was estimated as the square root of (ktg/ttrans), where ttrans is the total time in the trans state (i.e., 4 ns × N × ft). b Calculated separately for each lipid and averaged.
a headgroup distribution closer to the micelle center, and a decreased water overlap with the headgroups, supporting this hypothesis. Headgroup clusters on the surface of the R-isomer are larger (an average of 5.7 lipids per cluster (energy cutoff of -0.5 kcal/mol)) than for the β-isomer 27mer (an average of 3.9). Discussion Micelle structure has long been understood to result from hydrophobic aggregation of nonpolar tails and energetically favorable hydration of polar heads. Small micelles are assumed to be spherical, increasing in size with increasing aggregation until a maximum radius equal to the extended length of a lipid molecule is reached. At higher aggregation numbers, micelles would become rod-shaped (or ellipsoidal) with their radial axis limited by the extended lipid length. Micelles smaller than the maximum spherical size would be expected to have tails in less extended forms and, therefore, with an increased number of gauche torsion angles. The structures of the different size micelles described here (aggregation numbers of 10 to 75) do not conform to these simple assumptions. Instead of having a static shape, each fluctuates between spherical and various approximately ellipsoidal forms. The average eccentricity (≈0.6) and average of the ratios of the moments of inertia (1.1-1.7) does not show strong size dependence (Figure 6); i.e., the micelles show no sign of becoming more rod-shaped with increasing size. The average length of OG is 15.2 Å (independent of aggregation number), which is smaller than the micellular radius, RS, for all but the 10mer (eq 5). This unintuitive result arises because the micelles are not static spheres and their surfaces are very rough (Figures 3 and 5). Trans to gauche isomerization rates are mostly independent of size (Table 4), implying that the micelle interiors are similar. The average tail length, 8.2 Å, is also independent of micelle size (Table 3) and close to the fully extended length of 9.0 Å. Hence, there are no significant distortions of the tail dihedral angles as micelles change size. Rather, the average angle between tails and the radial vector increases as size decreases to resolve packing conflicts (Table 3 and Figure 10). Headgroup clustering and large exposed hydrophobic patches are evident in Figure 3, and were demonstrated to be nonrandom. The size dependence of these surface properties (Figures 8 and 9) is consistent with the area/volume increase of a smooth surface. Because the surface is not smooth, it may be surmised that the level of roughness is relatively constant. While there were very few instances of completely isolated water molecules in the micelle interior (Type III, as defined in the text), approximately 2 waters/lipid were found to associate with the headgroup from “underneath” or were otherwise partially buried (Table 2, and associated discussion of Type I and Type II).
The results were also compared with those from simulations of an OG bilayer and neat octane. The average tail lengths for the micelles and the bilayer are indistinguishable. The average angle of the tail with the radial vector for the 75mer is nearly the same as the average angle between tails in the bilayer and the vector normal to the surface (Table 3). One important difference between these two forms is the accessible surface area, where an OG in the bilayer has less than half the ASA of one in the 75mer, and less than 8% of its tail exposed (Figure 9). The dihedral distibutions for the micelle tails are comparable to neat octane, but the isomerization rates are significantly smaller. This result is similar to that of previous calculations for a sodium dodecyl sulfate (SDS) micelle36 where the tails had comparable dihedral distributions and decreased fluidity compared with neat decane. It differs from the results of simulations of a dipalmitoylphosphatidylcholine (DPPC) bilayer and hexadecane,61 where the fluidities of the bilayer interior and the neat alkane were almost identical. It is not surprising that OG (with its octane tail) should be more similar to SDS (a dodecane tail) than to DPPC (a hexadecane tail). Not all models, of course, are idealized. Tanford predicted a rough micelle surface,62 and Menger fit his experimental data to a model containing a rough surface, hydrophobic patches, transient water filled cavities, terminal methyls near the surface, and a nearly water free hydrophobic core,63-65 all features observed here. Nevertheless, it is interesting to compare the physical properties of OG micelles determined in this work with those of spherical grid based simulations.14-16 The micelle conformations created with this grid method are confined to a spherical geometry with headgroups at the same radial distance from the micelle center. These constraints are inconsistent with two important aspects of micelle structure observed here: that the micelles are nonspherical and that they have rough, uneven surfaces. Incorrectly modeled constraints such as these can lead to errors in other aspects of micelle structure. Dill and Flory14 find their interiors to be disordered on the outside and highly ordered at the center. As discussed in the Results Section, isomerization rates and dihedral distributions in the tail segments in the outer and inner regions of the hydrophobic core are similar. The different conclusions derived from the MD required allocating the large amount of computer resources required for atomic level simulations. The accuracy came at the cost of being limited to 4 ns and 20 000 atoms: insufficient for simulating a system consisting of multiple micelles interacting with a surrounding solution of free lipids (as discussed in the final subsection of the Methods Section). For example, in simulations of the 5mers, the micelles disassociated, even though they were run at a nominal concentration of 3 times the cmc. In a more realistic, though presently hypothetical, MD simulation, numerous micelles and a large number of free surfactants would be
5470 J. Phys. Chem. B, Vol. 104, No. 23, 2000 combined. Many 5mers would disassociate, but a few might grow to a stable size. Grid-based simulations and other simple models are much less computationally demanding, and thereby allow greater conformational sampling, albeit with simplifying assumptions. These assumptions could, however, be tested and systematically refined by comparison with MD. A combined approach might be effective for answering questions related to the free energy of micellular systems such as the optimal size (or distribution of sizes). In summary, the OG micelles studied here exhibit nonspherical shapes that fluctuate on the hundreds of ps time scale. The preference for extended tail conformations produces rough surfaced micelles. The structural aspects of micelles are strongly influenced by the structure of the amphiphile. The short tail of OG is more rigid than a longer tail would be. A longer tail may allow packing constraints to be satisfied more easily, and may thus exhibit less surface roughness. Other lipid types, especially those with charged headgroups, may imbue their micelles with different physical properties, somewhat limiting the ability to generalize about micelle behavior based solely on the OG simulations presented here. These simulations do reveal that many of the early assumptions of micelle properties need to be reexamined, and lay the foundation for a detailed understanding of octylglucoside/peptide complexes. Acknowledgment. We thank Bernard Brooks of the National Heart, Lung, and Blood Institute at the National Institutes of Health for the use of LoBoS, a Beowulf class supercomputer. Computations also utilized the IBM SP supercomputer at the Center for Information Technology, NIH, Bethesda, MD. References and Notes (1) Hantgan, R. R.; Braaten, J. V.; Rocco, M. Biochemistry 1993, 32, 3935. (2) Michel, H.; Oesterhelt, D. Proc. Natl. Acad. Sci. U.S.A. 1980, 77, 1283. (3) Garavito, R. M.; Rosenbusch, J. P. J. Cell. Biol. 1980, 86, 327. (4) Nilsson, F.; Soderman, O.; Johansson, I. Langmuir 1996, 12, 902. (5) Straathof, A. J. J.; van Bekkum, H.; Kieboom, A. P. G. Starch/ Starke 1988, 40, 438. (6) VanAken, T.; Foxall-VanAken, S.; Castleman, S.; Ferguson-Miller, S. Meth. Enzymol. 1986, 125, 27. (7) D’Aprano, A.; Giordano, R.; Jannelli, M. P.; Magazu, S.; Maisano, G.; Sesta, B. J. Molec. Struct. 1996, 383, 177. (8) Lasser, H. R.; Elias, H. H. Kolloid-Z. Z. Polymere 1972, 250, 58. (9) Roxby, R. W.; Mills, B. P. J. Phys. Chem. 1990, 94, 456. (10) Lorber, B.; Bishop, J. B.; DeLucas, L. J. Biochim. Biophys. Acta 1990, 1023, 254. (11) La Mesa, C.; Bonincontro, A.; Sesta, B. Colloid Polym. Sci. 1993, 271, 1165. (12) Kano, K.; Ishimura, T. J. Chem. Soc., Perkin Trans. 1995, 2, 1655. (13) Haan, S. W.; Pratt, L. R. Chem. Phys. Lett. 1981, 79, 436. (14) Dill, K. A.; Flory, P. J. Proc. Natl. Acad. Sci. U.S.A. 1981, 78, 676. (15) Dill, K. A.; Koppel, D. E.; Cantor, R. S.; Dill, J. D.; Bendedouch, D.; Chen, S. Nature 1984, 309, 42. (16) Dill, K. A.; Naghizadeh, J.; Marqusee, J. A. Annu. ReV. Phys. Chem. 1988, 39, 425. (17) Rodrigues, K.; Mattice, W. L. J. Chem. Phys. 1991, 95, 5341. (18) Rodrigues, K.; Kausch, C. M.; Kim, J. G.; Quirk, R. P.; Mattice, W. L. Polym. Bull. 1991, 26, 695. (19) Rodrigues, K.; Mattice, W. L. Langmuir 1992, 8, 456. (20) Wang, Y. M.; Mattice, W. L.; Napper, D. H. Langmuir 1993, 9, 66. (21) Wang, Y. M.; Mattice, W. L.; Napper, D. H. ACS Symposium Series 1993, 532, 45.
Bogusz et al. (22) Haliloglu, T.; Mattice, W. L. Chem. Eng. Sci. 1994, 49, 2851. (23) Haliloglu, T.; Mattice, W. L. Abstracts of Papers of the American Chemical Society 1994, 207, 187. (24) Adriani, P.; Wang, Y. M.; Mattice, W. L. J. Chem. Phys. 1994, 100, 7718. (25) Zhan, Y. J.; Mattice, W. L. Macromolecules 1994, 27, 677. (26) Nguyenmisra, M.; Mattice, W. L. Macromolecules 1995, 28, 1444. (27) Haliloglu, T.; Bahar, I.; Erman, B.; Mattice, W. L. Macromolecules 1996, 29, 4764. (28) Xing, L.; Mattice, W. L. Macromolecules 1997, 30, 1711. (29) Xing, L.; Mattice, W. L. Abstracts of Papers of the American Chemical Society 1997, 213, 137. (30) Xing, L.; Mattice, W. L. Macromolecular Theory and Simulations 1997, 6, 553. (31) Xing, L.; Mattice, W. L. Langmuir 1998, 14, 4074. (32) Karaborni, S.; O’Connell, J. P. Langmuir 1990, 6, 905. (33) Karaborni, S.; O’Connell, J. P. J. Phys. Chem. 1990, 94, 2624. (34) Wijmans, C. M.; Linse, P. Langmuir 1995, 11, 3748. (35) Shelley, J.; Watanabe, K.; Klein, M. L. Int. J. Quant. Chem. 1990, 103. (36) Mackerell, A. D. J. Phys. Chem. 1995, 99, 1846. (37) Bocker, J.; Brickmann, J.; Bopp, P. J. Phys. Chem. 1994, 98, 712. (38) Wendoloski, J. J.; Kimatian, S. J.; Schutt, C. E.; Salemme, F. R. Science 1989, 243, 636. (39) Jonsson, B.; Edholm, O.; Teleman, O. J. Chem. Phys. 1986, 85, 2259. (40) Watanabe, K.; Ferrario, M.; Klein, M. L. J. Phys. Chem. 1988, 92, 819. (41) Watanabe, K.; Klein, M. J. Phys. Chem. 1989, 93, 6897. (42) Shelley, J. C.; Sprik, M.; Klein, M. L. Langmuir 1993, 9, 916. (43) Brooks, B. R.; Bruccoleri, R. E.; Olafson, B. D.; States, D., J.; Swaminathan, S.; Karplus, M. J. Comput. Chem. 1983, 4, 187. (44) Ha, S. N.; Giammona, A.; Field, M.; Brady, J. W. Carbohydr. Res. 1988, 180, 207. (45) MacKerell, A. D.; Bashford, D.; Bellott, M.; Dunbrack, R. L.; Evanseck, J. D.; Field, M. J.; Fischer, S.; Gao, J.; Guo, H.; Ha, S.; JosephMcCarthy, D.; Kuchnir, L.; Kuczera, K.; Lau, F. T. K.; Mattos, C.; Michnick, S.; Ngo, T.; Nguyen, D. T.; Prodhom, B.; Reiher, W. E.; Roux, B.; Schlenkrich, M.; Smith, J. C.; Stote, R.; Straub, J.; Watanabe, M.; Wiorkiewicz-Kuczera, J.; Yin, D.; Karplus, M. J. Phys. Chem. B 1998, 102, 3586. (46) Reiling, S.; Schlenkrich, M.; Brinkmann, J. J. Comput. Chem. 1996, 17, 450. (47) Durell, S. R.; Brooks, B. R.; Ben-Naim, A. J. Phys. Chem. 1994, 98, 2198. (48) Essmann, U.; Perera, L.; Berkowitz, M. L.; Darden, T.; Lee, H.; Pedersen, L. G. J. Chem. Phys. 1995, 103, 8577. (49) Feller, S. E.; Pastor, R. W.; Rojnuckarin, A.; Bogusz, S.; Brooks, B. R. J. Phys. Chem. 1996, 100, 17 011. (50) Ryckaert, J. P.; Ciccotti, G.; Berendsen, H. J. C. J. Comput. Phys. 1977, 23, 327. (51) Lee, B.; Richards, F. M. J. Mol. Biol. 1971, 55, 379. (52) Friedberg, R.; Cameron, J. E. J. Chem. Phys. 1970, 52, 6049. (53) Allen, M. P.; Tildesley, D. J. Computer Simulation of Liquids; Oxford University Press: Oxford, U.K., 1987. (54) Pastor, R. W.; Venable, R. M. In Computer Simulation of Biomolecular Systems; van Gunsteren, W. F., Weiner, P. K., Wilkinson, A. J., Eds.; Escom: Leiden, 1993; Vol. 2, p 443. (55) Feller, S. E.; Zhang, Y.; Pastor, R. W.; Brooks, B. R. J. Chem. Phys. 1995, 103, 4613. (56) Hoover, W. G. Phys. ReV. A. 1985, 31, 1695. (57) Marquardt, D. A. J. Soc. Ind. Appl. Math 1963, 11, 431. (58) Bonincontro, A.; Briganti, G.; Daprano, A.; La Mesa, C.; Sesta, B. Langmuir 1996, 12, 3206. (59) Thiyagarajan, P.; Tiede, D. M. J. Phys. Chem. 1994, 98, 10 343. (60) Zhang, Y.; Pastor, R. W. Mol. Sim. 1994, 13, 25. (61) Venable, R. M.; Zhang, Y.; Hardy, B. J.; Pastor, R. W. Science 1993, 262, 223. (62) Tanford, C. The Hydrophobic Effect: Formation of Micelles and Biological Membranes; 2nd ed.; Wiley: New York, 1980. (63) Menger, F. M.; Doll, D. W. J. Am. Chem. Soc. 1984, 106, 1109. (64) Menger, F. M. In Surfactants in Solution; Mittal, K. L., Lindman, B., Eds.; Plenum: New York, 1984; Vol. 1. (65) Menger, F. M. Nature 1985, 313, 603.
