Surface Tension Parameterization in Molecular Dynamics Simulations

Jul 2, 2004 - Weill Medical College of Cornell UniVersity, New York, New York 10021 ... the most important issues in the simulation methodology of lip...
0 downloads 0 Views 628KB Size
Download PDF
11802

J. Phys. Chem. B 2004, 108, 11802-11811

Surface Tension Parameterization in Molecular Dynamics Simulations of a Phospholipid-bilayer Membrane: Calibration and Effects Ramasubbu Sankararamakrishnan* and Harel Weinstein* Department of Physiology and Biophysics, Mount Sinai School of Medicine, New York, New York 10029 and Weill Medical College of Cornell UniVersity, New York, New York 10021 ReceiVed: March 8, 2004; In Final Form: May 11, 2004

Molecular dynamics simulations of interfacial systems such as phospholipid bilayers require experimental estimates of either the surface area per lipid or the surface tension (γ). To investigate the mode of selection of γ values, we have carried out the molecular dynamics (MD) simulations described here for a neat DMPC bilayer comprising 90 lipids using NPNγT and NPNAT ensembles. Influence of γ values coupled with the effect of force truncation for both NPNγT and NPNAT ensembles was analyzed in multiple simulations. The membrane structural properties were studied comparatively using different γ values in NPNγT ensembles with a spherical cutoff of 18 Å and particle mesh Ewald method (PME) that takes long-range electrostatic interactions into account. In parallel, we tested the sensitivity of calculated γ value in simulations using NPNAT ensembles with four different spherical cutoffs and PME. Our studies suggest that when PME is used to calculate the nonbonded forces, the value of surface tension for the system of this size should be in the range of 25 to 30 dyn/cm (per interface). On the other hand, for simulations with a cutoff of 18 Å, shrinkage of box size of about 7-10% is observed for γ values ranging from 0 to 85 dyn/cm. The surface tension value has to be raised to 625 dyn/cm to keep the surface area/lipid close to the experimental estimate for neat DMPC bilayers. Molecular dynamics simulations were also carried out with the opioid peptide dynorphin A(1-13) embedded within the DMPC bilayers. Several γ values were tested with two electrostatic schemes. Our results demonstrated that for γ values 0, 45, and 65 dyn/cm, due to the shrinkage of box size, the R-helical structure in the N-terminal region of the peptide is distorted with an 18 Å cutoff. When the γ value is increased to 625 dyn/cm with an 18 Å cutoff or when PME is used with γ ) 25 dyn/cm, the peptide’s R-helical structure is maintained and its orientation is similar to that obtained from earlier simulations35 using NVE ensembles. The aromatic side chains in the peptide (Tyr-1 and Phe-4) interacted similarly in all the simulations except when γ)0 is used with a cutoff of 18 Å. All constant surface tension simulations were compared with constant area ensembles in which the surface area/lipid remains fixed. Our results suggest that neglect of long-ranged electrostatic interactions and/or incorrect surface tension value results in serious artifacts in membrane simulations and suggest that caution is required when constant surface tension simulations are carried out on neat bilayers. In simulations of membranes with an embedded solute, the choice of an appropriate γ value is more complex, as it will affect both the structure and the interactions of the solute within the bilayer.

Introduction The central role of membrane proteins in the function of biological systems has increased the interest in molecular dynamics simulation studies of such systems. A key element in such studies is the representation of the lipid bilayers in the biologically relevant liquid-crystalline state. Therefore, it was very gratifying to note that MD simulations are now capable of reproducing a number of important structural and dynamic parameters observed experimentally, such as bilayer thickness and order parameter profiles.1-9 Analysis showed that among the most important issues in the simulation methodology of lipid bilayers are the nature of the force field, the boundary conditions, * Correspondence: Dr. R. Sankararamakrishnan; Department of Biological Sciences & Bioengineering, Indian Institute of Technology, Kanpur 208016, India; Email: [email protected]; Tel: +91 512 2594014; Fax: +91 512 2594010. Dr. Harel Weinstein, Department of Physiology and Biophysics, Box 75, Weill Medical College of Cornell University, 1300 York Avenue, New York, NY 10021; E.mail: [email protected]; Tel: (212) 241 7018; Fax: (212) 860 3369.

and the treatment of nonbonded interactions. Thus, force fields such as CHARMM10,11 and GROMACS12-15 have been developed for both saturated and unsaturated phospholipids and the approaches are being improved continuously, for example, with the inclusion of the Particle Mesh Ewald (PME) algorithm16 in the calculation of Coulombic forces in lipid simulations. The PME algorithm has significant advantages over traditional spherical cutoff methods because it is efficient and effectively takes into account all interactions by calculating the longdistance nonbonded forces in reciprocal space. For boundary conditions, the majority of the membrane simulations have used Periodic Boundary Conditions (PBC),1,2,5,7 either with constant volume or constant pressure algorithms. In the constant volume ensembles NVE and NVT, the energy and temperature are kept constant, respectively, and the dimensions of the box are fixed throughout the simulations. However, for interfacial systems such as lipid bilayers, it is often advantageous to simulate a flexible cell. This is possible in constant-pressure ensembles,

10.1021/jp048969n CCC: $27.50 © 2004 American Chemical Society Published on Web 07/02/2004

