J. Phys. Chem. 1992, 96, 10497-10506 Bookbinder, D. C.; Wrighton, M. S.J. Electrochem. Soc. 1983, 130, 1080. (d) Nomura, K.; Hirayama, K.; Ohsaka, T.; Nakanishi, M.; Hatozaki, 0.; Oyama, N. J. Macromol. Sci.-Chem. 1989, A26, 593. (2) (a) Martigny, P.; Anson, F.C. J. Electroanal.Chem. 1982, 139, 383. (b) Oyama, N.; Oki, N.; Ohno, H.; Ohnuki, Y.; Matsuda, H.; Tsuchida, E. J. Phys. Chem. 1983,87,3642. (c) Coche, L.; Moutet, J.-C. J. Am. Chem. Soc. 1987, 109. 6887. (d) Coche, L.; Moutet, J . 4 . J. Electroanal. Chem. 1987,224, 1 1 1 . (e) Coche, L.; Moutet, J.-C. J. Electroanal. Chem. 1988,245, 313. ( 3 ) (a) Dominey, R. N.; Lewis, T. J.; Wighton, M. S.J. Phys. Chem. 1983,87,5345. (b) Lewis,T. J.; White, H. S.;Wrighton, M. S.J. Am. Chem. Soc. 1984,106,6947. (c) Shu, C.-F.; Wrighton, M. S.J. Phys. Chem. 1988, 92,5221. (d) Bidan, G.;Deronizier, A.; Moutet, J.-C. J. Chem. Soc., Chem. Commun. 1984,1185. (e) Downard, A. J.; Surridge, N. A.; Gould, S.;Meyer, T. J.; Deronizer, A.; Moutet, J. C. J. Phys. Chem. 1990, 94, 6754. (4) (a) Mortimer, R. J.; Anson, F. C. J. Elecrroanal. Chem. 1982, 138, 325. Ib) Smith. D. K.: Lane. G. A.: Wrinhton. M. S.J. phw. Chem. 1%. 92; 2616. (c) Smith, D. K.; Tender, L. M l Lane, G. A.; Licit, S.;Wrighton; M. S.J. Am. Chem. Soc. 1989, I l l , 1099. ( 5 ) (a) Bookbinder, D. C.; Wrighton, M. S.J. Am. Chem. Soc. 1980,102, 5123. (b) Abruna. H. D.: Bard. A. J. J. Am. Chem. Soc. 1981,103,6898. (c) Oyama, N.; Yamaguchi, S.;Kaneko, M.; Yamada, A. J. Electroanal. chem. 1982, 139, 215. (6) Dalton, E. F.; Murray, R. W. J. Phys. Chem. 1991, 95,6383. (7) (a) Ohsaka,T.; Nakanishi, M.; Hatozaki, 0.;Oyama, N. Electrochim. Acta 1990, 35, 63. (b) Ohsaka, T.; Oyama, N.; Sato, K.; Matsuda, H. J. Electrochem. Soc. 1985, 132, 187 1. (8) Bider, A. J.; Michaels, A. S.In Encyclopedia of Polymer Science and Technology;Mack, H. F., Ed.; Wiley: New York, 1969; Vol. 10, p 765. ( 9 ) (a) Murray, R. W. In Electroanalytical Chemistry;Bard, A. J., Ed.; Marcel Dekker: New York, 1984; Vol. 13, p 191. (b) Murray, R. W., Ed. Molecular Design of Electrode Surfaces;Wiley: New York, 1992. (c) Buttry, D. A. In Electrochemical IrerfacepModern Techniquesfor in situ Irerface Characterization;Abruna, H. D., Ed.; VCH Publishers: New York, 1991; p 529. (d) Oyama, N.; Ohsaka, T. Electron Transfer and Mass Transfer Kinetics in Electroactive Polymer Films; Pergamon Press, in press. (IO) Ohsaka, T.; Yamamoto, H.; Oyama, N. J. Phys. Chem. 1987, 91, 3775. (1 1 ) (a) Saveant, J. M. J. Electroanal. Chem. 1988,242,l. (b) Andriex, C. P.; Saveant, J. M. J. Phys. Chem. 1988, 92, 6761. (c) Buck, R. P. J. Electroanal. Chem. 1989. 258, 1. (12) Buttry, D. A. ElectroanalyticalChemistry;Bard, A. J., Ed.; Marcel
10497
Dekker: New York, 1989; Vol. 17, p 1. (13) Oyama, N.; Ohsaka, T.; Yamamoto, H.; Kaneko, M. J. Phys. Chem. 1986, 90, 3850. (14) (a) Oyama, N.; Anson, F. C. J. Electrochem. Soc. 1980,127,640. (b) Oyama, N.; Yamaguchi, S.;Nishiki, Y.; Tokuda, K.; Matsuda, H.; Anson, F. C. J. Electroanal. Chem. 1982,139,371. (c) Buttry, D. A.; Anson, F.C. J. Am. Chem. Soc. 1983,105,685. ( d ) Anson, F. C.; Buttry, D. A.; Saveant, J. M. J. Phys. Chem. 1984,88, 3086. ( 1 5 ) Bard,A. J.; Fauher, L. R. Electrochemical Methods, Fundamentals & Applications; Wiley: New Yo&, 1980. (16) Kielland, J. 1. Am. Chem. Soc. 1937, 59, 1675. (17) (a) Ikeda, T.; Schmehl, R.; Denisevich, P.; Willman, R.; Murray, R. W. J. Am. Chem. Soc. 1982, 104, 2683. (b) Anson, F. C.; Ohsaka, T.; Saveant, J. M. J. Phys. Chem. 1983,87,640. (c) Oyama, N.; Ohsaka, T.; Okajima, T.; Hirokawa, T.; Maruyama, T.; Ohnuki, Y. J. Electnxvurl.Chem. 1985, 187, 79. (d) Gottesfeld, S.;Raistrick, I. D.; Srinivasan, S.J. Electrochem. Soc. 1987, 134, 1455. (e) Ohsaka, T.; Hirokawa, T.; Miyamoto, T.; Oyama, N. Anal. Chem. 1987.59, 1758. ( 1 8 ) Bird, C. L.; Kuhn, A. T. Chem. Soc. Rev. 1981, IO, 49. (19) (a) Kaufman, F. B.; Engler, E. M. J. Am. Chem. Soc. 1979,101,547. (b) Kaufman, F.B.; Schroeder, A. H.; Engler, E. M.; Kramer, S.R.; Chambers. J. 0. J. Am. Chem. Soc. 1980.102.483. (20) fickup, P. G.;Kunter, W.; Leidder, C. R.; Murray, R. W. J. Am. Chem. Soc. 1984,106, 1991. (21) (a) Daum, P.; Lenhard, J. R.; Rolison, D.; Murray, R. W. J. Am. Chem. Soc. 1980,102,4649. (b) Degrand, C.; Miller, L. L. J. Electroanal. Chem. 1981. I l l . 267. (cl Demand. C.: Miller. L. L. J. Electroonu/. Chem. 1982, 132, 163. id) Denikvick P.;’Abruna, H. D.; Leidner, C.R.; Meyer, T. J.; Murray, R. W. Inorg. Chem. 1982, 21, 2153. (22) (a) Daum, P.; Murray, R. W. J. Phys. Chem. 1981.85, 389. (b) Buttry, D.A,; Anson, F. C. J. Electroanal. Chem. 1981,130,333. (c) Majda, M.; Faulkner, L. R. J. Electroanal. Chem. 1984, 169, 77. (d) Majda, M.; Faulkner, L. R. J. Electroanal. Chem. 1984, 169, 97. (23) Bock, C. R.; Connor, J. A.; Gutierrez, A.; Meyer, T. J.; Whitten, D. G.;Sullivan, B. P.; Nagle, J. K. Chem. Phys. Lett. 1979, 61, 552. (24) Leidner, C. R.; Murray, R. W. J. Am. Chem. Soc. 1985, 107, 551. (25) Winograd, N.; Kuwana, T. J. Am. Chem. Soc. 1970, 92, 224. (26) Facci, J.; Murray, R. W. J. Phys. Chem. 1981,85, 2870. (27) Facci, J.; Murray, R. W. J. Electroanal. Chem. 1981, 124, 339. ( 2 8 ) Tsou, Y.-M.; Anson, F. C. J. Phys. Chem. 1985, 89, 3818. (29) Buttry, D. A.; Saveant, J. M.; Anson, F.C. J . Phys. Chem. 1984,88, 3086.
Nonequlllbrlum Molecular Dynamics Slmulation of Oxygen Dlffuslon through Hexadecane Monolayers with Varying Concentrations of Cholesterol Stuart J. McKinnon,* Louisiana State University Eye Center, New Orleans, Louisiana 70112 Scott L.Whittenburg, Department of Chemistry, University of New Orleans, New Orleans, Louisiana 70148
and Bernie Brooks National Institutes of Health, Bethesda, Maryland {Received: April 2, 1992; In Final Form: September 14, 1992)
The effect of varying concentrationsof cholesterol on the static properties of a hexadecane monolayer and on the diffusion rate of oxygen across the monolayer has been studied by using a novel application of nonequilibrium molecular dynamics. The results show an increased ordering of the hexadecane chains in the region of the monolayer near the added cholesterol with a general increase in monolayer thickness with increasing concentration of cholesterol. The calculation of the diffusion coefficients from the work required to pull oxygen through the monolayer by using a forcing potential suggests an enhanced diffusion rate of oxygen with added cholesterol and a decrease in the diffusion coefficient at higher cholesterol concentrations.
IntNMlUCtiOll
Nothing is more important for respiratory function than lung surfactant. Lung surfactant is a complex mixture of phospholipids, proteins, and cholcateml, which is manufactured by type I1 alveolar pneumacytes. It first appears in human lung tissue between 24 and 28 weeks gestation. It serves to act as a surface-activeagent whose major function is to prevent lung collapse during the expiratory phase of respiration.’+ Lung surfactant acts to stabilize 0022-3654/92/2096- 10497$03.00/0
alveoli by reducing surface tension at the air-liquid interface2v4*s and may also form a rigid monolayer that keeps the alveoli open at low internal alveolar diameters.)~~ Surfactant was postulated to exist as early as 1929,’ but it was not until the early 1950s when Pattie* and Clements9 demonstrated a highly surface-active component in lung washings and also in extracts of whole lung tissue? Once the existence of surfactant was realized, it became the object of major research efforts. The 0 1992 American Chemical Society
10498 The Journal of Physical Chemistry, Vol. 96, No. 25, 1992
detrimental effects of a lack of this precious substance are manifest in a variety of human pathological conditions. Premature infants aged less than 28 weeks of gestation do not manufacture surfactant in significant enough quantities to maintain proper lung function. Thisdeficiency has a role in the pathog& of nspiratory distress syndrome (RDS),6Joa major cause of mortality in neonates.*" These infants develop atelectasis (collapse) within a few hours of birth followed either by the formation of a proteinamus hyaline membrane or by death or by recovery with the possibility of long-term lung pathology." Adults who are subject to traumatic injuries such as smoke or chemical inhalation or who suffer viral or bacterial respiratory infections are also subject to the development of a spectrum of respiratory pathology similar to that of the newborn. Surfactant is now known to contain a number of lipids and proteins. The primary lipid constituent is dipalmitoylphosphatidylcholine (DPPC). The surfactant phospholipids are zwitterionic or polar headgroups, with a negatively charged phosphate group and a positively charged amine group bonded to nonpolar, saturated, long-chain hydrocarbons (as in the case of DPPC). Thus the phospholipid has a hydrophobic end, which comsponds to the nonpolar "tail", and a hydrophilic end, which corresponds to the polar "headgroup". The phospholipids orient themselves in such a way as to form membranes and biological compartments in and around cells. The essential role of phospholipids in the process of life thus cannot be overstated. The phospholipid composition of lung surfactant is distinct from that of other tissues such as cell membranes, due to the relatively high concentrations of DPPC and PG and relatively low amcentrations of sphingomyelin and phosphatidylethanolamine. Other reported components include other phospholipids such as phosphatidylinositol, phosphatidylserine, and cardiolipin, as well as cholesterol and high molecular weight apoproteins. The goal of thii research project is to investigate the effect of naturally occurring substances on the rate of oxygen transport across lung membranes. As a result of this study nonantigenic surfactant additives may be identified. The tool chosen for this investigation is the molecular dynamics technique. In the results presented in this manuscript the additive chosen is cholesterol, a naturally occurring substance found in other membranes. Various concentrations of cholesterol will be substituted for long-chain hydrocarbons in a molecular monolayer, and the oxygen diffusion coefficients, hydrocarbon order parameters, and other properties will be determined to see if an optimal concentration of cholesterol exists that enhances oxygen transport across a model lipid surfactant. If cholesterol is found to give the model phospholipid the proper orientation for enhanced oxygen diffusion, then an improved artificial lung surfactant can be designed that is also nonantigenic. The information obtained in this simulation may provide new insights into the design of artificial lung surfactants uscd in conjunction with pulmonary assist devices in clinical settings such as neonatal intensive care and acute respiratory emergencies. Over the past several years, the usefulness of molecular dynamics and the empirical energy functions for investigating the physical properties of a wide variety of molecules has been demonstrated.1c21 Studies range from those concerned with the structures, vibrational frequencies, and energies of small molec u l e ~ to ~ those ~ - ~ dealing ~ with Monte Carlo and molecular dynamics simulations of pure liquids and s o l ~ t i o n s Other . ~ ~ studies deal with the analysis of conformational energies and fluctuations of large molecules such as proteins and nucleic acids in vacuum, solution, or crystal e n ~ i r o n m e n t s , exploration 3~~~ of the temporal behavior of macromolecular systems,'647 exploration of energy decay,'* solvent analy~is,'~ determination of three-dimensional structures using twedimmsional nuclear Overhauscr effect data from NMR experiments,50comparison of relative free energies of similar sysfcnt(L:1-~ and sampling of dormations for structure d~termination.~~"' Molecular dynamics simulationsof lipid systems have ranged from purely theoretical approaches to full molecular representations of micelles and bilayers with electrostatics and water sol-
McKinnon et al. vation. These simulations have resulted in mixed succeas. The first numerical treatment of a lipid monolayer was by MarceljaS8who applied the Maier-Saupe theory5gfor smectic liquid crystals to monolayers. He analyzed a representative 10band molecule and replaced the attractive potential due to all other molecules by a potential of mean force proportional to the order parameters. Marcelja found isotherms that showed a firstsrder transition that was in agreement with experimentallymeasured isotherms near the main transition."' CotterilP performed a molecular dynamics calculation on a monolayer model consisting of 240 rodlike molecules. The interactions were modeled as Lennard-Jones 6-12 potentials acting between the rod end points and a Coulombic potential acting between the rod midpoints. This system underwent a phase transition from an ordered state to an isotropic fluid state. T o x v a ~ modeled d ~ ~ 256 particles that represented point-centcrs of force with a parameterized interaction potential similar to the Lennard-Jones potential. The pressure was calculated for different densities and temperatures by using the virial theorem. Calculated values of the isothermal compressibility are at least an order of magnitude too small, due to the shortcomings of the model.6o Both Cotterill and Toxvaerd only considered models in which the individual molecules are d u c e d to purely two-dimensional objects that move in the plane of the monolayer. The influence of the steric interactions of the chains and their large number of configurations is neglected. These experiments simulated a gel phase with the lipid chains in the all-trans state. Also, the mechanism believed to be responsible for the main phase transition in lipids, the 'melting" of the tails, is absent in these models. Micellar repmentations of lipid interactions have recently been extensively simulated. An MD calculation by Haile and 0 ConnelP was performed on a simplified model of a spherical micelle composed of 40 13-member c h a i i with headgroups fmed on a spherical shell of 16-A radius. Explicit head-head and headdvent interactions were not included. Chain dormations and the distribution of gauche configurations were studied, and it was shown that the chains were clearly not in the all-trans configuration. The general trend was that chains were directed toward the micelle center or parallel to the outer shell. The motion of the tail groups was rapid, which the authors ascribed to the combined effects of group collisionsat high density plus high-speed rotation induced by gauchetrans conformational changes. The model micelle also contained an average of 29% dihedrals in the gauche conformation, which compared favorably to the ideal gas value of 34%. A further MD study of micelles was performed by Jonsson et alau on 15 sodium octanoate monomers with explicit solvation in 1094 water molecules for a total number of 3447 atoms. The interatomic interactions included a Lennard-Jones term and a Coulombic term. Due to deficiencies in the SPC water model, which is speculated to yield inadequate screening of the ionic interactions due to a low dielectric permittivity, the dynamics of the simulation were found to be too fast. These included rotational and translational diffusion of the micelle, reorientation rates of the hydrocarbon chains, and rotational and translational diffusion of the water molecules. The water molecules close to the micellar surface behaved in a similar fashion to that of bulk water, which contradicted the c(~ccptof "bound" water or water molecules that move more slowly in the region of large molecules. This study was significant in that a more accurate model was used which included electrostatic terms and water solvation. Another MD calculation of a sodium octanoate micelle in an aqueous environment was performed by Watanabe et al.65 The system represented consisted of 15 monomers and 1068 water molecules in a cubic box with a si& length of 34.2 & which yielded a micellar concentration roughly at the critical micelle concentration of 0.4 M.63 Initially, the am of the micelle (15 octanoate ions) was equilibrated in the absence of solvent by constraining the head groups to remainon the surface of a sphere. The radius of the sphere was gradually reduced until the am region achieved a density characteristic of liquid The solvent was prepared by equilibrating a system of 1331 water molecules in
Oxygen Diffusion through Hexadecane Monolayers a box with side length of 34.2A. Next, the micelle was inserted at the center of the solvent box, removing all overlapping water molecules. After a 50-ps equilibration, a 250-pssimulation was run. Interestingly, the percentage of gauche conformationswas found to be 22% which was smaller than that found in the study of Haile and O'C0nnell.6~This difference was explained as being due to the we of a different dihedral potential and is in accordance with s s c studies which show that micellization demases the proportion of gauche conformations in the alkyl chains.67dg A recent study of micellar dynamics was performed on a lysophosphatidylethanolamine(LPE) monomer conformation derived from the crystal structure of dilaurylphosphatidylethanolamine70 by Wendoloski et al." The model was assembled by aligning the extended hydrocarbon chains of the LPE molecules along radii of a sphere of radius 48 A, energy minimizing, and then hydrating with a two-layer shell of 1591 water molecules. This model differs from the above simulations in that the LPE monomers include zwitterionic head groups capable of forming a hydrogen-bonded surface network among themselves and with solvent molecules. As was noted experimentally,'2 the time-averaged radial probability functions show increasingly broad distributions for atoms near the end of the hydrophobic tail. This distribution is presumed to be due not only to trans-gauche conformationalstatistics of individual LPE hydrocarbon chains but also to relative motion of the LPE molecules and to flattening of the spheroid. Molecular dynamics calculations have also been carried out on planar lipid configurations, either monolayers or bilayers. A study of a planar hexane monolayer was done by K o x ~in~1980. Kox improved upon previous monolayer simulations by constraining the head group in the plane of the substrate, thus simulating the interaction between the headgroups and the molecules of the substrate. Several exclusions also exist in this study. Polar interactions between the head groups are not included, and the energy difference between trans and gauche conformations was not taken into account, implying no difference between the liquid-expanded and the liquid-condensed ~ h a m . 6The ~ simulation used Verlet's algorithm with two-dimensional periodic boundary conditions. The model showed a fmt-order phase transition from an ordered fluidlike state to a disordered gaslike state. Also, clustering of the head groups and formation of small pores were noted. These pores may play a part in membrane permeability. This study neglected dihedral angle potentials, which caused serious deviations from experimental data. A more realistic simulation of a lipid bilayer was performed by van der Ploeg and B e r ~ n d s e nin~ ~1982. The lipid bilayer consisted of two layers of 16 decanoate molecules that were periodic in two dimensions. Interactions included Lennard-Jones, bond angle, and the dihedral angle potential of Ryckaert and belle man^.^^ Bond lengths were constrained, and head groups were confined near the bilayer surface by harmonic forces r e p resenting their interaction with the water layer. Experimentally measured order parameters were perfectly reproduced, and an interesting cooperative tilt of the molecules, persisting over tens of picoseconds, was observed. This bilayer simulation involved realistic potentials but suffered from a cell size that was too small to exclude artifacts produced by the two-dimensional periodic boundary conditions. A structural correlation existed over the whole unit cell (of dimension 2 nm X 2 nm),which was ptulatad to have been enhanced by the periodic boundary conditions. The authors concluded that a larger sized simulation was required. The study of a larger system of a bilayer consisting of two layers of 64 decanes was done by van der Ploeg and Berend~en~~ in 1983. The simulation was in general agreement with that of the above smaller system, except for the observed cooperative tilt of the decane molecules. The spatial tilt correlation was found to extend over distances of about 1 nm. The percentage of gauche conformations was 21% which was smaller than that expected from Boltvnann statistics, which would yield a value of 38%. The fraction of gauche conformations was found to increase near the end of the chain, in agreement with the experimentally decreased order parameters. The transition dynamics were not discussed
The Journal of Physical Chemistry, Vol. 96, NO. 25, 1992 10499 within this article. Dynamic properties such as diffusion were calculated in this simulation by two different methods, one from the velocity autocorrelation function and the other from the lateral displacements. The lateral diffusion constant of the decanes was found to be an order of magnitude faster than the experimental value. The lack of inclusion of a water solvent may have eliminated a hydrodynamic damping effect, which may be the limiting factor in lateral diffusion. The simulation of a sodium decanoate-decanol-water system was reported in 1988 by Egberts and Berendsen." This system differs from biological membranes in the types of lipids (only single-chain lipids with ionic headgroups) and in the nearby presence of other bilayers, as the model system is a multibilayer system. This simulation was more accurate than the simulations described above in that the head groups of the amphiphilic molecules were created in full atomic detail, enabling the discrimination between decanoate ions and decanol molecules. partial atomic charges were included, aqueous degrees of freedom were explicitly treated, and sodium ions were included to constitute counterions for decanoate ions. An 80-ps MD simulation was performed by using a leapfrog Verlet algorithm with constrained bond lengths. The unit cell contained a bilayer of 52 decanoate ions and 76 decanol molecules, a water layer of 526 water molecules, and 52 sodium ions, for a total of 3166 atoms. Periodic boundary conditions were applied in three dimensions, and the system was weakly coupled to a temperature bath of 300 K by the method of velocity scaling. Lennard-Jones, bond angle, dihedral angle, and Coulombic potentials were utilized, and bond lengths were constrained. The total length of the simulation was 100 ps. The percentage of gauche conformations was found to be 21% for the decanoate ions and 22.5% for the decanol molecula, which was in good agreement with the previous simulation of van der Ploeg and B e r e n d ~ e n . ~ ~ The calculated order parameters were in reasonable agreement with experimentally derived order parameters as measured by Seelig and Niederberger?"* The diffusion constant of water was calculated by the method of lateral displacements. The diffusion constant of water is a factor of 2 larger than that measured experimentally in a palmitate system.83 The authors felt that the difference in diffusion constants could be attributed to the fact that the SPC model water diffuses more rapidly than real water. The effect of the addition of cholesterol on the model membrane has been previously studied by the Monte Carlo simulation technique. Scott used the Monte Carlo simulation method to study the effects of the inclusion of model proteins and also cholesterol 35 decane molecules into lipid monolayers. In a 1986 and one small model polypeptide or protein were simulated, and order parameters were calculated. The main conclusion drawn from this study is that a single polypeptide chain aligned perpendicular to the monolayer plane does not greatly perturb the equilibrium lipid chain states as measured by the order parameter profiles. The motional properties of a lipid chain in a dipalmitoylphosphatidylcholine (DPPC) bilayer were studied in a Brownian dynamics (BD) simulation by Pastor and Venable in 1985.85BD simulations are advantageous, as compared with molecular dynamics (MD), in that longer time-scale (10-100 ns) motions can be explored. The DPPC bilayer was modeled by a mean field derived from an extension of the Marcelja model, and the lipid chain motions were analyzed with and without the presence of the mean field. A comparison of the free-chain and membrane simulations found that the internal dynamics of the lipid chain in the DPPC bilayer were similar to that of neat alkane. A more interesting study was presented by Scott and Kalaskar in 1989 and pertained to the effects of one cholesterol molecule on a model monolayer of 90 hexadecane chains.86 A LennardJones 6-12 interaction was employed, and, as an approximation, the carbons on the rigid ring of the cholesterol were assigned methylene 6-12parameters. The interchain spacing was 6 A and the area per chain was 31.18 A*. The bottom carbon of the cholesterol and the bottom atoms of the hexadecanes were constrained in the same plane, and none of these molecules were
loso0 The Journal of Physical Chemistry, Vol. 96, No. 25, 1992 allowed to move perpendicular to the monolayer plane. Another simplification in this study was the lack of a glycerol backbone connecting pairs of chains. It was found in an earlier work by Scotts7 that chain ordering is affected mainly by the accessible free volume of the system, whether the chains are connected or not. Also, a 16-A cutoff was used in the calculations of nonbonded interactions. From analysis of the calculated order parameters, it was seen that cholesterol strongly hindered, but did not prevent, rotational changes in the lower segments of the hexadecane chains nearest the rigid ring moiety of the cholesterol molecule. The authors conjecture that each cholesterol molecule partially or wholly hinders the ability of about six neighboring lipid chains to change their rotameric state. This has implications for the mechanism by which cholesterol affects the lipid phase transition. By hindering the melting of the tails, the phase transition may occur at a lower temperature or be broadened over a larger range of temperatures. Oxygen diffusion has been measured experimentallyin phospholipid dispersions and erythrocyte plasma membranes through the use of quenched pyrene fluorescence by Fishkoff and Vand e r k ~ o in i ~1975. ~ The oxygen diffusion coefficient was quantitatively determined by the use of the diffusion-limited reaction between the paramagnetic species (oxygen) and the singlet excited state of a fluorescent probe molecule located in the membrane. Pyrene was chosen as a probe because its intrinsic fluorescence lifetime is long. Every collision between excited-state pyrene and oxygen causes quenching. Pyrene also has the added advantage of being a lipophilic hydrocarbon that readily partitions into the lipid phase of the membrane. Also the rate of diffusion of pyrene through membranes is approximately 2 orders of magnitude slower than that of oxygen, which should produce negligible effects on the measurement of the oxygen diffusion coeffcient.9* In previous studies, pyrene derivatives have been used to quantitate oxygen in proteins99and in rat liver cells." In the study of Fishkoff and V a n d e r k ~ o ithe ~ ~ fluorescent lifetimes for pyrene were measured in differing phospholipid vesicles or membranes as a function of oxygen concentration and temperature. The diffusion constant of O2in phospholipids agrees closely with that in red blood cell (RBC) membranes and is close to literature values in water or tissue. The diffusion barrier, or energy of activation (Ea),was found to be under 3 kcal/mol for oxygen diffusion in all membranes, which suggests that the activation barrier imposed by the membrane is not large.96 The experimental values of the oxygen diffusion constants are also shown in Table IV. A few interesting observations can be noted from this study. First, the diffusion constant of oxygen is approximately an order of magnitude smaller than that of oxygen in the gas phase. This may be explained by noting that, in the process of diffusing through a membrane, an oxygen molecule may experience areas of higher viscosity than that seen in a homogeneous gas phase. Second, this study postulates that oxygen diffusion in a membrane is isotropic with respect to direction. If O2had a lateral (perpendicular to the membrane normal) diffusion rate that differed from the transverse (parallel to the membrane normal) diffusion rate, the decay of pyrene fluorescence would be described by more than one exponential. However, only one exponential was seen, so the assumption that oxygen diffusion in a membrane is isotropic was supported. Last, an increase of the rate of diffusion, by a factor of about 3 or 4, was noted when the membranes were heated from the gel phase (designated LB) through the transition temperature to the liquid crystal phase (designated La). For a finely controlled biological system as that which exists in lung surfactant, it may be advantageous to keep the phospholipids in the more fluid crystal phase by the addition of cholesterol, which broadens the phase transition temperature.95,101-103 In a study by Windren and Plachy'O" in 1980, an ESR spin exchange measurement of the oxygen diffusion coefficient determined in phosphatidylcholine bilayers was found to be in good agreement with those obtained from fluorescence quenching measurements within the experimental error.96
McKinnon et al. Molecular dynamics (MD) simulations have been employed to simulate oxygen transportin the protein myoglobin, which serves to store oxygen in muscle tissue. The heme group of myoglobin is protected from water intrusion by a tightly packed globular protein, known as globin. Water is detrimental in that it would oxidize the heme iron from a ferrous to a ferric state and render it unable to bind oxygen. The packing is so dense that if the myoglobin proteins were fixed in the static configuration, as in that found by X-ray crystallographic determinations, the oxygen would be unable to get in or out of the myoglobin. Molecular dynamics studies have shown that atomic fluctuations lower the energy barriers, preventing oxygen from diffusing into or out of the heme pocket of myoglobin.46 This more realistic dynamic approach lowers the energy barrier from 100 kcal/mol (as derived from X-ray data) to about 14 kcal/mol. It was shown that the rotation of the side chains of three amino acids accounted for a significant part of the energy barrier lowering effect. In a study by Case and M c C a m m ~ nin ' ~ 1986, ~ oxygen diffusion in and out of myoglobin was studied by a free-energy method. In this study, the dynamics of two amino acid side chains were shown to dominate the nature of the kinetic energy barrier. A comparison of the free energy profile with the potential energy profile suggested that entropy effects are dominant. S i c e entropy can be related to the number of codigurational states accessible to the protein side chains, it was crucial to incorporate a simulation method that allowed the local motions of involved proteins. These simulations show the molecular dynamics method to be a useful tool in the analysis of oxygen transport through various macromolecules of interest.
Experimental Section Transverse diffusion constants can be calculated in a variety of ways from a molecular dynamics trajectory. The diffusion coefficient may be evaluated either from the mean square displacement of the particle along its trajectory or from the decay of the velocity autocorrelation function.'" Both of these methods would require a large number of simulations to arrive at a convergent value for the diffusion coefficient, and it would be prohibitively expensive in terms of computer time. Another difficulty with the straightforward approach is that the monolayer, particularly with the inclusion of cholesterol, is a nonhomogeneous environment. It is not expected that the diffusion rate will be identical in all regions of the system, and the actual rate of oxygen transport will be determined by a weighted average of the diffusion coefficient in all of the regions of the monolayer. In this study, we have chosen to use nonequilibrium molecular dynamics (NEMD)" to derive reasonable diffusion coefficientsfrom a finte number of simulations. The approach used a novel variant of a method used by Ciccotti and Jacuccilo7in which drift velocity is examined at steady state when influenced by a constant applied force. The method used here applies a harmonic restraint force where the minimum position changes at constant velocity and the work required to achieve this velocity is calculated and related to the diffusion constant. The Einstein equation relates the diffusion coefficient to the friction coefficient via D = keT/f where kBis Boltzmann's constant, Tis the absolute temperature, and f is the friction coefficient. The friction coefficient is related to the frictional force by the Langevin equation d(u)/dt = (Ffr) - f ( ~+) (W))
(2)
where ( u ) is the average velocity, (Ffr)is the frictional force averaged over the simulations,f is the friction coefficient, and ( k ( t ) )is the averaged random force due to the solvent. We can simplify this expression by taking the average value term-by-term. The average value of the random force is by definition zero. Also, noting for this simulation that the oxygen molecule has a zero average acceleration, eq 2 reduces to (Ffr) = f(0)
(3)
The Journal of Physical Chemistry, Vol. 96, No.25, 1992 10501
Oxygen Diffusion through Hexadecane Monolayers where we have dropped the vector notation since we are concemed with one-dimensionaldiffusion. Solving for f and substituting eq 3 into eq 1 yields
D = k ~ T ( u/)(Fr,)
(4)
The relations ( u ) = d / r and
(Ffr)= W/d
(5)
where d is the distance the diffusing particle travels in time t and W is the work required to force the particle along the trajectory can be substituted into eq 4, yielding a final equation for the diffusion coefficient using quantities derived from the molecular dynamics simulations:
D=-
kBTdz Wt
where D is the diffusion coefficient, kB is Boltzmann's constant, Tis the absolute temperature, d is the membrane thickness, W is the work done to force the oxygen through the monolayer (measured as the free-energy difference derived from the umbrella sampling procedure), and t is the time taken to traverse the monolayer (30 or 150 ps for all simulations). The use of small perturbations to determine self-diffusion coefficients has been described previously by Ciccotti and J a ~ u c c i . 'However, ~~ these derivations utilized correlation functions, as opposed to the free-energy method used in this work. To drive the oxygen across the membrane a simple harmonic restraining potential was used of the form
(7) where ku is the harmonic force constant, z(t) is the present z coordinate of the oxygen molecule as defined by the center of the bond between the two oxygen atoms at time t, and zco,(t) is the center of the restraining potential as it moves along the reaction coordinate, which is normal to the monolayer plane (z coordinate), given by
+ (Zfinal + zinitiaJt/T (8) where t is the current time, Tis the length of the simulation, zinitial is the initial restraint position above the monolayer plane, and zAnaI is the final restraint position below the monolayer plane. There were no restraints in directions parallel to the monolayer plane. To increase the probability of sampling configurations leading to diffusion of the oxygen moiety, the added potential, eq 7, is systematically shifted in a series of windows centered at different positions along the diffusion coordinate. These window samples were then corrected to remove the biasing effect of the umbrella potential.17J8 With the NEMD approach used in this study of diffusion, the simulation lengths are bounded. The speed of the oxygen in the monolayer must not be so fast as to be in the turbulent regime. The Reynolds number is the ratio of the internal to viscous forces and equals the velocity times the density times the characteristic length divided by the viscosity. For estimations based on a characteristic length comparable to the radius (2 A) and using standard values for the other quantities, it is seen that for a 3 0 9 simulation the Reynolds number is much less than one. Simulations where the oxygen is transported across the membrane in less than 5 ps should exhibit turbulent behavior. This is also consistent with the idea that the velocity relaxation time (l/y) 0.02 ps, implying that there is almost no inertial motion. The method also assumes that the applied velocity is greater than the terminal velocity so that U, ( u ) >> k ( t ) relative to eq 2. Assuming that (z2) = 2Dt, where z = 12 A and D = 2.3 X 1C5cm/s, the characteristic traverse time is roughly 300 ps. The simulation lengths of 30 ps were chosen after extensive experimentation of run times of 150 ps, where it was found that much of the motion was spontaneous and did not reflect motion in the terminal velocity regime. The coordinates for the molecules used in this work were determined either from the building routines of CHARMM (as in zcon(2)
-
Zinitial
the case of molecular oxygen and hexadecane) or from X-ray crystallographic data used in conjunction with the CHARMM building facility (as with cholesterol). The unit cell was then constructed with 19 hexadecane chains in the runs with 0% cholesterol. The samples with differing concentrations of cholesterol were obtained by replacing a given hexadecane chain with a cholesterol moiety. Finally, a monolayer was constructed by replicating the unit cell into six identical nearest neighbors in the (xy) plane with hexagonal symmetry. Molecular oxygen was initialized above the (x,y) plane with a random coordinate and velocity. High-energy van der Waals contacts were then removed by energy minimization. The bottom layer of atoms was constrained by a spherical harmonic potential and the molecular oxygen was rigidly constrained for the minimization. The system was then minimized with 200 steps of adopted basis Newton-Raphson minimization (ABNR). After the energy was minimized, the system was slowly heated to 300 K and the coordinates of the unit cell were written to a file to be used in subsequent dynamics runs. The heating was accomplished by velocity scaling with 10-deg temperature increments with 10-fstime intervals between each heating increment to allow for relaxation of the structure and distribution of the kinetic energy throughout the configuration. The energy was computed for a given configuration by using the standard CHARMM force field. Because the hexadecane chains were built from atoms of similar electronegativities and the partial charge on the atoms of the oxygen molecule are zero no electrostatic term was included in the force calculation. Also, no waters of hydration were included in the computation. Before each simulation in which properties of the system were computed the system was equilibrated for 20 ps. In the computation of the nonbonded interactions a 16-Acutoff was employed. A procedure developed especially for this project calculates the order parameters of the hexadecane molecules during the simulation. The coordinates of the three carbon bonds that describe the order parameter were taken from the MD trajectory, and the (x,y,z) vectors connecting the atoms were calculated. The plane containing the two hydrogens and the center carbon atom was calculated by simple vector algebra. Because the order parameter relevant to this simulation only involves the z coordinate, the calculation of the order parameter was simplified. The order parameter for each of the carbons, from carbons 2 to 15, was averaged and returned when the order parameter calculation was performed. The average fluctuation for each carbon and hydrogen atom was also computed. The simulated system is not large by present standards. This system includes 19 lipid molecules and the long dimension of the repeating hexagon is roughly 52 A. The nonbonded cutoff for electrostatics and van der Waals interactions is 16 A,which is large by current molecular dynamics simulation standards (but still small relative to the hexagon dimensions), so as not to include artificial ripples in the potential for long-range interactions that could adversely affect the results. Since the diffusion coefficient is dominated by short-range interactions, a significant increase in system size should not significantly alter the results. Size considerations are not likely to be a major factor. All simulations were performed at the National Institutes of Health, Division of Computer Research and Technology, Molecular Graphics and Simulation Section, by using a modified version of CHARMM.88 The simulation was run with a timestep of 1 fs, with atomic positions saved every 0.1 ps (100 timesteps) for analysis. The temperature was maintained at 300 K within a 15-Kwindow by periodic velocity scaling. The total computer time for this simulation project was 2450 cpu h on an Apollo DN-10040(roughly equivalent to 250 cray X/MP h). Results and Discussion chain Tilt. The chain tilts as measured by the average angular deviation from the membrane normal are tabulated in Table I and illustrated in Figure 1. There is a decrease in hexadecane chain tilt as a function of increasing cholesterol concentration. As previously discussed, the use of periodic boundary conditions
losO2 The Journal of Physical Chemistry, Vol. 96, No. 25, 1992 TABLE I: Hexadecane Chain Tilt as a Function of Cholesterol CoaeCnbrtiOa cholesterol hexadecane hexadecane cholesterol concn, mol % chain tilt, deg concn, mol % chain tilt, deg 30 25.2 36.5 0 10 21.8 40 23.1 20 30.2 50 12.6 TABLE Ik Monolayer Thickwas and Hexadecane Carbon Volume as a Funeti~aof Cholesterol Concentration cholesterol monolayer hexadecane carbon vol, concn, mol % thickness, A A'fcarbon 0 12.5 25.6 10 13.8 28.3 20 14.0 28.7 30 14.3 29.5 40 15.2 31.4 50 15.9 33.0
McKinnon et al. Hexadecane Length PI. Mole X Choleaterol
I
9 3
a-
B 3 H
"
I
'
I
I
I
16.0
15.0
& ; j: ] 0
10
20
30
40
50
Mole X Cholesterol
Figure 2. Hexadecane chain length (x) as a function of cholesterol concentration.
-
35
- ln
30';
8
:
-.- 25; ? I
.-c
2
-
-
f
20-
0
-
15
0
l o " ~ " " ~ " " " " ' l " " l " " l ' 0 10 20 30
40
50
Cholesterol Concentration (mole %)
Figure 1. Chain tilt (0)as a function of cholesterol concentration. The line is to guide the eye.
can increase the amount of tilt seen. However, one could argue that these effects are constant over all of the simulations and the increase in tilt is due primarily to the presence of cholesterol. The use of a small total number of molecules (19) and periodic boundary conditions may produce an artifact that increases chain tilt relative to the monolayer normal. The phenomenon of chain tilt with a small simulation has been noted before by van der Ploeg and Berendsen in 1982.74 With an increased chain tilt, the order parameters will be, on the average, decremd. One might question whether the diffusion of oxygen would be a f f d if the modeled structure was tilted too much. However, it is postulated that the diffusion of oxygen through diphosphatidylcholine (DPPC) is isotropic,96 and thus its diffusion rate would be independent of whether it moves laterally in the membrane or diffuses transversely through the membrane. Therefore, one might expect that the diffusion coefficient of oxygen would not be significantly affected with increased chain tilt. Effect of Cholesterol on Membrane Thickness. Membrane thickness increases in a roughly linear fashion as a function of cholesterol concentration, as shown in Figure 2. The fluctuations of the monolayer thickness over the time of the simulations are also shown in Figure 2. The effect that cholesterol exerts on the monolayer, namely, an overall increase in order parameters, is reflected in the increasd thickneas. Thii increased thickness has the effect of increasing the acccssible volume per hexadecane carbon in the monolayer, as tabulated below in Table 11, and as shown graphically in Figure 3. The average end-to-end length of the hexadecane chains in each of the simulations are compared to determine if increasing cholesterol concentration affects chain length (as opposed to monolayer thickness, which depends on chain tilt). The hexadecane chain lengths at 0 and 10 mol % cholaterol are 15.5 and 15.6 A, respectively. The chain length reaches a maximum of 16.5 A at 10 mol % cholesterol, which is approximately a 6% increase
Cholesterol Concentration (mole %)
Figure 3. Hexadecane carbon volume (0)as a function of cholesterol concentration. The line is to guide the eye.
in chain length over the 0 mol % chain length. Cholesterol will affect an increase in hexadecane chain length only if the average number of trans configurations is increased. The increase in hexadecane carbon volume with increasing cholesterol ooncentration may allow diffusion paths that have more accessible volume and thus lower the work required to diffuse through the monolayer. This may be important in the mechanism of increased diffusion coefficients of oxygen with increased cholesterol concentration. On the other hand, a thicker layer may increase the time required for oxygen to diffuse across the monolayer. The fluctuations shown in Figure 2 are small compared to the overall thickness (2%). These small deviations in length show that the density fluctuations are also small and that cholesterol exerts a strong and constant ordering effect on the monolayer. Order Panmeters. The order parameter is a measure of the anisotropy of the carbondeuterium bonds of the long-chain hydrocarbons being studied. The order parameter Sd is a direct measure of the deviation of the plane containing the carbon and both attached deuterium from perpendicularity to the lipid plane normal. The order parameter is defined as
sd=y 2 ( 3 ~ 2 e - i )
(9)
where 8 is the angle between the bilayer normal and the vector bisecting the D c - D bond, which lies in the plane defined by the carbon and its attached deuteriums, and the brackets denote the average over all configurations. Therefore, the order parameter can vary from in a fully ordered geometry (e = 90°) to +1 at the most disordered geometry (e = 180O). For example, if the hydrocarbon chain were in the fully extended all-trans configuration, the planes defined by the carbon and attached deuteriums would all be perpendicular to the bilayer normal, and all order parameters would have the identical value of -1/2. It is conventional to report the order parameters in terms of the methylene segment molecular order parameter Smol, which is
The Journal of Physical Chemistry, Vol. 96, No. 25, 1992 10503
Oxygen Diffusion through Hexadecane Monolayers geometrically related to the S, by the following equation:89
= -2s, (10) The order parameter can be easily derived computationally from a molecular dynamics trajectory, as explained in ref 28. Experimentally, the order parameter is obtained from 2H NMR measurements. The long-chain hydrocarbon is selectively deuterated at the desired carbon atom position, and the spectrum is recorded. The spectrum is characterized by two sharp lines corresponding to the two allowed NMR transitions of the deuterium Z = 1 nucleus. The separation between the two ftsonances is the midual deuterium quadrupole splitting, AuQ which provides information about the average orientation and fluctuation of the C 3 H bond vector. The quadrupole splitting is related to the order parameter through eq 589 Smol
AVQ= (3 /4)(e24Q/h)&
(11)
where (eZqQ/h)is the static quadrupole coupling constant, which ' bond.9l is 170 kHz for an aliphatic CH In this manner, direct comparisons between experimental and computational measurements of membrane structure can be performed. Order parameters are a sensitive measure of the structure of lipid aliphatic chains and have been measured for dipalmitoylphosphatidylcholine (DPPC) and various other lipids with and without the inclusion of cholesterol. Nonsonicated bilayers of DPPC were deuterated at nine different positions in a 2HNMR study by Seelig et al. in 1974.92 This study showed that the order parameters were constant for the first nine chain segments and then decreased over the remainder of the chain. The chain ordering was explained on the basis of the rotameric model for hydrocarbon chains. In the region of the constant order parameter, gauche conformations only occur in complementary pairs, which would leave the chains essentially parallel to each other. This leads to a well-ordered bilayer with disordered ~hains.9~ A 2H NMR study in 1978 by Oldfield et al. of a dimyristoylphosphatidylcholine and 30 mol 96 cholesterol mixture evaluated the order parameters describing motions of the bilayer.% The addition of cholesterol to the bilayer brought about an overall increase in the order parameters as compared to the above study. The authors showed that the chain lengths were very sensitive to the presence of cholesterol, and the addition of 30 mol 96 cholesterol increased the membrane thickness by 2196, which agrees with the increased values of the order parameters. Another 2H NMR study of dimyristoylphosphatidylcholine in concentrations of cholesterol varying from 0 to 50 mol 96 was performed by Jacobs and Oldfield in 197gSg5The quadrupole splittings (which are directly proportional to the order parameters) were found to increase in a roughly linear fashion with increasing cholesterol concentration, denoting an increase in order and membrane thickness with the addition of cholesterol. This study implies an increase in the number of trans conformations with increased cholesterol concentration. Seelig and Macdonald in 1987 showed an interesting order parameter profile in an *H NMR study of a 1:1 mixture of cholesterol/lysopalmitoylphosphatidylcholine,m Not only were the overall values increased over the pure DPPC values by almost a factor of 2, but the order parameters were markedly increased toward the middle of the long hydrocarbon chain. These results could be explained by the presence of the rigid ring of the cholesterol molecules, which strongly hinders chain flexibility in these regions. The moluxlar order parameters have been calculated for each of the simulations from 0 to 50 mol 46 cholesterol. The values of the molecular dynamics order parameters are tabulated along with the experimentally derived deuterium NMR order parameters in Table 111. Shown in Figure 4 are the MD order parameters in 30 mol 96 cholesterol with the deuterium NMR order parameters of 30 mol %.m The MD order parameters at all cholesterol concentrations in these simulations are in close agreement with the deuterium NMR order parameters. The order parameter of atom 2 is generally high most likely due to the added constraint of main-
TABLE IIk ExperiwaW and Molecular Dynrmics Order Parameters 2H C no. 0% 30% 50% 0% 10% 20% 30% 40% 2 0.47 0.538 0.367 0.289 0.330 0.393 0.461 0.326 3 0.43 0.436 0.381 0.499 0.510 0.378 0.734 4 0.47 0.463 0.718 0.385 0.440 0.512 0.498 0.445 5 0.46 0.481 0.382 0.417 0.511 0.510 0.520 6 0.498 0.811 0.308 0.374 0.496 0.504 0.581 7 0.499 0.271 0.376 0.503 0.510 0.684 8 0.500 0.227 0.410 0.443 0.475 0.674 9 0.45 0.493 0.764 0.236 0.395 0.422 0.467 0.668 10 0.43 0.486 0.267 0.423 0.382 0.452 0.616 11 0.446 0.275 0.357 0.348 0.434 0.588 12 0.35 0.406 0.577 0.285 0.313 0.347 0.393 0.556 13 0.337 0.280 0.345 0.347 0.378 0.510 14 0.30 0.267 0.258 0.334 0.330 0.330 0.459 15 0.21 0.214 0.307 0.278 0.182 0.406
50%
0.681 0.734 0.733 0.653 0.692 0.797 0.777 0.751 0.603 0.513 0.476 0.453 0.433 0.349
a 2HN M R order parameters of pure dipalmitoylphosphatidylcholine (DPPC) at the phase transition temperature of 314 K.92b2HNMR order parameters of a dimyristoylphosphatidylcholine/30mol % cholesterol mixture at 296 K." C2HNMR order parameters of a lysopalmitoylphosphatidylcholine/50 mol 5% cholesterol mixture at 3 18
K." Order Parametera of Hexadecane in 30 Mole 0.71 0.8
? 0.3 o.2
1
I
I
1
1
1
1
I
I
I
I
2
Cholesterol
I
calc. (14))
expt. (30X)
I
x-
I
--
0-
1 t t
0.0'
I
I
2
3
4
I
I
I
I
5
8
7
8
I
'
I
I
'
I
I
Q101112131415
atom number
Figure 4. Order parameter as a function of atom number along the hexadecane chain for 30 mol % cholesterol. (x) Denotes the MD values; (0) are the experimental values.
taining the head group in the plane of the membrane. The close agreement of the order parameters obtained from the molecular dynamics simulations with the experimentally derived deuterium NMR order parameters, especially from 20 mol 96 up to 50 mol 96 of cholesterol, represents an affirmation of the monolayer model used in these simulations. The accurate portrayal of the structural characteristics of the model monolayer lends credence to the MD derived results of the oxygen molecule diffusion coefficient. In Figure 5 , the molecular order parameters as a function of chain position and cholesterol concentration are shown. The order parameters of hexadecane carbon atoms 6 and 9 increase in a roughly linear fashion with increasing cholesterol concentration. This is easily explained by the presence of the rigid sterol ring of the cholesterols, which increases the number of trans configurations in the lower region of the monolayer, thereby lengthening the chains and increasing the monolayer thickness. The order parameters of atoms 12 and 15 do not increase as significantly with increasing cholesterol concentration. The aliphatic tail of the cholesterol may resemble the tail of the hexadecanes in this region. Cholesterol, with its unique structure composed of rigid and fluid components, not only adds order into the lower segments of the monolayer but also maintains a highly active tail region. The molecular order parameter fluctuations, labeled S,+, are shown in Figure 6. These values represent the changes in the order parameters averaged over the time of the simulation. The fluctuations do not increase as a function of cholesterol concen-
10504 The Journal of Physical Chemistry, Vol. 96, No. 25, 1992 Order Parametera
VI.
McKinnon et al.
Mole X Cholesterol
Cholenterol Plane
:8::
2
I
z
u
K
0.7
0.6 0.8
- Plane Correlations
-
I
I
I
I
l
1
I
I
I
I
I
I
I
I
I
I
I
I I
I I
I I
I I
I I
I I
I I
I I
I I
-
2
I
U
0
n I
2
V
1 /1
gni 0.2
W
t V
0.2
0
20 30 40 Mole X Cholesterol 10
I
::;-1 j 0.2
2
s!
I
n
g
0.1 0.1
-
0.3
-
0.2 0.2 0.1 0.1
I
I
VI. I
I
0.7 0.8
I
I
I
I
I
I
I
I
I
2
0.6
7
1
Mole 71 Cholesterol I
p
-
K - r ] I
1.0 r
.2 0.8
50
Figure 5. Order parameter for carbon atoms 6, 9, 12, and 15 of the hexadecane chain as a function of cholesterol concentration. Scd Fluctuations
0.4 0.3 0.6
2
3
4
5
8
7
8
8
10
Time in P I
Figure 7. Cholesterol plane-plane correlation as a function of time for several cholesterol concentrations.
I
20
I
::
~
I
I
I
I
I
I
I
I
I
I
I
I I
I I
I I
I I
I I
I I
-
-
0 Mole I
0
0 0
20 30 40 Mole X Cholesterol
10
50
Figure 6. Fluctuation in the order parameter as a function of cholesterol concentration.
tration but instead seem to increase only as a function of upward position on the hexadecane chain. These reorientations of the upper regions reflect the higher number of accessible geometries being sampled by the tails of the hecadecane chains. cbdcsterd PhnePlme Correlations. The cholesterol planeplane correlations are displayed in Figure 7 for cholesterol concentrations varying from 10 to 50 mol %. The two atoms selected are located at the edge of the unit cell and as close to the center of the unit cell as possible, and the identical ( x y , z ) vectors spanning the sterol rings are used for the correlation. In all but one of the concentrations (20 mol a), there is an extremely high correlation among the selected cholesterol planes. What is most interesting is the high correlation at 10 mol %, where only two cholesterol molecules are present in the unit cell. They are far enough apart that a direct steric interaction can be ruled out. A probable explanation for the high correlations is that the time necessary for rotational diffusion is much longer than the 30-ps time length of the simulation. Dynamic Properties. Free-brgy Calculations. Free-energy perturbation simulations using an umbrella sampling technique
20
50
20
50
Time la pa
Figure 8. Work required to drag the oxygen molecule through the monolayer as a function of time for each of the cholesterol concentrations employed in this study. Each of the trajectories used in the calculation of the diffusion coefficient is shown.
were performed with the restraint potential given by eq 7. The relative freeenergydifference plots for the simulations are shown in Figure 8. For each concentration of cholesterol, five or six simulations were performed. For plots of oxygen diffusion through no cholesterols, 19 hexadecanes represent the simulation at 0 mol %; for plots of oxygen diffusion through two cholesterols, 17 hexadecanes represent the simulation at 10 mol % for plots of oxygen diffusion through four cholesterols, 15 hexadecanes r e p resent the simulation at 20 mol %; for plots of oxygen diffusion through six cholesterols, 13 hexadecanes represent the simulation at 30 mol% for plots of oxygen diffusion through eight cholesterols, 11 hexadecanes represent the simulation at 40 mol 96; for plots of oxygen diffusion through ten cholesterols, nine hexadecanes represent the simulation at 50 mol %. In actuality, the exact concentrations of cholesterol are 0, 10.5,21.1,31.6,42.1, and 52.6 mol k, but in this manuscript they are referred to as 0, 10, 20,
Oxygen Diffusion through Hexadecane Monolayers
The Journal of Physical Chemistry, Vol. 96, No. 25, 1992 10505
TABLE Iv: Molecular Dynamics Oxygen Diffusion Coefficients Compared with Thoee Obtrined by Qwnched Pyrene Fluorescence molecular dynamics quenched pyrene fluorescence cholesterol diffusion cholesterol diffusion concn mol % coefficient, cm*/s concn mol 96 coefficient, cm*/s 2.3 (i0.2)x lo-’ 0 2.6 (t0.2) x 10-5 00 10 3.6 (t0.3) X IO-’ 20 3.3 (f0.3)X lo-’ 29 3.2 (f0.2) x lo-’ 30 4.1 (i0.9)X lo-’ 40 4.5 (i0.5)x lo-’ 50 3.9 (f0.7)x lo-’ 5oc 2.4 (~0.1)x 10-5 LI In dipalmitoylphosphatidylcholineat 298 K.% In 3:1 dimyrisotylphosphatidylcholine/cholesterol at 298 K.% In erythrocyte plasma membrane at 298 K?5
30, 40,and 50 mol 9%. The free energy represents a difference and is plotted relative to zero. The general pattern of the plots shows a low freeenergy bamer to diffusion during the initial 10-1 5 ps of the simulation. In reality, this monolayer would be part of a larger system, so the regions of plots that show a negative work are relative to the system. This cormponds to the easy movement of the oxygen molecule through the regions of the hydrocarbon tails, which have a high degree of disorder and molecular motion. This disorder is also reflected in the lower order parameters in the tail regions. Also, the attractive portion of the Lennard-Jones potential predominates at this point with few hard-sphere repulsions. As the oxygen molecule continues to diffuse through the monolayer, it encounters the lower half of the monolayer. The chains in this region are in a more rigid configuration, which is reflected by higher order parameters. In this region, the slope of the free-energy difference, which represents the force encountered by the oxygen molecule, is larger. The slope in this region is lower as the mole percent of cholesterol increases, up to 40 mol 96 (eight cholesterols, 11 hexadecanes). At 50 mol 96, the slopes are steeper on average. It is reasonable to assume that the presence of the rigid sterol ring of the cholesterol forces the lower segments of the hexadecanes into a more extended configuration, and this extended configuration reduces the work required to move the oxygen molecule through the monolayer. Diftuaon COetRcients. The oxygen diffusion coefficients have been derived from the nonequilibrium molecular dynamics procedure according to eq 6 for concentrations of cholesterolvarying from 0 to 50 mol 96. The derived values are compared to the oxygen diffusion coefficients for 0,25, and 50 mol 9% cholesterol obtained from quenched pyrene fluorescence108in Table IV. The MD obtained value of 2.6 X cm2/s for 0 mol 9% cholesterol compares quite favorably to the value of 2.3 X cm2/s obtained from quenched pyrene fluorescence.’Og This MD value represents only a 12% increase over the quenched fluoreucence value. When compared to previous simulations, where the diffusion coefficient of water was calculated to be 1009%faster this simulation method calculates than that seen the diffusion coefficient with remarkable accuracy and efficiency. The average values for the MD calculation show an enhancement of the diffusion coefficient over all cholesterol concentrations, with a broad maximum at roughly 40 mol %. This maximum represents an approximately 2-fold increase of the diffusion coefficient over the MD value (0 mol %) and quenched pyrene diffusion coefficient for pure dipalmitoylphosphatidylcholine (DPPC).A plot of the M D and quenched pyrene fluorescence derived diffusion coefficients is shown as a function of cholesterol concentration in Figure 9. The MD derived values are fit to a splined curve, which shows an inmasing diffusion coefficient with increased cholesterol concentration. The standard errors are calculated and also plotted. It is interesting to note the variation in the trajectories depicted in Figure 8 as a function of cholesterol concentration. At low concentration of cholesterol (
