A criterion to identify the equilibration time in lipid bilayer simulations

Rodolfo D. Porasso,¹ J.J. Lopez Cascales,²

 

With the aim of establishing a criterion for identifying when a lipid bilayer has reached steady state using the molecular dynamics simulation technique, lipid bilayers of different composition in their liquid crystalline phase were simulated in aqueous solution in pres-ence of CaCl₂ as electrolyte, at different concentration levels. In this regard, we used two different lipid bilayer systems: one composed by 288 DPPC (DiPalmitoylPhosphatidyl-Choline) and another constituted by 288 DPPS (DiPalmitoylPhosphatidylSerine). In this sense, for both type of lipid bilayers, we have studied the temporal evolution of some lipids properties, such as the surface area per lipid, the deuterium order parameter, the lipid hydration and the lipid-calcium coordination. From their analysis, it became evident how each property has a different time to achieve equilibrium. The following order was found, from faster property to slower property: coordination of ions deuterium order parameter > area per lipid hydration. Consequently, when the hydration of lipids or the mean area per lipid are stable, we can ensure that the lipid membrane has reached the steady state.

I. Introduction

Over the last few decades, different computational techniques have emerged in different fields of sci-ence, some of them being extensively implemented and used by a great number of scientists around the globe. Among others, the Molecular Dynamics (MD) simulation is a very popular computational technique, which is widely used to obtain insight with atomic detail of steady and dynamic proper-ties in the fields of biology, physics and chemistry. In this regard, a critical aspect that must be iden-tified in all the MD simulations is related to the re-quired equilibration time to achieve a steady state.

*E-mail: rporasso@unsl.edu.ar

Instituto de Matemática Aplicada San Luis (IMASL) -Departamento de Física, Universidad Nacional de San Luis/CONICET, D5700HHW, San Luis, Argentina.

Universidad Politécnica de Cartagena, Grupo de Bioin-formática y Macromoléculas (BioMac) Aulario II, Cam-pus de Alfonso XIII, 30203 Cartagena, Murcia, Spain.

This point is crucial in order to avoid simulation ar-tifacts that could lead to wrong conclusions. Cur-rently, with the increment of the computing power accessible to different investigation groups, much longer simulation trajectories are being carried out to obtain reliable information about the systems, with the purpose of approaching the time scale of the experimental phenomena. However, even when this fact is objectively desirable without further ob-jections, nowadays, much longer equilibration times are arbitrarily being required by certain reviewers during the revision process. From our viewpoint, this should be thoroughly revised due to the follow-ing two main reasons: first, because it results in a limiting factor in the use of this technique by other research groups which cannot access to very ex-pensive computing centers (assuming that authors provide enough evidence of the equilibration of the system). Second, to avoid wasting expensive com-puting time in the study of certain properties which do not require such long equilibration times, once the steady state of the system has been properly identified.

Phospholipid bilayers are of a high biological rele-vance, due to the fact that they play a crucial role in the control of the diffusion of small molecules, cell recognition, and signal transduction, among oth-ers. In our case, we have chosen the Phosphatidyl-Choline (PC) bilayer because it has been very well studied by MD simulations [1-7] and experimen-tally as well [8-14]. Furthermore, studies of the effects of different types of electrolytes on a PC bi-layer have also been studied, experimentally [15-24] and by simulation [25-33].

As mentioned above, the Molecular Dynamics (MD) simulations have emerged during the last decades as a powerful tool to obtain insight with atomic detail of the structure and dynamics of lipid bilayers [34-36]. Several MD simulations of mem-branes under the influence of different salt con-centrations have been carried out. One of the main obstacles related to these studies has been the time scale associated to the binding process of ions to the lipid bilayer. Considering the literature, a vast dispersion of equilibration times associated to the binding of ions to the membrane has been reported, where values ranging from 5 to 100 ns have been suggested for monovalent and divalent cations [25, 27-29, 32, 37, 38]. In this regard, we carried out four independent simulations of a lipid bilayer formed by 288 DPPC in aqueous solutions, for different concentrations of CaCl₂ to provide an overview of their equilibration times. Among other properties, the surface area per lipid, the deuterium order parameters, lipid hydration and lipid-calcium coordination were studied.