Surface Tension Parameterization in Simulations in which a pressure tensor P is specified. The commonly used approaches include the NPNAT and NPNγT ensembles, in which a constant normal pressure (PN) is applied with a fixed surface area (A) or a constant surface tension (γ). In both these ensembles, the length of the unit cell along the Z-axis (bilayer normal) is adjusted to maintain the constant normal pressure PN. In an NPNγT ensemble, the surface area is allowed to fluctuate in addition to the Z-axis, to keep a constant nonzero surface tension. In the NPT ensemble, the surface tension is zero and an isotropic pressure tensor is applied. In both NPNγT and NPT simulations, the simulation cell is fully flexible. In constant volume and NPNAT ensembles, an experimental estimate of the surface area per molecule is needed to set up the simulations. This parameter can be obtained from gravimetric data, X-ray diffraction, or NMR experiments.17,18 However, the surface area determined from different experiments over the years has an enormous range, and it has been shown that this quantity can vary as much as 25% for the biologically relevant fully hydrated liquid crystalline dipalmitoylphosphatidylcholine (DPPC) bilayers.17 On the other hand, in the NPNγT ensemble, it is the surface tension that must be specified. Surface tension is a macroscopic quantity and the precise value of γ is the subject of debate in the recent literature.19-22 It has been argued that the surface tension must be zero for an unstressed bilayer at its free energy minimum,20 and hence an NPT ensemble has been used in many simulations.6,9 Although direct determination of γ for bilayers is not available, based on the values obtained from analogous systems (oil/water and fatty acid/water) and also from monolayer data,23-27 a finite surface tension value has been suggested for microscopic membrane patches typically used in the simulations.2,19 The membrane patches in the simulations are only a few nanometers long, whereas the actual size of the membrane is ∼ 1 µm. Feller and Pastor21 argued that in order to produce macroscopic features, simulations of microscopic membrane patches should be carried out with a nonzero applied surface tension even if the macroscopic surface tension is zero or close to zero. A compromise has been implemented in several simulations by using the NPNAT ensemble in which the surface area of the lipid was taken from experimental values and only the Z-axis was adjusted against the atmospheric pressure.7,28 The results of these simulation studies obviously depend on the accuracy of the experimental estimate of the surface area. Simulations with a fully flexible cell should, in principle, yield the membrane dimensions (and hence the surface area per lipid) directly from the results. Allowing this flexibility in the parameters of the system is especially important since membranebound compounds such as peptides, proteins, and other molecules that are often included in the simulated systems are likely to influence their surface area and other membrane properties. Estimating the surface area per lipid in such complex heterogeneous systems will be complicated, and in these cases the use of an NPNγT ensemble seems appropriate. In the simulations of a pure dimyristoylphosphatidylcholine (DMPC) bilayer, Chiu et al.19 have used the NPNγT ensemble and applied a negative tangential pressure that corresponds to 28 dyn/cm per interface. Feller and Pastor2 simulated a membrane patch of 72 DPPC lipids and showed that the value of applied surface tension for this system should be 35-45 dyn/cm per interface. It should be noted that the choice of γ value is sensitive to the potential energy function in the force field, to the method of force truncation, and to system size.22 To investigate the mode of selection of γ values, we have carried out the molecular dynamics simulations described here for a neat DMPC bilayer comprising 90 lipids using NPNγT

J. Phys. Chem. B, Vol. 108, No. 31, 2004 11803 and NPNAT ensembles. We have analyzed the effect of force truncation for both NPNγT and NPNAT ensembles. The membrane structural properties were studied comparatively in NPNγT ensembles with a spherical cutoff of 18 Å and PME, using different γ values. In parallel, we tested the sensitivity of the calculated γ value in NPNAT ensembles for four different spherical cutoffs and compared the γ values with that obtained from the PME simulation. In addition to the neat bilayer simulations, simulations are also presented here for the membranebound opioid peptide dynorphin A(1-13), in which the structure and peptide orientation in the membrane were investigated with different γ values in NPNγT simulations using a spherical cutoff of 18 Å and PME. In total, 18 simulations of neat bilayers and 6 simulations of the peptide in the bilayers were carried out for periods of 2 ns to 4 ns. While the results of neat bilayer simulations using PME agree with the earlier simulation studies using a CHARMM force field2 utilizing similar γ values, but the same γ values with cutoff 18 Å in NPNγT ensembles yield significantly different membrane properties. Notably, the γ value obtained from the NPNAT ensemble using the 18 Å cutoff is more than 1 order of magnitude higher than that obtained using PME. Our results demonstrate that long-ranged electrostatic interactions with the proper choice of γ value are needed to produce membrane structural and physical properties that are in good agreement with experiment. Spherical cutoffs as long as 18 or 21 Å still result in serious artifacts. In the simulations of the membrane-bound peptide, the γ values are critical, because those that produce smaller surface area per lipid were found to give rise to distortions of the peptide structure, or to unlikely orientations of peptide side chains. Methods. The program CHARMM10 was used with the allatom PARM22b4b parameter set29 that includes parameters for phospholipids11 and TIP3P water potentials.30 The CPT (constant pressure and temperature) module of CHARMM was used for the constant pressure/temperature simulations. The simulations were carried out at 330 K (57 °C) and the temperature was maintained by means of the Hoover thermostat.31 The membrane normal was oriented along the Z axis with the center of the bilayer at Z ) 0. The SHAKE32 algorithm was used to fix the bond lengths involving all hydrogen atoms. A time step of 0.002 ps was used. The initial structure of pure bilayers and the membrane-bound peptide dynorphin A(1-13) were constructed using the protocol developed by Woolf and Roux.33 Neat Bilayers. NPNγT Ensembles. The periodic simulation cell contained 90 DMPC molecules (45 per monolayer) and ∼ 30 water molecules per lipid amounting to a full hydration of the bilayer (total: ∼ 18 700 atoms). In the NPNγT ensemble, a fully flexible simulation cell was employed with the Z-direction adjusted to maintain the normal pressure of 1 atm. The X and Y dimensions were varied to maintain the desired γ value. Pressure and surface tensions were maintained at their chosen values by means of the Langevin Piston algorithm that is a variant of the extended system formalism.34 A piston mass of 500 amu and a collision frequency of 5 ps-1 were used in the simulations. The initial area of the phospholipid was set at 63.1 Å2, a value close to the recent experimental estimate18 for DMPC bilayers at 50 °C. This estimate for the value at 50 °C was obtained by Nagle and Tristram-Nagle18 by using the area per DMPC molecule at 30 °C and the thermal expansion adjustment; notably, this value is quite close to the experimentally determined area per DPPC molecule at 50 °C. The periodic system consisted of a rectangular box with dimensions 53.3 × 53.3 × 62.4 Å3. Two schemes were used to calculate the nonbonded interactions, namely a spherical cutoff at 18 Å and PME.

11804 J. Phys. Chem. B, Vol. 108, No. 31, 2004

Sankararamakrishnan and Weinstein

TABLE 1: Summary of the Simulations Performed on Neat Bilayers (90 DMPC lipids) and Dyn A(1-13)-bilayer Systems system neat bilayers neat bilayers neat bilayers neat bilayers neat bilayers neat bilayers neat bilayers neat bilayers dyn A(1-13) - bilayers dyn A(1-13) - bilayers dyn A(1-13) - bilayers

cutoff/PME ensemble cutoff 18 Å PME PME-Ia cutoff 12 Å cutoff 15 Å cutoff 18 Å cutoff 21 Å PME cutoff 18 Å PME cutoff 18 Å

NPNγT NPNγT NPNγT NPNAT NPNAT NPNAT NPNAT NPNAT NPNγT NPNγT NPNAT

