Article Information

• PDF (352 kb)
• Full Text with Large Figures
• Cited by in Scopus (2)

PubMed

• Articles by Yoel Rodríguez
• Articles by Roman Osman

Related Articles

• Correction

  Correction Biophysical Journal, Volume 87, Issue 2, 1 August 2004, Pages 1398 Full Text | PDF (43 kb)

• Influence of Monolayer-Monolayer Coupling on the Phase Behav...

  Influence of Monolayer-Monolayer Coupling on the Phase Behavior of a Fluid Lipid Bilayer Biophysical Journal, Volume 93, Issue 12, 15 December 2007, Pages 4268-4277 Alexander J. Wagner, Stephan Loew and Sylvio May Abstract We suggest a minimal model for the coupling of the lateral phase behavior in an asymmetric lipid membrane across its two monolayers. Our model employs one single order parameter for each monolayer leaflet, namely its composition. Regular solution theory on the mean-field level is used to describe the free energy in each individual leaflet. Coupling between monolayers entails an energy penalty for any local compositional differences across the membrane. We calculate and analyze the phase behavior of this model. It predicts a range of possible scenarios. A monolayer with a propensity for phase separation is able to induce phase separation in the apposed monolayer. Conversely, a monolayer without this propensity is able to prevent phase separation in the apposed monolayer. If there is phase separation in the membrane, it may lead to either complete or partial registration of the monolayer domains across the membrane. The latter case which corresponds to a three-phase coexistence is only found below a critical coupling strength. We calculate that critical coupling strength. Above the critical coupling strength, the membrane adopts a uniform compositional difference between its two monolayers everywhere in the membrane, implying phase coexistence between only two phases and thus perfect spatial registration of all domains on the apposed membrane leafs. We use the lattice Boltzmann simulation method to also study the morphologies that form during phase separation within the three-phase coexistence region. Generally, domains in one monolayer diffuse but remain fully enclosed within domains in the other monolayer. Abstract | Full Text | PDF (618 kb)

• Concentration Effects of Volatile Anesthetics on the Propert...

  Concentration Effects of Volatile Anesthetics on the Properties of Model Membranes: A Coarse-Grain Approach Biophysical Journal, Volume 88, Issue 3, 1 March 2005, Pages 1524-1534 Mónica Pickholz, Leonor Saiz and Michael L. Klein Abstract To gain insights into the molecular level mechanism of drug action at the membrane site, we have carried out extensive molecular dynamics simulations of a model membrane in the presence of a volatile anesthetic using a coarse-grain model. Six different anesthetic (halothane)/lipid (dimyristoylphosphatidylcholine) ratios have been investigated, going beyond the low doses typical of medical applications. The volatile anesthetics were introduced into a preassembled fully hydrated 512-molecule lipid bilayer and each of the molecular dynamics simulations were carried out at ambient conditions, using the NPT ensemble. The area per lipid increases monotonically with the halothane concentration and the lamellar spacing decreases, whereas the lipid bilayer thickness shows no appreciable differences and only a slight increase upon addition of halothane. The density profiles of the anesthetic molecules display a bimodal distribution along the membrane normal with maxima located close to the lipid-water interface region. We have studied how halothane molecules between the two maxima of the bimodal distribution and we observed a different mechanism at low and high anesthetic concentrations. Through the investigation of the reorientational motions of the lipid tails, we found that the anesthetic molecules increase the segmental order of the lipids close to the membrane surface. Abstract | Full Text | PDF (429 kb)

• …more

Copyright © 2007 The Biophysical Society. All rights reserved.
Biophysical Journal, Volume 92, Issue 9, 3071-3080, 1 May 2007

doi:10.1529/biophysj.106.089078

Biophysical Theory and Modeling

Association Free Energy of Dipalmitoylphosphatidylserines in a Mixed Dipalmitoylphosphatidylcholine Membrane

Yoel Rodríguez, Mihaly Mezei and Roman Osman^,  

Department of Molecular Physiology and Biophysics, Mount Sinai School of Medicine, New York, New York

Address reprint requests to Roman Osman, Tel.: 1-212-241-5609; Fax: 1-212-860-3369.

