INTRODUCTION

TOP ABSTRACT INTRODUCTION METHODS RESULTS CONCLUSIONS APPENDIX REFERENCES

The main goal of the present article is to evaluate atomistic molecular dynamics (MD) simulations as tools for investigating monolayers of zwitterionic and anionic phospholipids. The understanding of the molecular level phenomena in lipid monolayers is important for providing insight into the behavior of biological membranes and for conferring the design certainty that is necessary to engineer molecular solutions for diagnosis and treatment of systems involving protein-lipids interactions. These systems are at the heart of central biological processes, being involved in such diverse phenomena as in sensory stimulus, the intracellular communication, viral fusion, regulation of nutrient and waste flow in and out of the cell, and sperm/egg fusion. Lipid monolayers at the air/water interface have been widely and successfully used as models to approximate the behavior of cellular membrane lipid bilayers with the underlying assumption that a bilayer can be considered as two weakly coupled monolayers (Brockman, 1999; MacDonald, 1996; Marsh, 1996; Feng, 1999). Besides their attractive reduced dimensionality, use of monolayers enables the direct control of the temperature and the surface pressure or its thermodynamic conjugate property, the surface area per molecule at the interface. Such control, accomplished using a Langmuir trough and not possible in bilayers, allows the unambiguous investigation of the interaction of membrane-active molecules, such as proteins, with the lipid matrix. For example, measurements of changes in the surface pressure at constant area or changes in the area of the monolayer at constant surface pressure, upon insertion of protein molecules in the monolayer provide a clear picture of the strength of protein-lipid interactions (Baszkin et al., 1999; Maget-Dana, 1999).

Moreover, there has been increased interest in the development and employment of sophisticated biosensor techniques based on lipid monolayers immobilized on solid supports for investigating biomolecular interactions. These techniques greatly facilitate the visualization and probing of association phenomena between biological molecules (Mrksich and Whitesides, 1996; Mozsolits et al., 1999; Nyquist et al., 2000).

In addition, of primary physiological importance is the lipid-protein mixture monolayer that lines the alveolar space of the lung (Clements and Avery, 1998; Johansson and Curstedt, 1997). This monolayer facilitates the proper respiratory function in vivo by reducing the surface tension at the interface between the air and the mucus layer that hydrates the lung surface epithelial cells (the alveoli).

Thus, it becomes evident that a clear picture of molecular level phenomena in lipid monolayers is important for improving our understanding of the underlying physicochemical principles determining the behavior of biological systems and for developing design rules for molecular products and processes of biological importance.

Computer simulations can provide such a picture. They have been extensively used for investigating lipid bilayer membranes and the interaction between lipids, proteins, ions, and small organic molecules (Mertz and Roux, 1996; Jakobsson, 1997; Tieleman et al., 1997; Forrest and Sansom, 2000). Recent MD simulations have modeled systems of hundreds of lipids and investigated time scales of tens of nanoseconds (Lindahl and Edholm, 2000), extending the simulated time and length scales by one order of magnitude. Such simulations shed light on mesoscopic phenomena and allow for the development of multiscale simulation methodologies.

However, in the vast majority of the efforts up-to-date, biological membranes have been modeled as bilayers of zwitterionic lipids, primarily dipalmitoylphosphatidylcholine (DPPC). Only a limited number of efforts have been attempted to simulate lipids that are ionized at physiological pH (Cascales et al., 1996; Bandyopadhyay et al., 1998), although numerous cell membranes contain considerable quantities of anionic and/or cationic lipids with important structural and functional roles (Yeagle, 1993).

Monolayers of lipid molecules have been traditionally simulated using simple coarse-grained models (Karaborni and Toxvaerd, 1992; Ahlstroem and Berendsen, 1993; Esselink et al., 1994; Mauk et al., 1998), and early atomistic simulations were limited to sampling the NVT ensemble, where the number of the atoms (N), the volume of the system (V), and the temperature (T) were kept constant (Alper et al., 1993). An alternative approach to simulating interfacial systems was recently proposed by Wong and Pettit (2000). These authors introduced a new boundary condition with the simulation box used being the asymmetric unit of space group Pb.

