Print version ISSN 0365-0375
An. Asoc. Quím. Argent. vol.94 no.1-3 Buenos Aires Jan./July 2006
A study by molecular dynamics simulation of the effect of the ionic strength on the properties of a model DPPC/DPPS asymmetric membrane.
López Cascales, J.J. *
Universidad Politécnica de Cartagena, Centro de Electroquímica y Materiales Inteligentes (CEMI), Aulario II, Campus de Alfonso XIII, 30203 Cartagena, Murcia, Spain.
FAX: +34-968-325932, E-Mail: email@example.com
Received April 28th, 2006. In final form June 5th, 2006
Dedicated to Prof. Imre G. Csizmadia on the occasion of his 75th birthday
In this work we focused our interest on the study of the effect of the ionic strength of the solution on the properties of the lipids that form an asymmetric membrane composed by DPPC and DPPS-. The Molecular Dynamics Simulation was the simulation technique chosen to carry out the study. The asymmetry in the lipid bilayer was generated by placing the charged lipids (DPPS-) in only one of the two leaflets of the membrane. In the opposite leaflet, the layer was composed only by neutral lipids (DPPC). The ionic strength of the aqueous solution was simulated at two different concentrations of NaCl, corresponding to 0.25 and 0.5 M. From the analysis of the simulated trajectories, information related to steady and dynamics properties of the system with atomic detail was obtained. Thus, the translational diffusion coefficient of the lipids, lipid head orientation, lipid hydration or the structure of the hydrocarbon zone of the membrane has been studied.
En el presente trabajo, nuestro interés ha estado centrado en el estudio del efecto de la fuerza iónica de una solución acuosa sobre las propiedades de los lípidos que constituyen un modelo de membrana asimétrica compuesta por DPPC y DPPS-, mediante la técnica de simulación de dinámica molecular. El modelo de membrana fue generada al distribuir asimétricamente las moléculas de DPPS- en una de las dos capas de lípidos que forman la membrana. El efecto de la fuerza iónica fue estudiada para dos concentraciones diferentes de NaCl (0.25 y 0.5 M). Del análisis de las trayectorias de simulación generadas se obtuvo información con resolución atómica de diferentes propiedades tanto estacionarias como dinámicas de los lípidos en ambos lados de la membrana, tales como el coeficiente de difusión translacional, la orientación de las cabezas lipídicas respecto a la normal a la membrana o la ordenación de la membrana en su zona hidrocarbonada.
Phospholipids play a crucial role in biological cells controlling the diffusion of small molecules between the interior and the exterior of the membrane, and providing a proper environment for other molecules embedded or anchored to the membrane , such as proteins, polysaccharides, and others. In addition, in eukaryotic cells, lipids are asymmetrically distributed between both leaflets of the membrane cell. Thus, the phosphatidylserine (PS) lipids are asymmetrically located in the inner leaflet of the membrane where the lost of this asymmetry seems to be a hallmark of the early stages of the cell apoptosis and thus, dying cells are removed from the blood stream .
In this work, we focused our interest on the study of the behavior of an asymmetric lipid membrane composed by DPPC and DPPS-, at different ionic strength in an aqueous solution. In this setting, Molecular Dynamics Simulation [3,4] is the computational technique used for obtaining information with atomic detail of the behavior of the different components forming the system. Thus, during the last decades, numerous articles have been published in the field of MD simulations of symmetric bilayers [5,6,7,8,9,10]. More recently, with the improving of program codes and the raising of computing power, longer trajectories and bigger systems have been achievable [11,12] to further studies. However, the common point of all these works is related to the fact that all of them correspond to studies of symmetric membranes were the same type of phospholipids were distributed in both leaflets of the membrane. However, the study of asymmetric bilayers is crucial due to its relevance to biological membranes. In this setting and on the basis of the recently work of López Cascales et al.  in which a model of asymmetrical phospholipid bilayer has been developed, in this work we focused our interest on the study of the effect of the ionic strength of an aqueous solutions on the properties of the membrane.
Among other properties, the lipid hydration, lipid-lipid charge bridge interactions, atomic density across the membrane, lipid head orientation and translational diffusion coefficient of lipids will be investigated.
Method and models
With the goal of studying the effect of the ionic strength on the lipid properties in the membrane, four different simulations were carried out:
1. DPPC bilayer in water, divided in two layers of 36 DPPC lipids each.
2. 54DPPC/12 DPPS-/12Na+ in water divided in two leaflets of 30 DPPC in the first leaflet and 12 DPPS plus 24 DPPC in the second one.
3. 54DPPC/12 DPPS-/12Na+, in a 0.25M NaCl aqueous solution, with lipid distributed on both leaflets like in case 2.
4. 54DPPC/12 DPPS-/12Na+, in a 0.5M NaCl aqueous solution, with lipid distributed on both leaflets like in case 2The DPPC bilayer was generated as follows: A three dimension periodical system was generated, starting from a single DPPC molecule, figure 1, placed perpendicularly to the water layer. This DPPC molecule was randomly rotated and copied 36 times in both leaflets of the bilayer. The existing gap above and below the lipid bilayer was filled with 2501 water molecules from an steady and well equilibrated computing box of water of the SPCE water model .
The system composed by 54DPPC/ 12DPPS- / 12Na+, 12 DPPC were randomly substituted by 12 DPPS-, figure 1, in one of the two leaflets. In the another leaflet (composed only by DPPC), 6 DPPC molecules were removed so as to reproduce the surface area of a pure DPPC bilayer. Next, 12 water molecules were substituted by 12 Na+ to maintain the charge-balance of the system.
Figure 1. Atomic numbering of a single DPPC and DPPS lipid, such as they will be referred to in this article
The other two systems corresponding to 0.25M and 0.5 M in NaCl in aqueous solution, 24 and 48 additional water molecules were substituted by 24 (12 Na+ + 12 Cl-) and 48 (24 Na+ + 24 Cl-) Na respectively. figure 2 depicts a snapshot of the systems generated following the procedure described above.
Figure 2. A snapshot of the DPPC+ DPPS- bilayer. Big dark beads correspond to sodium ions (Na+), medium size light beads to chloride ions (Cl-)
and small gray beads to water molecules. DPPS- are plotted by sticks and DPPC by wire-frames
The engine of the MD simulations described in this work was the GROMACS 3.1.4 package [15,16]. The force field parameters used to run the simulations were the previously proposed by Marrink et al.  and Lopez Cascales et al.  in previous simulations of pure bilayers of DPPC and DPPS-. The electrostatic interactions were evaluated using the particle mesh Ewald method [17,18], using a 0.9 nm cutoff and the reciprocal space interactions were evaluated on a 0.12 nm grid with a fourth order spline interpolation.
After a steepest descent energy minimization procedure was applied to the whole system, 500 ps of simulation was carried out at the temperature of 500 K to remove undesired interaction between neighboring atoms. In this case, the leap frog integrator algorithm was used, with a time step of 4 fs (femtoseconds). The temperature coupling time was of 0.1 ps, and an anisotropic pressure external bath of 1 atm, with a coupling constant of 0.5 ps  was used. Finally, 20 ns of trajectory length was carried out at the temperature of 350K, well above those of the transition temperature 314 and 326 K for bilayers of DPPC [20,21] and DPPS- [22,23] bilayers in their liquid crystalline respectively.
In this setting, the first 5000 ps of our generated trajectories were discarded, due to the fact that this part of the trajectory was required to equilibrate the systems, in a good agreement with our previous simulations . The data from the last 15 ns of the trajectories were used to extract the properties that are presented in this work. All the simulations were carried out on a parallel HP ES 45 cluster Alpha Server with four processors, 3 days of intensive computing was required for each simulation.
Results and discussion
From the surface area of the DPPC of the leaflet in absence of DPPS-, the following values were obtained: 0.703, 0.68, 0.64 and 0.65. Each of these values corresponds to each of the four cases studied in this work: DPPC bilayer, asymmetric DPPC/DPPS bilayer in absence of salt, asymmetric bilayer with a 0.25M and 0.5M in NaCl respectively. We observe how the surface area of the DPPC bilayer agrees with the experimental results of a DPPC bilayers in absence of salt, considering the spread of results that we can get from the bibliography . For the neutral DPPC/DPPS membrane in the absence of salt, we observe how we are able to reproduce the surface are of a DPPC bilayer, avoiding the presence of artifacts in our simulations. Finally, for the two simulations corresponding to 0.25 and 0.5M in NaCl, a shrinking of the lipid surface is measured compared with the two other cases, which are in a good agreement with the reported results of DPPC bilayers in presence of salt [11,12]. In this regard, we must remark how for the case of 0.25M, it seems as if the system would have reached a saturation point compared with the case of 0.5M in NaCl
Figure 3 displays the atomic density across the phospholipid bilayers for the system at different salt concentrations. This figure shows how the asymmetry of the system, due to the presence of DPPS- in only one of the two leaflets, introduces a certain dehydration and asymmetry in the membrane, compared with a symmetric DPPC bilayer. This dehydration is enhanced by increasing the concentration of NaCl. This effect appears to be due to sodium ions displacing solvating water molecules and lipid-lipid bridge interactions between neighboring lipids.
Figure 3. Atomic density across the asymmetric lipid bilayer. The density of the DPPC bilayer was overlapped
on the different figures to observe the respective changes.
The radial distribution function g(r) is often employed to identify close interactions between neighboring atoms. In this regard, the radial distribution function g(r) is defined as:
g(r) = N(r)/ 4pr2rdr (1)
where, N(r) is the number of atoms in a spherical shell at distance r and thickness dr from a reference atom, and r is the number density taken as the ratio of atoms to the volume of the computing box.
To analyze the hydration under different salt concentration in both leaflets of the membrane, we calculated the g(r) around one of the two oxygens of the phosphate oxygen group (atom number 10 in figure 1). From numerical integration of the area of the radial distribution function under the first peak, we are able to obtain the hydration numbers around the atom 10 (phosphate oxygen) of the DPPC. Thus the values of 2.06, 1.95 and 1.96 were obtained for the above mentioned atom of DPPC in the leaflet without DPPS- and 1.96, 1.90 and 1.91 for DPPC in the leaflet with DPPS-, corresponding to the three studied cases: in the absence of salt, 0.25M and 0.5M of sodium chloride (NaCl) respectively. Concerning the DPPS-, we obtained the hydration numbers of 2.13, 1.97 and 2.03 for the same atom 10 (phosphate oxygen) in the absence of salt, 0.25M and 0.5M respectively. The simulations indicate that in all the cases, the DPPS- is slightly more hydrated than DPPC and almost no dehydration was observed from our simulations for the oxygen carbonyl group with the salt concentration.
In addition, to explore on the charge bridge between phospholipid-phospholipid, the radial distribution function g(r) was also used. Thus, the NH3 of the DPPS- (atom 1 in figure 1) was taken as the reference atom to calculate the radial distribution function, and several oxygen atoms of the DPPC as the atoms of coordination. figure 4 depicts the radial distribution function of two different DPPC atoms around the NH3 group: phosphate oxygen (atom 10, figure 1), and carbonyl oxygen group (atom 16, figure 1) corresponding to one of the two lipid tails. In the both cases, the bridges between neighboring lipids are strongly affected by the salt in solution. Thus, the bridges between the NH3 of the DPPS- and the phosphate oxygen and the oxygen carbonyl group depend of the salt concentration in solution: higher salt concentration promote the bridge formation for the deepest atoms of the membrane (oxygen of the carbonyl groups) but for lower salt concentration, bridge charges with phosphate oxygen groups prevail over the carbonyl groups. This effect is strongly related to the fact that higher salt concentration promotes the cation penetration inside the membrane and thus, the bridges with the carbonyl groups, in good agreement with the results proposed by Berkowitz et al  from a mixed bilayer of DPPC and DPPS.
Figure 4. Radial distribution function of several DPPC oxygens around the NH3 of the DPPS-, for several concentrations of NaCl.
As a consequence of the lipid-lipid and lipid-ion interactions, figure 5 depicts the distribution of head group center-of-mass for the leaflet containing DPPS-, in the presence and in the absence of NaCl salt. It appears that the presence of the salt induces the formation of domains, presumably due to cross-linking of the anionic DPPS- head group by sodium cations. Hence and due to this inhomogeneous distribution of lipids on the leaflet, transient domains with different density, hydration and electrostatic charge emerge in the leaflet.
Figure 5. DPPC and DPPS head distribution on the leaflet 2 containing DPPC and DPPS without NaCl (left)
and with 0.25M NaCl (right) respectively. DPPC (O) and DPPS (•)
With the goal of gaining an insight related the structure of the hydrocarbon zone of the lipid membrane, the deuterium order parameter SCD has been used. Thus, the order parameter -SCD is an experimental measurement of the disorder of a membrane, measured from NMR experiments. This order parameter is associated with the orientation of the hydrogen of a CH2 group with respect to the normal to the lipid bilayer. Due to the fact that our simulations do not explicitly take into account the hydrogens of the CH2 groups, the order parameter -SCD on the i +1 CH2 was defined as the unitary vector normal to the vector defined from the i to the i+2 CH2 group and contained in the plane formed by the groups i, i+1, and i+2. As a consequence, the -SCD on the ith CH2 was defined as follows:
-SCD = (3 cos2(q) - 1)/2 (2)
Where q is the angle substended by the unitary vector defined above and the Z axis. Note that the -SCD can adopt any value in a range from -0.5 (which corresponds to parallel to the lipid/water interface) to 1 (oriented along the normal to the membrane surface).
In this regard, Seelig et al.  and Browning et al.  provided experimental data of the order parameters of DPPC and DPPS- bilayers respectively, obtained from NMR experiments. The deuteron quadrupole splitting is related to the order parameter -SCD as follows:
DvQ = - ( 6e2qQ/8h )SCD (3)
where eq is the electrostatic field gradient parallel to the C-D bond direction in the laboratory coordinate frame, e is the electric charge of the electron, Q is the nuclear quadrupole moment and h is the Planck constant. The deuteron quadrupole constant e2qQ/h has been found to be 170 kHz [20,25]. Considering this value, we were able to obtain the order parameters that are introduced below.
Figure 6 shows the results obtained from simulation compared with the experimental data  obtained from NMR experiments of a DPPC bilayer at the temperature of 350K. Thus, for the leaflet containing only DPPC, we observe a perfect agreement between simulation and experimental data in the DPPC bilayer and asymmetric DPPC/DPPS bilayer in the absence of NaCl.
Figure 6. DPPC order parameter for both lipid leaflets: (a) corresponds to leaflet 1 without DPPS lipids and (b) corresponds to leaflet 2 with DPPS. Circles refer to our simulation results of the DPPC bilayer and stars to the experimental results of a DPPC bilayer .
When the salt concentration increased, the order parameters of DPPC placed in the leaflet in absence of DPPS increased, in good agreement with the simulation results of Berkowitz et al.  for a DPPC bilayer in the presence of NaCl. On the other hand, the order parameters of DPPC placed in the leaflet containing DPPS-, the presence of DPPS rises up the order parameter of the DPPC, compared with the values from the other leaflet, where there were only DPPC molecules. In all the cases, we obtained a good correlation between the increasing in the order parameter of the lipid hydrocarbon tails and the shrinking of the area per lipid head obtained above. On the other hand, focusing our interest on the order parameters of the DPPS- with the salt concentration, figure 7 compares the variation of simulated and experimental order parameters of the DPPS hydrocarbon tails at different salt concentrations . In this regard, we must point out the good agreement found between the simulation results and the experimental data in the absence of salt. In this case, the DPPS- order parameters follow a similar trend to that obtained for the POPS-  in the presence of salt, where the order parameter increases with the salt concentration.
Figure 7. DPPS- order parameter. Symbol * refers to experimental results .
Another aspect that was explored in this work was related to the lipid headgroup orientation as a function of salt concentration. In this regard, an angle q was defined as the angle of P-N vector of the lipid-heads with the outward normal of the bilayer. Thus, figure 8 depicts the angle distribution function for the asymmetric membrane. The leaflet composed by DPPC in absence of DPPS was called as leaflet 1, and the leaflet containing DPPS as leaflet 2. Our results agree with previous results [11,12] where they measure that P-N angle remains parallel to the membrane surface independently of the salt concentration.
Figure 8. Distribution of the angle between the vector joining the P-N vector of DPPC and the outward normal of the bilayer on both asymmetric phospholipid leaflets. Leaflet 2 correspond to the lipid layer composed by DPPC+ DPPS- and the leaflet 1 to the formed only by DPPC
Table 1 shows the mean angle q of DPPC in both leaflets. In this regard, a slight reduction of less than 5 degrees in the tilt of the P-N vector was estimated compared with the angle measured for a symmetric DPPC membrane. In this sense, we must remark that even when the mean angle remained almost invariable with the salt concentration, the shape of the angular distribution function of the DPPC (in the leaflet in presence of DPPS-) changed with the variation of the salt concentration. However, for the leaflet in absence of DPPS-, neither the mean angle nor the shape of the angular distribution function of the lipid head orientation was affected by the salt concentration in the range of ionic strength studied.
Table 1. Mean angle between the vector joining the P-N of the DPPC lipid head group and the outward normal to the lipid bilayer.
Translational diffusion coefficient of phospholipids.
Concerning the study of the effect of salt concentration on the dynamic properties of the lipids, we focused our interest on the study of the translational diffusion coefficient Dt of lipids in the membrane plane. This property can be calculated from the mean of the free square displacement of the center of mass of the lipids as follows:
< x2 + y2> = 4Dtlat t, (4)
< z2 > = 2Dttranv t
Where the x and y axes correspond to the axes on the membrane plane and the z axis is perpendicular to the membrane plane. However, in the case of a bilayer, a free displacement can be only considered in the x-y plane. The translational diffusion coefficient along the z axis is not relevant in these simulations because it doesn't fit to the explicit condition imposed by equation 4. Thus, discarding the first 1000 ps of simulation where rapid motions of the lipids take place associated to vibrational motions, instead of diffusive processes, we are able to estimate the translational diffusion coefficient of DPPC in both leaflets of the membrane for different salt concentrations. table 2 depicts the calculated diffusion coefficients on both leaflets at different salt concentrations. These data show that, in the leaflet without DPPS-, the translational diffusion coefficient of DPPC remains unaltered with the presence of salt. On the other hand, the diffusion of DPPC in the leaflet with DPPS-, increases from 1.19x10-6 cm2/s in the absence of salt to 1.85x10-6 cm2/s (around a 55 %) at 0.5M in NaCl. This behavior maybe related to the breaking of DPPC/DPPS- charge bridge interactions between neighboring lipids such as we described above, increasing in this way the freedom of the lipids to move on the lipid plane. Related the diffusion of DPPS-, the diffusion coefficient remained constant with the salt concentration. Even when these values are typically one order of magnitude higher than the experimental ones Dxy ~ 10-7-10-8 cm2/s ), we must consider the difference in temperature at which we performed our simulations.
Table 2. Translational diffusion coefficient in cm2/s parallel to the membrane plane. Leaflet 1 and 2 corresponds to the leaflets without and with DPPS respectively.
An asymmetric lipid membrane was studied by Molecular Dynamics simulations from several equilibrated models with atomic detail. From the analysis of the simulated trajectories, we realized the effect of the ionic strength of the solution on the structure of a model biological membrane. Thus, from our simulations, we observed how each leaflet of the lipid membrane follows a different behavior with the variation of the salt concentration. Thus, the first conclusion of this work is related to the fact that the behavior of the DPPC molecules on both leaflets behaves almost independently of each other. This conclusion was obtained after comparison of the translational diffusion coefficient, lipid head orientation, lipid-lipid interactions, lipid hydration and order parameters of DPPC in both leaflets.
In this sense, the presence of a charged lipid DPPS- in one of the two leaflets modified noticeably the properties of the neutral DPPC (sharing the same leaflet) compared with the DPPC of the opposite leaflet in the absence of DPPS-. This different behavior between DPPC of opposite leaflets arises from the ionic strength of the solution for both, dynamic and steady properties of the lipids. In this regard, the presence of salt in solution diminished the lipid-lipid interaction between neighboring lipids and as a consequence, an increasing of the translational diffusion coefficient of the DPPC was measured. In addition, a dehydrating of the lipid heads was evidenced when the salt concentration increased.
In summary, an increasing of the salt concentration perturbed the properties of the lipid in each leaflets, and on the other hand, each leaflet follows an independent behavior, in that the changes in one leaflet didn't correlate with the lipid behavior of the other.
The author acknowledges the Spanish Ministry of Education and the Seneca Foundation of the Murcia Region for financial support, projects BQU-2001/0477 and AR-8-0271/FS/02, respectively.
 Yagle, P. Ed. The structure of Biological Membranes; CRC Press, Inc.,1991. [ Links ]
 Fadok, V.; Bratton, D.; Rose, D.; Pearson, A.; Ezekewitz, R.; Henson, P. Nature 2000, 405, 85. [ Links ]
 van Gunsteren, W.F.; Berendsen, H.J.C., Angew. Chem Int. (Ed. Engl.), 1990, 29, 992. [ Links ]
 Frenkel, D.; Smit, B., Understanding Molecular Simulations; Academic Press, 2002. [ Links ]
 Bassolino-Klimas, D.; Alper, H.; Stouch, T., Biochemistry 1993, 32, 12624. [ Links ]
 Heller, H., Schaefer, M.; Schulten, K., J. Phys. Chem. 2003, 97, 8343. [ Links ]
 Alper, H.; E. Bassolino-Klimas, D.; Stouch, T.R., J. Chem. Phys. 1993, 99, 5554. [ Links ]
 Egberts, E.; Marrink, S.J.; Berendsen, H.J.C., Eur. Biophys. J. 1994, 22, 423. [ Links ]
 Tieleman, D.; Marrink, S.J.; Berendsen, H.J.C., Biochemica et Biophysica Acta (BBA) 1994, 1331, 235. [ Links ]
 López Cascales, J.J.; García de la Torre, J.; Marrink, S.J.; Berendsen, H.J.C., J. Chem. Phys. 1996, 104, 2713. [ Links ]
 Pandit, S.; Bostick, D.; Berkowitz, M., Biophysical Journal 2003, 84, 3743. [ Links ]
 Sachs, J.; Nanda, H.; Petrache, H.; Woolf, T., Biophysical Journal 1980, 86, 3772. [ Links ]
 López Cascales, J.J.; Otero, T.F.; Smith, B.; Gonzalez, C.; Marquez, M., J. Phys. Chem. B 2006, 110, 2358. [ Links ]
 Berendsen, H.J.C.; Grigera, J.; Straatsma, T., J. Phys. Chem. 1987, 91, 6269. [ Links ]
 Lindahl, E.; Hess, B.; van der Spoel, D., J. Mol. Mod. 2001, 7, 306. [ Links ]
 Berendsen, H.J.C.; van der Spoel, D.; van Drunen, R., Comp. Phys. Comm. 1995, 91, 43. [ Links ]
 Darden, T.; York, D.; Pedersen, L., J. Chem. Phys. 1993, 98, 10089. [ Links ]
 Essmann, U.; Perea, L.; Berkowitz, M.; Darden, T.; Lee, H.; Pedersen, L., J. Chem. Phys. 1995, 103, 8577. [ Links ]
 Berendsen, H.J.C.; Postma, J.P.M.; van Gunsteren, W.F.; DiNola, A.; Haak, J.R., J. Chem. Phys. 1984, 8, 3684. [ Links ]
 Seelig, A. Seelig, J. Biochemistry 1974, 23, 4839. [ Links ]
 Young, L.R.D.; Dill, K.A., Biochemistry 1993, 27, 5281. [ Links ]
 Cevc, G.; Watts, A.; Marsh, D., Biochemistry 1981, 20, 4955. [ Links ]
 Hauser, H.; Paltauf, F.; Shipley, G.G., Biochemistry 1982, 21, 1061. [ Links ]
 Nagle, J.; Tristam-Nagle, S., Biochimica et Biophysica Acta 2000, 1469, 159. [ Links ]
 Browning, J. L.; Seelig,, J. Biochemistry 1980, 19, 1262. [ Links ]
 Brown, M.F., J. Phys. Chem. 1982, 3, 1576. [ Links ]
 Mukhopadhyay, L.M.P.; Tieleman, D., Biophysical Journal 2004, 86, 1601. [ Links ]