Molecular Dynamics Simulations of Liquid Condensed to Liquid

J. Phys. Chem. B 2010, 114, 1325–1335

1325

Molecular Dynamics Simulations of Liquid Condensed to Liquid Expanded Transitions in DPPC Monolayers Delara Mohammad-Aghaie,†,‡ Emilie Mace´,† Charles A. Sennoga,† John M. Seddon,*,† and Fernando Bresme*,† Department of Chemistry, Imperial College London, SW7 2AZ London, United Kingdom, and Department of Chemistry, Shiraz UniVersity of Technology, Shiraz 71555-313, Iran ReceiVed: June 30, 2009; ReVised Manuscript ReceiVed: October 1, 2009

We have investigated the phase behavior of DPPC (dipalmitoylphosphatidylcholine) monolayers at the water-air interface using molecular dynamics simulations, where the phospholipids and the water molecules are modeled atomistically. We report pressure-area isotherms in the interval of 273-310 K. Our results show evidence for a liquid condensed (LC) to liquid expanded (LE) phase transition and indicate that ordered condensed phases can nucleate from a starting disordered phase on a time scale of approximately 50 ns. The existence of the phase transition is confirmed with structural analyses of the phospholipid pair correlation functions and of the monolayer thickness. We find that the change in the monolayer thickness associated with the LC-LE transition is largely due to a shortening of the hydrocarbon chains, with little modification in the average tilt angle of the choline head group. This result is compatible with recent sum frequency spectroscopy experiments, which concluded that the transition occurs without major changes in the orientation of the head group with respect to the monolayer plane. The dependence of the simulated pressure-area isotherms on temperature, in particular, the reduction in width of the coexistence plateau with increasing temperature, is consistent with published experimental pressure-area isotherms. Introduction Investigation of the phase behavior of phospholipid monolayers is attracting considerable interest. Phospholipid monolayers play a major role in the biological processes concerned with the regulation of the surface tension of the air-alveolar interface, which is essential for breathing, and it has been suggested that the ability of these self-assembled structures to undergo order-disorder transitions (solid-liquid) may influence signal transduction1 and protein transport. Moreover, phospholipid self-assembled monolayers are being used as stabilizers in the formulation of emulsions and the production of gas-filled microbubbles for ultrasound imaging, drug delivery, and blood substitutes.2 Phosphatidylcholine is one of the major phospholipids in most animal cell membranes, and dipalmitoylphosphatidylcholine (DPPC) accounts for approximately 50% of the surfactants in the alveoli of the lungs. Hence, DPPC has been the subject of many investigations, both experimental and theoretical. It is well-known that hydrated DPPC monolayers exhibit a liquid condensed (LC) to liquid expanded (LE) transition close to physiological temperatures upon cooling below 41.5 °C. The bulk phospholipid transforms from a fluid lamellar LR (liquid crystal) to an ordered lamellar structure (gel phase), where the phospholipid chains arrange themselves with a two-dimensional hexagonal packing.3 Very recent experiments of DPPC monolayers using sum frequency generation spectroscopy indicate that the LE-LC transition occurs via dehydration of the phospholipid head group, without major changes in the orientation of the head group with respect to the monolayer plane.4 * To whom correspondence should be addressed. E-mail: j.seddon@ imperial.ac.uk (J.S.); [email protected] (F.B.). † Imperial College London. ‡ Shiraz University of Technology.

This result seems to agree with earlier experiments using X-ray studies of DPPC, where it was concluded that the orientation of the head group is not sensitive to the packing of the DPPC chains.5 On the other hand, neutron reflectivity investigations indicated significant differences in the orientation of the DPPC head group of the LE and LC phases.6 The LC-LE transition is very sensitive to temperature and sample composition. The coexistence plateau of the transition shifts to higher surface pressures with increasing temperature.7 Experiments have not yet demonstrated whether this transition ends in a critical point or whether this is preceded by a collapse of the phospholipid monolayer. Head group charge, chain unsaturation, and degree of deuteration can affect the transition and may even tend to inhibit it.6,8,9 Experiments also indicate that the phase behavior of calf lung surfactant extracts, which consist of a mixture of phospholipids and proteins, can be markedly different from that of pure DPPC,10 suggesting that such a composition may enable supercompression of the film without inducing the transition. Similarly, molecular gases alter the packing of phospholipid membranes11 and consequently the onset of the transition.12 These previous studies emphasize the relevance of composition and thermodynamic conditions in tuning the phase behavior of phospholipid films. In addition to composition and temperature, the phase behavior of DPPC monolayers may be influenced by dynamic factors, such as the compression rate.13 Such an effect has been observed in experimental studies of DPPC14 and was interpreted in terms of the formation of nonequilibrium structures. Horie and Hildebrant already indicated in an early paper that it is difficult to deduce static behavior from dynamic data and vice versa.15 These articles emphasize the importance of dynamic effects on the mechanical response of the film. Such dynamical effects are expected to be significant in alveolar mechanics, where the alveoli may experience significant and rapid changes

10.1021/jp9061303  2010 American Chemical Society Published on Web 12/28/2009

1326

J. Phys. Chem. B, Vol. 114, No. 3, 2010