γ (dyn/cm) 0, 10, 20, 45, 65, 85, 625 0, 20, 25, 30, 45 45

0, 45, 65, 625 25

a A neighbor list, used for calculating the LJ potential and the real space portion of the Ewald sum, was kept to 18 Å using an atombased cutoff and updated every 10 fs. LJ potential was switched smoothly to zero over the region from 14 to 17 Å. In all other PME simulations, a neighbor list was kept to 14 Å using an atom-based cutoff and the LJ potential was switched to zero from 10 to 12 Å.

For “spherical cutoff 18 Å”, the nonbonded list was generated using a group-based cutoff of 18 Å. Both the electrostatic and van der Waals interactions were smoothly switched to zero from 14.0 to 17.0 Å. In the PME scheme, the Lennard-Jones (LJ) potential was switched smoothly to zero over the region from 10 to 12 Å. PME was used to calculate the electrostatic interactions with a B-spline order of 4 and a FFT grid of 64 × 64 × 81. A realspace Gaussian-width kappa of 0.32 Å-1 was used. A neighbor list, used for calculating the LJ potential and the real space portion of the Ewald sum, was kept to 14 Å using an atombased cutoff and updated every 10 fs. For a spherical cutoff of 18 Å, we used a series of γ values: 0, 10, 20, 45, 65, and 85 dyn/cm. A higher γ value of 625 dyn/ cm (see below) was also used with a spherical cutoff of 18 Å. Surface tensions of 0, 20, 25, 30, and 45 dyn/cm were used in the PME simulations (Table 1). The values of surface tension are given on a per-interface basis throughout this paper. To investigate the influence of van der Waals (vdW) contributions on membrane properties, an additional simulation was carried out in the PME scheme with the LJ potential switched smoothly to zero over the region 14.0 to 17.0 Å as was done in the spherical cutoff simulations. A γ value of 45 dyn/cm was used in this simulation (PME-I in Table 1). NPNAT Ensembles. Five NPNAT simulations were carried out using PME and four different spherical cutoffs. The lateral dimensions of the simulation cell, X and Y, were fixed in these simulations so that the surface area per lipid remains close to the recent experimental estimate of 63.1 Å218 throughout the simulations. Cutoff values used to calculate the nonbonded interactions were 12, 15, 18, and 21 Å. The nonbonded lists were generated using a group-based cutoff and they were smoothly switched from 8 to 11, 11 to 14, 14 to 17, and 17 to 20 Å for cutoffs 12, 15, 18, and 21 Å respectively (Table 1). The details of PME simulations are the same as those in the NPNγT simulations above. The surface tension γ (per interface) was calculated for each of the NPNAT simulations from the pressure tensor using the equation

γ ) 1/2 {}

(1)

where Lz denotes the length of the unit cell along the bilayer normal and Pzz is the component of the pressure tensor normal to the membrane surface. Pxx and Pyy are tangential components of the pressure tensor.

Dynorphin A(1-13) in DMPC Bilayers. In the periodic system, the top layer contained 41 DMPC lipids and the opioid peptide dynorphin A(1-13) [Dyn A(1-13)], and the bottom layer contained 45 lipids. The construction of the dynorphin peptide-lipid complex has been described in detail.35,36 The amino acid sequence of this peptide is YGGFLRRIRPKLK. The NMR conformation of dynorphin in DPC micelles37 was used as a structural template for the peptide in which the first nine residues assumed an R-helical conformation. Previous simulations of this peptide were carried out using an NVE ensemble.35 Here, effects of the flexible cell in influencing the peptide structure and orientation were investigated using NPNγT ensembles. As in the case of neat bilayers, two different schemes were used to calculate the nonbonded interactions, namely the spherical cutoff of 18 Å and PME. When PME was used in the simulations, 5 chloride ions (Cl-) were added near the site of the basic residues so that the entire system was electrically neutral. Surface tension values used in the 18 Å cutoff simulations were 0, 45, 65, and 625 dyn/cm. The γ value used in the PME simulation is 25 dyn/cm. For comparison purpose, a simulation using the NPNAT ensemble was also carried out with the spherical cutoff of 18 Å. The details of simulations using the spherical cutoff of 18 Å and PME are the same as those used in the neat bilayer simulations. The summary of all the simulations performed for neat bilayers and the membrane-bound peptide is given in Table 1. Simulation Details. Initial structures of neat bilayers were first coupled to a heat bath at 330 K, and Langevin dynamics simulations of 100 ps were carried out at constant volume. A planar harmonic restraint was applied on the center of mass (CM) of the lipid atoms to maintain the planarity of the membrane. A further 150 ps of simulation were carried out using an NVE ensemble during which time the restraints were gradually removed. In the case of the Dyn A(1-13)-lipid system, in addition to the restraints on the lipid atoms, positional harmonic restraints were applied on the peptide atoms during the 100 ps Langevin dynamics at constant volume. During the subsequent 150 ps simulations in the NVE ensemble, a cylindrical harmonic restraint was applied on the CM of the R-helical residues 1-9 and, as in the neat bilayers, a planar harmonic restraint was applied on the CM of lipid atoms. At the end of the 150 ps simulations, all of the restraints were removed. The systems were then simulated using NPNγT or NPNAT ensembles for a period of 2 to 4 ns. The first nanosecond was considered the equilibration phase and hence all the analyses were done on the structures saved after 1 ns. To facilitate comparison, the results from the analyses are presented from 1 to 2 ns for all the simulations. Results Neat Bilayers: Cutoff 18 Å vs PME in NPNγT Ensembles. Surface Area. Knowledge of surface area per lipid can be correlated with many of the membrane structural and physical properties (e.g., order parameters, bilayer thickness, conformation of acyl chains, etc.) and hence this is considered to be the most central structural quantity. Molecular dynamics trajectories of surface area per lipid are plotted for simulations using a spherical cutoff of 18 Å (Figure 1a) and PME (Figure 1b). Use of γ ) 0 results in contraction of the box in both the cutoff 18 Å and PME simulations. In PME, this effect is much more dramatic, as observed in earlier studies.2 The nonzero γ values

Surface Tension Parameterization in Simulations