In a remarkable contribution, Feller and co-workers successfully modified the Andersen extended system approach (Andersen, 1980) and developed a MD algorithm that enables the sampling of lipid monolayer statistical ensembles at constant surface pressure or constant surface density (Feller et al., 1995a,b; Zhang et al., 1995). This allows for the direct comparison between simulations and surface pressure/tension isotherm measurement experiments.

The drawback of the approach developed by Feller and co-workers is the chosen geometry for simulating lipid monolayers at the air/liquid interfaces. Limited by periodic boundary conditions considerations, Feller and co-workers chose the geometry shown in Fig. 1 A, essentially simulating two lipid monolayers separated by a layer of water. The simulation cell is effectively periodic in the lateral x and y directions, because in the z direction the distance between the lipid tails of the one leaflet and the periodic image of the second leaflet can be defined to be arbitrarily larger than the range of interatomic forces. In addition, care has to be exercised in choosing the thickness of the water layer so that the lipid head groups of the opposing monolayers do not interact. Feller and co-workers used 36 DPPC molecules for each monolayer and 2511 water molecules. The layer thus formed was ~20 Å, and little interaction was anticipated between the two surfaces.

View larger version (21K): [in this window] [in a new window]   FIGURE 1   Schematic representation of simulation geometries employed by Feller and co-workers (A) and used in the present work (B).

However, this choice of geometry renders computationally very expensive all simulations that require systems with a large lateral area. Such systems are useful when long-range interactions are important or when it is desired to investigate large-scale phenomena. For example, when interactions between protein molecules and the monolayer are the research objective, it is important for the lipid layer matrix to accommodate the protein molecule(s) in a way that prohibits direct interactions between the protein and its periodic images. In such instances where large numbers of lipids are required per monolayer, simulation of two independent monolayers might be unnecessarily cumbersome and computationally expensive.

Moreover, when using Feller's method to simulate small proteins in lipid monolayers for long enough times to properly capture the nature of protein-lipid interactions, there is the danger of the protein molecules of the opposing monolayers diffusing and interacting with one another in the bulk. To avoid such instances, an unnecessarily thick layer of water should be used. Recent simulations of ours have resulted in a protein with 25 amino acids diffusing up to 10 Å away from the interface in the course of a 3-ns simulation, having started with the protein positioned at the interface. Also, simulations that investigate the diffusion of small proteins to the lipid/water interface and do start with the protein molecule in the bulk become problematic. Inordinate amounts of water molecules would be required to assure the protein molecules do not interact. The methodology that Wong and Pettit (2000) introduced has the same drawbacks, because protein molecules might interact with their mirror images in the direction normal to the monolayer's interface.

In this study we report on the results of simulating lipid monolayers at the air/water interface using an approach that is largely based on the one developed by Feller and co-workers but avoids the drawbacks of their system geometry. Simulation of the same number of lipids and water molecules, hence similar computational cost, results in doubling of the simulated area with the proposed methodology. This is achieved by using the geometry depicted in Fig. 1 B. Instead of two opposing monolayers, a wall potential is applied that disallows the diffusion of water molecules in the z direction. A similar wall was also used by Alper and co-workers in their simulations of DPPC monolayers (Alper et al., 1993). The simulation cell is periodic in the x and y directions. In the z direction the distance between the wall and the periodic image lipid tails is arbitrarily chosen to be 200 Å, significantly larger than the effective range of interatomic interactions, in a fashion similar to the one adopted by Feller and co-workers. We are interested in creating a solid framework for simulating lipid monolayers that facilitates the inclusion of biological molecules such as proteins in the lipid matrix and the investigation of numerous interesting biological phenomena.

We have simulated monolayers of DPPC and dipalmitoylphosphatidylglycerol (DPPG), an anionic lipid with the same acyl chains. DPPG lipids have attracted considerable interest, because they are important components of pulmonary surfactants (Nag et al., 1994; Veldhuizen et al., 1998), and of the thylakoid membrane of the chloroplast (Webb and Green, 1990; Fragata et al., 1997), among other systems. Pure systems of DPPC and DPPG are simulated at four different areas per lipid. As discussed in the Results section, simulation of systems at different surface densities enables comparisons to be made between the simulations and experimental measurements of surface pressure isotherms.