in surface area. Similarly, microbubbles currently under investigation for use as ultrasound contrast agents to image internal organs16 undergo large changes in surface area in response to the pressure field generated by the ultrasound.17 These situations pose important questions regarding the state of the monolayer at high compression rates, particularly whether the LC phase nucleates under fast compression conditions. Despite the considerable amount of experimental work on the surface-pressure isotherms of DPPC, there is still a certain degree of uncertainty regarding the true isotherms. Duncan and Larson18 have reviewed this situation in an extensive article. As discussed by these authors, the pressure-area isotherm depends to a certain extent on the experimental conditions used; low pH, salt conditions, leakage in Langmuir troughs, different compression rates, impurities, and, particularly, the spreading solvent used in the experiments can significantly affect the measured isotherms. Similarly, there is some indication that the experimental method used can also affect the results. Differences between the isotherms obtained with the Langmuir trough approach and by axisymmetric drop shape analysis have been reported.13 On the other hand, good agreement between the isotherms obtained using the captive bubble and micropipet approaches and the Langmuir trough technique has been reported.7,19 Lee and co-workers19 have suggested that the discrepancies in surface pressure values reported in the literature using different techniques can be ascribed to the long equilibration times required in Langmuir trough experiments. The development of accurate force fields has advanced considerably the application of computer simulations to investigate phospholipid films. Many simulation studies have focused on DPPC phospholipid bilayers.18,20-31 By comparison, only a small number of simulation studies have been devoted to phospholipid monolayers,18,32-38 and only a few of these have addressed the computation of the pressure-area isotherm. Earlier simulations32 reported pressure-area curves using short simulation times (subnanosecond). More recent work indicates that very long simulation times, at least tens of nanoseconds, are necessary to equilibrate the monolayer and enable the formation of the liquid condensed phase. Very long simulation times, on the order of microseconds, have been achieved using coarse-grained models, where several groups of atomic sites in the phospholipid chain and the solvent are represented as pseudoatoms. Baoukina and co-workers34 showed that such a model can reproduce the LE-LC transition of DPPC. This has been confirmed in more recent simulations by different authors.18 These isotherms show a reasonable overall agreement with some experimental results, although, in general, they tend to overestimate the area compressibility modulus. In a recent study,35 some evidence has been provided for the formation of the liquid condensed phase, but the pressure-area isotherm was not investigated. The effect of additives, cholesterol, on the phase stability of the DPPC monolayer has been discussed recently, concluding that high concentrations of cholesterol tend to stabilize the LC phase.38 Many of the isotherms reported to date correspond invariably to the liquid expanded phase (see, for instance, ref 18 for an extensive discussion of simulation results). Atomistic models, where the electrostatics and the dielectric properties of the solvent are considered explicitly, may help to understand the general phase behavior of phospholipid films. For example, the miscibility of DPPC/DPPG (phosphatidylglycerol) mixtures is sensitive to the presence of Ca2+ ions.39 Water can bind very strongly to charged groups, and the dielectric properties of interfacial water may differ significantly from those of bulk water, affecting the screening of the

Mohammad-Aghaie et al. Coulombic interactions. We have previously reported that water confined in surfactant thin films exhibits an anomalous dielectric behavior.40-42 Moreover, small changes in the phospholipid chain architecture, for example, DPPC versus POPC (palmitoyloleoylphosphatidylcholine), lead to large modifications of the phase diagram. To address these effects, particularly those concerned with electrostatics, it is necessary to include the effect of the solvent explicitly. In this paper, we report pressure-area isotherms of a DPPC monolayer adsorbed at the water surface using an atomistic model. The simulated isotherms show the existence of the LC-LE transition within this model. The temperature dependence of the surface pressure of the transition, as well as the jump in surface area, agrees with the experimental observations. Structural analyses of the monolayers provide further evidence for the transition. Computational Details We have performed molecular dynamics simulations of DPPC monolayers adsorbed at the water-air interface. The phospholipids were adsorbed at one of the water film-air surfaces, whereas the other water surface was left free of phospholipids. See Figure 1 for a representative snapshot of the simulation setup and Figure 2 for a molecular model of the DPPC phospholipid. The phospholipids were modeled using the Berger et al. force field.21 This force field uses a united atom approach to represent the CH2 and CH3 groups. All of the interatomic bonds are rigid, and the angles and improper dihedrals are modeled through the GROMOS force field. The torsional interactions of the hydrocarbon chains were modeled through the Ryckaert-Bellemans potential.43 Nonbonded interactions were handled through a combination of Lennard-Jones and Coulombic terms

Uij(r) ) 4εij

[( ) ( ) ] σij r

12

-

σij r

6

+

1 qiqj 4πε0 r

(1)

where r is the distance between atoms of species i and j. In this force field, the effective atom diameters σij and energies εij can be derived using the combination rules, σij ) (σiσj)1/2 and εij ) (εiεj)1/2. The appropriate parameters for the atom partial charges, energies (ε), atom diameters (σ), and bonding, bending, and torsional constants can be found in ref 21. Most simulations of phospholipids reported so far, including the original parametrization of Berger et al.,21 have considered the single-point charge (SPC) model of water.44 This model and its variant, the SPC/E,45 predict surface tensions that are too low in comparison with the actual surface tension of water. This has an impact in the surface pressure computations, shifting the pressure isotherms to lower values as compared with those in the experimental case. We note that similar problems arise in other models, for example, TIP3P, which is also widely employed in the literature.29 In the present work, we have considered the TIP4P/200546 model. This model provides a more accurate representation of the surface tension of water.47 In general, it provides a good account of the solid-liquid and liquid-vapor phase diagrams of water, as well as liquid densities and critical properties. The model represents water as a rigid four-site molecule, where the atoms interact through a combination of Lennard-Jones and Coulombic forces (cf. eq 1). The effective atomic diameters and energies needed to calculate the phospholipid-water interactions were obtained using the combination rules discussed above.

Liquid Condensed to Liquid Expanded Transitions in DPPC

J. Phys. Chem. B, Vol. 114, No. 3, 2010 1327

Figure 1. DPPC monolayer adsorbed at the water-air interface. T ) 310 K, ADPPC ) 0.639 nm2, and Π ) 19.1 mN/m.

Figure 2. Molecular model of dipalmitoylphosphatidylcholine (DPPC). The numbering of the atoms used in the simulations is indicated (see ref 21).

