15942
J. Phys. Chem. 1996, 100, 15942-15946
Comparison of the Structures of Dimyristoylphosphatidylcholine in the Presence and Absence of Cholesterol by Molecular Dynamics Simulations Razif R. Gabdoulline,† Garret Vanderkooi, and Chong Zheng* Department of Chemistry, Northern Illinois UniVersity, DeKalb, Illinois 60115 ReceiVed: May 17, 1996; In Final Form: July 31, 1996X
Molecular dynamics simulations were carried out on bilayers consisting of dimyristoylphosphatidylcholine dihydrate (DMPC) and DMPC-cholesterol at a 1:1 mole ratio. The calculations were carried out under NPT conditions (constant number of molecules, pressure and temperature, with volume and bilayer area per molecule free to change). Raising the simulation temperature from 200 to 325 K caused the surface area per DMPC molecule to increase from 39 to 50 Å2 in the pure bilayer; the area per heterodimer of DMPC and cholesterol increased from 79 to 91 Å2 for the same temperature range. The properties of the high-temperature state approximates that of a liquid crystalline array, although the area per molecule in the pure DMPC bilayer remains smaller than the experimentally determined value of 65-69 Å2. The number of gauche bonds, the P-N vector angle distribution, and the internal dipole potential are also reported for the simulated liquid crystalline state in the presence and absence of cholesterol.
Introduction The effects of cholesterol on the physical properties of lipid bilayers and membranes are of considerable biomedical interest and have been studied extensively by a variety of experimental methods, but the interpretation of the results obtained has been hampered by a lack of detailed structural information on these mixed bilayers. X-ray diffraction analysis has been successfully employed to elucidate the three-dimensional structures of several pure phospholipid bilayers, including that of dimyristoylphosphatidylcholine dihydrate (DMPC),1,2 but long-range disorder in cholesterol-containing bilayers has prevented the use of X-ray diffraction for the structural analysis of these systems. For this reason, it is particularly appropriate to apply computational structure determination methods to cholesterol-containing bilayers. One of us3 has recently presented two detailed structural models for representing the short-range order in DMPCcholesterol mixed bilayers at a 1:1 mole ratio, as obtained by energy minimization of a two-dimensional bilayer crystal lattice structures. In both of these low-energy structures, a hydrogen bond is found between the cholesterol hydroxyl and the SN-2 carbonyl oxygen of a neighboring DMPC molecule. The structures were shown to be consistent with the available experimental data. In the present paper, we report the application of molecular dynamic simulations to the DMPC bilayer crystal structure1,4 and to the proposed DMPC-cholesterol crystal structure.3 Calculations were done at three temperatures, viz. 200, 300, and 325 K. The first of these is well below the thermal phase transition temperature (296 K), the second is slightly above the transition temperature, and the third is well above the transition temperature. These calculations were carried out under the NPT conditions, i.e., a fixed number of atoms, constant pressure and temperature, but with the volume and hence the surface area per molecule being free to vary. This approach differs from, and is more realistic than, all previous calculations on lipid † Permanent address: Institute of Mathematical Problems of Biology of the Russian Academy of Sciences, Pushchino, Moscow Region, 142292 Russia; present address: EMBL, Meyerhofstr. 1, 69012 Heidelberg, Germany. * To whom correspondence should be addressed. X Abstract published in AdVance ACS Abstracts, September 1, 1996.
S0022-3654(96)01445-1 CCC: $12.00
bilayers of which we are aware, which held the molecular area fixed at a value derived experimentally for liquid crystalline bilayers. In the present work, all simulations were started from the crystal structure, with lattice expansion occurring dynamically as the temperature was raised. Several papers have appeared in the past few years on the application of molecular dynamic methods to the study of lipid bilayer properties. Pastor and co-workers5,6 have used a mean field approach to model the Brownian dynamics of the interior of a lipid bilayer in the liquid crystalline state. Scott7 used a Monte Carlo approach to study the effects of cholesterol on the acyl chain dynamics in a simulated array, at constant area, and Edhom and Nyberg8 also studied the effects of cholesterol on the acyl chains of phosphatides using molecular dynamics. Both of the latter studies emphasized hydrocarbon chain dynamics, with the details of the group interactions being omitted. The dynamical properties of bilayers consisting of dilauroylphosphatidylethanolamine have been studied by Damodaran and co-workers9 using molecular dynamics. This calculation used 48 lipid molecules and 553 water molecules, with periodic boundary conditions. The lipids were initially placed in the conformation obtained by X-ray diffraction,10 but with the crystal structures being artificially dilated to give the area per molecule (51 Å2) expected for the liquid crystalline state, on the basis of experimental data. Stouch11 has reported molecular dynamics calculations on lipid bilayers of DMPC at 320 K, using periodic boundary conditions and a bilayer segment consisting of 36 DMPC molecules solvated with 481 water molecules. The bilayer segment had lateral dimensions of 34.5 × 34.5 Å, with a fixed surface area of 66 Å2 per DMPC. The latter value was based upon experimental determinations of area per molecule in the liquid crystalline state. Materials and Methods Atomic coordinates for the starting bilayer conformations were taken from the results of the energy minimization studies of one of us.3,4 For DMPC bilayers, the coordinates of structure II were employed,4 whereas for DMPC-cholesterol (1:1 mole ratio), the coordinates of structure B were used.3 For the molecular dynamics simulations, a segment of 36 DMPC lipid molecules surrounded by 621 water molecules was used as one © 1996 American Chemical Society
Structures of Dimyristoylphosphatidylcholine cell of periodic dynamics. The total number of atoms is 3519 and the cell dimensions are 29.5 × 31.0 × 80.0 Å3 at the starting configuration. For the mixed bilayer, 16 DMPC and 16 cholesterol molecules together with 664 water molecules were used in the simulation. The total number of atoms is 3192 and the cell dimensions are 35.3 × 17.8 × 78.0 Å3 at the start of the simulation. Water molecules were represented using the SPC/E (extended single point charge) model.12 The GROMOS molecular dynamics simulation program13 was used to simulate the dynamics of the water-bilayer system. An adjustment of the GROMOS force field parameters was undertaken to match the previous energy-minimized structure of the bilayers.3,4 In particular, the Lennard-Jones radii of the carbon groups CH2 of the acyl chains were increased by 5.2% of the GROMOS standard values. NPT conditions were applied, i.e., the number of atoms (N) in the simulation box was kept constant, pressure was set equal to 1 atm, and the temperature was fixed at a specified value. Accordingly, the size of the simulation box is free to deviate from the starting size, which is determined by the gel-state conformation of the lipid bilayer. The temperature was kept constant by coupling to a thermal bath with a relaxation time of 0.05 ps for the lipid molecules and 0.1 ps for the water. The dynamics simulations were carried out by the time steps of 4 fs for DMPC and 2 fs for the mixed lipid bilayer. The total simulation time for each of bilayers was 0.6 ns and is composed of the following: first, the minimized bilayer was subject to a 0.1 ns simulation at 200 K; within a subsequent 0.3 ns simulation time, apparent equilibration of the bilayers at 300 K was achieved; and an additional 0.2 ns simulation was carried at 325 K. The quantities described below were calculated using the last 100 ps of the simulations at different temperatures. Average values were computed from the trajectories for several physical parameters, including bilayer area per molecule as a function of temperature, fraction of gauche bonds as a function of temperature, the P-N vector angle distribution, and the transbilayer electrostatic potential. The potential profile across bilayers was calculated using the time-average charge density F(z) along the coordinate z perpendicular to the bilayer, as described by, for example, Chiu and co-workers.14 The potential dependence V(z) can be obtained using integration on z:15
V(z) )
1 z ∫ dz′∫-∞z′ dz′′ F(z′′) -∞
where is the permittivity of the medium. Simulations were carried out in a periodic box without any fixed center; therefore we have recentered the simulation box in accordance with the current position of the bilayer molecules. Results Figure 1 shows a cross section of a DMPC bilayer in the energy-minimized gel state (Figure 1a) and after molecular dynamic equilibration at 300 K. The acyl chains show tilt in the gel state but appear to be perpendicular to the bilayer plane on the average in the high-temperature liquid crystalline state. It is visually apparent from comparison of Figure 1a,b (which are drawn to the same scale) that the thickness of the bilayer decreases upon raising the temperature and hence also that the area per molecule increases. The volume of the simulation box, on the contrary, was constant up to a few percent. The molecular structure shown in Figure 1b was obtained as a result of molecular dynamic simulation applied to the regular structure shown in Figure 1a; only the temperature was changed and no artificial lattice dilation was applied to the gel-state structure to obtain the disordered liquid crystalline structure.
J. Phys. Chem., Vol. 100, No. 39, 1996 15943 Figures 2a,b shows the DMPC-cholesterol bilayer in its energy-minimized gel state (Figure 2a) and in the hightemperature state, following molecular dynamic equilibration at 300 K. Here also the bilayer thickness is seen to decrease as a result of raising the temperature, with a concomitant increase in area per molecule. The locations of the cholesterol molecules in the bilayer became disordered, as expected, but remain aligned in the bilayer plane along with the DMPC molecules, with the long axes of both molecules being on the average parallel to the bilayer normal. Area per Molecule. The bilayer area per DMPC molecule is presented as a scatter diagram in Figure 3 as a function of the fraction of gauche bonds in the acyl chains of the DMPC molecules. Results are given for 200, 300, and 325 K. The structure at 200 K approximates that of the energy-minimized state. The area per molecule in the energy-minimized lattice is 38.6 Å2 and its fraction of gauche bonds is 0.04.4 The latter value is nonzero because the B4 bond in the crystal structure is gauche. Both the fraction of gauche bonds and the area per molecule increase as the temperature is raised, and the range of scatter increases. It may be noted that both of these parameters also increase upon going from 300 K (which is slightly above the DMPC transition temperature) to 325 K (considerably above the transition temperature). The centers of the distribution correspond to approximately 47 and 50 Å2, and 0.19 and 0.22 fraction gauche, at 300 and 325 K, respectively. It is evident that the equilibrated area per molecule at 325 K is less than that observed experimentally for the liquid crystalline state. Lis and co-workers16 reported a value of 71.2 Å2 and Albon obtained 65-69 Å2; both of these values are for dipalmitoylphosphatidylcholine in excess water above the phase transition temperature.17 The number of gauche bonds found here is also at the low end of the range of 0.2-0.3 determined by Mendelsohn and co-workers18 for the liquid crystalline state using infrared spectroscopy. The reason for this apparent discrepancy between the computed and experimental results will not be known until further calculations are carried out using a larger number of temperatures. Judging from the data at hand, the problem does not appear to be one of insufficient equilibration time. Having only three temperature data points as yet, it is not possible to say whether, or how faithfully, our computational approach will reproduce the experimental phase transition. Considering the relatively small size of the periodic unit cell, which is the cooperative unit for phase transition purposes, it may be that no discrete phase transition can be reported, but rather that the change in area will change smoothly with temperature over a fairly broad temperature range. Should that be the case, we would not expect the computed molecular area to approach that of the experimental value until the temperature has been raised to a considerably higher value. At this point, it is perhaps sufficient to say that it has proven to be possible to carry out the bilayer calculation using area and volume as free variables and to observe that the area does indeed spontaneously increase in the simulations as the temperature is raised. Figure 4 gives the area per heterodimer of DMPC and cholesterol as a function of DMPC acyl chain gauche bonds, for the same three temperatures as were reported in Figure 3 for DMPC alone. At 200 K, the area per dimer is similar to that found in the energy-minimized structure, which is 78.6 Å2.3 When the temperature is raised to 300 or 325 K, the number of gauche bonds also increases, but only to 0.12 as compared to 0.19 or 0.22 in the absence of cholesterol. This reflects the chain ordering tendency of cholesterol on DMPC which has been observed and discussed by previous authors.7,8 There is
15944 J. Phys. Chem., Vol. 100, No. 39, 1996
Gabdoulline et al.
Figure 1. Comparison of the cross-sectional appearance of a DMPC bilayer in the crystalline state (a) and in the fluid state (b), following molecular dynamic equilibration at 300 K, as described in the text.
no further increase in gauche bonds on going from 300 to 325 K, as occurs for pure DMPC. The area per dimer increases to 87.5 or 91.5 Å2 at 300 and 325 K, respectively. The actual change in area upon going from 200 to 300 or 325 K is about the same for DMPC and DMPC-cholesterol, but since the area per heterodimer is twice that of a DMPC monomer, the percentage change in area is only one-half as great in the mixed bilayer. This implies that the cholesterol contributes very little to the area change. P-N Vector Angle. The polar headgroup of phosphatides lies roughly parallel to the bilayer plane in most lipid crystal structures which have been determined.2 Now it is of interest to determine the degree of motional freedom of the headgroup which may be expected in the liquid crystalline state. The P-N vector is a convenient parameter to use for this purpose, since it lies approximately in the same direction as the headgroup dipole. Figure 5 gives the distribution of P-N vector orientations obtained for DMPC and DMPC-cholesterol. The distribution is expressed as a fraction of cos(φ), where φ is the angle between the outwardly directed bilayer normal and the P-N vector. The maximum of the distribution lies at φ ) 82° for DMPC, which implies that the P-N vector makes an angle of
8° away from the bilayer plane. This may be compared with the values of 9.4° and 26.7° in the two nonequivalent DMPC molecules in the energy-minimized crystal structure.4 The range observed for φ in Figure 5 runs from 37° to 126°, however, showing that for part of the time the N atom lies closer than the P atom to the bilayer midplane (negative values of cos(φ) in Figure 5). A similar result has been reported by Stouch.11 The maximum of the distribution function of the P-N vector angle in the presence of cholesterol is shifted to a slightly higher value of cos(φ) as compared to DMPC alone and corresponds to an angle of 13° with respect to the bilayer plane. The major difference between the two is, however, that there is a greater weighting of values at negative cos(φ); this implies that with cholesterol present there is a somewhat greater tendency for the headgroup of DMPC to be pointing down toward the bilayer than is the case in the absence of cholesterol. The cholesterol acts as spacers between the DMPC molecules and provides a greater latitude of freedom in which the DMPC headgroups can move. Internal Dipole Potential. The internal potential of a lipid bilayer is the difference between the potential at the bilayer midplane and the potential in the surrounding aqueous environ-
Structures of Dimyristoylphosphatidylcholine
J. Phys. Chem., Vol. 100, No. 39, 1996 15945
Figure 2. Cross-sectional views of a 1:1 DMPC-cholesterol bilayer, in the crystalline state (a) and in the fluid state (b), obtained by molecular dynamic equilibration at 300 K.
Figure 3. Relationship between the area per DMPC molecule and the fraction of acyl bonds in the gauche configuration, following equilibration at three temperatures. The plotted symbols represent discrete points along the molecular dynamic trajectory. For comparison, the area per molecule in the energy-minimized structure is 38.6 Å2 and the fraction of gauche bonds is 0.04.
ment. Experimental measurements on the differential permeability of otherwise equivalent positive and negative hydropho-
Figure 4. Plot of the area per heterodimer of cholesterol and DMPC (1:1 mole ratio) as a function of the fraction of gauche bonds in the acyl chains of the DMPC molecules. In the energy-minimized structure, the area per heterodimer is 78.6 Å2 and the fraction of gauche bonds is 0.04.
bic ions across phosphatidylcholine bilayers has led to the conclusion that the potential at the bilayer midplane is positive relative to water, with a value in the range 100-200 mV (see
15946 J. Phys. Chem., Vol. 100, No. 39, 1996
Gabdoulline et al. difference between the bilayer midplane and water is about 180 mV, which is somewhat larger than the value of 138 mV previously reported for the gel-state bilayer.21 The dashed line in Figure 6 gives the potential for the mixed DMPC-cholesterol bilayer. Surprisingly, the cholesterol significantly decreases the magnitude of the potential difference to about 60 mV. We are not aware of any experimental data with which this value may be compared, but it would indicate that cholesterol-containing bilayers are less discriminatory between positive and negative ions. This point deserves experimental verification. Discussion
Figure 5. Distribution of polar group orientations in DMPC and DMPC-cholesterol bilayers equilibrated at 300 K. The headgroup orientation is expressed as the angle φ between the outwardly directed bilayer normal and the vector from the P to the N atoms of the phosphocholine group. The solid line is for DMPC, and the dashed line is for DMPC-cholesterol.
The results presented here show the utility of using the molecular dynamic approach under NPT conditions for studying the properties of pure and mixed lipid bilayers. By this means, both gel-state and liquid crystalline bilayers can be modeled in a self-consistent manner. A stringent test of the method and the parameters employed will be to attempt to reproduce the thermal phase transition by raising the simulation temperature in small increments. It will be interesting to see if doing so actually displays a cooperative melting transition or simply a smoothed curve running through what would experimentally be the phase change discontinuity. Finally, it will be important to calculate the surface tension and compare the values with the experimentally measured ones. However, one difficulty would be the accurate calculation of the surface area by simulation methods. We hope to address these problems in future studies. Acknowledgment. We thank the National Research Council for a CAST grant that made R.R.G.’s visit at Northern Illinois University possible. This work was also supported in part by NSF through the Presidential Young Investigator grant (CHE91-57717) awarded to C.Z. The computations were carried out on a Silicon Graphics workstation at NIU.
Figure 6. Comparison of the electrostatic potential for DMPC and DMPC-cholesterol as a function of position across the bilayer, so calculated following molecular dynamic equilibration at 300 K. The curves have been arbitrarily aligned with the zero point corresponding to the potential at the bilayer midplane.
refs 18 and 19 for a review of the experimental background on this subject). It has been of interest to attempt to determine the molecular origin of this potential. In our first set of calculations,19 we simply calculated the potential for DMPC and dilauroylphosphatidylethanolamine, using the crystal structure coordinates. We concluded from these calculations that the dipole potential created by the charge distribution in the lipid molecules themselves could not account for the experimentally observed positive value of the internal potential. We inferred from these results that the water of hydration must be responsible for the internal potential. The same conclusion was also arrived at by Gawrisch and co-workers20 on the basis of experimental work. Thereafter we carried out further calculations in which surrounding water was included in the calculations, being modeled as in this present paper using molecular dynamics, but with the lipid molecules being held fixed in their energy-minimized structures.21 From that work we realized that the internal potential of bilayers only appears to be positive because the potential of bulk water computes as being negatiVe relative to a vacuum or a hydrocarbon medium. The key to the problem is hence the electrical properties of water. Now we compute the internal potential of DMPC and DMPC-cholesterol for the simulated liquid crystalline state. The results are given in Figure 6. Since we are only interested in potential differences, the plots are given with the potential at the bilayer midplane arbitrarily set to zero. The potential
References and Notes (1) Pearson, R. H.; Pascher, I. Nature (London) 1979, 281, 499-501. (2) Hauser, H.; Pasher, I.; Pearson, R. H.; Sudell, S. Biochim. Biophys. Acta 1981, 650, 21-51. (3) Vanderkooi, G. Biophys. J. 1994, 66, 1457-1468. (4) Vanderkooi, G. Biochemistry 1991, 30, 10760-10768. (5) Paster, R. W.; Verable, R. M.; Karplus, M. J. Chem. Phys. 1988, 89, 1112-1127. (6) Paster, R. W.; Verable, R. M.; Karplus, M. Proc. Natl. Acad. Sci. U.S.A. 1991, 88, 892-896. (7) Scott, H. L. Biophys. J. 1991, 59, 445-455. (8) Edholm, O.; Nyberg, A. M. Biophys. J. 1992, 63, 1081-1089. (9) Damodaran, K. V.; Merz, K. M., Jr.; Braber, B. P. Biochemistry 1992, 31, 7656-7664. (10) Hitchcock, P. B.; Mason, R.; Thomas, K. M.; Shipley, G. G. Proc. Natl. Acad. Sci. U.S.A. 1974, 71, 3036-3040. (11) Stouch, T. R. Mol. Simul. 1993, 10, 335-362. (12) Berendsen, H. J. C.; Grigera, J. R.; Straatsma, T. P. J. Phys. Chem. 1987, 91, 6269-6271. (13) van Gunsteren, W. F.; Berendsen, H. J. C.; Hermans, J.; Hol, W. G. J.; Postma, J. P. M. Proc. Natl. Acad. Sci. U.S.A. 1983, 80, 4315-4319. (14) Chiu, S.-W.; Gulukota, K.; Jacobsson, E. In Membrane Proteins: Structures, Interactions and Models; Pullman, A., et al., Eds.; Kluwer Academic Publishers: The Netherlands, 1992. (15) Wilson, M. A.; Pohorille, A.; Pratt, L. R. J. Chem. Phys. 1988, 88, 3281. (16) Lis, L. J.; McAlister, M.; Fuller, N.; Rand, R. P.; Parsegian, V. A. Biophys. J. 1982, 37, 657-666. (17) Albon, N. J. Phys. Chem. 1985, 89, 3147-3151. (18) Mendelsohn, R.; Davies, M. A.; Brauner, J. W.; Schuster, H. F.; Dluhy, R. A. Biochemistry 1989, 28, 8934-8939. (19) Zheng, C.; Vanderkooi, G. Biophys. J. 1992, 63, 935-941. (20) Gawrisch, K.; Ruston, D.; Zimmerberg, J.; Parsegian, V. A.; Rand, R. P.; Fuller, N. Biophys. J. 1992, 63, 1213-1223. (21) Gabdoulline, R. R.; Zheng, C.; Vanderkooi, G. Unpublished results.
JP9614456