We also simulated mixtures of DPPC and DPPG lipids at a molar ratio of 7:3, which is the ratio between nonionic and anionic lipid species in the lung surfactant. Experimental measurements of the excess free energies of mixing suggest that under certain conditions of temperature and ionic atmosphere DPPC and DPPG mix nearly ideally (Williams et al., 1995). Even under the assumption of near-ideal mixing, we did not expect the lateral distributions to reach equilibrium due to the very long time scales of lateral diffusion. Diffusivities of phospholipids in monolayers are of the order of 10⁶ cm² s¹ or ~10 Ų ns¹ (Adamson and Gast, 1997). Our simulations reach 1 ns of equilibrated run, thus we do not expect to observe many lipids interchanging their positions on the xy plane. Nevertheless, we felt that it was methodologically useful to attempt such simulations starting with random lateral distributions of the DPPG lipids in the matrix of DPPC lipids as discussed in the Methods section.

Moreover, of interest is the structural and functional influence of cations on anionic lipids (Watts et al., 1981; Flach et al., 1993). This influence is investigated by simulating the lipid monolayers in the presence of 150 mM NaCl. In some systems, 50 mM CaCl₂ is added to specifically characterize the influence of divalent cations on the monolayer structure.

The simulation results of 24 systems allow us to construct a clear picture of the structural properties of zwitterionic and anionic phospholipid monolayers in various density and ionic environments. Our goal is to demonstrate the stability of the simulations for long time scales and provide a clear picture of the nature of the lipid/water interface depending on the lipid and the electrolyte solution.

In the next section, details of the simulation methods are presented. The simulation results are discussed in the Results section. The Conclusion section provides an assessment of the methods used as a means of investigating lipid monolayer systems.

	METHODS

TOP ABSTRACT INTRODUCTION METHODS RESULTS CONCLUSIONS APPENDIX REFERENCES

The simulated system consists of a periodically replicated cell containing 40 lipids and ~2800 water molecules. Using the geometry depicted in Fig. 1, the system is effectively periodic in the lateral x and y dimensions. The program Chemistry at Harvard Molecular Mechanics (CHARMM) is used for the simulations with the PARAM26b2 parameter set (Brooks et al., 1983; Schlenkrich et al., 1996). The water molecules are simulated using the TIP3P model (Jorgensen et al., 1983).

The lipids are randomly chosen from a library of 2000 preequilibrated and prehydrated DPPC molecules (Venable et al., 1993), generated by Monte Carlo simulations with a Marcelja mean field (Hardy and Pastor, 1994). The monolayer is constructed following largely the protocol described by Woolf and Roux (1994). The phosphorus atoms of the lipid heads are positioned randomly on the xy plane at z = 0 with the acyl tails stretching in the positive z direction. Bad contacts between the lipids are reduced by systematic translations in the xy plane and rotations around the z axis. When simulating lipid mixtures, 12 of the lipids are randomly chosen and changed from DPPC to DPPG modifying their head atoms. In the remainder of this article we use the designation "PCPG" to refer to the mixtures. All of the DPPC molecules are modified to DPPG molecules for the pure DPPG systems.

In all of the cases, four different surface densities are simulated at 55, 60, 65, and 70 Å²/lipid. The size of the simulation cell in the x and y directions is adjusted accordingly to accommodate 40 lipids. Subsequently, water is added in the negative z direction. The number of water molecules ranges from 2600 to 3000 depending on the size of the simulation cell in the x and y directions. Subsequently, 0.15 M NaCl is added with each ion taking the place of a randomly chosen water molecule. In addition, for the systems containing the negatively charged DPPG molecules, sodium counterions are added to neutralize the simulated system. In an additional set of 12 systems (pure DPPC, pure DPPG, and mixture of DPPC/DPPG at four different surface densities) 50 mM CaCl₂ was added with the 150 mM NaCl. Altogether, each simulated system contains ~12,000 to 14,000 atoms, depending on the surface density of the monolayer.