Negatively charged phospholipids such as phosphatidylserines (PS) in a lipid membrane composed of neutral phosphatidylcholines (PC) provide an anchor for vitamin K-dependent zymogens and thus play an essential role in activating of the blood coagulation cascade 1,2,3,4. Upon vascular tissue damage, PSs are exposed on the surface of the epithelial cells, and in a Ca²⁺-dependent process blood-borne proteases—factors VII, IX, X, and prothrombin—are anchored to the negatively charged phospholipids through the N-terminal γ-carboxyglutamic acid containing domains (Gla domains). The importance of such an event has been amply demonstrated (e.g., 5) and it serves as a basis for anticoagulant therapy through the use of inhibitors that prevent the vitamin K-dependent posttranslational modification of glutamates to Gla 6.

The molecular details of the interaction between PS and Gla domains have been recently illustrated in the crystal structure of the Gla domain of prothrombin fragment 1 (PT1) in complex with a lysoPS 7. The structure shows that the anchoring of the negatively charged Gla domain to a negatively charged PS is mediated via Ca²⁺ ions. The high conservation of the Gla domains in other proteins suggests a similar mechanism of anchoring of proteins such as zymogene II, factor Xa, and cofactor Va 8,9. In addition to anchoring proteins via their Gla domains, PSs also regulate allosterically both factors Xa and Va 9,10,11, and have been proposed to act as a second messenger because they link platelet activation (IIa, collagen) to thrombin generation 5.

Electrostatic interactions are the driving force for binding peripheral proteins to negatively charged membranes 12. The thermodynamics of this process has been studied both experimentally and theoretically 13,14,15,16. The distribution of the negatively charged lipids in the membrane before protein association has to be taken into consideration to properly account for the entropic contribution of demixing to the binding free energy and absorption isotherms 17,18,19,20. The model of Denisov et al. 17 postulates that electrostatic interactions with small basic peptides produce lateral membrane domains enriched in acidic lipids. The model by Heimburg et al. 18, describes the influence of lipid redistribution upon protein adsorption on mixed lipid membranes. This model, with appropriate interaction parameters, predicts the experimental adsorption isotherm of cytochrome c on mixed dioleoyl phosphatidylglycerol/dioleoyl phosphatidylcholine bilayer membrane. The model of May et al. 19 starts with a mixture of negatively charged (e.g., PS) and zwitterionic (e.g., PC) lipids distributed homogenously, and introduces the contribution of lipid demixing to the free energy of protein-membrane association through an entropic term. The minimization of the free energy functional establishes the relationship between electrostatic and entropic contributions in the binding process between charged proteins and lipids. The critical assumption of the model is the homogenous distribution of the lipids at initial conditions, despite the fact that the PS/PC mixtures are nonideal 21. Therefore, a better estimation of the initial conditions would be required to properly estimate the binding free energy of proteins to charged lipids.

In mixed bilayers, such as those formed by PC and PS molecules, lateral phase separation occurs in the presence of Ca²⁺ and Mg²⁺ due to lipid-ion interactions as well as due to lipid-lipid interactions 22,23,24,25. Previous works have dealt with the problem of the effect of Ca²⁺ ions in leading to aggregation in membranes containing negatively charged lipids 26,27 demonstrating that in liposomes made from PS, addition of Ca²⁺ leads to aggregation, followed by vesicle fusion and leakage. It has been recently shown that the effectiveness of cations in inducing aggregation and fusion in N-palmitoylphosphatidylethanol-amine and N-palmitoylphosphatidylserine liposomes is Ca²⁺>Mg²⁺ ≈ Na⁺28. Likewise, a calcium-mediated association between the carbohydrate headgroups of the galactosylceramide-I³-sulfate and galactosylceramide has been demonstrated by vesicle aggregation and electrospray ionization mass spectroscopy 29,30,31. Thus, the involvement of Ca²⁺ in lipid aggregation plays an important and as yet not fully understood role.

We have designed a lipid membrane system to study the energetics of lipid-lipid association by using molecular dynamics simulation, free energy perturbation (FEP), and thermodynamic integration (TI) 32,33. The aim of this work is to elucidate the association of two dipalmitoylphosphatidylserines (DPPS) in the environment of a dipalmitoylphosphatidylcholine (DPPC) membrane and estimate the role of Ca²⁺ in this process. The first part of the article describes the methods and the models used in our studies. The free energy thermodynamic integration formalism and the free energy perturbation method are briefly reviewed. We then present the results of simulations of the lipid membrane with and without Ca²⁺ ions. We present and discuss the free energy calculations of the association between two DPPSs in the environment of DPPCs and use the results to estimate the effect of lipid-lipid interaction on the nonideality of the mixed membrane. We find that the association between DPPSs in the environment of DPPCs is favorable only in the presence of Ca²⁺ ions.