Figure 1. Molecular dynamics trajectories of surface area per lipid for the neat bilayers using NPNγT ensembles. (A) Spherical cutoff of 18 Å and (B) PME. γ ) 0 results in contraction of the box length in both a cutoff of 18 Å and PME simulations. The results show that value of γ for a system of 90 lipids should be in the range of 25 to 30 dyn/cm when PME is used in the simulations. In the simulations using a cutoff of 18 Å, the γ value has to be increased to 625 dyn/cm to keep the surface area per lipid close to the experimental estimate of 63 Å2. This is more than an order of magnitude higher than that used in the PME simulations. In the case of NPNAT ensembles, the surface area per lipid was taken close to the experimental estimate for the neat DMPC bilayers.

yield significantly different results in PME compared to the cutoff 18 Å simulations. In PME, the surface area becomes larger (> 70 Å2) at the end of the 2 ns simulation for γ ) 30. For γ ) 45, the area expansion is sudden and very large right from the beginning of the simulation, reaching a value of more than 75 Å2 within the first 500 ps of the simulation. The average value of the surface area becomes much larger (∼ 110 Å2) at the end of the 2 ns simulation indicating that the system no longer represents the liquid crystalline state of a DMPC bilayer. For γ ) 20, the area stays reasonably close to the experimental estimate of ∼ 63 Å2

J. Phys. Chem. B, Vol. 108, No. 31, 2004 11805 during the first 1 ns of the simulation, but later drifts to ∼ 55 Å2. For γ ) 25, the area remains close to the experimental value for the major part of the simulation. Even here, at the end of the 2 ns simulation, a drop is observed in the value of the surface area. With all the γ values tested, γ ) 25 appears to give the most reasonable value for the surface area. This agrees well with similar studies carried out for a system of 72 DPPC lipids.2 For that system, the correct surface tension to apply was recommended in the range of 35-45 dyn/cm. The present system of 90 lipids is larger and hence a smaller value of surface tension is anticipated.21 With the results from the reasonably long simulations of 2 ns, our studies suggest that the Value of surface tension for the system containing 90 lipids should be in the range of 25 to 30 dyn/cm when PME is used to calculate the nonbonded forces. In the simulations using a cutoff of 18 Å, the γ values ranging from 0 to 85 dyn/cm yield properties that are significantly different from those observed in PME simulations. Even with γ ) 85, the simulation results in a smaller box size in cutoff 18 Å simulations. The average surface area at the end of these simulations lies between ∼ 57 to ∼ 59 Å2, a decrease of 4-6 Å2 from the experimental estimate of 63 Å2.18 When the simulation for γ ) 45 dyn/cm was extended for an additional 2 ns, no change was observed in the area that remained close to 57 Å2 (data not shown). Surface tension values were calculated from NPNAT simulations using Equation 1 (see Methods). Trajectories of surface tension are plotted for a cutoff of 18 Å and PME simulations in Figure 2. Using 50 ps block averages to estimate γ values, PME simulation yields values between 15 and 35 dyn/cm during the last 500 ps of simulation. For the same period, the average γ values calculated for cut off 18 Å simulations is close to 625 dyn/cm, more than 1 order of magnitude higher than that calculated from PME simulations. When the γ ) 625 dyn/cm value is used with 18 Å cutoff, the surface area per lipid remains close to the experimental value (Figure 1a). The average surface area and the tangential pressure calculated from the last 500 ps of simulations are plotted against their respective surface tensions for both the spherical cutoff of 18 Å (Figure 3a) and PME (Figure 3b). The curves were obtained using a polynomial curve fitting procedure.38 Clearly, only a slight increase in the surface area is observed in the 18 Å cutoff when the surface tension was increased from 0 to 85 dyn/cm; the decrease of the tangential pressure is more pronounced. In PME, both the increase in surface area and a decline in tangential pressure were steep when γ was raised from 20 to 45 dyn/cm. The tangential pressure was close to -80 atm for γ ) 25 dyn/cm (PME) and ∼ -1950 atm for γ ) 625 dyn/cm (cutoff of 18 Å). It should be noted that a tangential pressure value of -100 atm for a system of 100 DMPC lipids with the hydrated lipid thickness of 55 Å was derived using monolayer data.19 The γ value of 25 to 30 dyn/cm obtained from the PME simulations for the system of 90 lipids is similar to the surface tension value derived from MD simulations for a system of 72 lipids.2 Bilayer Thickness. The average bilayer and hydrophobic thickness calculated during the last 1 ns production runs are given in Table 2 for all the neat bilayer simulations. It is clear that for cutoff 18 Å simulations the bilayer thickness is about 2 Å higher when the value of γ is between 0 and 85 dyn/ cm. For PME simulations, when γ ) 0 (or 45), the thickness is too high (or low), indicating that the system no longer represents a liquid-crystalline state. The bilayer and hydrophobic thickness

11806 J. Phys. Chem. B, Vol. 108, No. 31, 2004

Figure 2. MD trajectories of surface tension calculated for the neat bilayer simulations using NPNAT ensembles. (A) Cutoff of 18 Å and (B) PME. Both 1 ps (black) and 50 ps (white) block averages are plotted.

evaluated for γ ) 625 dyn/cm with a cutoff of 18 Å and for the two NPNAT simulations are very close to the values obtained from experiments and previous molecular dynamics simulations for DMPC bilayers.5,39-41 Analysis for the PME simulations again confirms that the appropriate γ value for simulations is likely to be between 25 and 30 dyn/cm for the system considered in the present study. Order Parameters. The values of the order parameters for the sn-2 chain are given in Figure 4 for the simulations using γ ) 25 and 30 (PME) and 625 (cutoff of 18 Å). Data from the two NPNAT simulations are also plotted in the same figure. It is clear that the general features observed experimentally are present in the NPNAT simulations and for the analyzed γ values. For PME, γ ) 30 dyn/cm shows the best agreement with the experimentally observed profiles.42

Sankararamakrishnan and Weinstein

Figure 3. Average surface area (b) and the tangential pressure (9) calculated for the last 500 ps are plotted against the surface tension for the neat bilayer simulations using NPNγT ensembles. (A) Spherical cutoff of 18 Å and (B) PME method. A polynomial curve fitting available in the KaleidaGraph package was used to obtain the curves.

Long-range Electrostatics or Van der Waals? The significant differences observed for the same γ value in the two sets of simulations (spherical cutoff of 18 Å and PME) point clearly to the fact that the neglect of long-range interactions lead to serious artifacts. In simulations where long-range electrostatic interactions are explicitly included, a 12 Å cutoff was used to calculate the van der Waals (vdW) interactions. In spherical cutoff simulations, vdW interactions were truncated at 17 Å. Beyond the contact distance, vdW interactions are always negative and this may give rise to a speculation that the increased vdW cutoff in spherical cutoff simulations might be the cause of the shrinkage of area for a given γ value. To investigate the contribution of vdW interactions, we carried out an additional PME simulation with γ ) 45 dyn/cm. In this simulation, a neighbor list to calculate the LJ potential and the real space portion of the Ewald sum was kept to 18 Å using an atom-