The initial configuration of the DPPC monolayer was prepared by replicating a single DPPC molecule. The sn-1 and sn-2 chains in the DPPC molecule were set in the all-trans conformation. Sixty-four DPPC molecules were arranged in a square twodimensional lattice, and the head groups were partially immersed in one of the water surfaces of a pre-equilibrated water film containing 4 × 103 water molecules. The water-DPPC system was set parallel to the (x,y) plane, and two vacuum regions were left above and below the film in the z direction. We employed full periodic boundary conditions in all of the simulations. The short-range Lennard-Jones interactions were truncated at 1.7 nm (about 5.5 times the water-oxygen diameter), and the longrange Coulombic interactions were computed using the Particle Mesh Ewald summation method.48,49 Because the system is asymmetric in the z direction, in order to eliminate spurious effects arising from correlations between the dipoles of the simulation box images, we applied the correction proposed by Yeh and Berkowitz.50 Following ref 50, we ensured that the vacuum region above the monolayer was larger than the box length in the parallel direction of the simulation box. This large vacuum region was much larger that the cutoff employed for the dispersion interactions, hence ensuring the absence of spurious effects connected to the interactions between the molecules in the central box and their images. We did not consider any long-range corrections for the Lennard-Jones term

since for the cutoff value employed here, such corrections are expected to be very small. We note that a straight cutoff of the Coulombic interactions results in a large force discontinuity, which can strongly influence the ordering of water molecules. This issue has been extensively discussed in the simulation literature, and it has been investigated in the case of DPPC bilayers too.51 With this simulation setup, the surface pressure was obtained from Π(A, T) ) γw(T) - γm(A, T), where A is the area per phospholipid, γw is the surface tension of the TIP4P-2005 water at temperature T, and γm is the surface tension of the monolayer. The latter surface tension was extracted from the simulation as follows, γm(A, T) ) γt(A, T) - γw(T), where γt is the total surface tension imposed on the DPPC-water layer. The constant surface tension and temperature simulations were performed using the Berendsen barostat/thermostat with coupling constants of 1 and 0.1 ps, respectively. The compressibility for the surface tension computations was set to 4.5 × 10-5 bar-1. The crosssectional area of the interface employed in the simulations was large enough to avoid the periodic error associated with the computation of the surface tension.52 The simulations involved a pre-equilibration of the initial configuration. The area per phospholipid for the initial configuration was set to 0.55 nm2. This configuration was equilibrated using the Berendsen thermostat at T ) 310 K for tens of nanoseconds, before performing the simulations under constant surface tension and constant temperature conditions. The initial configuration employed for all of the simulations corresponded to a liquid expanded phase. One of the objectives of our work was to investigate whether the model investigated here is able to nucleate the liquid condensed phase in the simulation times considered. As we will see below, this phase can be obtained within 50 ns simulation times. Such simulation times are needed to obtain reliable results for the pressure-area isotherm. We used a time step of 2 fs for all of the simulations. All of the calculations were performed with the code GROMACS 3.3.53,54

1328


Figure 3. Dependence of the surface pressure on the dispersion interaction’s cutoff and on the water model. The simulations correspond to a liquid expanded phase at 310 K. The error bars are representative of the error associated with the areas per phospholipid. The lines joining the points are just a guide to the eye.

Results and Discussion Prior to the production runs, we performed a number of tests to assess the sensitivity of the surface pressure to the force field details. Our cutoff for the dispersion interactions is significantly longer than the one used in previous studies. It has been shown in computer simulations of simple liquids that the cutoff used to truncate the dispersion interactions can have a significant impact on the phase diagram obtained. In general, decreasing the cutoff results in a reduction of the critical temperature of the liquid-vapor phase diagram.55,56 To address the effect of the cutoff on the surface pressure isotherms, we performed several simulations using different cutoffs. The impact of the water model was also addressed by computing the isotherms with the SPC/E and TIP4P-2005 models. Figure 3 shows the results of our analysis for a liquid expanded phase at T ) 310 K. Our results clearly show that long cutoffs result in a significant lowering of the surface pressure. This effect mainly follows from the increasing interactions between the surfactants which render a lower pressure. This effect is similar to that observed in simple liquids. The surface pressure isotherm is also sensitive to the water model employed, although in comparison, this dependence is weaker than the cutoff. Hence, to minimize cutoff effects, we chose a long cutoff, 1.7 nm, for the dispersion interactions. Figure 4 shows the dependence of the simulation box length with time for two representative temperatures. At 310 K, thermodynamically stable states reach equilibration within 15-20 ns. However, we find that metastable states can survive for very long times. These states can be easily identified by inspecting the surface pressure curves (c.f. Figure 6), with Π > γw for small areas and Π < 0 for large areas. Near the liquid-vapor instability, simulation times on the order of 40 ns are needed before the system can undergo the transition to the vapor phase (see Figure 4). On the other hand, states near the collapse region, (Π > γw) can survive for even longer times. These states are not stable, and according to the experiments, the monolayer should buckle and form a three-dimensional structure. Clearly, such transformation requires times much longer than the ones considered here (60 ns). At lower temperatures, 293 K (cf. Figure 4), the equilibration time slows down. We find that the metastable states corresponding to high compressions, Π > γw, show a steady decrease for the whole duration of the simulation (up to 60 ns). On the other hand, the