The construction of the lipid membrane followed a standard protocol 34. The general strategy is to randomly select phospholipids from a preequilibrated and prehydrated library of DPPC generated by Monte Carlo simulations in the presence of a mean field 35,36,37. The lipids are placed in a periodic system and the number of core-core overlaps between heavy atoms is reduced through systematic rotations around the z axis and translations in the xy-plane. All waters that overlap with the hydrocarbon interior of the bilayer, between ±12Å in the z-direction, were deleted. The initial configuration of the model comprised 48 DPPCs and 2009 water molecules, which represent a hydration level of 52%. Two DPPCs in each layer were replaced by two DPPSs in such a way as to generate in both layers of the membrane two PSs separated by one PC in the middle (Scheme I). The radial distribution function of the phosphorous atoms of DPPC yields a distance of ∼6Å between the P atoms, which was used to set up the initial configuration. To prepare the solvated DPPS in the DPPC environment, we removed all the water and performed grand canonical ensemble Monte Carlo simulations (GCMC) to solvate the system again. At the end, the system consisted of 44 DPPC, four DPPS, and 1960 water molecules, which amounts to a total of 12,092 atoms. To neutralize the system, we substituted two water molecules by two calcium ions (Ca²⁺). In the new initial configuration the calcium ions were placed approximately in the middle between two DPPS. The initial distances between the phosphorus atoms of each DPPS and the Ca²⁺ in the upper layer (lower layer) was 12.54Å (9.76Å) and 11.06Å (10.30Å) and between the DPPS was 11.63Å (11.62Å). We decided on this initial placement of the Ca²⁺ ions after having seen in preliminary simulations that Ca²⁺ ions placed far from the DPPSs invariably moved into their vicinity. At this level of treatment we also decided against using added salt to represent physiological ionic strength since at this system size it would only add 5–10 ion pairs, resulting in increased fluctuations and thus the need for excessively long runs.

Display large version of this figure

Scheme I

Interconversion between DPPS and DPPC form the PS/PC/PS to a PS/PS/PC configuration in the environment of PC. The most probable distance, obtained from the radial distribution function of the phosphorous atoms of DPPC, is ∼6Å.

Hexagonal periodic boundary conditions were applied in the xy-plane to maximize the distance between two periodic images 38. The edge of the hexagon and the length of the prism were 24.3Å and 75.0Å, respectively, with a center-to-center distance between neighboring cells of 42.1Å and the center of the bilayer membrane located at Z=0. The all-atom CHARMM27 force field for phospholipids 39 and the TIP3P water model 40 was used in all calculations. Long-range interactions were treated by a group-based spherical cutoff at 12Å. The van der Waals and electrostatic interactions were smoothly switched over a distance of 8.0Å. The equilibration was carried out using Langevin dynamics with decreasing planar harmonic constraints of 10, 5, 1, 0.5, 0 kcal/mol/Ų on water and lipids for 25ps in each stage, so that by the end of 125ps, the full system was completely unrestrained. The GCMC simulation to obtain the solvation of the mixed system was then performed for 9.5million MC steps. The chemical potential value used in the (μ, V, T) ensemble was obtained from a set of preliminary test runs, in which the water density in a 10-Å layer farthest from the lipids was monitored while the excess chemical potential was varied through the choice or the parameter B as expressed in the equation where is the Boltzmann constant, T the absolute temperature of the system, and is the average number of waters. The final value of produced an average density of 0.99719 g/ml. This was followed with a constant pressure, temperature, and area protocol (CPTA) production molecular dynamics for 5ns. The temperature of the system was maintained at 330K, which is above the gel-liquid crystal phase transition of DPPC. In the production stage the temperature was maintained using the Nose-Hoover scheme. The length of all bonds involving hydrogen atoms was kept fixed with the SHAKE algorithm 41. The equations of motion were integrated with a time step of 2 fs. All simulations were performed using the CHARMM program 42 and the Metropolis Monte Carlo program (MMC) 43.