Surface Tension Parameterization in Simulations

J. Phys. Chem. B, Vol. 108, No. 31, 2004 11807 TABLE 3: Average (Standard Error) Surface Tension (γ), Tangential (Ptangential), and Normal (Pnormal) Pressure Calculated from the NPNAT Simulationsa cutoff

γ (dyn/cm)

Ptangential (atm)

Pnormal (atm)

12 Å 15 Å 18 Å 21 Å PME

297.4 (3.7) 622.9 (3.9) 625.6 (3.6) 486.2 (4.1) 28.6 (3.3)

-942.2 (7.9) -1997.6 (8.3) -1981.7 (8.1) -1568.7 (8.3) -92.8 (6.9)

-0.9 (10.3) 0.1 (10.1) 12.6 (10.2) -13.6 (11.0) -2.9 (10.0)

a The parameters were calculated for the last 500 ps of 2 ns production runs.

Figure 4. Order parameter profiles for the sn-2 chain from the neat bilayer simulations using NPNγT and NPNAT ensembles. Only those NPNγT simulations in which the γ value keeps the surface area per lipid close to the experimental value are included (625 dyn/cm with a cutoff of 18 Å; 25 and 30 dyn/cm with PME). For NPNAT simulations, data are shown both for a cutoff of 18 Å and PME.

TABLE 2: Average and Root-mean-squared Fluctuations of Bilayer and Hydrophobic Thickness. The Parameters Are Calculated for the Last 1 ns Production Runs for All Neat Bilayer Simulations γ (dyn/cm) cutoff of 18 Å

bilayer thickness (Å)

hydrocarbon thickness (Å)

0 10 20 45 65 85 625 NPNAT γ (dyn/cm) - PME 0 20 25 30 45 NPNAT

37.7 (0.2) 37.8 (0.2) 37.1 (0.3) 37.3 (0.2) 37.2 (0.2) 37.0 (0.2) 35.1 (0.2) 34.8 (0.3)

26.6 (0.2) 27.0 (0.2) 25.8 (0.2) 26.5 (0.2) 26.7 (0.2) 26.2 (0.2) 24.4 (0.2) 24.0 (0.2)

40.8 (0.2) 38.5 (0.5) 37.0 (0.5) 33.0 (0.8) 25.5 (1.9) 35.2 (0.3)

28.3 (0.2) 26.4 (0.4) 24.9 (0.5) 21.2 (0.8) 14.5 (1.8) 23.7 (0.2)

based cutoff and the LJ potential was switched smoothly to zero over the region from 14 to 17 Å, as was done in spherical cutoff simulations. All other parameters were same as that of the other PME simulation with the same γ value. Comparison of MD trajectories of surface-area per lipid (Figure 1b) shows that the two PME simulations with γ ) 45 dyn/cm are almost similar, although the vdW cutoff to calculate the LJ potential is increased by 4 Å. If the neglected long-ranged vdW forces are responsible for the shrinkage of area, then the PME simulation with a longer real-space cutoff (PME-I) would have resulted in a smaller surface area per lipid. This clearly demonstrates that the shrinkage of area observed in the simulation using an 18 Å cutoff with γ ) 45 dyn/cm is NOT due to the increased vdW interactions and is due to the neglect of long-ranged electrostatic interactions.

Clearly, the surface tension of 625 dyn/cm that can maintain the experimental area of the lipid in the cutoff 18 Å simulations is too high, and this must be an artifact of the spherical cutoff method. When PME is used for a system of 90 lipids, a value of γ between 25 and 30 dyn/cm is expected to keep the surface area per lipid, other membrane structural properties, and lateral pressure. Thus the above set of simulations demonstrates that the neglect of long-range electrostatic interactions will have serious shortcomings in lipid simulations. NPNAT Ensembles: Cutoffs 12, 15, 18, and 21 Å, and PME. Influence of force-truncation methods in determining the surface tension values was investigated in simulations with NPNAT ensembles. Surface area (A) in these simulations was maintained close to the experimentally determined value.18 Average (( standard error) values of surface tension, tangential, and normal pressure were analyzed (Table 3). The values from the cutoff simulations were compared with that of PME simulation using the same NPNAT boundary conditions. Calculated surface tension values for different NPNAT ensembles range from ∼ 300 (for 12 Å) to 625 (for 15 and 18 Å) dyn/cm with no clear trend available. This is true also for average tangential and normal pressure values. Even when the cutoff is increased to 21 Å, the surface tension value is still high, and higher than that obtained from a smaller cutoff of 12 Å. However, the calculated surface tension for NPNAT simulation with PME is 28.6 dyn/cm and the tangential pressure is -92.8 atm. For a system of 90 lipids, these values appears very reasonable.2,19 Also from the set of NPNγT simulations with different γ values, when long-ranged electrostatic interactions were explicitly included, a γ value between 25 and 30 dyn/cm is expected to maintain the membrane fluid-phase properties for the system considered here (see above). Dyn A(1-13) in DMPC Bilayers NPNγT Ensembles. The opioid peptide Dyn A(1-13), like other peptide hormones, is believed to bind the membrane first before reaching the target receptors.43,44 We have previously simulated this peptide in the NVE ensemble and studied the membrane-binding properties of this peptide.36 The results of these simulations have shown that the R-helical region of the peptide is stabilized within the bilayers by special interactions that are due to the nature and position of the aromatic and the basic residues in the peptide. The long side chains of the arginine residues interact with all components of the lipids, giving rise to a type of interactions identified as the “snorkel-model”. The aromatic residues Tyr-1 and Phe-4 are key for the positioning of the peptide in the bilayer,36 but behave differently. Thus, the simulations showed Tyr-1 interacting with the lipid headgroups, while the side chain of Phe-4 prefers to be close to the more hydrophobic membrane interior. The N-terminal R-helical structure of the peptide remained stable from residues 4 to 9 throughout the simulations. To evaluate the effect of surface tension on the calculated properties of the membrane-bound peptide, we simulated Dyn A(1-13) in DMPC bilayers using NPNγT ensembles. A fully

11808 J. Phys. Chem. B, Vol. 108, No. 31, 2004

Sankararamakrishnan and Weinstein

Figure 5. Molecular dynamics trajectories of surface area per lipid for the bilayers in which Dyn A(1-13) is imbedded. Data are for four different γ values (0, 45, 65, and 625 dyn/cm) that used a cutoff of 18 Å, and γ ) 25 dyn/cm with PME using NPNγT. Results from simulations using an NPNAT ensemble in which the surface area is close to the experimental value for the neat bilayer are shown for comparison.

