Langmuir 1998, 14, 5255-5266
5255
Structural Characterization of Self-Assembled Lipid Monolayers by NπT Simulation Andrew W. Mauk,† Elliot L. Chaikof,†,‡ and Peter J. Ludovice*,† School of Chemical Engineering, Georgia Institute of Technology, Atlanta, Georgia 30332-0100, and Department of Surgery, Emory School of Medicine, Atlanta, Georgia 30322 Received September 2, 1997. In Final Form: June 19, 1998 Molecular dynamics simulations in which the number of molecules (N), the two-dimensional pressure (π), and the temperature (T) are held fixed (NπT) were performed on phospholipid monolayers in twodimensional periodic boundary conditions. Four different phospholipids were simulated: dipalmitoylphosphatidylcholine (DPPC), disteroylphosphatidylethanolamine (DSPE), 1-palmitoyl-2-oleoylphosphatidylcholine (POPC), and 1-steroyl-2-arachidonylphosphatidylethanolamine (SAPE) in both low- and highdensity states. The simulations reproduced pressure-area isotherms measured on a Langmuir-Blodgett trough. The alkyl chain alignment in the low-density simulations also compared well to results of H2 NMR investigations from the literature. Simulations in which the area of the lipid monolayer was fixed (NAT) were also performed. Simulations in the NπT ensemble were superior to the NAT simulations. Inaccurate equilibrium pressures and alkyl chain order resulted from the NAT simulations. Experimental results indicated the presence of a phase transition in the DPPC isotherm. The simulations indicate this phase transition results from a desolvation of the lipid headgroups as the monolayer is pressurized. However, these simulations have not completely equilibrated. This makes any conclusions drawn from these or the many similar simulations in the literature suspect until independent validation is obtained.
Introduction Computer simulations of lipid-water systems are of great interest because of the simulation’s ability to obtain a molecular level insight into the structure and dynamics of lipid assemblies. This insight is useful in understanding both cell membranes and synthetic liposomes used in applications such as drug delivery. Additionally, accurate simulation models are required to accurately simulate transmembrane protein channels and receptors. Numerous simulation protocols have been used to simulate lipid assemblies. Due to their computational efficiency, Monte Carlo simulations were initially used. Scott used the Monte Carlo technique to simulate membrane phenomena such as hydrocarbon chain orientation,1 the effects of cholesterol and membrane proteins in the lipid bilayer,2-4 surface tension of the lipid-water interface,5 and density profiles of water interacting with phospholipid headgroups.6 Other studies included using a lattice model to simulate lipid membranes,7,8 lateral diffusion in lipid assemblies,9,10 and many other Monte Carlo simulations of lipid structure.11-18 * To whom correspondence should be addressed. † Georgia Institute of Technology. ‡ Emory School of Medicine. (1) Scott, H. L. Biochim. Biophys. Acta 1977, 469, 264-271. (2) Scott, H. L.; Kalaskar. S. Biochemistry 1989, 28, 3687-3691. (3) Scott, H. L.; Cherng, S. L. Biochim. Biophys. Acta 1978, 510, 209-215. (4) Scott, H. L. Biochemistry 1986, 25, 6122-6126. (5) Scott, H. L.; Lee, C. Y. J. Chem. Phys. 1980, 73, 5851-5353. (6) Scott, H. L. Biophys. J. 1991, 59, 445-455. (7) Dill, K. A.; Flory, P. J. Proc. Natl. Acad. Sci. U.S.A. 1981, 78, 678-680. (8) Pink, D. A. Can. J. Biochem. Cell Biol. 1984, 62, 760-777. (9) Saxton, M. J. Biophys. J. 1992, 61, 119-128. (10) Galle, J.; Volke, F. Biophys. Chem. 1995, 54, 109-117. (11) Pastor, R. W.; Venable, R. M.; Karplus, M. Proc. Natl. Acad. Sci U.S.A. 1991, 88, 892-896. (12) Scheringer, M.; Hilfer, R.; Binder, K. J. Chem. Phys. 1992, 96, 2269-2276. (13) Larson, R. G. J. Chem. Phys. 1996, 11, 7904-7917. (14) Kramer, D.; Ben-Shaul, A.; Chen, Z.-Y.; Gelbart, W. M. J. Chem. Phys. 1992, 96, 2236-2252.
With the recent increase in the computational efficiency of computers, molecular dynamics (MD) studies of lipids have become more prevalent. MD simulations can access the dynamic behavior of these systems, and they do not require the choice of a set of Monte Carlo moves to run the simulation. Initial MD studies focused only on the hydrocarbon region of lipid systems.19-21 Later studies have expanded to include much more detailed and complex systems. Damodaran and co-workers simulated 48 molecules of dilauroylphosphatidylethanolamine (DLPE) with 553 water molecules, making it the first large-scale phospholipid-water system studied.22 Schulten and coworkers constructed very large simulations consisting of up to 200 phospholipid molecules with thousands of water molecules.23,24 Jonsson et al. and Watanabe et al. performed explicit molecular dynamics simulations of a sodium octanoate micelle in aqueous solution.25,26 Other recent simulations include those by Berkowitz and Raghavan,27 Huang and co-workers,28 Stouch and coworkers,29 and Tu and co-workers.30,31 (15) Aloisi, G.; Foresti, M. L.; Guidelli, R.; Barnes, P. J. Chem. Phys. 1989, 91, 5592-5596. (16) Kroll, D. M.; Gompper, G. Science 1991, 255, 968-971. (17) Pastor, R. W.; Karplus, M. J. Chem. Phys. 1989, 1, 211. (18) Englemann, A. R.; Llanos, C. M.; Nyholm P. G.; Tapia, O.; Pascher, I. J. Mol. Struct. (THEOCHEM) 1987, 151, 81-102. (19) van der Ploeg, P.; Berendsen, H. J. C. J. Chem. Phys. 1982, 76, 3271-3276. (20) Edholm, O.; Berendsen, H. J. C.; van der Ploeg, P. Mol. Phys. 1983, 48, 379-388. (21) Applegate, K. R.; Glomset, J. A. J. Lipid Res. 1986, 27, 658680. (22) Damodaran, K. V.; Merz, K. M.; Gaber, B. P. Biochemistry 1992, 31, 1. (23) Heller, H.; Schaefer, M.; Schulten, K. J. Phys. Chem. 1993, 97, 8343-8360. (24) Zhou, F.; Schulten, K. J. Phys. Chem. 1995, 99, 2194-2207. (25) Jonsson, B.; Edholm, O.; Teleman, O. J. Chem. Phys. 1986, 85, 2259-2271. (26) Watanabe, K.; Ferrario, M.; Klein, M. L. J. Phys. Chem. 1988, 92, 819-821. (27) Berkowitz, M. L.; Raghavan, K. Langmuir 1991, 7, 1042-1044. (28) Huang, P.; Perez, J. J.; Loew, G. H. J. Biomol. Struct. Dynam. 1994, 11, 927-956.
S0743-7463(97)00987-6 CCC: $15.00 © 1998 American Chemical Society Published on Web 08/08/1998
5256 Langmuir, Vol. 14, No. 18, 1998
Mauk et al.
While the aforementioned simulations concentrated mainly on lipid systems in a bilayer configuration, relatively few Monte Carlo32-34 or molecular dynamics35-39 studies have been performed on lipid systems in monolayer form to reproduce the behavior of lipid films assembled at a liquid-air interface known as Langmuir monolayers. However, recent advances in the possible applications of Langmuir monolayers in fields such as microelectronics and medicine have made such molecular dynamic simulations very worthwhile in order to study the thermodynamics and packing of said monolayers. Furthermore, the packing behavior of phospholipid monolayers is directly measured on a Langmuir-Blodgett trough as opposed to phospholipid bilayers in a vesicle configuration. This study simulated four different types of phospholipids in monolayer form in contact with water at two different densities. The four phospholipids were chosen because of our interest in characterizing a supported bilayer for biomedical applications; the four phospholipids are the most common in arterial endothelial cell membranes. Pressure-area isotherms of all four phospholipids were determined experimentally on a Langmuir-Blodgett trough in order to compare the experimentally observed packing to that predicted by simulation. Experimental Section Pressure-Area Isotherm Formation. The four phospholipids used in isotherm formation were (Figure 1) dipalmitoylphosphatidylcholine (DPPC), 1-palmitoyl-2-oleoylphosphatidylcholine (POPC), disteroylphosphatidylethanolamine (DSPE), and 1-steroyl-2-arachidonylphosphatidylethanolamine (SAPE). All phospholipids were obtained from Avanti Polar Lipids Inc. and were used without further purification. Solutions of all four types of lipids were created for deposition onto the water surface of the Langmuir-Blodgett trough. Both phosphatidylcholine solutions were prepared with pure chloroform, and both phosphatidylethanolamine solutions were prepared with a chloroform/ methanol mixture in a 2:1 ratio. The solutions were all created at a concentration of 0.5 mg/mL. A Langmuir-Blodgett trough, constructed by Nima Technology Ltd.,40 was used for isotherm measurements in all cases. The trough was filled with pure, double distilled water obtained from a Millipore Modulab Analytical system (Continental Water Systems Corp.). The isotherms were collected at 21 °C and a compression barrier speed of 5 cm2/min. Compression was continued until collapse of the monolayer occurred, evidenced by a sharp discontinuity in pressure. Computational Details. All models were constructed and all simulations were run using the programs QUANTA version 4.042 and CERIUS2 version 3.5,42 both developed by Molecular Simulations Inc. The computers used in the simulations included a Silicon Graphics Indigo2 XZ Graphics workstation and a Silicon (29) Stouch, T. R.; Alper, H. E.; Bassolino-Klimas, D. Supercomput. Appl. High Performance Comput. 1994, 8, 6-23. (30) Tu, K.; Tobias, D. J.; Blasie, J. K.; Klein, M. L. Biophys. J. 1996, 70, 595-608. (31) Tu, K.; Tobias, D. J.; Klein, M. L. Biophys. J. 1995, 69, 25582562. (32) Harris, J.; Rice, S. A. J. Chem. Phys. 1988, 88, 1298. (33) Millik, M.; Kolinski, A.; Skolnick, J. J. Chem. Phys. 1990, 93, 4440. (34) Makovsky, N. N. Mol. Phys. 1991, 72, 235. (35) Alper, H. E.; Bassolino, D.; Stouch, T. R. J. Phys. Chem. 1993, 98, 9798-9807. (36) Alper, H. E.; Bassolino-Klimas, D.; Stouch, T. R. J. Phys. Chem. 1993, 99, 5547-5559. (37) Karaborni, S.; Toxvaerd, S. J. Chem. Phys. 1992, 96, 55055515. (38) Northrup, S. H. J. Phys. Chem. 1984, 88, 3441-3446. (39) Northrup, S. H.; Curvin, M. S. J. Phys. Chem. 1985, 89, 47074713. (40) Grunfeld, F. Rev. Sci. Instrum. 1993, 64, 548-555. (41) Voet, D.; Voet, J. G. Biochemistry; John-Wiley & Sons: New York, 1990. (42) Molecular Simulations Inc., San Diego, CA, 1994.
Figure 1. Chemical structures of (a) DPPC, (b) POPC, (c) DSPE, and (d) SAPE. The upper alkyl chain is the sn-1 chain and almost consists of completely saturated carbons. The lower chain is the sn-2 chain, and it is the chain that contains unsaturated carbons, if any. The sn designation stands for stereospecific numbering, which assigns the number one position to the chain occupying the pro-S position of a prochiral center.41 Graphics Power Challenge Series M workstation. All initial monomeric structures were derived from the X-ray crystal structure of DLPE.43,44 Cis double bonds and extra CH2 units were added as necessary to achieve the desired hydrocarbon tails in the all trans geometry of the DLPE for the phospholipid in question. Construction of the phosphatidylcholine group was accomplished by replacing the hydrogens in the phosphatidylethanolamine amine group with methyl groups. All methylene units on the hydrocarbon chains were treated as unified atoms, where the methylene group was represented by a single effective atom. The new phospholipid structure was then reproduced to form a monolayer of 25 phospholipids in a 5-by-5 lattice with P1 symmetry, making the computational model part of a twodimensional periodic phospholipid monolayer-water system as opposed to an isolated complex. For each of the four phospholipids in question, a low-density and a high-density phase were constructed, determined by their lateral packing densities. A slab of water molecules approximately 16 Å in depth, whose potential energy was already minimized, was placed adjacent to each hydrophilic headgroup region. These slabs were generated by randomly packing water at a density of 1 g/cm3 in a cell of the initial dimensions of the two-dimensional lipid simulations with a thickness of 16 Å. The energy of this periodic cell of water was minimized before joining it to the hydrophilic side of the initial lipid cell. The TIP3P model was used for the water.45 The resulting periodic unit cell was then extended to 250 Å in the direction normal to the surface to eliminate potential energy interactions in the normal direction. Given that the potential was truncated, this prevented the water molecules on one side of the periodic cell from interacting with the alkane chains on the other side of the crystal. This effectively converted the simulation to a two-dimensional periodic system. The parameters used for the lipids in all cases were obtained from the CHARMM22.0 force field with the potential energy cutoff using (43) Hitchcock, P. B.; Mason, R.; Thomas, K. M.; Shipley, G. G. Proc. Natl. Acad. Sci. U.S.A. 1974, 71, 3036-3040. (44) Elder, M.; Hitchcock, P.; Mason, R.; Shipley, G. G. Proc. R. Soc. London 1977, 354, 157-170. (45) Jorgensen, W. J. Am. Chem. Soc. 1981, 103, 335-342.
Structures of Self-Assembled Lipid Monolayers
Langmuir, Vol. 14, No. 18, 1998 5257
Table 1. Atom Type Charges Obtained from CHARMM22.0 Templates atom type
chargea
description
C CH1E CH2E CH3E CT CUA1 HA
0.700 0.275 0.026 0.027 0.125 -0.200 0.128
HC HT NT O OAC OC OS1 OS2 OW PO4
0.350 0.417 -0.300 -0.550 -0.550 -0.570 -0.670 -0.400 -0.834 1.110
carbonyl carbon carbon containing one virtual hydrogen carbon containing two virtual hydrogens carbon containing three virtual hydrogens tetrahedral carbon connected to four other atoms carbon that is part of a double carbon-carbon bond in alkyl chains hydrogen that is part of methyl group on headgroup of phosphatidylcholines or hydrogen connected to a double bonded carbon in alkyl chains hydrogen connected to nitrogen in headgroup of a phosphatidylethanolamine hydrogen attached to oxygen in a water molecule nitrogen in headgroup of phospholipids one of two oxygens connected to atom type PO4 carbonyl oxygen one of two oxygens connected to atom type PO4 ester oxygen in phospholipids ester oxygen in phospholipids oxygen in water molecules phosphate in headgroup of phospholipids
a
In fraction of an electron charge.
a spline function between 11.0 and 14.0 Å.46,47 This potential cutoff is slightly larger than the commonly used spline range of 8-10 Å based on recent work that suggests using cutoffs in the 10-15 Å range.48,49 However, more recent work suggests that the orientational structure of water adjacent to self-assembled lipids is affected significantly for effective cutoffs below 30 Å.50 CHARMM22.0 is a classical force field using quadratic force constants and was specifically created for biomolecules. The force field consists of bond stretch, bond angle, torsion angle, and inversion terms. Charges for all atoms were obtained from the CHARMM22.0 templates and the charges on the carbon atoms smoothed to produce a net zero molecular charge using the algorithm in CHARMM. The resulting charges are listed in Table 1. Because the smoothing algorithm treats all atoms of the same atom type identically, this results in a small positive charge on the united atoms on the alkyl chains (atom types CH2E and CH3E). This is somewhat unrealistic, as this positive charge would normally be distributed in a nonuniform fashion with most of it localized near the headgroup. Given the relatively small size of the charge we are assuming that this will not have a significant effect on the resulting structure. After construction of the model, the entire crystal was minimized using the conjugate gradient algorithm.51 Minimization was continued until an root mean square average of the potential energy gradient of 0.1 kcal/ (mol Å), or less, was achieved. After minimization was completed, MD simulations in which the number of molecules (N), the two-dimensional pressure (π), and the temperature (T) were held constant (NπT ensemble) simulations were performed. The Parrinello and Rahman52 method, which allows not only changes in the volume of the cell but also its shape, was used with the Verlet algorithm,53 using a boundary weight parameter of 0.1 times the total weight in a periodic cell. To use the pressure algorithm, the threedimensional equivalent of the two-dimensional pressure (γ) obtained from a Langmuir-Blodgett trough was used. This is facilitated by taking the desired two-dimensional pressure, in units of force/distance and dividing this value by the length of the crystal in question, effectively producing a three-dimensional pressure with dimensions of force per unit area. The two(46) Brooks, B. R. J. Comput. Chem. 1983, 4, 187. (47) Quanta CHARMM Parameter Handbook; Molecular Simulations Inc.: Burlington, MA, 1993. (48) Loncharich, R. J.; Brooks, B. R. Proteins: Struct., Funct., Genet. 1989, 6, 32. (49) Kitson, D. H.; Hagler, A. T. Biochemistry 1988, 27, 5246. (50) Alper, H. E.; Bassolino, D.; Stouch, T. R. J. Chem. Phys. 1993, 98, 9798-9807. (51) Jacoby, S. L. S.; Kowalik, J. S.; Pizzo, J. T. Iterative Methods for Nonlinear Optimization Problems; Prentice Hall: Englewood Cliffs, NJ, 1972. (52) Parrinello, M.; Rahman, A. Phys. Rev. Lett. 1981, 52, 71827190. (53) Verlet, L. Phys. Rev. 1967, 159, 98-103.
Table 2. Surface Pressures and Initial Lateral Stresses for Simulations syst
surf press (mN/m)
XX and YY stresses (MPa)a
DPPC low density DPPC high density DSPE low density DSPE high density POPC low density POPC high density SAPE low density SAPE high density
-66.75 -57.75 -67.75 -50.75 -70.75 -42.75 -70.75 -44.75
0.267 2.30 2.71 2.03 2.83 1.71 2.83 1.79
a The lengthwise dimension corresponding to the phospholipid’s lengthwise axis of the crystal was changing, albeit very slowly, with time at a constant velocity. Therefore a pressure rescaling technique was utilized that reformulated the stresses based on the new measure of the lengthwise cell parameter. This reformulation of the external stress tensor was executed every 7 ps. This effectively kept the stresses at, or very near to, the desired value at all times. Also, in the sign convention of Cerius2, negative pressures and positive stresses both correspond to the system being in tension.
dimensional pressure is related to the pressure obtained from the Langmuir-Blodgett trough by using eq 1,
π ) γ′ - γ
(1)
where γ′ is the surface tension of water (72.75 dyn/cm), γ is the surface tension of the phospholipid in question, and π is the two-dimensional compressive stress placed on the film by the trough barrier. Equation 1 simply states that the actual tension on the film is the surface tension of water less the compressive stress due to the barrier. As is typical of the literature, the terms in eq 1 are written as scalar quantities, when in fact they are vectors. It is important to note that the surface tension terms (γ and γ′) are two-dimensional tensional stresses and π is a twodimensional compressive stress. This three-dimensional pressure value was then converted into the YY and ZZ components of the stress tensor, but the XX component, which is the axis normal to the phospholipid surface, was kept at zero. This exerts force on the sides of the phospholipids without exerting any force lengthwise, approximating the conditions of the LangmuirBlodgett trough by assuming that the air on top of the lipid tails does not significantly perturb the lipid structure. Table 2 shows the two pressures, along with the corresponding initial stresses at which each of the systems was simulated. The thermostat used in all cases was the temperature damping algorithm developed by Berendsen et al.;54 even though the Berendsen thermostat does not reproduce the canonical distribution, it was (54) Berendsen, H. J. C.; Postma, J. P. M.; van Gunsteren, W. F.; Di Nola, A.; Haak, J. R. J. Phys. Chem. 1984, 81 (8), 3684-3690.
5258 Langmuir, Vol. 14, No. 18, 1998
Figure 2. Average internal energy of the simulation vs time (system: DSPE low density).
Figure 3. Average pressure of the simulation vs time (system: DSPE low density). chosen because of the lack of stability of the Hoover thermostat when applied to our simulations. For all simulations, the temperature was set at 294.15 K and the time step used was 1 fs. The temperature relaxation time factor for the Berendsen algorithm was set to 0.1 ps to facilitate quicker equilibration times. All simulations were run until equilibrium was approached, determined by viewing the approach of the system’s average energy, average temperature, and average pressure to a constant value over time. Simulations in which the number of molecules (N), the twodimensional area (A), and the temperature (T) were held fixed (NAT ensemble) were also performed with DPPC and DSPE in the low- and high-density phases. These were used to determine if the equilibrium pressure produced by the NAT simulation was correct compared to experimental data. In these cases the correct area per molecule for the system was constructed for the crystal before the simulation was started. The thermostat used in these simulations was also the temperature damping thermostat55 with a temperature factor of 0.1 ps. The time step and temperature were the same as in the NπT simulations.
Results and Discussion Simulations. All simulations were run until equilibration was approached, which was between 100 and 140 ps. This was determined by viewing the respective system’s running average energy, running average pressure, and running average temperature, with time. The running average is calculated by taking all previous points at a specific time and averaging them. In all cases, these properties changed significantly at the initial stages of the simulations but eventually approached asymptotic values. Figures 2-4 give results for the DSPE low-density system, but these are typical of all the systems studied. Also viewed was the instantaneous area per molecule versus time (Figure 5). After the phospholipid systems reached a near-equilibrium state, the properties of interest were obtained. Like most self-assembled lipid simulations, these simulations produce a near-equilibrium structure because
Mauk et al.
Figure 4. Average temperature of the simulation vs time (system: DSPE low density).
Figure 5. Instantaneous area per molecule vs time (system: DSPE low density).
Figure 6. Pressure-area isotherm of DPPC with simulated points.
they neglect certain long-time motions that cannot be addressed in the nanosecond range of typical simulations. The nature of this near-equilibrium state is discussed below. Consequently, any conclusions regarding the structure of the lipid molecules in these and other simulations are tempered by the knowledge that some long-time motions are not included in the simulation. Pressure-Area Isotherms. Figures 6-9 depict the DPPC, DSPE, POPC, and SAPE isotherms with the simulated points in the low- and high-density phases. The hyperbolic shapes of the DSPE, POPC, and SAPE Langmuir-Blodgett isotherms indicate that these monolayer systems are above the critical temperature for liquidphase-solid-phase transitions. A phase transition would appear as a flat region in the isotherm such as the one observed in Figure 6 which suggests that a liquid-phasesolid-phase transition may be occurring in the DPPC monolayer at the measurement temperature. The short span of this flat region suggests that the phase transition
Structures of Self-Assembled Lipid Monolayers
Langmuir, Vol. 14, No. 18, 1998 5259 Table 3. Difference in Initial and Final Areas per Molecule
system DPPC low density DPPC high density DSPE low density DSPE high density POPC low density POPC high density SAPE low density SAPE high density
Figure 7. Pressure-area isotherm of DSPE with simulated points.
Figure 8. Pressure-area isotherm of POPC with simulated points.
Figure 9. Pressure-area isotherm of SAPE with simulated points.
is spanning the upper end, or near-critical, area of the two-phase region. The high-density DPPC simulation occurs in the solid region, whereas the low-density simulation appears to be very near the liquid-solid saturated point in the isotherm curve (Figure 6). Since size limitation effects in NPT simulations typically preclude them from containing coexisting phases,55 the low-density simulation should reflect the liquid structure of the film. The area fluctuations in the simulation, discussed below, do qualitatively support this supposition. Table 3 shows the cross-sectional area per molecule before and after each molecular dynamics simulation. In all cases, the simulated points compare very favorably to the experimentally determined isotherms. The maximum deviation in the cross-sectional area from any one of the simulations compared to its pressure-area isotherm was 1.59% (POPC in the high-density phase). The 95% confidence interval for the surface pressure was calculated for all cases and was on the order of 10-2 mN/m which is (55) Strandburg, K. J. Rev. Mod. Phys. 1988, 60 (1), 161.
area per molecule area per molecule before simulation after simulation percent (Å2/molecule) (Å2/molecule) change 70.00 50.00 70.06 60.00 94.00 61.00 116.00 70.38
69.85 49.95 69.77 59.53 93.57 60.15 115.76 70.07
0.21 0.094 0.41 0.78 0.46 1.39 0.21 0.44
within the size of the circles in Figures 6-9. The simulations reproduce the packing behavior of the phospholipids at high, as well as low, surface pressures. The good reproduction of the experimental packing densities is due in part to the fortuitous choice of the initial volume, which was based on experimental data. NπT simulations similar to the ones presented here for the four lipid systems suggest that these simulations will not deviate significantly from the initial starting point.56 These simulations are identical except the two-dimensional pressure in the simulation has been shifted significantly from the experimentally measured pressure. Despite this change in pressure there is no significant change in the area produced after 120 ps of equilibration. This implies that an accurate initial guess based on experimental data is critical to the performance of the simulations. The failure of the simulations to properly equilibrate is currently under investigation. Other simulations of constant area and temperature (NAT ensemble) by Stouch and co-workers have observed large-scale relaxations that occur well beyond the 1 ns range.29 However, the exact impact of these long-time equilibrations on the structure of the equilibrated lipid layer was not obvious. The small degree to which the area changes in our simulations (see Table 3) and the preliminary results of our perturbation studies56 indicate that a good initial guess for the area density of the lipid layer is required. One possible source of error that may be responsible for the slow equilibration of the cross-sectional area is the failure of the united atom model in crystalline systems. To investigate this, energy minimizations using the conjugate gradient method52 were run on the crystal form of DLPE that was used to generate the initial structures with and without explicit hydrogens. When the crystal structure of DLPE containing all the explicit hydrogens was minimized, the initial and ending structures showed very little difference in conformation. However, energy minimization with the united atom force field resulted in the unrealistic crystal structure pictured in Figure 10. This failure of united atom force fields to reproduce accurate crystal structures is typical in alkane structures. The symmetry of the alkane hydrogens is required to reproduce stable crystal structures. Despite this potential drawback, the united atom force field did reproduce the experimentally measured alignment of the lipid tails fairly well (see discussion below). This occurs because the hydrocarbon tails of the lipids are not crystalline, even in the high-density simulations. This can be observed directly by viewing the simulations. Typical behavior is seen in the end of the DPPC high-density simulation in Figure 11. The alkyl chains on high-density-phase phospholipids are more ordered than the tails viewed in the low-density state. Although aligned, the hydrocarbon tails are not so ordered, so as to be considered crystalline, (56) Mauk, A.; Ludovice, P. Manuscript in preparation.
5260 Langmuir, Vol. 14, No. 18, 1998
Mauk et al.
tally, more so than those witnessed in the constant pressure simulations presented here. This large discrepancy between the NAT and NπT simulations is also observed in the recent results by Tieleman and Berendsen.59 The structure of the hydrocarbon chains was not reproduced well when compared to experiment as discussed below. Order Parameters. The order parameter is one of the most widely used methods in determining the structural features of hydrocarbon chains in phospholipids. Experimentally, order parameters are obtained from selectively deuterating different methylene chains in the alkyl chains of the phospholipids and analyzing the quadrupole splitting in the 2H NMR spectrum while the phospholipids are in the form of vesicles.60,61 The quadrupole splitting is dependent on the time-averaged orientation of the C-D bond vector with respect to the bilayer normal. This time-averaged orientation is usually :60,61 defined in terms of an order parameter, Snmr j
Snmr ) 1/2〈3 cos2 θi - 1〉 i Figure 10. Structure of DLPE (a) on the basis of the crystal structure from X-ray diffraction and (b) after energy minimization using the united atom approximaton. Table 4. Calculated and Experimental Pressure from NAT Simulations system
calcd press (N/m)
exptl press (N/m)
DPPC low density DPPC high density DSPE low density DSPE high density
147.7 ( 0.53 122.9 ( 0.61 -489.0 ( 0.18 -426.6 ( 0.37
66.75 57.75 67.75 50.75
as opposed to simulations involving alkanethiols on gold where the chains are crystalline in arrangement and the inclusion of explicit hydrogen would be important.57,58 Because the lipid tails are not crystalline, the united atom model is not a significant source of error in the structures resulting from the simulations. Constant Area Simulations. Constant area simulations were performed on DSPE and DPPC as described in the Experimental Section, and the final pressure was calculated. The results of these NAT simulations, as well as the pressures determined experimentally from Langmuir trough experiments, are given in Table 4. The range for the simulated pressure is the 95% confidence interval based on the last 50 ps of the NAT simulation. As displayed in Table 5, the pressure determined from the NAT simulations is off by an order of magnitude or more in all cases. Furthermore, the calculated pressure for the DPPC simulations is positive, indicating the force on the lipids is pushing them together, as opposed to pulling them apart as with the case of the DSPE simulations. The above results indicate the inability of constant volume-temperature simulations to adequately simulate the properties of these phospholipid systems. It appears that the lack of fluctuations in the periodic cell restricts the phospholipids from assuming conformations which would allow transition to lower energy structures to occur. This inability to assume lower, more realistically representative, structures is seen in the constant area simulation of the DLPE bilayer by Damodaran and co-workers.22 In this study the order parameters of the hydrocarbon tails are noticeably higher than those calculated experimen(57) Mar, W.; Klein, M. L. Langmuir 1994, 10, 188-196. (58) Shnidman, Y.; Eilers, J. E.; Ulman, A.; Sellers, H. Mater. Res. Soc. Symp. Proc. 1993, 291, 211-216.
(2)
where θi is the angle between the bilayer normal and the C-D vectors at the carbon atom j. Since united-atom carbons were used in all simulations, the order parameter cannot be determined explicitly in this manner. Therefore, a new parameter, Smol i , needs to be defined. It can be can be defined and related to the above shown that Smol i equation by60,61
) 1/2〈3 cos2 βi - 1〉 Smol i
(3)
) Smol 2Snmr i i
(4)
where βi is the angle between the bilayer normal and the vector joining the midpoints of the carbons j - 1 and j + 1. The relationship between the measured and calculated order parameter of eq 4 comes from assuming an axially symmetric movement of the CH2 segment about the long molecular axis. The order parameter can vary between 1 for chains that are perfectly normal to the surface to -1/2 for chains that are parallel to the surface. Typical Smol values range from near 1, in an all trans conformai tion, to zero in the amorphous disordered state. Also, since the vector joining the midpoints of carbons j - 1 and j + 1 is always perpendicular to the Cj-D bond vector, the conceivable problem in identifying the segmental rotation axes due to hydrocarbon chains with cis double bonds is not involved in the calculation. All order parameter plots are done in terms of Smol i . Figure 12 compares the order parameter of the simulated sn-1 and sn-2 alkyl chains of DPPC in the low-density phase with that determined experimentally. The error bars indicate the standard deviation of the order parameter of each carbon in the alkyl chains over the time the simulation was at equilibrium. The carbons near the headgroups of both chains are slightly higher than the experimentally determined ones, while the later carbons match up with the experimental results fairly well. The decrease in the order parameter as the carbon atom number increases is due to the enhanced freedom of the carbons as distance away from the headgroup is increased. The “kink” in carbon 3 is generally attributed to the (59) Tieleman, D. P.; Berendsen, H. J. C. J. Chem. Phys. 1996, 105, 4871-4880. (60) Seelig, J. Q. Rev. Biophys. 1977, 10, 353-418. (61) Seelig, J.; Seelig, A. Q. Rev. Biophys. 1980, 13, 18-61.
Structures of Self-Assembled Lipid Monolayers
Langmuir, Vol. 14, No. 18, 1998 5261
Figure 11. Snapshot of the DPPC solid-phase monolayer. Snapshot was taken near the completion of the simulation.
Figure 12. Comparison of molecular order parameters of simulated DPPC in the low-density phase with experimentally determined order parameters of DPPC. The experiment was carried out at a temperature of 41 °C in the vescicle bilayer state60 vs 21 °C and a Langmuir monolayer in the simulation.
Figure 13. Comparison of the order parameter of DSPE simulated in the low-density phase with the order parameter determined experimentally from DMPE, carried out at a temperature of 50 °C.63
conformation the carbon needs to take in the sn-2 chain in order for the rest of the chain to proceed roughly parallel to the bilayer normal. Both the simulated and experimental systems are believed to be in the low-density phase. The experiment was carried out at 41 °C,60 which is believed to be in the liquid phase for liposomes based on calorimetry data.62 The slight discrepancies between the results in Figure 12 are probably due to the simulation being a flat monolayer while the experiment was made on a bilayer in a liposome. No experimental order parameters of DSPE were available. However, the order parameter has been measured for dimyristorylphosphatidylethanolamine
(DMPE), the only phosphatidylethanolamine to have its order parameter resolved experimentally. By comparing with DSPE, some attributes, especially at the carbons near the headgroup, may be elucidated with DMPE. The DMPE experiment was carried out on DMPE liposomes at 50°63 at which point the DMPE should be in its liquid phase. With the exception of some long-range correlations three-dimensional systems above the critical temperature have molecular structure that is very liquidlike. Therefore we assume that the supercritical DSPE low-density simulation should have similar structure to low-densityphase DSPE liposomes. Figure 13 illustrates the qualitative similarity of the DSPE simulation with the liposomal DMPE. Both systems exhibit a flat region in the order parameter that
(62) Marsh, D. The Handbook of Lipid Bilayers; CRC Press: Boca Raton, FL, 1990.
(63) Marsh, D.; Watts, A.; Smith, I. C. P. Biochemistry 1983, 22, 3023-3026.
5262 Langmuir, Vol. 14, No. 18, 1998
Mauk et al.
Figure 14. Comparison of the order parameters of DSPE simulated in the low- and high-density phases. Figure 17. Comparison of molecular order parameters of simulated DPPC in the low-density phase obtained from an NAT simulation with experimentally determined order parameters of DPPC.60
Figure 15. Comparison of the order parameter of POPC simulated in the low-density phase with the order parameter of POPC determined experimentally at a temperature of 27 °C.64 Sn-1 chain. Figure 18. Comparison of the order parameter of DSPE simulated in the low-density phase obtained from an NAT simulation with the order parameter determined experimentally from DMPE.63
Figure 16. Comparison of the order parameter of POPC simulated in the low-density phase with the order parameter of POPC determined experimentally at a temperature of 27 °C.64 Sn-2 chain.
indicates a high degree of alignment near the headgroup. Figure 14 compares the high- and low-density DSPE simulations, suggesting similar qualitative behavior, but with a higher degree of alignment in the high-density simulation. As was expected all the simulated order parameters had a higher degree of alignment in the highdensity phase. Overall the simulations have a higher degree of alignment than that of the liposomal DMPE. This could be due to inherent differences between the liposomal bilayer phase and the monolayer. Given the longer hydrocarbon chains in the DSPE, a greater degree of alignment is expected for these systems. Figures 15 and 16 show the differences of simulated high-density POPC with that determined experimentally for the sn-1 and sn-2 chains, respectively. The POPC experiment was run at 27 °C which is in the solid phase for liposomal POPC. Although the simulated and ex-
perimental results are qualitatively similar, the simulated order parameter is lower, suggesting less alignment than indicated in the experimental system.64 Again, these differences could be due to structural differences between the liposomal bilayer and monolayer configurations, suggesting that the mobility of the hydrocarbon chains might be reduced in bilayer form. In bilayer vesicles the end of the hydrocarbon chains interface with the hydrocarbon chains of the opposing half of the bilayer, whereas in the simulation these are exposed to vacuum. The cis double bond on the sn-2 chain of POPC allows the sn-1 chain to become much more mobile because the packing of the phospholipids becomes much looser as opposed to a sn-2 chain with no double bonds. The rigid double bond effectively reduces tight packing. In Figure 16, the principal feature is the ability of the simulation to capture the behavior of the cis double bond on the sn-2 chain. The double bond, located between carbons 9 and 10, produces a sharp decrease in the order parameter. The NAT simulations did not reproduce the experimentally characterized order parameter as well as the NπT simulations. The order parameters for the lowdensity simulations of DPPC and DSPE are seen in Figures 17 and 18, respectively. These figures clearly illustrate the failure of the NAT simulations in this regard. The high degree of order in the sn-2 chain of DPPC in Figure 17 suggests that this chain is in a very highly ordered, crystalline alignment, while the sn-1 chain has lower than (64) Seelig, A.; Seelig, J. Biochemistry 1977, 16, 45-50.
Structures of Self-Assembled Lipid Monolayers
Langmuir, Vol. 14, No. 18, 1998 5263
expected order parameters. Figure 18, depicting the DSPE alkyl chains, shows these same types of trends expressed in the DPPC NAT low-density simulation. Volume Fluctuations. The more accurate alkane chain structure produced by the NπT simulation may be due to the area fluctuations permitted by this simulation algorithm. These fluctuations could be responsible for the improved equilibration of the NπT algorithm over the NAT algorithm. The levels of these fluctuations are determined directly by the boundary weight parameter used in the Parrinello and Rahman constant pressure algorithm.53 A boundary weight value of 0.1 times the unit cell weight was used in all the simulations. Whether or not the fluctuations are correct can be determined by comparing the mean squared area fluctuations to the Langmuir isotherm experiments. The mean squared area fluctuations 〈δA2〉 for the NπT ensemble are defined as follows.
〈δA2〉NπT ) 〈(A - 〈A〉)2〉 ) 〈A2〉 - 〈A〉2
(5)
The averages are computed by summing the area for a particular state times that probability associated with that state (Pj)
〈δA2〉NπT )
∑j Aj2Pj - 〈A〉2
(6)
where
Pj ) e-πAj/kT/Z
(7)
and the state density function for the NπT ensemble (Z) is given by
Z)
∑j Q(N,Aj,T)e-πA /kT j
(8)
The term Q(N,Aj,T) is the state density function for the constant area (A) and temperature (T) ensemble (NAT), which is the two-dimensional equivalent of the NVT partition function. Substituting eq 6 into eq 5, we obtain
1
Aj2e-πA/kT - 〈A〉2 ∑ Z j
〈δA2〉NπT )
(9)
Equation 9 can be further modified as follows
〈δA2〉NπT )
-kT ∂ Z
〈δA2〉NπT )
Aj2e-πA/kT - 〈A〉2 ∑ ∂T j
-kT ∂ (〈A〉Z) - 〈A〉2 Z ∂π
〈δA2〉NπT ) -kT
∂〈A〉 kT〈A〉 ∂Z - 〈A〉2 ∂π Z ∂π
(10)
(11) (12)
Realizing
1 ∂Z -〈A〉 ) Z ∂π kT
(13)
we can substitute eq 13 into eq 12 to obtain the following result.
〈δA2〉NπT ) -kT
∂〈A〉 ∂〈A〉 + 〈A〉2 - 〈A〉2 ) -kT ∂π ∂π
(14)
Table 5. Comparison Between Theoretical and Simulated Area Fluctuations simulationa
exptl (Å4)
simuld (Å4)
factor off
DPPC low density DPPC high density DSPE low density DSPE high density POPC low density POPC high density SAPE low density SAPE high density
3318.781 338.835 178.164 154.690 1211.538 248.389 1472.151 367.471
1.370 463.495 5146.903 1628.487 6851.287 2343.410 6822.831 2741.612
0.000413 1.368 28.888 10.527 5.655 9.435 4.635 7.461
The derivative of the area with respect to the twodimensional pressure (π) can be estimated from the slope of the experimental Langmuir-Blodgett isotherms. This can be compared to the mean squared area fluctuations from the simulation through eq 14 above. The mean squared area fluctuations can be extracted directly from the simulation using eq 5 and then normalized by dividing by 25 (the number of lipids) to produce the specific area in units of area per molecule. The results of this comparison are listed in Table 5. With the exception of the high-density DPPC simulation, the fluctuations are within 1 order of magnitude of what they should be estimated from the Langmuir-Blodgett isotherms. This would suggest that a slightly higher value of the boundary weight parameter might be used to better reproduce the fluctuations. With the exception of DPPC, the fluctuations are lower for the higher density systems in both the simulation and experimental results. This is consistent with the suppressed fluctuations one would expect at the higher density. The large discrepancy between the simulation and the estimated fluctuations for the high-density DPPC case may be due to the phase transition apparent in Figure 6. If the high-density point in the DPPC isotherm is in the two-phase region, as it appears in Figure 6, then the above derivation of the mean squared fluctuations is invalid. The above derivation suggests that when the slope of the Langmuir-Blodgett isotherm becomes flat, the fluctuations approach infinity. This is exactly what happens at the critical point when the isotherm contains an inflection point. However, the fluctuations estimated from lowdensity measurement of DPPC are incorrect because the above derivation used to extract them from the data is valid for a single-phase system only. The low fluctuations observed in the simulation suggest that the simulation is qualitatively more like a two-phase region than a nearcritical fluid. However, two-phase regions are rarely observed in NPT simulations of this size.56 Most likely the low-density DPPC simulation is typical of the condensed liquid phase. Although it is difficult to ascertain whether the flat region in Figure 6 is a two-phase region or a critical point, the simulations suggest that a phase separation does occur. Radial Distribution Functions. The aqueous environment in the immediate vicinity of the lipid headgroups can be better understood through the use of the radial distribution function, or gij(r), which is proportional to the probability of finding atom i at a distance r from atom j. The function gij(r) is defined as
gij(r) )
Vtot. dNij Nij,tot. dV
(15)
where Vtot. is the total volume of the periodic cell, Nij,tot. is the total number of ij atom pairs, dV(r) is the volume differential between a sphere of radius r and r + δr, and dN(r) is the number of pairs contained within that volume
5264 Langmuir, Vol. 14, No. 18, 1998
Figure 19. gij(r) of the water oxygens (OW) with respect to the amide nitrogens (NT) in the DPPC headgroups. Shown are DPPC in the low- and high-density phases.
element. Given the two-dimensional nature of these simulations, the differential volume element in eq 15 is no longer spherical for r greater than half the thickness of the two-dimensional film. The differential volume was adjusted in this region of r similar to previous calculations on three-dimensional systems.65 The value for Vtot. used in the normalization of eq 14 was based on an estimated average of the thickness of the two-dimensional system. Figure 19 is the pair distribution function for wateroxygen and lipid-nitrogen pairs and depicts the structure of water around the nitrogen atom in the headgroup for DPPC in both phases. The peaks in the low-density phase indicate solvation shells around the nitrogen atom contained in the headgroup. The position of this first peak is consistent with simulations done by Alper et al.36 It should also be noted that in Alper’s simulations, the phospholipids were minimized and subsequently fixed while the water was able to move freely, making the simulation not entirely realistic. These solvation shells are more pronounced at longer ranges than those observed in either experiment or other simulations. One possibility is the difference in water models used. The simulations presented here used the TIP3P model while other simulated results used the SPC water model. However, on the basis of separate molecular dynamic simulations of water, comparisons between the two models66 do not depict any major structural differences, leaving the explanation of the long-range coordination shells unclear. Bulk water NPT simulations using both TIP3P and SPC water models produced unusually strong solvation peaks. This suggests that the long-range order solvation shells in the water may be an artifact of the NπT simulation carried out in this study. The first solvation peak in the DPPC high-density phase is almost negligible. This comparison in Figure 19 suggests that low-density DPPC is significantly more hydrated than high-density DPPC. Very different behavior is observed for the other lipids (DSPE, POPC, and SAPE) whose water pair distribution functions were similar for both the high- and low-density phases. Representative behavior can be observed for DSPE in Figure 20. The position of the first solvation shell is consistent with DLPE bilayer simulations by Damodaran and co-workers.22 The simulations predict a drastic change in the solvation of DPPC when the packing density is increased, which is not seen in the other lipid systems. This could be due to the fact that the DPPC is undergoing a phase change and the other systems are not. This phase (65) Theodorou, D. N.; Suter, U. W. J. Chem. Phys. 1985, 82, 955. (66) Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D. J. Chem. Phys. 1983, 79, 4926-4935.
Mauk et al.
Figure 20. gij(r) of the water oxygens (OW) with respect to the amide nitrogens (NT) in the DSPE headgroups. Shown are DSPE in the low- and high-density phases.
change is observed in the Langmuir-Blodgett isotherms for DPPC but not for the other systems. If this is the case, then the simulations indicated that a change in headgroup solvation accompanies the low-density-high-density phase transition in DPPC. However, this indication is suspect because of the failure of these systems to completely equilibrate. Both the volume fluctuations and the pair distribution functions above indicated that a phase transition may be occurring in the DPPC system. In an attempt to discern whether a phase transition had occurred in the DPPC simulations, the two-dimensional pair distribution functions of the lipid headgroups were calculated. This is the two-dimensional analogue of eq 15 that considers interactions only in the plane parallel to the surface of the monolayer. This function characterizes any changes in the packing of the lipid headgroups that might occur during a phase change. Although all of these functions showed a moderate difference in the order between the high- and low-density simulations, this difference was not significantly greater for the DPPC case, suggesting that a classic phase change is not occurring. The lack of a large difference in two-dimensional order between the two DPPC simulations does not preclude a phase change. This can occur because the solid phase of a two-dimensional system can be incommensurate with the other component of its interface. This occurs in monolayer adsorption of gas molecules on solid substrates.56 Although a phase transition is observed thermodynamically, the resulting higher density phase is not highly ordered relative to the low-density phase because its dense packed geometry is incommensurate with the lattice spacing of the atoms in the solid substrate. A similar phenomenon is possible in these systems because the lipid headgroups may pack in a fashion that is incommensurate with the geometry of the order induced in the water from hydrogen bonding. What results are two distinct phases with similar levels of structural order. This is similar to the transition from a melt to a glass, except that this is kinetically stabilized, whereas the two-dimensional version is stabilized by the incommensurate geometry between the substances at the interface. Given the simulation protocol used, it is possible that the poor solvation in the DPPC high-density phase is an artifact of simulation. Previous simulations position water molecules around the headgroups to ensure that they are solvated.22,27-31,35 The simulation approach used here placed an equilibrated layer of water in contact with the lipid layer and allowed solvation to occur during the simulation. This solvation did occur for seven out of eight simulations, suggesting that this is not energetically
Structures of Self-Assembled Lipid Monolayers
Figure 21. Dihedral distribution of the low-density DPPC sn-1 chain.
Figure 22. Dihedral distribution of the high-density DPPC sn-1 chain.
prohibited during these relatively short simulations. Given the structures of these four lipids, it is not unreasonable that only DPPC would expel water from its hydration shell upon compression. The large size of the phosphatidylcholine headgroup compared to phosphatidylethanolamine suggests the headgroup repulsions would dominate as the density of these two-dimensional layers is increased. Since only the headgroups are solvated, water would be expelled as the headgroups approach a close-packed density. The POPC system is less likely than DPPC to allow headgroup repulsions to dominate because the effective area occupied by its tail groups is much larger. This occurs because the unsaturation produces a kink in the hydrocarbon tail that causes the tail to occupy additional volume. Therefore, of these four lipids, DPPC would be the first to be dominated by headgroup repulsion and therefore expel waters from its headgroup hydration shell. Dihedral Angle Distribution. The dihedral angle distributions were calculated for all systems to further evaluate the differences between different types of phospholipids and phases. The probability of finding a dihedral angle of the hydrocarbon chains in the trans (180 ( 30°), anticlinal (+120 ( 30 and -120 ( 30°), gauche (+60 ( 30 and -60 ( 30°), and syn (0 ( 30°) were plotted for selected cases of DPPC, POPC, and SAPE, shown in Figures 2124. Figures 21 and 22 depict the sn-1 chain of DPPC in both phases. The results for the sn-1 chain for both phases are as expected. The low-density phase shows a high amount of trans probability throughout the entire chain. The high-density phase shows the chain as completely trans for almost its entire length, which is expected because the tight packing required in the high-density phase necessitates the need for the trans conformation for optimal packing. Despite the indication that essentially all dihedral angles are within (30° of 180°, visual examination of the DPPC solid phase confirms that this is not a crystalline phase. Both the variation of the
Langmuir, Vol. 14, No. 18, 1998 5265
Figure 23. Dihedral distribution of the low-density POPC sn-2 chain.
Figure 24. Dihedral distribution of the low-density SAPE sn-2 chain.
dihedral angle about 180° and the rotational and orientational variation of the alkane chain axis prevent crystallization. The sn-2 chain of DPPC gives results slightly different from those seen for the sn-1 chain (not shown). The latter part of the chain depicts the trans state in both phases as expected, but the section near the head shows strong percentages of gauche and anticlinal states. This is due to the sn-2 chain’s need to bend near the headgroup to achieve a structure in which the rest of the chain can optimally pack with neighboring chains. Also, as described earlier, the sn-2 chain protrudes from the glycerol backbone at an angle closer to parallel to the bilayer normal than perpendicular to it, further supporting the case for initial chain bending. The most prevalent feature of Figure 23 is the exclusive syn conformation at the position carbon 10. The syn conformation is the only one found, and it is located at the cis double bond in the chain. Furthermore, the usually unfavored anticlinal conformation is preferred over all other conformations for the bonds located on the other side of the double bond. After the double bond, the chain starts to assume the more normally preferred trans conformations. The phenomena of a preferred anticlinal conformation on either side of a double bond is no more apparent than with the sn-2 chain of SAPE (Figure 24). A high anticlinal percentage is apparent throughout the whole chain due to the four cis double bonds. Self-Diffusion Coefficients. The self-diffusion coefficients can be determined using the Einstein relation for self-diffusion:
D)
1 〈(r (t) - ri(0))2〉 6∆t i
(16)
where 〈(ri(t) - ri(0))2〉 is the mean square displacement
5266 Langmuir, Vol. 14, No. 18, 1998
of the specified atoms at position vector ri(t) compared to the initial position vector ri(0). A linear fit of the mean squared displacement of the lipid headgroups increased more rapidly with time for the low-density systems, indicating that they have a higher diffusivity. This result is expected, but the fluctuations in this quantity preclude the results from being statistically significant. In fact the diffusivity is essentially zero because we do not sample enough time to allow noticeable diffusion of the phospholipid headgroup. Equation 16 can be used to illustrate that only a near-equilibrium state is obtained during typical simulations. A typical value for the translational diffusion constant of phospholipids is of the order of 10-7 to 10-8 cm2/s.58 Given the two-dimensional density obtained from Figures 6-9, the typical distance traveled between a lipid that swaps positions with its neighbor is of order 8 Å. Substituting the square of this value and the approximate diffusivity above into eq 16, we can calculate an approximation of the time it would take for each lipid to move one lattice site in the two-dimensional plane. This value represents a minimum value for equilibration of these two-dimensional liquids and is approximately 10-100 ns, which is much longer than the simulations done to date. The higher diffusivity of water allows the calculation of the its diffusivity, which can then be compared to experiment to determine the validity of the dynamic behavior of the water model used. The coefficient was determined for all simulations for the last 20 ps. Then results for each simulation were averaged together to obtain the final result. The calculated average diffusion constant of water from the simulations at 21 °C was 2.70 × 10-5 cm2 s-1. This result compares favorably with the experimental result of 2.43 × 10-5 cm2 s-1 67 measured at 25 °C, as well as the simulated value of 2.59 × 10-5 cm2 s-1 for the SPC/E model of water68 at 42 °C. It is expected that the results from the phospholipid simulations would be lower than normal instead of higher, as was indicated in all but two of the phospholipid simulations, because of the hydration shells around the phospholipids. The waters in this area have restricted mobility due to the bonding of the waters with various parts of the headgroups, which would lower the overall self-diffusion coefficient as compared to experiment and the SPC/E water simulation which used only “neat” water for their calculations. While some simulations result in lowered water diffusivity,22 Stouch and co-workers suggest that with longer range energy cutoffs the water adjacent to the lipid layer has a local diffusivity that is higher than the bulk value.51 Given the approximate nature of the water model employed, the simulated diffusivity is in good agreement with the experimental value.
Mauk et al.
Conclusions The NπT molecular dynamics simulations of monolayers of four different phospholipids provided a reasonably accurate picture of their atomic structure. The simulated pressure-area behavior at both low and high densities compared well with Langmuir-Blodgett trough measurements. The simulated order in the alkane chains at low density compared well with results of H2 NMR experiments. One of the four phospholipids (DPPC) exhibited a flat region in the Langmuir-Blodgett isotherm indicative of a phase change. The behavior of the simulation was qualitatively consistent with a phase change. The simulation indicates that this phase change involves a dehydration of the lipid headgroups of DPPC and not a change in the order or structure of the packing of the DPPC lipids. This result is important in the design of lipid-based biomimetic materials. Deposition of lipids on substrates at a high density may produce poor biomimetic materials because the are poorly solvated. Recent atomic force microscopy measurements indicate that DPPC is deposited at this high density on hydrophobic substrates.69 Simulation of lipid monolayers in the NπT ensemble is superior to simulation in the NAT ensemble. Consistent with previous work, NAT simulations failed to equilibrate to the appropriate pressure, while NπT simulations equilibrated to areas close to their experimentally measured values. The NAT simulations failed to reproduce the experimentally measured structure of the lipid alkyl chains, while the NπT simulations accurately modeled this structure. However, the NπT simulations have not fully equilibrated, as indicated by the negligible diffusion of the lipid headgroups and preliminary results from perturbation studies currently underway.57 Therefore, NπT simulations depend on an accurate initial density for the simulation. The failure of the NπT simulations to reach a true equilibrium means that any conclusions drawn from these and the many similar simulations published in the literature are suspect. We speculate that the rate limiting step in the equilibration of these simulated systems is the equilibration of the water around the headgroups of the lipids. This phenomenon is currently under investigation. Acknowledgment. The authors gratefully acknowledge the financial support of the Whitaker Foundation Biomedical Engineering Research Grant Program, the Georgia Tech Molecular Design Institute, and the EmoryGeorgia Tech Biomedical Technology Research Center. We also are grateful for the insightful discussions with Dr. Jeff Morris of Georgia Tech. LA970987R
(67) Skelland, H. P. Diffusional Mass Transfer; Krieger Publishing Co.: Malabar, FL, 1985. (68) Berendsen, H. J. C.; Grigera, J. R.; Straatsma, T. P. J. Phys. Chem. 1987, 91, 6269-6271.
(69) Stephens, S. Structual Characterization of Model Membranes using Infrafed Spectroscopy and Scanning Probe Microscopy, Ph.D. Thesis, University of Georgia, Athens, GA, 1996.