Finally, in order to generalize our results, a bi-layer formed by 288 DPPS in its liquid crystalline phase, in presence of CaCl₂ at 0.25N, was simulated as well .

II. Methodology

Different molecular Dynamics (MD) simulations of lipid bilayer formed by 288 DPPC were carried out in aqueous solutions for different concentrations of CaCl₂, from 0, up to 0.50 N. Furthermore, with the aim of generalizing our results, a bilayer of 288 DPPS in presence of CaCl₂ at 0.25 N was simulated as well. Note that the concentration of CaCl₂ in terms of normality is defined as:

 

 

 

To build up the original system, a single DPPC lipid molecule, or DPPS lipid (Fig. 1), was placed with its molecular axis perpendicular to the mem-brane surface (xy plane). Next, each DPPC, or DPPS, was randomly rotated and copied 144 times on each leaflet of the bilayer. Finally, the gaps ex-isting in the computational box (above and below the phospholipid bilayer) were filled using an equi-librated box containing 216 water molecules of the extended simple point charge (SPC/E) [39] water model.

Thus, the starting point of the first system of Table 1 was formed by 288 DPPC in absence of CaCl₂. Once this first system was generated, the whole system was subjected to the steepest descent minimization process to remove any excess of strain associated with overlaps between neigh-boring atoms of the system. Thereby, the following DPPC systems in presence of CaCl₂ were gener-ated as follows: to obtain a [CaCl₂]=0.06 N, 15 water molecules were randomly substituted by 5 Ca²⁺ and 10 Cl^-. An analogous procedure was ap-plied to the rest of the systems, where 36, 69 and 138 water molecules were substituted by 12, 23 and 46 Ca²⁺ and 24, 46 and 92 Cl^-, to obtain a [CaCl₂] concentration of 0.13 N, 0.25 N and 0.50 N, respec-tively. Finally, the DPPS bilayer was generated fol-lowing the same procedure described above for the DPPC, starting from a single DPPS molecule and once the lipid bilayer in presence of water passed the minimization process, 324 water molecules were substituted by 204 Ca²⁺ and 120 Cl^- (note that 144 of the 204 calcium ions were added to balance the negative charge associated with the DPPS).

GROMACS 3.3.3 package [40,41] was used in the simulations, and the properties showed in this work were obtained using our own code. The force field proposed by Egberts et al. [2] was used for the lipids, and a time step of 2 fs was used as integra-tion time in all the simulations. A cut-off of 1.0 nm was used for calculating the Lennard-Jones in-teractions. The electrostatic interaction was evalu-ated using the particle mesh Ewald method [42,43]. The real space interaction was evaluated using a 0.9 nm cut-off, and the reciprocal space interaction using a 0.12 nm grid with a fourth-order spline in-terpolation. A semi-isotropic coupling pressure was used for the coupling pressure bath, with a refer-ence pressure of 1 atm which allowed the fluctu-ation of each axis of the computer box indepen-dently. For the DPPC bilayer, each component of the system (i.e., lipids, ions and water) was coupled to an external temperature coupling bath at 330 K, which is well above the transition temperature of 314 K [44,45]. For DPPS bilayer, each component of the system was coupled to an external temper-ature coupling bath at 350 K, which is above the transition temperature [46, 47]. All the MD sim-ulations were carried out using periodic boundary conditions. The total trajectory length of each sim-ulated system was of 80 ns of MD simulation, where the coordinates of the system were recorded every 5 ps for their appropriate analysis.

Finally, in order to study the effect of the temper-ature, only the case corresponding to 0.25 N CaCl₂ was investigated at two additional temperatures, 340 K and 350 K.

III. Results and discussion

i. Effect of the CaCl2 concentration

a. Surface área per lipid