The CHARMM potential contains terms for bond lengths, bond angles, torsional angles, and improper torsional angles. Nonbonded atoms interact with a 6-12 Lennard-Jones potential. Atoms are also assigned point charges positioned at their center of mass. In the simulations, the nonbonded van der Waals interactions were smoothly switched off over a distance of 3.0 Å, between 9 and 12 Å. The electrostatic interactions are simulated using the particle mesh Ewald summation (Essman et al., 1995) with no truncation. A real space Gaussian width of 0.42 Å¹, a B-spline order of 6, and a fast Fourier transform grid of approximately one point per Å (60 × 60 × 200) are used. The SHAKE algorithm (Ryckaert et al., 1977) is used to set constraints on the water hydrogen bond lengths.

A weak harmonic potential applied at z  30 Å prevents the bulk water molecules from diffusing in the negative z direction. We found that a wall potential U_wall = 5/2(z_p  3)² kcal/mol, in which z_p is the distance from the wall (z_p = 0 Å at the wall), never allowed water molecules to diffuse, without significantly distorting the structure of the bulk water. The potential is applied only on the molecules that diffuse below the wall. The exact position of the wall potential is fine tuned so that the density of the bulk portion of waters has a value of 0.0333 molecules per ų. We found that the wall effects do not propagate for more than 10 Å, allowing for an adequately large volume of water with bulk properties. In principle, there is the possibility of water molecules diffusing in the positive z direction through the lipid monolayer to the vacuum space, but in none of the simulations was this observed. The wall potential can be considered as an effective finite representation of an infinite bulk system. Beglov and Roux (1994) developed a rigorous theoretical framework for a solvent boundary potential that incorporates the influence of the bulk on the finite simulated system. Specifically, their formulation is the result of an exact separation of the multidimensional configurational Boltzmann integral in terms of the solvent nearest the solutes and the remaining bulk solvent molecules. They developed a generalized solvent boundary for a spherical geometry, and the exact form was the result of an approximate semiempirical process. We preferred to design the wall potential in a purely empirical manner so that bulk properties would be recovered at short distances from the wall.

After constructing the simulation cell, the energy is minimized using the adopted-basis Newton Raphson algorithm for 1000 steps. Minimization is followed by a 100-ps MD heating period. The final temperature is 323.15 K. An equilibration period of 300 ps follows. All the systems are then simulated for 1000 ps. The calculations were performed simulating the canonical ensemble, i.e., the number of the atoms (N), the volume (V), and the temperature (T) were kept constant. The temperature was set at 323.15 K using the Hoover temperature control with a mass of 1500 kcal ps² for the thermal piston (Hoover, 1985). Although the volume of the simulation cell is kept constant, the presence of explicit vacuum in z direction results in the density of the system adjusting to its equilibrium value.

The equations of motion are integrated using a leap-frog integrator (2-fs time step). In practice, we used the constant pressure-temperature module of CHARMM, and by setting all the components of the piston mass array to zero, the constant volume-temperature equations of motions described by Hoover are recovered.

The system's height does relax in the z direction so that the normal pressure component fluctuates ~0 atm at equilibrium. The lateral direction pressure components obtain negative equilibrium values, as discussed in the Results section, because of the inhomogeneity of the system in the normal direction. This, in turn, gives rise to a positive monolayer surface tension, which is one of the most important properties calculated during the simulations. The surface tension calculation allows the direct comparison of the simulations with surface pressure isotherm experimental measurements. A detailed description of the surface tension calculation is given in the Appendix.

	RESULTS

TOP ABSTRACT INTRODUCTION METHODS RESULTS CONCLUSIONS APPENDIX REFERENCES

The parsed structures of DPPC and DPPG in functional groups are shown in Fig. 2. The difference is that in DPPG instead of a positively charged choline group, we find a glycerol. Thus, at physiological pH, DPPG is acidic due to the negatively charged phosphate group, whereas DPPC is zwitterionic with two groups being ionized at physiological pH with opposite charges. We should note here that it is a simplifying assumption that all the DPPG molecules are ionized in the monolayer. Indeed, Brewster angle microscopy experiments (Grigoriev et al., 1999) suggest that DPPG is not fully deprotonated, although a substantial degree of dissociation is observed. For the purposes of the simulations, the assumption of full deprotonation is warranted.