flexible simulation cell will allow the lipid bilayer to expand or contract when a peptide is added so that the peptide can assume an optimum orientation inside the bilayers. Here again, an estimate of surface tension value is necessary to use NPNγT ensembles. This becomes more complicated due to the fact that addition of compounds in the membrane is expected to influence the surface tension of the bilayer. As in the neat bilayer simulations, we have carried out simulations using NPNγT ensembles for different γ values. Two schemes (cutoff of 18 Å and PME) to calculate the nonbonded forces were used in the simulations as described above. For comparison purposes, Dyn A(1-13) in the DMPC bilayers was simulated using the NPNAT ensemble with a cutoff of 18 Å, and the surface area/lipid used in NPNAT simulations was taken from the experimental estimate for the neat bilayers.18 The influence of γ on the structure and orientation of the peptide is studied in comparison to the neat membranes. The simulation results of NPNγT ensembles are also compared with simulations using NPNAT and NVE35,36 ensembles. The molecular dynamics trajectories of surface area are plotted for simulations with different γ values in Figure 5. As observed in the neat bilayers, the area decreased rapidly before the first 500 ps simulations for γ ) 0, 45, and 65 dyn/cm when a spherical truncation was used. The average area for the last 500 ps simulations is in the range of 57-59 Å2. This is about 4-6 Å2 less than the experimental value estimated for the neat bilayers.18 For γ ) 625 dyn/cm (a value obtained from NPNAT neat bilayer simulations with a cutoff of 18 Å) the area is maintained close to the experimental estimate for neat bilayers. For the system of 90 DMPC lipids studied here, γ ) 25 dyn/ cm, and the use of PME yields a surface area close to the experi-

mental estimate for neat bilayers. It is also clear that the fluctuation in the surface area is larger in PME simulation. In all these simulations of the peptide-bilayer system, the surface area/ lipid is similar to that observed in neat bilayers in NPNγT ensembles. Dyn A(1-13) Peptide Structure in Bilayers. The structures of Dyn A(1-13) peptides at the end of 2 ns production runs from all the simulations are plotted in Figure 6. For comparison, Figure 6 includes the previously published initial structure, and the structure saved at the end of 8 ns production run from NVE simulation.36 It is clear that the N-terminal helical region of the peptide remained within the bilayers in all the simulations. In NVE simulations, residues 4 to 9 retained an R-helical structure at the end of 8 ns production run.36 However, closer examination reveals that the helix is distorted for some of the γ values used in the simulations. The average values of backbone dihedral angles φ and ψ for residues 4 to 9 are plotted for all simulations in Figure 7. Average φ and ψ values calculated from the first 500 ps of the production runs (Figure 7a and 7b) are compared to those from the last 500 ps of 2 ns productions runs (Figure 7c and 7d) in all the simulations. For γ ) 0, 45, and 65 that resulted in shrinkage of the box size, the φ, and ψ values deviated significantly from the R-helical structure (Figure 7c). The γ values that maintained the lipid surface area (625 for a cutoff of 18 Å and 25 for PME), or the NPNAT simulations in which the area is fixed to the experimental value, maintained the helical structure (Figure 7d). Even a distortion in the C-terminal end of the helix observed in the first 500 ps of the production run in the γ ) 625 dyn/cm with an 18 Å cutoff simulation (Figure 7b) disappeared during the last 500 ps of the 2 ns production run (Figure 7d). It appears that the smaller

Surface Tension Parameterization in Simulations

J. Phys. Chem. B, Vol. 108, No. 31, 2004 11809

Figure 6. Structures of Dyn A(1-13) at the end of 2 ns production runs are plotted for all the peptide-bilayer simulations. For comparison, the initial structure and the structure saved at the end of the 8 ns production run from the previous NVE simulation are also shown. The γ values used in the simulations, and the use of PME, are indicated in each plot as appropriate (all other simulations use the cutoff of 18 Å). For details of NVE simulation, see.35,36 For clarity, lipid molecules are not shown in this figure.

Figure 7. Average backbone dihedral angles φ (blue) and ψ (red) calculated for the first (Figure 7a and 7b) and last (Figure 7c and 7d) 500 ps of production runs. Simulations using NPNγT ensembles with a cutoff of 18 Å are plotted for γ values 0 (b), 45 (9), and 65(2) in the left. Data for NPNγT ensembles with a cutoff of 18 Å and PME are plotted for γ values 625 dyn/cm (b) and 25 dyn/cm (2). For comparison, results from NPNAT ensembles using a cutoff of 18 Å (9) are also shown.

surface area leaves little breathing space for the peptide, and as a result a structural distortion is observed for γ ) 0, 45, and 65 dyn/cm. In previous NVE simulations,35,36 the two aromatic residues Tyr-1 and Phe-4 showed different behavior that reflected the chemical nature of the side chains. The side chains of Tyr-1

and Phe-4 pointed toward the membrane-water interface and the bilayer interior, respectively. When γ ) 0 dyn/cm was used with a cutoff of 18 Å, the smaller box size appears to have restricted the movement of the peptide helix(Figure 5). As a result, instead of moving upward toward the lipid headgroup region, the side chain of the Tyr-1 residue moved in this