Mohammad-Aghaie et al. existence of an instability at larger areas is more evident at this temperature. For surface pressures that are below the limit of mechanical stability of the monolayer, the box length shows a dramatic increase with time (see Figure 4; Π ) -15.96 mN/ m). Similar to the higher temperature case, metastable states can survive for the whole duration of the simulation considered here (see Figure 4; Π ) -5.95 mN/m). For highly compressed monolayers (see Figure 4; Π ) 84.1 mN/m), above the stability limit defined by Π ) γw, we find again a steady slight decrease of the simulation box length in the whole simulation interval. Clearly, the convergence of the simulation box length with time is not a good criterion to decide on the stability of the systems simulated. This information has to be combined with an analysis of the surface pressure. Hence, care has to be exercised when combining the surface tensions of the simulations with the experimental surface tension of water, as has been done in previous studies. We believe that this approach might introduce some problems, particularly if the surface tension of the water model employed in the simulations is significantly lower than the experimental surface tension. In these cases, thermodynamic states that are nonstable according to the surface pressure (Π < 0) would appear stable (Π > 0) when using the experimental surface tension of water. This concern has prompted us to use a water model that provides a better description of the interfacial properties of water. Experimental studies have shown that DPPC monolayers at temperatures below the melting temperature of the bilayer (314.2 K) exhibit a liquid expanded (LE) to liquid condensed (LC) transition.57 Computer simulation studies have shown the existence of this transition but only using a coarse-grained approach.18,34 This transition has not been investigated in detail using an atomistic description, although there is some evidence that such a model can lead to the formation of small phospholipid domains with the characteristic hexagonal close-packing arrangement of the LC phase.35 One question that we want to address in the present paper is whether such a transition can be observed using an atomistic approach and therefore whether the liquid condensed phase can nucleate within acceptable simulation times. To answer this question, we have computed the pressure-area isotherms of DPPC monolayers. In order to ensure that the LC phase that nucleated during the simulation was stable, we have performed a time-dependent analysis of the pressure-area isotherm and the monolayer structure. Figure 5 shows the dependence of the surface pressure with time for T ) 293 K. At this temperature, the diffusion is slower than that at the higher temperatures, and hence, it represents a more stringent test for equilibration. We find that after 30 ns, the surface pressure does not show any drift and fluctuates around a well-defined average. The fluctuations for the system at lower surface pressure, 4.05 mN/m, are, as expected, larger. In addition, we have performed an analysis of the time dependence of the structure of the LC phase for the latter surface pressure. The radial distribution function does not show any particular trends after 20 ns and exhibits the longrange order characteristic of a liquid condensed phase. The structure of this phase will be discussed below in detail. We ensured with similar tests that all of the thermodynamic states investigated in this work had reached equilibration within the time scale of the simulation. Figure 6 shows one of the main results of our paper, and Table 1 compiles the corresponding numerical data. All of the results were obtained from subaverages over 10 ns blocks. The pressure-area isotherms in the interval (293 and 310 K) show a characteristic change in the slope at a specific surface pressure,



Figure 4. (Top) Evolution of the simulation box length parallel to the monolayer with time at 310 K. The lines represent results for different surface pressures. From bottom to top, 79.1, 49.1, 29.1, 19.1, 9.1, -0.88, and -10.9 mN/m. (Bottom) The same as (Top) for T ) 293 K. From bottom to top, 84.1, 54.1, 34.1, 14.1, -5.95, and -15.9 mN/m.

indicating the possible existence of a liquid condensed to a liquid expanded phase transition. Large changes in area, about 25%, are obvious at low temperatures, 293 K, where the transition would happen at a surface pressure very close to zero. We expect a similar behavior at even lower temperatures, 273 K, although within the simulation time scale, we did not observe a transition from the LC to the LE phase.Our results indicate that the area change associated with the LC-LE transition decreases as the temperature increases (see Figure 6). This observation agrees with experimental results.7 Such a dependence of the area change with temperature was not observed in previous molecular dynamics investigations using coarse-grained models.34 The LC-LE transition is a first-order transition; therefore, hysteresis effects are expected to be present in the simulations. We have addressed this point by investigating the impact that the initial configuration has on the final equilibrium area of the monolayer. To do this, we started either from a LE phase or a

LC phase, which were simulated at the same surface pressure using the constant surface tension coupling. At 310 K, the hysteresis is found to be small. Above the discontinuity in the surface pressure, we find that a starting LE configuration converts into a LC one in about 70 ps (see Π ≈ 60 mN/m in Figure 6). Similarly, below the pressure discontinuity, a LC phase relaxes toward the equilibrium area obtained for the LE phase. Only at surface pressures very near the discontinuity do we find some signature of hysteresis. The hysteresis effects are enhanced at lower temperatures. We performed in this case one simulation at a surface pressure above the discontinuity using as a starting configuration the LE phase. In this case, the system gets trapped in the LE phase for the whole duration of the simulation, nearly 80 ns. The presence of hysteresis supports the existence of the LC-LE first-order transition. To investigate further the occurrence of a LC to LE transition, we have performed a structural analysis of the monolayers at

1330


Mohammad-Aghaie et al.

Figure 5. (Top) Dependence of the surface pressure with time for the monolayer in the liquid condensed phase at 293 K. Each set of points corresponds to different surface pressures, 4.05, 14.1, and 24.1 mN/m. Each point corresponds to a subaverage over 10 ns. (Bottom) Time evolution of the radial distribution function of the carbon atoms in the phospholipid chains. The correlation functions were obtained using three carbon atoms in the middle of each phospholipid chain, (21, 22, 23) in the sn-2 chain and (40, 41, 42) in the sn-1 chain (see Figure 2 for definitions of the carbon atoms).

Figure 6. Pressure-area isotherms of DPPC monolayers at the water-air interface. Surface pressure in mN/m and area in nm2 per phospholipid. The horizontal lines represents the limits of stability of the monolayers, Π ) γw and 0. The lines joining the points are just a guide to the eye. The open symbols at T ) 310 and 293 K correspond to the hysteresis tests discussed in the main text. These tests involve the analysis of trajectories 80 ns long. Circles: initial state; squares: intermediate state; and diamonds: final state. The symbols correspond to subaverages from 0 to 10 ns for the initial state, and 30-40 and 70-80 ns for the intermediate and final states, respectively. The arrows indicate states referred to in the text and in Figures 5 and 6.

different surface pressures, looking for structural signatures of the LC and LE phases. Figure 7 shows the two-dimensional pair correlation function (in the interface plane) of the hydrocarbon chains. To compute this radial distribution function, we used three carbon atoms in the middle of each phospholipid chain. The liquid condensed phase is characterized by the existence of long-range order in the correlation function. Our results show long-range order for those surface pressures above the suggested transition (see the arrows in Figure 6 and points in Figure 7). We find a clear loss of long-range order in a surface pressure interval that is consistent with the change of slope observed in the pressure-area isotherms. The formation of the liquid condensed phase can be further confirmed by inspecting