A conventional approach to the evaluation of free energy of association in fluid mixtures relies on the calculation of the free energy profile of the selected species as a function of the distance that separates them. For lipid bilayers such a calculation presents formidable difficulties due to the extremely slow lateral diffusion of lipid molecules. We therefore decided to perform simultaneous exchanges of lipid headgroups from a (PS/PC/PS)_PC configuration into a (PS/PS/PC)_PC configuration, as shown in Scheme I. This results in the conversion of a PC-separated pair of PSs into a near-neighbor adjacent pair. In addition to eliminating the calculation of a distance-dependent free energy profile, this approach has the additional advantage that the mutations involve only the headgroups since the hydrocarbon chains of DPPC and DPPS are identical. Because the lipid portions of PC and PS are the same we used a dual topology of the PS-PC hybrid that involves only the headgroups (see Fig. 1). In this approach, the parts of the system, which are not the same in the initial and the final states, coexist at all times as the free energy simulation is carried out. They interact with the environment but not with each other 44. Thus, in the potential energy of this system, expressed as a function of a parameter λ that describes such a transformation, only the phosphocholine of the PC and phosphoserine of the PS are weighted by λ:

(1)

Display large version of this figure

Figure 1

Local dual topology of the PS-PC hybrid residue. The conformation shown is the starting geometry used in the dynamics. The atoms of the choline headgroup (N, C₁₂, C₁₁, blue) are superimposed on the corresponding atoms in the serine headgroup (N, C_A, C₄, green). The atoms were slightly shifted to make them visible.

The free energy calculations for two of the three runs were performed in two steps. When the calculation was performed in two steps, step 1 represented the interconversion of PS/PC/PS to PS/PS/PS and step 2 the interconversion of PS/PS/PS to PS/PS/PC. The potential energy of the systems is expressed in Eq. (1) with the steps represented by the superscripts. When the calculation is performed in one step, i.e., PS/PC/PS to PS/PS/PC, the subscripts a and b represent both residues; that is, a=PC/PS and b=PS/PC. and are the contributions of PS, PC, and the rest of the system, respectively 33,45. Phosphates of the two lipids as well as the first carbon of the glycerol with its hydrogens are included in the dual topology because they belong to the same group and have different partial charges in the CHARMM force field 39 (see Fig. 1).

The free energy simulations at different λ-values were performed using the BLOCK module in the CHARMM program 46,47. Three blocks were defined: one for the reactant, one for the product, and another one for the rest of the system. Blocks 2 and 3 consisted of the choline and phosphate groups of PC and serine and phosphate groups of PS, respectively, as described in Fig. 1. The rest of the system forms block 1. The interaction energy between block 2 and block 3 was set to zero to eliminate unphysical interaction terms.

Both the free energy perturbation and the thermodynamic integration methods were used to compute the free energy for each λ-window on the same runs. The free energy difference for each step using FEP was computed with the following equation:

(2)

When TI is used, the total free energy difference between λ=0 and λ=1 is:

(3)

For linear λ-dependence, the derivative of the potential energy with respect to λ for each perturbation is:

(4)

Thus, the change in the free energy is given by:

(5)

Both methods have been shown to reproduce experimental values of the free energy differences in several systems (e.g., 48,49,50).

Three different free energy simulation runs were performed. In the first two the change in one of the membrane layers was calculated: in run 1 the lower-layer system was used and in run 2 the perturbation was in the upper-layer system. The free energy calculations for both runs were carried out in two steps. Run 3 corresponds to the lower-layer system as well, but with different initial conditions and the calculations of the free energy were executed in only one step.

The value of λ was incremented from 0.1 to 0.9 in steps of 0.2 in FEP simulations. In TI, two sets of λ-values were used: the same values as in FEP for the system with calcium and those dictated by a five-point Gaussian quadrature for the system without calcium. The latter values were λ=0.046910, 0.230765, 0.5, 0.769235, 0.953089 51. Initially, the system was heated and equilibrated for 100ps at λ=0.1. At each λ-value, the system was reequilibrated for 50ps and data were collected for another 100ps, during which the trajectory was recorded every 50 fs, producing 2000 frames for each λ-value.