11810 J. Phys. Chem. B, Vol. 108, No. 31, 2004 simulation to the opposite side, toward the other half of the bilayer, and attracted water molecules from that side. Discussion The importance of achieving an accurate representation of the phospholipid membrane environment in simulations of cellsurface systems is well established. Among the parameters that determine the quality of the representation in such simulations, the surface tension is being recognized as playing a key role. Thus, while experimentally based surface area/lipid knowledge has improved over the years and the value obtained for this quantity is more reliable from improved experimental techniques and MD simulations, the surface tension value is still debated and estimated indirectly. As demonstrated in this study and in previous work on such systems,2,21 it depends sensitively on various factors. The results show that the influence of forcetruncation methods on the surface tension is serious and cannot be overlooked. To evaluate the specifics of the dependence of simulation results on this key parameter, we have carried out 13 simulations of neat bilayers of 90 DMPC lipids using NPNγT ensembles with a cutoff of 18 Å and PME. Another 5 simulations of the same system were carried out using NPNAT ensembles with PME and four different truncation schemes. The results of these simulations show significant differences between a cutoff of 18 Å and PME. Although 18 Å is generally considered as a sufficiently long cutoff, the calculated surface tension value, lateral pressure (in NPNAT simulation), surface area per lipid, and other membrane structural properties (in NPNγT simulations) are significantly different from the experimental results18,27 or from other computational studies.2,19 Comparison of PME simulations with different γ values and different real space cutoffs clearly demonstrates the role of long-ranged electrostatic interactions in membrane simulations. In PME simulations, the surface tension value that maintains the properties of the bilayer of this size (90 lipids) is expected to be between 25 and 30 dyn/cm. On the other hand, surface tension value of more than an order of magnitude is required for 18 Å cutoff simulations to maintain the membrane fluid-phase properties. Differences between cutoff and PME simulations have been observed in earlier simulations also using the CHARMM force field.45 Simulations of DPPC bilayers have shown that the γ value calculated from a cutoff of 12 Å is about 20 dyn/cm larger than that from Ewald potentials. This is a very small difference compared to the present studies and may be due to the larger surface area considered during NPNAT ensembles. In general, the results obtained from PME simulations for the system of 90 DMPC lipids (with γ ) 25-30 dyn/cm) agreed with the experimental results and earlier simulation studies on 72 DPPC lipids.2 It should also be noted that the present simulations have been carried out for a longer duration of 2 ns with the first 1 nanosecond considered as an equilibration phase. Our simulation results with the CHARMM force field have been compared here to previous simulation studies from Feller and Pastor.2,21 Comparison of simulations using CHARMM and GROMACS force fields and between different GROMACS simulations has been much more difficult due to various factors. First, the CHARMM simulations include all atoms, while GROMACS uses a united atoms model for nonpolar hydrogens. In GROMACS simulations, Chiu et al.19 have used a nonzero surface tension for neat bilayers while Lindahl and Edholm6 advocated an NPT ensemble with zero surface tension. In addition to different van der Waals parameters used in the latter simulations, the 1,4 electrostatic interactions and the van der

Sankararamakrishnan and Weinstein Waals interactions were reduced by factors of 2 and 8, respectively. Two recent papers have studied the influence of truncating electrostatic interactions in DPPC bilayers using the GROMACS force field.46,47 They support the conclusion from the present studies that long-ranged electrostatic interactions are needed to maintain the fluidlike properties of bilayer membranes. Both studies used an NPT ensemble with a surface tension value of 0. A lateral and normal pressure of 1 bar was applied. Patra et al.46 have studied DPPC bilayers over a time scale of 20 ns using three different truncation schemes (18, 20, and 25 Å) and PME. In the simulations that used truncation schemes, the observations included small surface areas per lipid, enhanced ordering of acyl chains, and artificial ordering of lipids in the plane of the membrane; the lipid bilayer was no longer in a truly fluidlike state. Comparison with PME simulations led Patra et al.46 to conclude that truncation of electrostatic properties can have serious consequences for the structural and electrostatic properties of lipid bilayer systems. Studies of Anezo et al.47 reached the similar conclusion that for the simulation of fluidphase bilayers Ewald techniques are more appropriate than straight cutoff methods. In this series of simulations, increasing the electrostatic cutoff resulted in further decrease in the surface area per lipid. The area contraction observed in the straight cutoff techniques is attributed to the small lateral dipole moments of lipid headgroups. Notably, the significant difference between the studies of Patra et al.46 and Anezo et al.47 is in the treatment of surface tension value. While simulation studies with GROMACS parameters used a surface tension value of 0,46,47 we investigated here the influence of truncation on this crucial parameter with a nonzero value in membrane simulations (see Introduction). Previously suggested surface tension values for finite size membrane patch simulations typically lie in the range of 20 to 45 dyn/cm.2,19 A range of nonzero surface tension values was used in the present simulations coupled with different schemes to calculate the nonbonded interactions. It is evident that neglect of long-ranged electrostatic interactions results in serious artifacts for γ values less than 100 dyn/cm. With the truncation method, only a nonzero value of several hundred dyn/cm can maintain the structural and physical properties of the membrane. If longranged electrostatic interactions are taken fully into account, a surface tension value between 25 and 30 dyn/cm is expected to maintain the experimentally observed membrane properties for a system size of 90 lipids. This is also supported by NPNAT simulations with PME in which the calculated γ value is 28.6 dyn/cm. The lateral pressure calculated from these simulations is also close to the value derived from monolayer experimental data.19,27 The present studies also reveal a steep surface tension dependence if PME is used in combination with the NPNγT ensemble in bilayer simulations. Our simulations of Dyn A(1-13) in DMPC bilayers show that surface tension values that result in a smaller area/lipid are likely to affect the structure of the peptide embedded in the membrane. With zero surface tension, the residues of the peptide did not appear to have enough freedom to move due to the shrinkage of the box. All other simulations with nonzero surface tension show the characteristic interactions of aromatic residues. Chiu et al.41 have simulated the gramicidin channel in a DMPC bilayer using an NPNγT ensemble with a nonzero surface tension. Simulations of bacterial membrane proteins, KcSA potassium channel, and FhuA show less drift in the initial structure if PME is used to treat the long-range electrostatic interactions.48 Recent simulation studies on a potassium channel,49 aquaporin,50 and glycophorin A51 have all been performed