snapshots of the configurations at different surface pressures (cf. Figure 8). We have represented snapshots corresponding to LC and LE phases (see the arrows in Figure 6 and temperatures 310 and 293 K). The liquid condensed phase shows a characteristic hexagonal packing, and the chains are straight and tilted (31° with respect to the monolayer normal). Upon the transition, the structure becomes disordered; this is particularly clear in the loss of orientational order of the phospholipids. We find that at the higher temperatures, 310 K, some ordered (LC) nuclei are still visible in the LE phase. Overall, our results confirm the existence of the transition and show that simulation times on the order of 50 ns are long enough to nucleate the liquid condensed phase using atomistic simulations.



TABLE 1: Surface Pressure of DPPC Monolayers As a Function of Interface Area and Temperaturea 310 K

300 K

293.15 K

273.15 K

A

δA

Π

A

δA

Π

A

δA

Π

A

δA

Π

0.462 0.472 0.486 0.553 0.591 0.607 0.639 0.690 0.763 0.900

0.004 0.003 0.006 0.007 0.007 0.009 0.011 0.016 0.023 0.032

79.1 69.1 59.1 49.1 39.1 29.1 19.1 9.10 -0.88 -10.9

0.466 0.495 0.523 0.529 0.557 0.620 0.653 0.746 0.833

0.004 0.005 0.004 0.011 0.003 0.012 0.017 0.015 0.016

81.8 71.7 51.7 41.7 31.7 21.7 11.7 1.79 -8.30

0.487 0.508 0.523 0.540 0.550 0.554 0.571 0.609 0.753 0.783

0.004 0.004 0.006 0.006 0.003 0.003 0.006 0.004 0.020 0.012

84.1 54.1 44.1 34.1 24.1 14.1 4.05 0.05 -0.95 -5.95

0.519 0.528 0.536 0.538 0.547 0.556 0.561 0.589

0.002 0.002 0.003 0.003 0.005 0.004 0.004 0.004

79.8 59.8 49.8 39.8 29.8 19.8 9.75 -0.25

a

The surface tensions of water (TIP4P-2005) are 69.9 (273 K), 67.0 (293 K), 65.9 (300 K), and 64.6 (310 K). Areas are in nm2; δA is the error associated with the area, one standard deviation from the mean, and Π is the surface pressure in mN/m.

Figure 7. Radial distribution functions (RDF) of the carbon atoms in the phospholipid chains. The data correspond to the points represented in Figure 6, except for T ) 300 K, where the results for Π ≈ 70 mN/m have not been included. The dots correspond to the thermodynamic states signaled with an arrow in Figure 6 and indicate the onset of long-range oscillatory behavior in the RDF, characteristic of the liquid condensed phase. Successive RDF have been shifted vertically by 0.2 units for greater clarity. Bottom RDF, smallest Π value; top RDF, largest Π value. The correlation functions were obtained using three carbon atoms in the middle of each phospholipid chain, (21, 22, 23) in the sn-2 chain and (40, 41, 42) in the sn-1 chain (see Figure 2 for definitions of the carbon atoms).

Previous investigations of DPPC monolayers using the Berger force field have reported lateral separation of lipids, indicating the coexistence of gel-like liquid condensed and fluid-like liquid expanded phases.35 Our simulations indicate that at 293 K, the liquid condensed phase is stable for areas per phospholipid ≈ 0.60 nm2. This area is close to that investigated in ref 35. We note that simulations tend to overestimate the equilibrium area measured in the experiments, which is in the range of 0.42 nm2 in the Lβ (untilted gel) and 0.48 nm2 in the Lβ′ (tilted gel) phases. In ref 35, it was observed that the monolayer ruptures at 0.97 nm2/molecule, a process that indicates the existence of a liquid-vapor transition. In our work, we find that the surface pressure becomes negative at ∼0.75 nm2; therefore, 0.97 nm2 would be beyond the limit of stability in our case. Nonetheless, because the liquid-vapor transition is first order, metastable states can survive for very long simulation times (see, for instance, Figure 1) before undergoing the transition. In this situation, the surface pressure can be used to decide whether the system investigated is under equilibrium conditions. Vibrational sum frequency generation experiments have indicated the existence of a sharp transition at surface areas ≈ 1.10 nm2.58 This observation is most likely connected to the liquid-vapor transition. This would indicate that the Berger force field

underestimates the coexistence area per lipid corresponding to the LE branch. The snapshots reported in Figure 8 clearly indicate an order-disorder transition in the DPPC chains in going from the LC to the LE phase. This transition can be quantified by computing the length of the hydrocarbon chains as a function of the surface pressure (cf. Figure 9). We find that the sn-2 chain is consistently shorter than the sn-1. We note that such shortening has been observed in mixtures containing cholesterol.59 It may also be present in the gel phase, but it is difficult to establish experimentally. In general, the changes in the length of phospholipid chains correlate well with the pressure-area curves presented in Figure 6, with larger changes around surface pressure values where the transition is expected to appear. At low temperatures, T ) 293 K, we find a very abrupt change in the chain length, about 40% shorter in the liquid expanded phase. Earlier experimental studies of the dependence of the orientation of the DPPC head group with packing density have reported conflicting results. Hauser et al. performed X-ray experiments, concluding that the orientation of the head group is insensitive to packing. Brumm et al.6 used neutron reflectivity to investigate the conformational changes of the DPPC monolayer at the air-water interface as a function of the surface pressure. These

1332



Figure 8. Snapshots of liquid condensed and liquid expanded phases at two temperatures, T ) 293 K (top four snapshots) and 310 K (bottom four snapshots), with the liquid condensed phase above the liquid expanded phase in each group of four. Snapshots on the left correspond to a view from the top of the monolayer, and those on the right to a view on the plane of the monolayer. The snapshots correspond to the thermodynamic states indicated with arrows in the left two panels of Figure 6.