View larger version (27K): [in this window] [in a new window]   FIGURE 2   Parsing scheme for DPPC and DPPG.

The electron density profiles of the groups in the z direction are plotted in Fig. 3. The average electronic density of each group was calculated by assuming a Gaussian distribution of electrons centered at each atomic center of mass with a variance equal to the van der Waals radius. The profile was calculated by time averaging the number of electrons in slabs with height 0.1 Å.

View larger version (20K): [in this window] [in a new window]   FIGURE 3   Electron density profiles of molecular groups for a monolayer with DPPC and DPPG molecules at a surface density of 60 Å2/lipid.

As evident in Fig. 3, the hydrophilic groups of the lipids mix well with the water subphase, whereas the methyls and methylenes of the lipid tails stretch in the positive z direction. The electron density of water increases smoothly from zero at the wall to ~0.333 e/ų in less than 10 Å. This leaves a water layer of thickness ~20 Å with the appropriate bulk density. We can therefore safely assume that the wall does not significantly affect the water/lipid interfacial properties. The water electron density decreases to zero at the interface with the lipids. We can define the interfacial thickness as the distance between the point on the hydrocarbon side of the interface where the electron density of the water becomes less than 0.001 e/ų and the point on the water side where only a lipid head group has an electron density as large as 0.001 e/ų and all others have smaller densities. The thickness of the interface is reported in Table 1 for all the simulated systems.

For the DPPC monolayers the thickness of the interface increases with increasing area per lipid. Calcium appears to have no significant impact on the interfacial thickness of DPPC monolayers. Interestingly, for monolayers containing DPPG, the thickness of the monolayer decreases with increasing area per lipid. The effect is slightly more pronounced in the presence of calcium in the subphase. This behavior can be explained in terms of the repulsive electrostatic forces between the DPPG heads. A clear picture emerges when we plot the electron density profiles in the z direction of the phosphate groups (unless otherwise stated the subphase does not contain calcium). In Fig. 4 A, the trend of a more diffuse DPPC lipid head distribution with increasing area per lipid is clear. As suggested by experiments and simulations (Ho et al., 1995; Pasenkiewicz-Gierula et al., 1999) there is a mutual attraction between phosphocholine molecules, and this attraction inhibits the diffusive motion in the z direction.

View larger version (23K): [in this window] [in a new window]   FIGURE 4   Electron density profiles of DPPC (A) and DPPG (B) phosphate groups at different surface densities.

The attraction between zwitterionic phospholipids has been explained in terms of the formation of water cross-bridges between negatively charged groups in which a water molecule is hydrogen-bonding to two phosphocholine molecules and the weak electrostatic interactions between the choline and phosphate groups of neighboring lipids (Pasenkiewicz-Gierula et al., 1999). In Fig. 5, the intermolecular lipid head phosphorus-phosphorus and phosphorus-nitrogen radial distributions functions are plotted, confirming the existence of a strong attraction between the choline and the phosphate groups. When the area per lipid is increased, more water molecules penetrate into the freed volume spaces and screen the attractive electrostatic interactions, enabling the entropy-driven dispersion of the lipids in the z direction. The interfacial thickness increases accordingly. Calcium ions do not affect this behavior of the DPPC monolayers.

View larger version (12K): [in this window] [in a new window]   FIGURE 5   Intermolecular DPPC head phosphorus-phosphorus (P-P) and phosphorus-nitrogen (P-N) radial distributions functions.

In the case of DPPG monolayers the opposite trend is observed for the phosphate group electron density changing with the area per lipid (Fig. 4 B). This behavior can be explained in terms of the repulsive electrostatic interaction between the DPPG heads. This repulsion results in the lipid head position distribution in the z direction being broadened at small areas per lipid. As in the case of DPPC, increasing the area per lipid allows water molecules to penetrate into the lipid head group area of DPPG, effectively screening the electrostatic interactions. Other forces, such as hydrophobic, bring the lipid heads closer to each other along the z direction. Eventually, the equilibrium distribution of the DPPG phosphate groups at 70 Å²/lipid is similar with that for DPPC molecules at the same area per lipid. Calcium ions in the subphase apparently contribute to the screening of the repulsive forces between the DPPG heads, resulting in a more pronounced decrease of the interfacial thickness with increasing area per lipid as shown in Table 1.