The trapezoidal rule was used to evaluate the integral from the discrete values of the free energy derivative between λ=0.1 and λ=0.9 and between λ=0.046910 and λ=0.953089. The values of the five-point quadrature coefficients, c_i, were 0.118463, 0.239314, 0.284444, 0.239314, and 0.118463 51. The coefficient for the first set of λ-values is constant and equals 0.2.

To improve the evaluation of the integral (Eq. (5)), three different schemes to obtain the end-point contributions (regions near λ=0 and λ=1) were applied. Scheme 1 assumes that the free energy derivative is constant and does not change in these intervals and the second derivative of the free energy is zero. Scheme 2 calculates the end-point contributions by extrapolation assuming that the third derivative of the free energy is constant. In this scheme it is assumed that the end-point values exist and are finite. Scheme 3 is based on the fact that when the Lennard-Jones interaction energy is scaled linearly, the free energy in three dimensions behaves like 33,52. Thus, the free energy derivative is proportional to The end-point contributions and are calculated by:

(6)

The density profile of the different components in both layers is nearly symmetrical (Fig. 2). The number density of water outside the lipid is ∼0.0334Å⁻³, corresponding to a bulk solvent density, and it approaches zero in the hydrocarbon core region. The P and N atoms of the PC and PS headgroups are located approximately at the boundary between the lipid and the aqueous environment. The distribution of the angle between the phosphorous-nitrogen vector, and the outwardly directed bilayer normal shows that the headgroups are approximately parallel to the membrane plane. The angle is 83°±23° for DPPC and 74°±10° for DPPS. The density of the hydrocarbon chain is reduced near the center of the membrane in agreement with other experiments and models of the bilayer system. The general features of the lipid density profile are similar to those observed for pure DPPC bilayers 53,54.

Display large version of this figure

Figure 2

Density profile of the main components for the lipid membrane without calcium ions along the z axis. The distribution of waters (red), hydrocarbon chains (green), P (yellow), N (pink) atoms of PC, and P (cyan), N (orange), and carboxyl group (gray) of PS are shown. In panel B the distribution of the components of the PC and PS headgroups represented in panel A are magnified.

The values of the free energy difference ΔG_{PS/PC/PS→PS/PS/PC} obtained with both methods FEP and TI are summarized in Table 1. Runs 1 and 2 represent independent free energy evaluations in the lower and upper layers, respectively. The different columns of the TI method summarize the results from the three schemes used to evaluate the integrand at the boundaries when λ→0 or λ→1 (see Methods and Models). All the results show that the free energy of PS association is positive indicating an unfavorable process. Further examination of the TI results shows that although the values are somewhat different from each other, the standard deviations make them statistically indistinguishable. The FEP yields the smallest free energy difference and the largest standard deviation mostly because of the large difference between the free energy values for the upper and lower layers. The presence of the exponential function in Eq. (2) results in a significant increase in the fluctuations of the calculated averages. This can lead to poorer convergence of the calculated free energy difference. Its mean value is also included in the energy interval obtained with TI. The contributions of the end points to the free energy difference are also listed in Table 1. The sum of both end-point contributions to the free energy is positive for all runs of both methods, FEP and TI.

Table 1 Free energy differences for the exchange PS/PC/PS → PS/PS/PC in a DPPC membrane in the absence of calcium ions

		ΔG / (kcal/mol) [ΔGfirst end point-0.0 / ΔG1.0-second end point]
		TI method
	Run No.*	Scheme 1	Scheme 2	Scheme 3	FEP method
	1	19.34 [49.88 / −44.62]†	19.94 [55.29 / −49.42]†	23.94 [91.02 / −81.15]†	9.52 [78.76 / −67.69]‡
	2	25.21 [45.67 / −37.45]†	26.04 [50.46 / −41.41]†	31.66 [82.20 / −67.54]†	24.80 [80.69 / −68.44]‡
	Mean±SD§	22.28 (4.15)	22.99 (4.31)	27.80 (5.46)	17.16 (10.80)

Five-point Gaussian quadrature (λ=0.04691, 0.230765, 0.5, 0.769235, 0.953089) was used to calculate the free energy differences using the linear TI method. In the FEP method five points (λ=0.1, 0.3, 0.5, 0.7, 0.9) were used to calculate the free energy differences.