experiments were performed at 293 K. They found a large change in the monolayer thickness, 0.55 nm, in going from the LC to the LE phase. The experimental results were rationalized in terms of the modification of the average orientation of the P-N dipole of the DPPC head group. It was suggested that this dipole is parallel to the monolayer in the LE phase, whereas it is perpendicular in the condensed phase. It was not clear whether the perpendicular orientation was attained in the LC phase or was restricted to a highly packed solid phase. The orientation of the choline group in DPPC has also been investigated experimentally using sum frequency generation

spectroscopy.4 It was concluded that the choline head group orientation is not significantly different in the LE and LC phases. Our simulations provide microscopic insight into the orientational changes undergone by the DPPC head group in the LC-LE transition. At 293 K, we find a “jump” in the phospholipid chain length of about 0.7 nm (c.f. Figure 9). We note that this change is due only to changes in the conformation of the phospholipid chains, and it does not include the whole monolayer thickness. To analyze this point further, we have computed the average P-N dipole orientation of the phospholipid chains for both the LC and LE phases. The corresponding


Figure 9. Length of the hydrocarbon chains measured as the distance between carbons C15 and C31 in the sn-2 and C34 and C50 in the sn-1 chain (see Figure 2 for definitions). The ovals correspond to the thermodynamic states indicated with arrows in Figure 6. Circles, T ) 273 K; squares, T ) 293 K; diamonds, T ) 300 K; and triangles, T ) 310 K. Open symbols: sn-1 chain; filled symbols: sn-2 chain.

Figure 10. P-N orientation with respect to the monolayer normal obtained from an ensemble average for the system at 293 K and for representative LC and LE phases. The full line corresponds to the LC phase and the dashed line to the LE phase (see arrows in Figure 6). The P-N orientational distribution for a representative LC phase at 310 K (left arrow in Figure 6) is also shown (dash-dotted line).

histogram for T ) 293 K is represented in Figure 10. In agreement with previous studies, the orientation of the head group in the LE phase shows a single maximum centered around 90°. These results agree with other simulations of DPPC monolayers reported recently18,33 and confirm that in the LE phase, the head group is lying parallel to the water surface. The LC phase shows two distinctive maxima, one at 60° and the second one at about 100°. This bimodal distribution is similar to that observed in computer simulations of 1,2-dilignoceroylphosphatidylcholine (DLGPC) monolayers.60 The probability distribution is sensitive to temperature. Our results indicate that at higher temperature, the bimodal distribution disappears, and we recover a single maximum centered around 90° (cf. Figure 10). The image emerging from our computations supports the notion that the average orientation in the LC phase is similar to that of the LE phase, ∼90°, but the probability distribution depends on the temperature, low temperatures favoring a bimodal distribution. We do not find evidence for a configuration where the head group is perpendicular to the water surface as suggested in ref 6.


Figure 11. Simulated DPPC monolayer surface tensions using atomistic (AT) and coarse-grained models (CG). T ) 293 (filled circles, AT, this work), 300 (open squares, CG34), 310 (filled squares, AT, this work), and 323 K (open circles, AT;18 open diamonds, CG36).

Comparison with Previous Experimental and Simulation Results. In the following, we compare our pressure-area isotherms with previous simulation and experimental data. In simulation studies, it has been a common practice to use the experimental surface tension of water or a fitted surface tension in order to construct the pressure-area isotherm and provide agreement with the experiments. In our view, this practice can be problematic. As discussed above, metastable states can appear to be stable. A direct comparison of the monolayer surface tensions removes the effect of an adjustable solvent surface tension. Figure 11 compares the DPPC monolayer surface tensions obtained in this work with results in other simulations. We have selected only a few results, which are representative of atomic and coarse-grained (CG) models. Interestingly, we find that atomistic and CG models provide essentially the same result for the isotherm. In particular, our atomistic results for T ) 310 K coincide with the coarse-grained simulations performed at T ) 300 K.34 These coarse-grained simulations involved 4096 DPPC molecules and simulation times on the order of 400 ns. It is remarkable that similar results can be obtained from atomistic simulations involving 64 DPPC phospholipids and up to 60 ns trajectories. We note that the signature of a LC-LE transition in the isotherm has been shown before only using coarse-grained models. Our work shows for the first time that this can be investigated using an atomistic approach, within reasonable simulation times. The surface-pressure isotherm of DPPC has also been investigated using atomistic models and the same DPPC force field at higher temperatures, T ) 323 K.18 These authors did not find evidence of a phase transition at this temperature. This isotherm was shifted to larger surface areas. This can be explained as a combined effect of (a) higher temperature and (b) a shorter van der Waals cutoff, namely, 1.2 nm, versus 1.7 nm employed by us. Other investigations of this isotherm, T ) 323 K, using coarse-grained models offer contrasting results36 (c.f. Figure 11). These authors reported an isotherm that appears at significantly smaller areas than the ones found in most simulations. This isotherm shows good agreement with one set of experimental data, but it is not consistent with the results reported by different authors using the same coarsegrained force field. It has been suggested18 that these results are affected by a misuse of the periodic boundary conditions. Our area compressibility moduli

1334


KA ) -A

∂π ( ∂A )

T

( )

)A

∂γm ∂A

T


(2)