For systems with mixed lipid components (PCPG) the trends are similar with the ones in DPPG monolayers, albeit less pronounced presumably due to the screening effect the DPPC lipids have on the DPPG repulsion. In Fig. 6, the radial distribution function g(r) is plotted for the phosphate groups of DPPC and DPPG lipids. The peaks in the DPPG monolayer are more attenuated because of the more diffuse nature of the lipid head distribution. Note, however, that the first peak occurs at smaller distances in DPPG than in DPPC. Because the DPPG lipids are more dispersed in the z direction than the DPPC lipids, it appears that besides the electrostatic repulsion there is a source of attraction in the lateral directions for the DPPG molecules. Infrared spectroscopy studies suggest that intermolecular hydrogen bonding occurs between the glycerol hydroxyl and the phosphate or carbonyl groups of the phospholipid (Dicko et al., 1998). In Fig. 7, we plot the radial distribution function between the glycerol hydroxyl hydrogens and the phosphate and carbonyl oxygens. Whereas there is apparently no specific interaction between the carbonyl oxygens and the glycerol hydroxyls, there is a clear peak at ~2.6 Å in the radial distribution function of the hydroxyl-phosphate interaction confirming the experimental suggestion.

View larger version (14K): [in this window] [in a new window]   FIGURE 6   Radial distribution function g(r) between the phosphate group phosphorus atoms of DPPC and DPPG lipids.

View larger version (16K): [in this window] [in a new window]   FIGURE 7   Radial distribution function between DPPG glycerol hydroxyls and the phosphate oxygens.

Hence, a more unambiguous picture emerges for the distribution of the acidic lipid heads in the normal and lateral directions. Strong repulsive electrostatic interactions result in a broader distribution in the direction normal to the monolayer's interface, and this distribution is further stabilized through intermolecular hydrogen bonds between functional groups in the lateral directions. This is possible due to the offset between these groups along the z axis of the molecule.

In Fig. 8, snapshots at the end of simulations of DPPG monolayers at two different areas per lipid are shown. It is evident that there is a phase transition from a liquid condensed to a liquid disordered phase upon increase of the area. This transition is accompanied by a marked increase in the conformational disorder of the lipid tails. The lipid tail conformational order can be quantified by calculating the order parameter S_CD.

(1)

with each of the components calculated as

(2)

The brackets indicate an ensemble average and  is the angle between the carbon-hydrogen bond of the methylene and methyl groups and the monolayer normal. The subscript CD stands for carbon-deuterium and has been traditionally used to denote the bonds monitored experimentally by nuclear magnetic resonance techniques. Although in the simulations we use carbon-hydrogen bonds, we retain the nomenclature for consistency with the literature. In Fig. 9, the order parameters are plotted for all the carbon atoms of the lipid tails for DPPC and DPPG monolayers at different areas per lipid. Visual inspection reveals a general trend of reduced order parameters with increasing area per lipid. The decrease is more continuous for DPPG monolayers, although it is not linear, in agreement with infrared reflection-absorption spectroscopy measurements (Dicko et al., 1998). For the DPPC monolayers there appears to be a more abrupt transition from 60 Å²/lipid to 65 Å²/lipid. Interestingly, the values of the order parameter do not follow a decreasing trend with the number of the carbon atom, as is the case for lipid bilayers. In bilayers the order parameters have similar values for carbon atoms 2 through 10 and decrease in an almost linear fashion toward the end of the lipid tail (Seelig and Seelig, 1974). For the monolayers, the simulations suggest that only for carbon atoms 14 through 16 do order parameter values begin to decrease from the plateau values. Apparently, the absence of another lipid layer allows the lipid tails to order themselves without competition from the tails of the opposing leaflet of the bilayer. Experimentally, it is challenging to quantify the conformational order in monolayers, because nuclear magnetic resonance and infrared spectroscopic techniques, which are the methods traditionally used for bilayers, lack the necessary sensitivity (Gericke et al., 1996). Gericke and co-workers used infrared reflection-absorption spectroscopy and qualitatively suggested that for DPPC monolayers the region near the head group is more ordered than the region at the tail end and that compression results in increased ordering. The simulation results are in agreement with the experimental observations.