Surface área per lipid {A) is a property of lipid bilayers which has been accurately measured from experiments [48]. The calculation of mean área per lipid can be determined from the MD simulation

{A) =

x . y

(2)

where x and y represent the box sizes in the direc-tion x and y (perpendicular to the membrane sur-face) over the simulation, and N is the number of lipids contained in one leaflet, in our case N = 144. Focusing on the study of the time evolution of the area per lipid, Figure 2 depicts the running

 

Figure 2: Running area per lipid at T = 330 K in presence of [CaCl₂] at (A) 0.06 N, (B) 0.13 N, (C) 0.25 N, (D) 0.50 N and (E) 0.25 N (in this case, T = 350 K). Solid lines represent the mean area obtained from the last 70 ns of the simulated trajectories (see text for further explanation). The type of lipid is indicated in the legends.

surface área per lipid for different concentrations of CaCl2 and type of lipid. In general, for tlie 5 bilayers formed by DPPC or DPPS, tlie área per lipid achieved a steady state after 10 ns of simula-tion, being this equilibration time almost indepen-dent of tlie concentration of CaCl2 and type of lipid which composed the membrane.

In absence of salt, an average área per lipid of (A) = 0.663 ± 0.008 nm² was calculated from the last 70 ns of the simulated trajectory, discarding the first 10 ns corresponding to the equilibration time. This valué agrees with experimental data, where valúes in a range from 0.55 to 0.72 nm have been measured [10,11,48-51]. Table 2 shows the mean surface área per lipid (again, after discarding the equilibration time of 10 ns) with their corresponding error bar.

From the simulation results, a shrinking in the surface área per lipid with the increment of the ionic strength of the solution is observed. This shrinking is expected and attributed to the com-plexation of lipid molecules by calcium, such as it has been pointed out in previous studies [28,29,52].

b. Deuterium order parameter

The deuterium order parameter, Scd, is measured from ²H-NMR experiments. This parameter provides relevant information related to the disorder of the hydrocarbon región in the interior of the lipid bilayers by measuring the orientation of the hydro-gen dipole of the methylene groups with respect to the perpendicular axis to the lipid bilayer. Due to the fact that hydrogens of the lipid methylene groups (CH2) have not been taken into account (in an explicit way) in our simulations, the order parameter —Scd on the i + 1 methylene group was defined as the normal unitary vector to the vector defined from the i to the i + 2 CH2 group and con-tained in the plañe formed by the methylene groups i, i + 1 and i + 2. Thus, the deuterium order parameter —Scd on the i —th oí the CH2 group can be estimated by Molecular Dynamics simulations as follows:

—Scd = —(3 eos (9) — 1}                (3)

where 9 is the angle formed between the unitary vector defined above and the z axis. The expres-sion in brackets (...) denotes an average over all