are similar to those reported in previous simulation works.18 For the LE phase, we find, KA(310 K, 0.665 nm2) ) 130 mN/ m, and for the LC phase, KA (310 K, 0.479 nm2) ) 360 mN/m. The high LC value is on the same order as that reported in experiments conducted at high compression rates61 and using the captive bubble apparatus.7 The comparison of simulated and experimental pressure-area isotherms is not straightforward. The main problem lies in the lack of a single pressure-area isotherm. This situation has been reviewed very recently by Duncan and Larson.18 As these authors point out, the experimental conditions can strongly influence the pressure-area isotherm. Hence, trying to fit a force field to reproduce a specific set of results may be an ill-defined problem. Figure 12 shows a comparison of our results with selected experimental data. We have chosen the experimental data of Crane et al.7 as representative of high-quality measurements. These data were obtained using the captive bubble apparatus, which minimizes some of the problems found with the Langmuir trough, particularly at high surface pressure and temperatures. In addition, two sets of data obtained using a film balance apparatus are reported.57,62 In ref 57, the measurements were performed at the air-water interface with pH 7.4 Tris buffer with 150 mM NaCl. These two sets of data give an idea of the typical spread found in the experimental measurements of the isotherms. Our simulation results are close to the experiments reported by Crane.7 The general behavior of the isotherms, namely, the dependence of the surface pressure with temperatures, and the reduction of the width of the plateau with increasing temperature are in line with the experimental observations. Conclusions We have reported molecular dynamics simulations of the pressure-area isotherms of DPPC monolayers at the water surface using an atomistic force field. The force field consists of two terms, the Berger potential to model the phospholipids and the TIP4P-2005 model for water. It is shown that this force field provides a good representation of the properties of the DPPC monolayers and their associated order-disorder transitions. Our computations show that long simulation times, typically tens of nanoseconds, are necessary to equilibrate the monolayers, confirming previous observations. We show that these atomistic models exhibit a liquid expanded (LE) to liquid condensed (LC) transition. The transition is observed in the pressure-area isotherms as a distinctive discontinuity in the pressure. We have confirmed the existence of the transition by investigating different structural parameters. The two-dimensional phospholipid-phospholipid pair correlation function shows a clear loss of long-range order in going from the LC to the LE phase, with a concomitant shortening (“melting”) of the phospholipid chains. We find that the choline head group average tilt angle undergoes little change during the transition. This result supports recent experiments of DPPC monolayers using sum frequency generation spectroscopy4 and suggests that the large changes in tilt angle reported using neutron reflectivity experiments may correspond to a more ordered solid phase.6 The existence of the LE-LC first-order transition was further confirmed by the existence of hysteresis in the simulated pressure-area isotherms. Our results also indicate the existence of liquid expanded to vapor phase transitions in the DPPC monolayers. Considering

Figure 12. Simulated and experimental pressure-area isotherms of DPPC monolayers. Experimental isotherms: dashed line, T ) 298 K;62 full lines, T ) 297 (lower curve) and 310 K (upper curve);7 dash-dotted line, T ) 310 K.57 Computer simulations from this work: T ) 293 (full symbols) and 310 K (open symbols).

the simulated surface pressure isotherms, this transition would appear at approximately 0.70-0.80 nm2 per phospholipid, which is smaller than the area inferred from recent vibrational sum frequency generation experiments, ∼1.10 nm2.58 Given the considerable spread of experimental data of the DPPC pressure-area isotherm, it is not possible to compare the simulation results with a single set of data. Nonetheless, we find that the dependence of the surface pressure with temperature as well as the reduction of the width of the plateau with increasing temperature is consistent with the experiments reported by Crane et al.7 using the captive bubble apparatus. We find that the high compressibility of the LC phase observed in our simulations is compatible with experimental measurements conducted at high compression rates.61 An atomistic representation of water and phospholipids provides a route to analyze the influence of electrostatic effects on the monolayer phase transition, from a microscopic approach. Further work in this direction is under way. Acknowledgment. The simulations were performed at the Imperial College High Performance Computing Service. This work was supported by EPSRC Platform Grant EP/G00465X. References and Notes (1) Simons, K.; Toomre, D. Nat. ReV. Mol. Cell Biol. 2000, 1, 31–41. (2) Unger, E. C.; Porter, T.; Culp, W.; Labell, R.; Matsunaga, T.; Zutshi, R. AdV. Drug DeliVery ReV. 2004, 56, 1291–1314. (3) Tardieu, A.; Luzzati, V.; Reman, F. J. Mol. Biol. 1973, 75, 711– 733. (4) Ma, G.; Allen, H. C. Langmuir 2006, 22, 5341–5349. (5) Hauser, H.; Pascher, I.; Pearson, R.; Sundell, S. Biochim. Biophys. Acta 1981, 650, 21–51. (6) Brumm, T.; Naumann, C.; Sackmann, E.; Rennie, A. R.; Thomas, R. K.; Kanellas, D.; Penfold, J.; Bayed, T. M. Eur. Biophys. J. 1994, 23, 289–295. (7) Crane, J. M.; Putz, G.; Hall, S. B. Biophys. J. 1999, 77, 3134– 3143. (8) Lee, A. Biochim. Biophys. Acta 1977, 472, 285–344. (9) Evans, R.; Williams, M.; Tinoco, J. Lipids 1980, 15, 524–533. (10) Yan, W.; Biswas, S. C.; Laderas, T. G.; Hall, S. B. J. Appl. Physiol. 2007, 102, 1739–1745. (11) Venegas, B.; Wolfson, M. R.; Cooke, P. H.; Chong, P. L.-G. Biophys. J. 2008, 95, 4737–4747. (12) Gerber, F.; Krafft, M. P.; Vandamme, T. F.; Goldmann, M.; Fontaine, P. Biophys. J. 2006, 90, 3184–3192. (13) Li, J.; Miller, R.; Vollhardt, D.; Weidemann, G.; Mo¨hwald, H. Colloid Polym. Sci. 1996, 274, 995–999.