View larger version (71K): [in this window] [in a new window]   FIGURE 8   Snaphsots of DPPG monolayers at 60 Å2/lipid (A) and 70 Å2/lipid (B). Lipids are shown as wireframe, water molecules as cyan cpk, sodium ions as yellow cpk, and chloride ions as red cpk.

View larger version (29K): [in this window] [in a new window]   FIGURE 9   Order parameter profiles for DPPC (A) and DPPG (B) monolayers at the simulated surface densities.

In Fig. 9, it is also evident that the DPPC tails are generally more conformationally disordered that the DPPG ones at the same surface density. We have also found (results not shown) that the nature of the ionic environment does not significantly influence the shape of the order parameter curves.

We turn now to the position of the ions and their interactions with functional groups of the lipids. In Fig. 10 A, the electron density profile of the sodium and chloride ions is plotted in monolayers of various lipid contents at the same surface density. It is clear that the sodium ions preferentially diffuse in the lipid head group area of the DPPG monolayers. On the other hand the chloride ions are repelled from the DPPG monolayer interface. In DPPC neither the sodium ions nor the chloride ions show preferential interactions with the lipids, distributing themselves homogeneously in the water subphase. The areas below the profile curves are a measure of the number of ions in the simulated systems. The difference in the sodium ions is due to the addition of counterions to balance the charge of DPPG. Calcium ions also diffuse preferentially in the DPPG lipid head area as shown in Fig. 10 B.

View larger version (27K): [in this window] [in a new window]   FIGURE 10   Electron density profiles of ions and their relative position to the phosphate group and the water subphase: (A) sodium and chloride ions and (B) calcium ions.

To identify possible specificity in the interactions, the radial distribution functions of the ions with the phosphate groups of the lipids in the pure lipid component systems are plotted in Fig. 11. The differences depending on the nature of the lipid head are apparent. The DPPG lipids repel chloride ions. Sodium and calcium ions interact strongly with the DPPG phosphate group. Interestingly, there is another calcium-phosphate peak at ~3 Å., besides the one at ~6 Å. This peak reveals the highly specific nature of interaction between the divalent cations and the DPPG phosphate groups. Experimental studies of the effect of calcium on the nuclear magnetic resonance spectra of phosphoglycerol molecules have clearly demonstrated that phosphate is the group that calcium binds tightly to (Lau et al., 1981). The simulations confirm this finding.

View larger version (19K): [in this window] [in a new window]   FIGURE 11   Radial distribution functions between the ions in the simulation and the center of mass of the lipid phosphate groups.

A very interesting property calculated from the trajectories is the surface tension of the monolayer. The surface tension, as calculated by Eq. A2 in the Appendix, using the components of the pressure tensor, is a measure of the anisotropy in the pressure due to the inhomogeneity of the system in the normal direction. In Fig. 12 the normal and one of the tangential components of the pressure tensor, P_zz and P_yy respectively, are plotted for the first 1200 ps of a simulation. The first 200 ps are part of the equilibration process, followed by 1 ns of production run simulation. The calculated pressures fluctuate significantly because of the finite size of the system. Nevertheless, it is evident that these values are equilibrated for the 1 ns of production run. The normal pressure fluctuates around a value of 0 atm and the tangential component equilibrates to a negative pressure value. The other tangential component, P_xx, equilibrates to the same value like P_yy as expected due to the lateral dimension system symmetry. The negative tangential pressures suggest that there is a free energy cost for expanding the surface and give rise to positive surface tension values.

View larger version (44K): [in this window] [in a new window]   FIGURE 12   Pressure tensor components during equilibration and production run.