* Runs 1 and 2 correspond to the calculation of the free energy in the lower-layer and the upper-layer systems, respectively. The lipid membrane is symmetric. † The values in square brackets are the and end-point contributions to the free energy differences, respectively. ‡ The values in square brackets are the and end-point contributions to the free energy differences, respectively. § The standard deviation was computed as and is given in parentheses; n is the number of runs.

The total TI integrands (i.e., the contribution of the two steps, from PS/PC/PS to PS/PS/PC) of the TI method for both runs as function of λ are represented in Fig. 3. The values at the end points were obtained by extrapolation using Scheme 2 of integration (see Methods and Models). The dependence of the integrand on λ is very smooth, indicating that the TI method is a good choice for evaluating the free energy.

Display large version of this figure

Figure 3

TI integrand for the transition PS/PC/PS → PS/PS/PC. Five-point Gaussian quadrature fit polynomials were used to calculate the free energy differences in a lipid membrane without calcium ions. The values at the end points were obtained by extrapolation (see text). Run 1 (♦, lower layer) and run 2 (●, upper layer).

As expected, these results demonstrate that the association of two PS in the environment of PC without counterions is unfavorable. This binary unfavorable PS-PS association is clearly due to the electrostatic repulsion between the headgroups. However, it does not offer an opportunity to examine the forces that may lead to clustering in an ensemble of lipids.

The general density profiles of the water, the hydrocarbon chain, and the P and N atoms of the PC headgroup are similar to those in the lipid membrane without calcium ions. However, the distribution of the P and N atoms, and the carboxyl group of the PS headgroup in both layers is quite different. In both layers the distribution of the N atoms, the P atoms, and the carboxyl groups have two populations (Fig. 4). This is because the individual distributions represent different PSs, which are separated vertically. In both layers the distributions of calcium ions have only one peak, centered at 17.5Å in the upper layer and at 20Å in the lower layer. The distribution of calcium ion in the lower layer is located between the two peaks of the P atom distribution, whereas in the upper layer it coincides with one of the P atom distributions closer to the membrane. The distribution of the of PC with respect to the bilayer normal has a value of 82°±23° (average over all PC). The PS headgroups of the lower layer prefer an orientation more outwardly directed ( vector ). DPPS most likely prefers this conformation due to the negative charge of the serine carboxyl group, which prefers to point out of the membrane. Nevertheless, the PS headgroups of the upper layer have the same orientation as the PC headgroups (). Further study of the orientation of the lipid headgroups may be required 54,55. Due to the structural differences observed in the two layers we calculated the free energy difference in both layers separately. The results are presented in the next section.

Display large version of this figure

Figure 4

Density profile of the main components for the lipid membrane with calcium ions along the z axis. The distribution of the waters (red), hydrocarbon chains (green), P (yellow), N (pink) atoms of PC, and P (cyan), N (orange), and carboxyl group atoms of PS (gray), and Ca²⁺ (dark blue) are shown. In panel B the distribution of the components of the PC and PS headgroups represented in panel A are magnified.

The potential energy derivative, for each λ in run 1 (i.e., the upper layer) as a function of simulation time is shown in Fig. 5. These curves represent the cumulative average of for all λ-values. The equilibration period during the initial 50ps shows some fluctuations in the energy, which relax to the new λ-value and stabilize at the appropriate value of the integrand. This behavior is similar for other runs.

Display large version of this figure

Figure 5

TI integrand, as a function of simulation time (ps) for all λ-values, run 1.

The presence of Ca²⁺ ions changes the free energy difference of the association between two PSs from a repulsive interaction to an attractive one (Table 2). Runs 1 and 2 represent independent free energy evaluations in the lower and upper layers carried out in two steps, respectively, and run 3 corresponds to the lower-layer system performed in only one step. The three schemes used to evaluate the integrand in the TI method show that the free energy of PS association in the presence of Ca²⁺ ions is negative, indicating a favorable process. The FEP method also yields negative values of free energy differences. Even though the mean values are not very different from each other for each method (TI and FEP), the standard deviations also make them statistically equivalent.