Liquid Condensed to Liquid Expanded Transitions in DPPC (14) Tabak, S. A.; Notter, R. H.; Ultman, J. S. Biophys. J. 1976, 60, 117–125. (15) Horie, T.; Hildebrandt, J. J. Appl. Physiol. 1971, 31, 423–430. (16) Blomley, M. J. K.; Cooke, J. C.; Unger, E. C.; Monaghan, M. J.; Cosgrove, D. O. Br. Med. J. 2001, 322, 1222–1225. (17) Chetty, K.; Stride, E.; Sennoga, C. A.; Hajnal, J. V.; Eckersley, R. J. IEEE Trans. Ultrason., Ferroelectr., Freq. Control 2008, 55, 1333– 1342. (18) Duncan, S. L.; Larson, R. G. Biophys. J. 2008, 94, 2965–2986. (19) Lee, S.; Kim, D. H.; Needham, D. Langmuir 2001, 17, 5544–5550. (20) Tu, K.; Tobias, D. J.; Blasie, J. K.; Klein, M. L. Biophys. J. 1996, 70, 595–608. (21) Berger, O.; Edholm, O.; Ja¨hnig, F. Biophys. J. 1997, 72, 2002– 2013. (22) Feller, S. E. J. Phys. Chem. B 2000, 104, 7510–7515. (23) Lindahl, E.; Edholm, O. Biophys. J. 2000, 79, 426–433. (24) Venable, R. M.; Brooks, B. R.; Pastor, R. W. J. Chem. Phys. 2000, 112, 4822–4832. (25) Jensen, M.; Mouritsen, O. G.; Petersz, G. H. Biophys. J. 2004, 86, 3556–3575. (26) Marrink, S. J.; de Vries, A. H.; Mark, A. E. J. Phys. Chem. B 2004, 108, 750–760. (27) Marrink, S.; Risselada, J.; Mark, A. Chem. Phys. Lipids 2005, 135, 223–244. (28) de Vries, A. H.; Yefimov, S.; Mark, A. E.; Marrink, S. J. Proc. Natl. Acad. Sci. U.S.A. 2005, 102, 5392–5396. (29) Klauda, J. B.; Wu, X.; Pastor, R. W.; Brooks, B. R. J. Phys. Chem. B 2007, 111, 4393–4400. (30) Leekumjorn, S.; Sum, A. K. Biochim. Biophys. Acta 2007, 1768, 354–365. (31) Sonne, J.; Jensen, M.; Hansen, F. Y.; Hemmingsen, L.; Peters, G. H. Biophys. J. 2007, 92, 4157–4167. (32) Mauk, A. W.; Chaikof, E. L.; Ludovice, P. J. Langmuir 1998, 14, 5255–5266. (33) Dominguez, H.; Smondyrev, A. M.; Berkowitz, M. L. J. Phys. Chem. B 1999, 103, 9582–9588. (34) Baoukina, S.; Monticelli, L.; Marrink, S. J.; Tieleman, D. P. Langmuir 2007, 23, 12617–12623. (35) Knecht, V.; Mu¨ller, M.; Bonn, M.; Marrink, S. J.; Mark, A. E. J. Chem. Phys. 2005, 122, 024704. (36) Adhangale, P. S.; Gaver, D. P. Mol. Phys. 2006, 104, 3011–3019.

J. Phys. Chem. B, Vol. 114, No. 3, 2010 1335 (37) Rose, D.; Rendell, J.; Lee, D.; Nag, K.; Booth, V. Biophys. Chem. 2008, 138, 67–77. (38) Laing, C.; Baoukina, S.; Tieleman, D. Phys. Chem. Chem. Phys. 2009, 11, 1916–1922. (39) Williams, A. D.; Wilkin, J. M.; Dluhy, R. A. Biophys. J. 1995, 102, 231–245. (40) Bresme, F.; Faraudo, J. Langmuir 2004, 20, 5127–5137. (41) Faraudo, J.; Bresme, F. Phys. ReV. Lett. 2004, 92, 236102. (42) Faraudo, J.; Bresme, F. Phys. ReV. Lett. 2005, 94, 077802. (43) Ryckaert, J.; Bellemans, A. Chem. Phys. Lett. 1975, 30, 123–125. (44) Berendsen, H. J. C.; Postma, J. P. M.; van Gunsteren, W. F.; Hermans, J. Intermolecular Forces; Reidel: Dordrecht, The Netherlands, 1981. (45) Berendsen, H. J. C.; Grigera, J. R.; Straatsma, T. P. J. Phys. Chem. 1987, 91, 6269–6271. (46) Abascal, J. L. F.; Vega, C. J. Chem. Phys. 2005, 123, 234505. (47) Vega, C.; de Miguel, E. J. Chem. Phys. 2007, 126, 154707. (48) Darden, T.; York, D.; Pedersen, L. J. Chem. Phys. 1993, 98, 10089– 10092. (49) Essmann, U.; Perera, L.; Berkowitz, M.; Darden, T.; Lee, H.; Pedersen, L. J. Chem. Phys. 1995, 103, 8577–8592. (50) Yeh, I.; Berkowitz, M. J. Chem. Phys. 1999, 111, 3155–3162. (51) Ane´zo, C.; de Vries, A. H.; Ho¨ltje, H.-D.; Tieleman, D. P.; Marrink, S. J. J. Phys. Chem. B 2003, 107, 9424–9433. (52) Gonza´lez-Melchor, M.; Bresme, F.; Alejandre, J. J. Chem. Phys. 2005, 122, 104710. (53) Berendsen, H.; van der Spoel, D.; van Drunen, R. Comput. Phys. Commun. 1995, 91, 43–56. (54) Lindahl, E.; Hess, B.; van der Spoel, D. J. Mol. Model. 2001, 7, 306–317. (55) Smit, B. J. Chem. Phys. 1992, 96, 8639–8640. (56) Trokhymchuk, A.; Alejandre, J. J. Chem. Phys. 1999, 111, 8510– 8523. (57) Mansour, H. M.; Zografi, G. Langmuir 2007, 23, 3809–3819. (58) Roke, S.; Schins, J.; Mu¨ller, M.; Bonn, M. Phys. ReV. Lett. 2003, 90, 128101. (59) Sankaram, M.; Thompson, T. Biochemistry 1990, 29, 10676–10684. (60) Sun, F. Biophys. J. 2002, 82, 2511–2519. (61) Ahuja, R. C.; Moebius, D. Langmuir 1992, 8, 1136–1144. (62) Lee, Y.-L.; Lin, J.-Y.; Chang, C.-H. J. Colloid Interface Sci. 2006, 296, 647–654.

JP9061303

Molecular Dynamics Simulations of Liquid Condensed to Liquid

Recommend Documents