From the components of the pressure tensor we calculate the surface pressure, , as described in the Appendix. In Fig. 13, the surface pressure isotherms are plotted for the simulated systems. In Fig. 13 A, the values of  are plotted for systems without Ca⁺² ions, the influence of which is shown in Fig. 13 B. The pressure isotherms show a decreasing trend with decreasing surface density capturing qualitatively the behavior observed experimentally. For all the areas per lipid and all the ionic environments, the DPPG monolayers exhibit higher surface pressures than the DPPC ones. This, again, is in excellent agreement with experimental Langmuir trough measurements (Nag et al., 1994). For the DPPC monolayers, addition of divalent cations results in insignificant changes in the surface pressure except for the system at 55 Å²/lipid. Contrary to this result, experimental measurements show that there is no difference in surface pressure of DPPC at 55 Å²/lipid between systems with and without calcium, as expected by the minimal interaction between calcium ions and DPPC molecules. The disparity between the simulation and the experiment points to the limitations of the simulations to accurately predict the surface pressure of monolayers. These limitations are also pointed by the occasional unrealistic surface pressure increase with increasing area per lipid and the unphysical instance of  > _o for the DPPG monolayer system at 55 Å²/lipid, in which _o, is the air/water surface tension. Nevertheless, the general trends are encouraging enough to suggest that simulations can be used to investigate the interfacial behavior of phospholipid monolayers.

View larger version (20K): [in this window] [in a new window]   FIGURE 13   Simulation surface pressure isotherms: (A) comparison between systems of various lipid content without calcium ions, and (B) comparison between systems with and without calcium ions.

For DPPG monolayers, insertion of calcium results in a significant shift of the surface pressure isotherms to lower values. This response of the surface pressure to the change of ionic environment was also captured by Langmuir trough experiments (Nag et al., 1994). Nag and co-workers explained this change in terms of the expected significant dehydration of the DPPG headgroup region upon binding of calcium (Garidel and Blume, 1999). In Fig. 14, the radial distribution function between the center of mass of the phosphate group and the oxygen of water molecules is plotted.

View larger version (15K): [in this window] [in a new window]   FIGURE 14   Radial distribution function between lipid phosphate groups and water molecules.

It is observed that insertion of calcium results in a clear, albeit small, decrease in the phosphate group hydration confirming the suggestion of Nag et al. (1994). Interestingly, there is a stronger attraction between the hydration shell and the charged DPPG headgroup than the zwitterionic DPPC headgroup. A similar finding emerged from external reflection-absorption infrared studies of DPPC and dipalmitoylphosphatidylserine, another acidic phospholipid (Flach et al., 1993).

Finally, a comment is deserved on the diffusivities of the lipids. As discussed in the Introduction, phospholipids in monolayers diffuse ~10 Ų ns¹. During the 1 ns of our simulations, we observed no jump-type motions where lipid pairs interchange their position on the xy plane. We extended the simulation of the mixed DPPC/DPPG system at 70 Å²/lipid to 7 ns, and four instances of jump-type motions were observed. It thus becomes clear that is not yet possible to investigate miscibility phenomena in lipid monolayers using atomistic simulations. Nevertheless, the ever-increasing power of computers encourages optimism, because an order of magnitude increase in the performance of computers will allow for a thorough miscibility study in a timely and tractable manner.

	CONCLUSIONS

TOP ABSTRACT INTRODUCTION METHODS RESULTS CONCLUSIONS APPENDIX REFERENCES

We have presented results of atomistic MD simulations of DPPC and DPPG monolayers at the air/water interface. The thickness of the interface was calculated from electron density profiles, and its dependence on the type of the lipids, the surface density, and the ionic environments was clarified. The distributions of the lipids in the lateral and normal to the interface directions were investigated using atomistically detailed pair distribution functions and electron density profiles. The nonlinear response of the order parameter to surface pressure changes was also studied. Radial distribution functions clearly show that calcium ions bind with considerable strength and specificity on the phosphate groups of the anionic lipids. Surface pressure isotherms were calculated with encouraging qualitative agreement with experiments. These results clearly demonstrate the strength of simulations in investigating complex, multicomponent systems of biological importance.

We demonstrated that the chosen geometry for the simulation cell, where the simulation of a single monolayer is facilitated by the use of a wall, allows for a realistic representation of phospholipid monolayer interfaces and greatly facilitates the investigation of the interactions between membrane active molecules and lipids. Simulations of small lung surfactant proteins with DPPC and DPPG monolayers are under way with promising preliminary results.