Table 2 Free energy differences for the exchange PS/PC/PS → PS/PS/PC in a DPPC membrane in the presence of calcium ions.

		ΔG / (kcal/mol) [ΔG0.1–0.0 / ΔG1.0–0.9]
		TI method
	Run No.*	Scheme 1	Scheme 2	Scheme 3	FEP method
	1	−7.78 [47.39/ −49.59]	−7.52 [57.39/ −59.34]	−6.95 [120.71/ −122.09]	−5.01 [69.67/ −70.65]
	2	−7.26 [43.95/ −44.83]	−7.82 [51.75/ −53.20]	−9.78 [103.29/ −106.70]	−5.67 [73.04/ −72.10]
	3	−8.01 [39.19/ −41.40]	−8.82 [47.17/ −50.19]	−12.82 [98.50/ −105.53]	−9.80 [55.67/ −54.51]
	Mean±SD†	−7.68 (0.38)	−8.05 (0.68)	−9.85 (2.96)	−6.83 (2.60)

Five points (λ=0.1, 0.3, 0.5, 0.7, 0.9) were used to calculate the free energy differences using the FEP and TI methods. The values in square brackets are the and end-point contributions to the free energy differences, respectively.

* Runs 1 and 2 correspond to the calculation of the free energy in the lower-layer and the upper-layer systems, respectively. Run 3 corresponds to the lower-layer system as well, but different initial conditions were set up for the dynamics. The lipid membrane is symmetric. † The standard deviation was computed as and it is given in parentheses; n is the number of runs.

The picture emerging from the analysis of the PS headgroup orientation (see above distribution of the ) shows that the PSs in the two layers were sampling two neighboring substates of the system. However, the calculated ΔG values for both layers (Runs 1 and 2 in Table 2) are essentially the same within the error limit, reinforcing the overall conclusion that the association of two PSs in the presence of Ca²⁺ ions is favorable. Although this study did not explore extensively the dependence of the free energy on the relative configurations of the lipids, it appears that it is not strongly dependent on the particular headgroup orientation. This is in clear distinction from the simulations in the absence of Ca²⁺, where the free energy shows a much stronger dependence on the particular arrangement of PS headgroups in the lipid membrane (Table 1).

The free energy derivative and the cumulative free energy change as a function of λ are shown in Figure 6AB. The free energy resulting from the different schemes of TI and FEP varies between −6.8 and −9.9 kcal/mol. These mean value extremes correspond to the FEP method and Scheme 2 of the TI method, respectively. Similarly to the simulations of the membrane without calcium, the value of the free energy difference obtained from FEP is also included in the interval where the values of the free energy using TI are found, that is, kcal/mol. Table 2 also lists the end-point contributions to the free energy difference at and The sum of both end-point contributions to the free energy is negative for all runs and schemes of TI.

Display large version of this figure

Figure 6

(A) TI integrand for the transition PS/PC/PS → PS/PS/PC. Five points (λ=0.1, 0.3, 0.5, 0.7, 0.9) were used to calculate the free energy differences in a lipid membrane with calcium ions. Diamonds, circles, and triangles correspond to run 1 (lower layer), 2 (upper layer), and 3 (lower layer with different initial conditions), respectively. (B) Free energy differences showing the similarity on the final values (see Table 2).

The decomposition of the free energy into its bonded (bond, angle, Urey-Bradley, dihedral and improper angles) and nonbonded terms (van der Waals and electrostatics) shows that the contributions of the bonded terms to the free energy difference are small with values ranging between −0.43 and 0.47 kcal/mol. Bond terms and improper torsions contribute negative values and the remaining terms are positive. The main contribution to the free energy difference comes from the nonbonded terms, dominated by electrostatics (ΔG_elect=−5.66 kcal/mol). This is not surprising because the free energy is dominated by the interaction between the two negative charges of the PS headgroups and the Ca²⁺ ions. Since decompositions of the free energy are path dependent the results above are specific to the path chosen for the free energy calculation. However, the difference among the bonded and nonbonded (van der Waals and electrostatic) terms are large enough to have general significance. They also agree with the intuitive picture of headgroup interactions in the presence of Ca²⁺ ions.

The ensemble generated in our simulation can serve in assessing whether lipid-lipid interactions contribute to nonideal mixing and cluster formation. The theoretical model of Huang et al. 21 writes the total potential energy of a triangular lattice of PC/PS membrane as made up of two types of terms: the long range electrostatic interactions between PS molecules, and the short-range non-electrostatic interactions between the lipids. Thus, within the lattice approximation the total potential energy is:

(7)

(8)