Surface Tension Parameterization in Simulations with an NPNAT ensemble and a fixed surface area. Only the Z-direction (bilayer normal) was allowed to vary to maintain the constant normal pressure. It is clear from this work that inclusion of peptides and other compounds complicates the situation, and that the inserts have a direct bearing on the membrane properties such as surface area and surface tension. Results from the simulations of the membrane-bound compounds in NPNγT ensembles have to be treated with caution and preferably compared with simulations using NPNAT ensembles. Acknowledgment. This work was supported by NIH grants P01 DA-11470, DA-12923, and K05 DA-00060. Computational support was provided by the National Science Foundation Terascale Computing System at the Pittsburgh Supercomputing Center. References and Notes (1) Essmann, U.; Berkowitz, M. L. Biophys. J. 1999, 76. (2) Feller, S. E.; Pastor, R. W. J. Chem. Phys. 1999, 111, 1281. (3) Feller, S. E. Curr. Opin. Colloid Interface. Sci. 2000, 5, 217. (4) Feller, S. E.; Mackerell, A. D. J. J. Phys. Chem. B 2000, 104, 7510. (5) Ha Duong, T.; Mehler, E.; Weinstein, H. J. Comput. Phys. 1999, 151, 358. (6) Lindahl, E.; Edholm, O. Biophys. J. 2000, 79, 426. (7) Mashl, R. J.; Scott, H. L.; Subramaniam, S.; Jakobsson, E. Biophys. J. 2001, 81, 3005. (8) Moore, P. B.; Lopez, C. F.; Klein, M. L. Biophys. J. 2001, 81, 2484. (9) Shinoda, W.; Okazaki, S. J. Chem. Phys. 1998, 109, 1517. (10) Brooks, B. R.; Bruccoleri, R. E.; Olafson, B. D.; States, D. J.; Swaminathan, S.; Karplus, M. J. Comput. Chem. 1983, 4, 187. (11) Schlenkrich, M.; Brickmann, J.; Mackerell, A. D., Jr.; Karplus, M. An empirical potential energy function for phospholipids: criteria for parameter optimization and applications. In Biological membranes: A molecular persepctiVe from computation and experiment; Merz, K. M., Roux, B., Eds.; Birkhauser: Boston, 1996; p 31. (12) Berendsen, H. J. C.; van der Spoel, D.; van Drunen, R. Comput. Phys. Commun. 1995, 91, 43. (13) Berger, O.; Edholm, O.; Jahnig, F. Biophys. J. 1997, 72, 2002. (14) Marrink, S. J.; Berger, O.; Tieleman, D. P.; Jahnig, F. Biophys. J. 1998, 74, 931. (15) Chiu, S.-W.; Clark, M. M.; Jakobsson, E.; Subramaniam, S.; Scott, H. L. J. Phys. Chem. B 1999, 103, 6323. (16) Darden, T. A.; York, D.; Pedersen, L. G. J. Chem. Phys. 1993, 98, 10089. (17) Nagle, J. F.; Tristram-Nagle, S. Curr. Opin. Struct. Biol. 2000, 10. (18) Nagle, J. F.; Tristram-Nagle, S. Biochim. Biophys. Acta 2000, 1469, 159.

J. Phys. Chem. B, Vol. 108, No. 31, 2004 11811 (19) Chiu, S.-W.; Clark, M.; Balaji, V.; Subramaniam, S.; Scott, H. L.; Jakobsson, E. Biophys. J. 1995, 69, 1230. (20) Jahnig, F. Biophys. J. 1996, 71, 1348. (21) Feller, S. E.; Pastor, R. W. Biophys. J. 1996, 71, 1350. (22) Marrink, S. J.; Mark, A. E. J. Phys. Chem. B 2001, 105, 6122. (23) Parsegian, V. A. Trans. Faraday Soc. 1966, 62, 848. (24) Nagle, J. F. Annu. ReV. Phys. Chem. 1980, 31, 157. (25) Evans, E. A.; Waugh, R. J. Colloid. Int. Sci. 1977, 60, 286. (26) MacDonald, R. C.; Simon, S. A. Proc. Natl. Acad. Sci. U.S.A. 1987, 84, 4089. (27) Albrecht, O.; Gruler, H.; Sackman, E. J. Physique 1978, 39, 301. (28) Feller, S. E.; Veneable, R. M.; Pastor, R. W. Langmuir 1997, 13, 6555. (29) Mackerell, A. D., Jr.; Bashford, D.; Bellot, 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, T. K.; Mattos, C.; Michnick, S.; Ngo, T.; Nguyen, D. T.; Prodhom, B.; Reiher, W. E.; Roux, B.; Schlenkrich, B.; Smith, J.; Stote, R.; Straub, J.; Watanabe, M.; WiorkiewiczKuczera, J.; Karplus, M. J. Phys. Chem. B 1998, 102, 3586. (30) Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L. J. Chem. Phys. 1983, 79, 926. (31) Hoover, W. G. Phys. ReV. A 1985, 31, 1695. (32) Ryckaert, J. P.; Ciccotti, G.; Berendsen, H. J. C. J. Comput. Phys. 1977, 23, 327. (33) Woolf, T. B.; Roux, B. Proteins: Struct., Funct., Genet. 1996, 24, 92. (34) Feller, S. E.; Zhang, Y. H.; Pastor, R. W.; Brooks, B. R. J. Chem. Phys. 1995, 103, 4613. (35) Sankararamakrishnan, R.; Weinstein, H. Biophys. J. 2000, 79, 2331. (36) Sankararamakrishnan, R.; Weinstein, H. J. Phys. Chem. B 2002, 106, 209. (37) Tessmer, M. R.; Kallick, D. A. Biochemistry 1997, 36, 1971. (38) KaleidaGraph; Synergy software: Reading, PA. (39) Buldt, G.; Gally, H. U.; Seelig, J.; Zaccai, G. J. Mol. Biol. 1979, 134, 673. (40) Shen, L.; Bassolino, D.; Stouch, T. Biophys. J. 1997, 73, 3. (41) Chiu, S.-W.; Subramaniam, S.; Jakobsson, E. Biophys. J. 1999, 76, 1929. (42) Douliez, J.-P.; Leonard, A.; Dufourc, E. J. Biophys. J. 1995, 68, 1727. (43) Schwyzer, R. Biopolymers 1991, 31, 785. (44) Schwyzer, R. Biopolymers (Peptide Science) 1995, 37, 5. (45) Feller, S. E.; Pastor, R. W.; Rojnuckarin, A.; Bogusz, S.; Brooks, B. R. J. Phys. Chem. 1996, 100, 17011. (46) Patra, M.; Karttuen, M.; Hyvonen, M. T.; Falck, E.; Lindqvist, P.; Vattulainen. Biophys. J. 2003, 84, 3636. (47) Anezo, C.; de Vries, A. H.; Holtje, H. D.; Tieleman, D. P.; Marrink, S. J. J. Phys. Chem. B 2003, 107, 9424. (48) Faraldo-Gomez, J. D.; Smith, G. R.; Sansom, M. S. P. Eur. Biophys. J. Biophys. Lett. 2002, 31, 217. (49) Berneche, S.; Roux, B. Nature 2001, 414, 73. (50) de Groot, B. L.; Grubmuller, H. Science 2001, 294, 2353. (51) Petrache, H. I.; Grossfield, A.; MackKenzie, K. R.; Engelman, D. M.; Woolf, T. B. J. Mol. Biol. 2000, 302, 727.