Type of	[CaCy N	(A) (nm )	Hydration
Lipid			Number
DPPC	0	0.663 ±0.008	1.758 ±0.009
DPPC	0.06	0.658 ±0.008	1.740 ±0.009
DPPC	0.13	0.651 ±0.007	1.719 ±0.010
DPPC	0.25	0.641 ±0.009	1.680 ±0.015
DPPC	0.50	0.628 ±0.010	1.610 ±0.015
DPPS	0.25	0.522 ±0.007	2.552 ±0.010

Table 2: Área per lipid and lipid hydration number as a function of salt concentration (see text for fur-ther explanation). Note that the salt concentration is given in normal units. Error bars were calculated for each system separately from subtrajectories of 10 ns length. Simulation temperature = 330 K.

the lipids and time. Henee, note that the —Scd can adopt any valué between -0.5 (corresponding to a parallel orientation to the lipid/water inter-face) and 1 (oriented along the normal axis to the lipid bilayer).

Figure 3 shows the running —Scd for different carbons of the DPPC and DPPS tails and salt concentrations. Only the carbons which correspond to the initial (hydrocarbons 2 and 6), the middle (hydrocarbon 10) and final (hydrocarbons 13 and 15) methylene groups of the lipid tails were depicted in this figure. Each point of the figure represents the average valúes of —Scd over 5 ns of subtrajec-tory length, and the lines represent the mean valúes calculated from the last 70 ns of the trajectories simulated. From this figure, it is observed how in all the cases, the required equilibration time is less than 10 ns of simulation time, independently of the salt concentration and the type of lipid. Finally, it is noted that Figure 3 exhibits an increase in the deuterium order parameters with the salt concentration, consistent with the shrinking of the área per lipid described above.

c. Lipid hydration

To analyze the lipid hydration, the radial distri-bution function g[r) of water around one of the oxygens of the phosphate group (atom number 10 in Fig. 1 for DPPC and DPPS) was calculated. The radial distribution function g (r) is defined as follows:

Figure 3: Running deuterium order parameter, —Scd, in presence of [CaCy at (A) 0.06 N, (B) 0.13 N, (C) 0.25 N, (D) 0.50 N and (E) 0.25 N. DPPC simulations were performed at 330 K and DPPS simulation temperature was 350 K. Solid lines represent mean valúes —Scd obtained írom the last 70 ns oí the simulated trajectories. The type oí lipid is indicated in the legends. Symbols: o hydrocarbon 2; o hydrocarbon 6; < hydrocarbon 10; + hydrocarbon 13 and x hydrocarbon 15. Note that the error bars have the same size as symbol.

g(r) =

N(r)

(4)

4"7rr² pSr

where N (r) is the number of atoms in a spherical shell at distance r and thickness Sr from a reference atom. p is the density number taken as the ratio of atoms to the volume of the total computing box.

From numerical integration of the first peak of the radial distribution function, the hydration num-bers can be estimated for different atoms of the DPPC or DPPS. Figure 4 depicts the hydration number of phosphate oxygen (atom 10 in Fig. 1 for DPPC and DPPS) in presence of CaCl2, where each point represents the average of 5 ns subtrajec-tory length. These results show how this property reached a steady state for the cases (A), (B) and (E), after 10 ns of simulation. However, for the cases (C) and (D), 5 ns of extra simulation trajec-tory were required to reach a steady state. Table 2

 

Figure 4: Hydration number of the phosphate oxy-gen (atom 10 in Fig. 1) along the simulated tra-jectories in presence of [CaCl₂] at (A) 0.06 N, (B) 0.13 N, (C) 0.25 N, (D) 0.50 N (for DPPC T = 330 K) and (E) 0.25 N. In this case, T = 350 K. Solid lines represent the mean value of the hydra-tion number calculated from the last 70 ns of the simulated trajectories. The type of lipid is indi-cated in the legends. Note that the error bars have the same size as symbol.

shows the hydration numbers for the last 70 ns of the trajectory length, corresponding to four concen-trations of CaCl₂ and both types of lipids, DPPC and DPPS. In this regard, from Fig. 4, the signif-icant lipid dehydration with the increment of the ionic strength of the solution is evident, in good accordance with previous results [52].

d. Phospholipid-calcium coordination

Some authors have reported how the lipid coordi-nation by divalent cations widely varies . Thus, on the one hand, some authors [25] have reported that this is a very slow process, which requires about 85 ns of simulation time, but, on the other hand, other authors [26] have suggested that this process results much more rapid, taking less than 1 ns. In this sense, the coordination of DPPC-Ca²⁺ was studied by monitoring the oxygen-calcium coordination of the carbonyl oxygens (atoms 16 and 35 in Fig. 1) and phosphate oxygens (atoms 9 and 10 of DPPC

 

Figure 5: Left column represents the number of Ca²⁺ coordinated to lipids in presence of [CaCl₂] at (A) 0.06 N, (B) 0.13 N, (C) 0.25 N and (D) 0.50 N, for T = 330 K. Right column shows the quantity of calcium ions coordinated to lipids expressed in percentage, along the simulated trajectory.

in Fig. 1), as a function of time. The left column in Fig. 5 represents the oxygen coordination num-ber, while the right one depicts the percentage of calcium ions involved in the coordination process with respect to the total number of calcium ions present in the aqueous solution. Figure 5 shows how the DPPC coordination by calcium is a quick process, taking less than 5 ns of simulation time to achieve a steady state. The kinetic of this pro-cess appears to be related to the ratio between cal-cium/lipid. Thus, after the first 5 ns of simulation time, the Ca-lipid coordination presents some fluc-tuation along the rest of the simulated trajectory. However, in Fig. 5 (A) and (B) (for the cases of lower concentration), it is observed how the per-centage of coordination fluctuates between a 60% and a 100%. We consider that this broad fluctu-ation is related to the limited sample size of our simulations that introduces a certain noise in our results.

ii. Effect of temperature

This section focuses on the study of the role played by temperature on the equilibration process. In this regard, only the system corresponding to a

330	0.642 ±0.009	1.680 ± 0.010
340	0.650 ±0.007	1.683 ± 0.020
350	0.666 ±0.008	1.689 ± 0.015

Table 3: Area per lipid and the lipid hydra-tion number as a function of temperature. Error bars were calculated from subtrajectories of 10 ns length. DPPC bilayer in presence of [CaCl₂] = 0.25 N.

concentration of 0.25 N in CaCl₂ was studied, for a range of temperatures from 330 K to 350 K (all of them above the transition temperature of 314 K [44, 45] for the DPPC).

Figure 6 shows the running area along the trajec-tory. In this case, it was noticed how the systems achieve a steady state after a trajectory length of roughly 10 ns, where Table 3 shows the mean area per lipid calculated from the last 70 ns of simu-lation time. Figure 7 shows the deuterium order parameter of the methylene groups along the lipid tails, calculated from Eq. (3). Figure 7, on the one hand, clearly shows that for the three temper-atures the systems have reached the steady state before the first 10 ns of simulation. On the other hand, it shows an increase in the disorder of the lipid tails with temperature, which is closely re-lated with the increase of the area per lipid, such as it was pointed out above. Figure 8 depicts the results of the hydration numbers of DPPC for the three temperatures studied, where the equilibrated state was achieved after 10 ns of simulation time. Table 3 provides the calculated hydration numbers in the equilibrium, showing how the lipid hydra-tion remained invariable with the rising of the tem-perature. Concerning the lipid-calcium coordina-tion, Fig. 9 represents the lipid-calcium coordina-tion number, and the right column represents the calcium that participates in the coordination ex-pressed in percentage respect the total of calcium ions in solution. From simulation, it becomes ev-ident how calcium ions required less than 5 ns to achieve an equilibrated state for the three temper-atures studied. In summary, for all the properties studied in this section, a slight decrease in the equi-libration time with the increasing temperature was appreciated.

 

Figure 6: Running area per lipid for [CaCl₂] = 0.25 N at different temperatures, (A) T = 330 K, (B) T = 340 K and (C) T = 350 K. Solid lines represent the mean values obtained from the last 70 ns of simulation.

 

Figure 7: Deuterium Order Parameter, -S_CD, along the simulated trajectory for a concentration of [CaCl₂] = 0.25 N, for the following temperatures: (A) T = 330 K, (B) T = 340 K and (C) T = 350 K. Solid lines represent the average order parame-ter for the last 70 ns of simulation. Note that the error bars have the same size as symbol.

IV. Conclusions

The present work deals with the simulation time re-quired to achieve the steady state for a lipid bilayer system in presence of CaCl₂. In this regard, we studied two different systems: one with DPPC and another one with DPPS bilayer; both systems in

Figure 8: Hydration number of phosphate oxygen (atom 10 in Fig. 1) along the simulated trajectories for a [CaCl₂] = 0.25 N at different temperatures: (A) T = 330 K, (B) T = 340 K and (C) T = 350 K. Solid lines represent the average hydration number for the last 70 ns of simulation. Note that the error bars have the same size as symbol.

presence of CaCl₂ (at different level concentration). The salt free case was also studied, as control. The analysis of various lipid properties studied here in-dicates that some properties reach the steady state more quickly than others. In this sense, we found that the area per lipid and the hydration number are slower than the deuterium order parameter and the coordination of cations. Consequently, to en-sure that a system composed by a lipid bilayer has reached a steady state, the criterion that we propose is to show that the area per lipid or the hydration number have reached the equilibrium.

From our results, two important aspects should be remarked:

1.  The equilibration time is strongly depen-dent on the starting conformation of the sys-tem. Wrong starting conformations will re-quire much longer equilibration times, even of one order of magnitude higher than the re-quested from a more refined starting conformation.

2.  Temperature is a critical parameter for reduc-ing the equilibration time in our simulations, due to the fact that higher temperatures in-crease the kinetic processes, i.e., the sampling of the configurational space of the system.

Acknowledgements - Authors wish to thank the assistance of the Computing Center of the Universidad Polit´ecnica de Cartagena (SAIT), Spain. RDP is member of ‘Carrera del Investigador’, CON-ICET, Argentine.

[1] O Berger, O Edholm, F Jahnig, Molecular dy-namics simulations of a fluid bilayer of di-palmitoylphosphatidylcholine at full hydration, constant pressure, and constant temperature, Biophys. J. 72, 2002 (1997).

[2] E Egberts, S J Marrink, H J C Berendsen, Molecular dynamics simulation of a phospho-lipid membrane, Eur. Biophys. J. 22, 423 (1994).

[3] U Essmann, L Perera, M L Berkowitz, The origin of the hydration interaction of lipid bi-layers from MD simulation of dipalmitoylphos-phatidylcholine membranes in gel and crys-talline phases, Langmuir 11, 4519 (1995).

[4] S E Feller, Y Zhang, R W Pastor, R B Brooks, Constant pressure molecular dynamics simula-tions: The Langevin piston method, J. Chem. Phys. 103, 4613 (1995).

[5] W Shinoda, T Fukada, S Okazaki, I Okada, Molecular dynamics simulation of the dipalmi-toylphosphatidylcholine (DPPC) lipid bilayer in the fluid phase using the Nosr-Parrinello-Rahman NPT ensemble, Chem. Phys. Lett. 232, 308 (1995).

[6] D P Tieleman, H J C Berendsen, Molecular dy-namics simulations of a fully hydrated dipalmi-toylphosphatidylcholine bilayer with different macroscopic boundary conditions and param-eters, J. Chem. Phys. 105, 4871 (1996).

[7] K Tu, D J Tobias, M L Klein, Constant pressure and temperature molecular dynamics simulation of a fully hydrated liquid crystal phase dipalmitoylphosphatidylcholine bilayer, Biophys. J. 69, 2558 (1995).

[8] M F Brown, Theory of spin-lattice relaxation in lipid bilayers and biological membranes. Dipolar relaxation, J. Chem. Phys. 80, 2808 (1984).

[9] M F Brown, Theory of spin-lattice relaxation in lipid bilayers and biological membranes. ²H and ¹⁴N quadrupolar relaxation, J. Phys. Chem. 77, 1576 (1982).

[10] J F Nagle, R Zang, S Tristam-Nagle, W S Sun, H I Petrache, R M Suter, X-ray structure de-termination of fully hydrated L. phase dipalmi-toylphosphatidylcholine bilayers, Biophys. J. 70, 1419 (1996).

[11] R P Rand, V A Parsegian, Hydration forces be-tween phospholipid bilayers, Biochim. Biophys. Acta 988, 351 (1989).

[12] J Seelig, Deuterium magnetic resonance: The-ory and application to lipid membranes, Q. Rev. Biophys. 10, 353 (1977).

[13] J Seelig, A Seelig, Lipid conformation in model membranes and biological systems, Q. Rev. Biophys. 13, 19 (1980).

[14] W J Sun, R M Suter, M A Knewtson, C R Worthington, S Tristram-Nagle, R Zhang, J F Nagle, Order and disorder in fully hydrated un-oriented bilayers of gel phase dipalmitoylphos-phatidylcholine, Phys. Rev. E. 49, 4665 (1994).

[15] H Akutsu, J Seelig, Interaction of metal ions with phosphatidylcholine bilayer membranes, Biochemistry 20, 7366 (1981).

[16] M G Ganesan, D L Schwinke, N Weiner, Effect of Ca²⁺ on thermotropic properties of satu-rated phosphatidylcholine liposomes, Biochim. Biophys. Acta 686, 245 (1982).

[17] L Herbette, C A Napolitano, R V McDaniel, Direct determination of the calcium profile structure for dipalmitoyllecithin multilayers using neutron diffraction, Biophys. J. 46, 677 (1984).

[18] D Huster, K Arnold, K Gawrisch, Strength of Ca²⁺ binding to retinal lipid membrane: Con-sequences for lipid organization, Biophys. J. 78, 3011 (2000).

[19] Y Inoko, T Yamaguchi, K Furuya, T Mitsui, Effects of cations on dipalmitoyl phosphatidyl-choline/cholesterol/water systems, Biochim. Biophys. Acta 413, 24 (1975).

[20] R Lehrmann, J J Seelig, Adsorption of Ca²⁺ and La³⁺ to bilayer membranes: Measurement of the adsorption enthalpy and binding con-stant with titration calorimetry, Biochim. Bio-phys. Acta 1189, 89 (1994).

[21] L J Lis, W T Lis, V A Parsegian, R P Rand, Adsorption of divalent cations to a variety of phosphatidylcholine bilayers, Biochemistry 20, 1771 (1981).

[22] L J Lis, V A Parsegian, R P Rand, Binding of divalent cations to dipalmitoylphosphatidyl-choline bilayers and its effect on bilayer inter-action, Biochemistry 20, 1761 (1981).

[23] T Shibata, Pulse NMR study of the interaction of calcium ion with dipalmitoylphosphatidyl-choline lamellae, Chem. Phys. Lipids. 53, 47 (1990).

[24] S A Tatulian, V I Gordeliy, A E Sokolova, A G Syrykh, A neutron diffraction study of the in-fluence of ions on phospholipid membrane in-teractions, Biochim. Biophys. Acta 1070, 143 (1991).

[25] R A B¨ockmann, H Grubmu¨ller, Multistep binding of divalent cations to phospholipid bi-layers: A molecular dynamics study, Ange-wandte Chemie 43, 1021 (2004).

[26] J Faraudo, A Travesset, Phosphatidic acid do-mains in membranes: Effect of divalent coun-terions, Biophys. J. 92, 2806 (2007).

[27] A A Gurtovenko, Asymmetry of lipid bilay-ers induced by monovalent salt: Atomistic molecular-dynamics study, J. Chem. Phys. 122, 244902 (2005).

[28] P Mukhopadhyay, L Monticelli, D P Tieleman, Molecular dynamics simulation of a palmitoyl-oleoyl phosphatidylserine bilayer with Na⁺ counterions and NaCl, Biophys. J. 86, 1601 (2004).

[29] S A Pandit, D Bostick, M L Berkowitz, Molecular dynamics simulation of a dipalmitoylphos-phatidylcholine bilayer with NaCl, Biophys. J. 84, 3743 (2003).

[30] U R Pedersen, C Laidy, P Westh, G H Pe-ters, The effect of calcium on the properties of charged phospholipid bilayers, Biochim. Bio-phys. Acta 1758, 573 (2006).

[31] J N Sachs, H Nanda, H I Petrache, T B Woolf, Changes in phosphatidylcholine headgroup tilt and water order induced by monovalent salts: Molecular dynamics simulations, Biophys. J. 86, 3772 (2004).

[32] K Shinoda, W Shinoda, M Mikami, Molecular dynamics simulation of an archeal lipid bilayer whit sodium chloride, Phys. Chem. Chem. Phys. 9, 643 (2007).

[33] N L Yamada, H Seto, T Takeda, M Nagao, Y Kawabata, K Inoue, SAXS, SANS and NSE studies on unbound state in DPPC/water/CaCl₂ system, J. Phys. Soc. Jpn. 74, 2853 (2005).

[34] D Frenkel, B Smit, Understanding molecular simulations, Academic Press, New York (2002).

[35] J J L´opez Cascales, J Garc´ıa de la Torre, S J Marrink, H J C Berendsen, Molecular dynam-ics simulation of a charged biological mem-brane, J. Chem. Phys. 104, 2713 (1996).

[36] W F van Gunsteren, H J C Berendsen, Computer simulations of molecular dynamics: Methodology, applications and perspectives in chemistry, Angew. Chem Int. Ed. Engl. 29, 992 (1990).

[37] R A B¨ockmann, A Hac, T Heimburg, H Grubmu¨ller, Effect of sodium chloride on a lipid bilayer, Biophys. J. 85, 1647 (2003).

[38] A A Gurtovenko, I Vattulainen, Effect of NaCl and KCl on phosphatidylcholine and phos-phatidylethanolamine lipid membranes: In-sight from atomic-scale simulations for un-derstanding salt-induced effects in the plasma membrane, J. Phys. Chem. B. 112, 1953 (2008).

[39] H J C Berendsen, J R Grigera, T P Straatsma, The missing term in effective pair potentials, J. Phys. Chem. 91, 6269 (1987).

[40] H J C Berendsen, D van der Spoel, R van Drunen, A message-passing parallel molecular dynamics implementation, Comp. Phys. Comm. 91, 43 (1995).

[41] E Lindahl, B Hess, D van der Spoel, GRO-MACS 3.0: A package for molecular simula-tion and trajectory analysis, J. Mol. Mod. 7, 306 (2001).

[42] T Darden, D York, L Pedersen, Particle mesh Ewald: An N.log(N) method for Ewald sums in large systems, J. Chem. Phys. 98, 10089 (1993).

[43] U Essmann, L Perea, M L Berkowitz, T Dar-den, H Lee, L G Pedersen, A smooth particle mesh Ewald method, J. Chem. Phys. 103, 8577 (1995).

[44] L R De Young, K A Dill, Solute partition-ing into lipid bilayer-membranes, Biochem-istry 27, 5281 (1988).

[45] A Seelig, J Seelig, The dynamic structure of fatty acyl chains in a phospholipid bilayer mea-sured by deuterium magnetic resonance, Bio-chemistry 13, 4839 (1974).

[46] G Cevc, A Watts, D Marsh, Titration of the phase transition of phosphatidilserine bi-layer membranes. Effect of pH, surface electro-statics, ion binding and head-group hydration, Biochemistry 20, 4955 (1981).

[47] H Hauser, F Paltauf, G G Shipley, Structure and thermotropic behavior of phosphatidylser-ine bilayer membranes, Biochemistry 21, 1061 (1982).

[48] J F Nagle, S Tristam-Nagle, Structure of lipid bilayers, Biochim. Biophy. Acta 1469, 159 (2000).

[49] B A Lewis, D M Engelman, Lipid bilayer thick-ness varies linearly with acyl chain length in fluid phosphatidylcholine vesicles, J. Mol. Biol. 166, 211 (1983).

[50] R J Pace, S I Cham, Molecular motions in lipid bilayer. I. Statistical mechanical model of acyl chains motion, J. Chem. Phys. 76, 4217 (1982).

[51] R L Thurmond, S W Dodd, M F Brown, Molecular areas of phospholipids as deter-mined by 2H NMR spectroscopy, Biphys. J. 59, 108 (1991).

[52] R D Porasso, J J L´opez Cascales, Study of the effect of Na⁺ and Ca²⁺ ion concentration on the structure of an asymmetric DPPC/DPPS + DPPS lipid bilayer by molecular dynamics simulation, Coll. and Surf. B. Bioint. 73, 42 (2009).