The first term in Eq. (7) is constant so it does not contribute to the nonideal mixing. The mixing behavior thus depends entirely on the last two terms of which the ΔE_m controls the mixing behavior at constant If ΔE_m is positive the interaction between like lipids is more attractive reducing the number of PS-PC contacts and leading to the formation of clusters. In contrast, for ΔE_m negative the two types of lipids prefer to interact with each other and the system will mix uniformly to increase the number of PS-PC contacts. Sufficiently large ΔE_m will overcome the electrostatic repulsion leading to cluster formation. Huang et al. 21 show that at =0.4 nonideal mixing is observed at ΔE_m of 0.6 kT and that at ΔE_m=0.4 kT nonideal mixing appears throughout the entire range of χ_PS. Thus, an estimation of ΔE_m from the simulation results can indicate whether the PS/PC ensemble explored here shows nonideal mixing.

We calculated ΔE_m by averaging the PS-PS, PC-PC, and PC-PS interaction energies over the simulation. In our calculations we only include the headgroups because both lipids have the same hydrocarbon tails. The potential energy terms for PC-PC and PC-PS were averaged over pairwise combinations present in the system. The PS-PS potential energy term was calculated from a simulation of the membrane with a Ca²⁺ to provide a close contact PS-PS pair. Since our PS/PC membrane is not a two-dimensional ideal lattice the value of ΔE_m() was calculated as a function of the number of neighbors. To evaluate the proper number of neighbors for averaging the interaction energies, we have calculated the radial distribution function (rdf) of PC headgroups (Fig. 7) The integration of the rdf shows a cluster of two PCs separated at 6Å for the P–P distance. The major peak in the rdf is spread between 12 and 14Å and consists of six to eight neighbors. We thus calculated the average U_PC-PC using progressively increasing number of neighbors starting with six.

Display large version of this figure

Figure 7

Radial distribution function for the P atoms of the PC headgroup lipids and the number of neighbors (PC) as a function of the phosphorous atom distance (P–P distance).

The results in Table 3 show that ΔE_m is always positive and increases from 0.30 to 0.63 kT for an increasing number of neighbors. This clearly reflects a tendency of cluster formation, which confirms the small cluster of two PCs observed in the rdf.

Table 3 Nonelectrostatic excess mixing energy, ΔEm

		* (kT)
		No. of neighbors
	Lipid pairs	6	7	8	9
	PC-PC	−1.59	−1.43	−1.21	−1.07
	PC-PS	−1.82	−1.57	−1.38	−1.22
	PS-PS†	−2.64	−2.64	−2.64	−2.64
	ΔEm / (kT)	0.30	0.47	0.55	0.63

* is the pair average van der Waals interaction energy calculated from a 500-ps trajectory. † The fixed value was obtained from a simulation of the PS/PC system in the presence of Ca2+.

It is important to comment on the value of U_PS-PS that we have used to estimate ΔE_m. This system contains two PSs embedded in an ensemble of 22 PCs, which allows for a good estimate of the U_PC-PC and U_PC-PS terms. The relatively low prevents a similar investigation of the dependence of U_PS-PS on the number of members. However, the distribution of PS-PS distances obtained from the simulation extends all the way to 12Å. Thus, we have used the average U_PS-PS to estimate the ΔE_m. Equation (8) and the results in Table 3 clearly set the limit of U_PS-PS that is required for the appearance of nonideal mixing as a function of the number of neighbors. The value that we estimate from our simulations always satisfies this limit. Thus, it appears that nonideal mixing is a realistic possibility in PS/PC mixed lipid systems.

We have simulated a DPPS/DPPC bilayer mixture with the aim of calculating the free energy of association between two PSs in the environment of PCs. The importance of the behavior of mixed lipid systems is that divalent cations such as Ca²⁺ or proteins that associate with the lipid lead to lateral demixing of PS and PC. Thus, the free energy of peripheral protein association with charged membranes composed of zwitterionic and acidic lipids mixture may depend on the nonideal nature of the lipid system presented at the initial conditions.

Our results suggest that the association between two PSs in the environment of PCs in the presence of Ca²⁺ ions is thermodynamically favorable, which agrees with previous studies on domain separation mediated by Ca²⁺ ions. Furthermore, the nonelectrostatic interactions between the lipid headgroups lead to clustering. A careful description of the initial conditions of the mixed lipids is therefore essential for a proper evaluation of the thermodynamics of protein-lipid association.