PMCCPMCCPMCC

Search tips
Search criteria 

Advanced

 
Logo of membranesMDPIHomeThis articleThis journalInstructions for authorsSubscribeMembranes
 
Membranes (Basel). 2017 March; 7(1): 5.
Published online 2017 January 25. doi:  10.3390/membranes7010005
PMCID: PMC5371966

Effect of Sodium and Chloride Binding on a Lecithin Bilayer. A Molecular Dynamics Study

Hsueh-Chia Chang, Academic Editor

Abstract

The effect of ion binding on the structural, mechanical, dynamic and electrostatic properties of a 1-palmitoyl-2-oleoyl-sn-glycero-3-phosphocholine (POPC) bilayer in a 0.5 M aqueous NaCl solution is investigated using classical atomistic molecular dynamics simulation with different force-field descriptions for ion-ion and ion-lipid interactions. Most importantly, the repulsive Lennard–Jones parameters for the latter were modified, such that approximately similar binding of cations and anions to the lipid membrane is achieved. This was done to qualitatively improve the apparent ion-lipid binding constants obtained from simulations with the original force field (Berger lipids and GROMOS87 ions in combination with the SPC water model) in comparison to experimental data. Furthermore, various parameters characterizing membrane structure, elasticity, order and dynamics are analyzed. It is found that ion binding as observed in simulations involving the modified in comparison to the original force-field description leads to: (i) a smaller salt-induced change in the area per lipid, which is in closer agreement with the experiment; (ii) a decrease in the area compressibility and bilayer thickness to values comparable to a bilayer in pure water; (iii) lipid deuterium order parameters and lipid diffusion coefficients on nanosecond timescales that are very similar to the values for a membrane in pure water. In general, salt effects on the structural properties of a POPC bilayer in an aqueous sodium-chloride solution appear to be reproduced reasonably well by the new force-field description. An analysis of membrane-membrane disjoining pressure suggests that the smaller salt-induced change in area per lipid induced by the new force-field description is not due to the alteration of membrane-associated net charge, but must rather be understood as a consequence of ion-specific effects on the arrangement of lipid molecules.

Keywords: molecular dynamics, POPC bilayer, salt effects, lipid force field, ion force field, sodium chloride

1. Introduction

The presence of ions crucially influences membranes in terms of structural and mechanical properties, dynamics, thermodynamic stability and biological function [1,2,3,4,5,6,7]. Understanding the molecular mechanisms underlying membrane structure and mechanics is important not only because of the biological relevance of membranes in aqueous electrolyte solutions, but also because the structural and mechanical characteristics of membranes will, in the future, be of increasing interest in the development of synthetic materials or engineered tissues [8,9].

The strength of ion-membrane binding is governed by the equilibrium between solvation and the membrane-interaction of the ions. Intuitively, electrostatic interactions are a key determinant of this equilibrium, which is why, at first glance, considerable differences in ion-membrane binding may be expected depending on the ion valency and the membrane charge. Divalent ions show stronger binding than monovalent ones, e.g., the (intrinsic) binding constants of Na+, Cl and Ca2+ to zwitterionic 1,2-dipalmitoyl-sn-glycero-3-phosphatidylcholine (DPPC) liposomes were measured at 25 °C as 0.25 M1, 0.28 M1 and 37 M1, respectively [10]. For comparison, the intrinsic binding constant at 25 °C of Na+ to anionic phosphatidylserine-containing vesicles was suggested to be 0.6 M1 [11]. This illustrates the seemingly important role of electrostatic interactions. Considering the release of tightly-bound solvation-shell water molecules that accompanies the desolvation of a membrane-associated ion, entropy gain may also be seen as a main driving force for the binding of ions to membranes [5]. The two mechanisms (the role of electrostatic interactions and entropy gain of released water) are, however, intimately connected, because the Coulomb interaction energy (between two charged species in water, expressing the solvent degrees of freedom via the dielectric constant of the solvent) can, due to the temperature dependence of the solvent relative dielectric permittivity, be interpreted as a free energy [12], with entropic and enthalpic components per unit charge of TΔS=1.3 V and ΔH=0.3 V, respectively. Moreover, one may expect that the ionic strength of the electrolyte solution or the presence of other salts can affect ion binding constants, because other salts can, for instance, decrease the hydration of an ion in the electrolyte and, thus, enhance its association with the membrane.

Physical insight into the thermodynamics of ion-membrane association can be obtained experimentally via radiotracing of radioactive ion isotopes [13], measurement of membrane electrophoretic mobilities [10,11,14], isothermal titration calorimetry [5,14] or via computer experiments, e.g., through molecular dynamics (MD) simulation [14,15]. Due to the heterogeneous nature of the membrane (biological membranes are composed of a variety of lipids and proteins) and its environment (electrolyte solution), as well as the physical complexity of the forces at play (long-range electrostatic, short-range van der Waals interactions), neither a theoretical understanding, nor a predictive theoretical description of this equilibrium are straightforward. It was pointed out before [16] that current models of specific ion effects are not sufficient to explain the interaction of chaotropic anions and kosmotropic cations with lipid membranes, i.e., the fact that anions and cations show opposite Hofmeister-series behavior, which hints at different binding mechanisms for the two types of ions. For zwitterionic lecithin (phosphatidylcholine) lipids, ion-lipid association can, to a certain extent, be explained in terms of the empirical law of matching water affinities [17], based on which Garcia-Celma et al. [16] interpreted cation binding to phosphate groups (both of kosmotropic nature) and anion binding to choline groups (both of chaotropic nature). However, this interpretation fails for other lipids, e.g., anionic surfactants [16]. Leontidis and Aroti [18] suggested that the different interaction behavior of cations and anions with the lipid membrane is drawn back to cations interacting directly with lipid headgroups, thus giving rise to density inhomogeneities through a clustering of lipids, whereas anions invade the created space of low lipid density. They also stressed the potential MD simulation has in unraveling ion-membrane interactions.

Numerous MD studies on ion-membrane interactions were performed in the past, e.g., focusing on the kinetics, stoichiometry and structural aspects of ion-lipid binding [19,20,21,22,23] and the specific ion effects thereof [24,25,26,27]. Partly, the results have been found to depend strongly on the employed force field [22], and it was emphasized that careful force-field optimization might be required to achieve simulation results in agreement with experimental findings [14]. However, such a force-field recalibration is highly intricate as its aims, namely an accurate representation of both ion-ion, the ion-membrane functional group and ion-water interactions, might be irreconcilable with the approximations made by current “simple” force-field descriptions. In other words, it might be required to circumvent combination rules for Lennard–Jones parameters or/and account (explicitly) for polarization effects. Nevertheless, MD simulation is a valuable tool to investigate ion-membrane systems. Besides providing information about the thermodynamics of ion-membrane binding, it also offers atom-level insight concerning the structure, mechanics, order and dynamics of membranes.

Phospholipid membranes can exist in various alternative phases. The biologically most relevant phase is the liquid-crystalline phase, which is fluid-like, presenting the lipid chains in rather disordered configurations. The flexible nature of the membrane renders it susceptible to several types of deformation that may occur at room temperature and atmospheric pressure. The four most common types of membrane deformation include [28] changes of membrane curvature, area, thickness and deformation due to shear forces. For a phospholipid bilayer, the most important deformation at room temperature is a change of curvature, i.e., bending. The associated bending rigidity is around 10–20 kBT, where kB is Boltzmann’s constant and T the absolute temperature. The energetic cost of area changes is higher, Hooke’s law force constant for an area change being around 55–70 kBT·nm2. Stretching or compressing the membrane along its normal has a similar Hooke’s law force constant of around 60 kBT·nm2. Finally, shearing is not of relevance for fluid-like membranes, such as a membrane in the liquid-crystal phase, but may occur when membranes are attached to solid structures.

When lipid bilayers are immersed in aqueous electrolyte solutions in comparison to pure water, MD simulations evidence a reduction of membrane fluidity (i.e., of lateral lipid diffusion), a reduction of the area per lipid and an increase in membrane bilayer thickness [14,19,23,29]. Similar findings were obtained experimentally [4], although the influence of salt concentration on some structural membrane properties appears to be more complex in the case of certain salts. For example, in the case of NaCl, the bilayer separation (thickness of the inter-membrane water layer) increases monotonously with salt concentration, whereas it decreases with increasing CaCl2 concentration [4]. Zimmermann et al. [7] analyzed the effect of different salts on membrane fluidity in terms of a relationship between membrane fluidity and membrane charge. It was observed that, although chaotropic anions and kosmotropic cations can influence fluidity in a similar fashion, they have a different impact on membrane charge. A simple explanation for the effect of ions on membrane fluidity might be that the decrease of fluidity is caused by dehydration of lipid molecules due to ions near the membrane-water interface [29].

In general, lipid order [4] and area compressibility [1] are found to increase upon the addition of salt. Concomitant with the increase of order, the elasticity of the bilayer decreases [4]. A number of experimental studies on the bending stiffness of lipid membranes were performed [30,31,32], illustrating, e.g., the influence of the electrostatic characteristics of the membrane-water interface [30,32] or membrane-active compounds [31] on membrane mechanical properties. Membrane-water interfacial electrostatic properties are commonly described by the zeta potential, and it was found [30] that the relation between zeta potential and bending stiffness is qualitatively captured by the theoretical model of Helfrich and Winterhalter [33,34]. Thus, one can decipher how ion binding affects membrane mechanics via altering surface charge. Ion binding does not only influence the properties of the membrane-water interface, but also intramembrane electrostatic properties, such as the membrane dipole potential. However, supposedly, there is no immediate connection between the membrane dipole potential and membrane mechanical properties [30].

In MD simulations, surface charge is expected to be reflected in both physical (membrane stiffness) and artificial phenomena. As MD simulations of membrane-solvent systems are commonly performed under periodic boundary conditions, i.e., the separation of periodic bilayer copies (along the bilayer normal) is finite rather than infinite, membrane interfacial electrostatic properties translate into artificial repulsive forces between bilayer copies. The magnitude of the resulting disjoining pressure [35] is proportional to the membrane area compressibility and the inverse of the solvent layer thickness. Thus, it may be difficult to separate physical and artificial effects of the membrane surface charge induced by ion-membrane association. An intuitive reasoning suggests that surface charge correlates with inter-bilayer repulsion, hence increasing inter-bilayer separation (i.e., the thickness of the solvent layer) and (due to solvent incompressibility) decreasing area per lipid.

In the present study, MD simulations are performed to investigate the above reasoning about salt effects on membranes. In particular, an existing ion-lipid force field is recalibrated to achieve similar binding affinities for Na+ and Cl ions to a 1-palmitoyl-2-oleoyl-sn-glycero-3-phosphocholine (POPC) bilayer in an aqueous approximately 0.5 M NaCl solution. It is studied how membrane properties, as well as electrostatic interfacial properties are altered. The goal is to investigate whether a correct description of ion-membrane binding behavior reproduces experimental trends concerning salt effects on membrane structure, mechanics, order and dynamics. The paper is organized as follows. Section 2 gives a brief outline of the employed biophysical theory about membrane-ion interactions, membrane-solvent interfacial electrostatics and inter-bilayer electrostatic forces. Section 4 provides computational details about the MD simulations. Corresponding results are reported and discussed in Section 3. Finally, Section 5 summarizes the main conclusions.

2. Theory

2.1. Ion-Lipid Binding Constants

Assuming a membrane in a monovalent salt and a 1:1-stoichiometric association of an ionic species I with lipid molecules, one can define an apparent binding constant Kapp(I) as:

Kapp(I)=Nb(I)(NLNb(I))c¯blk=(c¯blk)1α(I)1α(I),
(1)

where the law of mass action was used, with α(I)=Nb(I)NL denoting the fraction of the number of ions bound to lipid molecules (Nb) and the total number of lipids (NL) and c¯blk denoting the mean salt concentration in the bulk [15,36]. Equation (1) captures an association equilibrium via counting concentrations of complexes and free components, i.e., straightforward application of the law of mass action. It implicitly assumes that the ratio of free and complexed species is representative of a macroscopic system in the sense that it does not suffer from the microscopic-system artifact encountered when determining a binding constant based on association frequency [37].

The apparent binding constant characterizes the attraction of an ion to the lipid membrane due to the ion-membrane affinity intrinsic to the nature of the formed ion-lipid complex, but also due to electrostatic forces that arise solely on the basis of the membrane zeta potential, Δζ. The latter contribution can be removed through the corresponding Boltzmann factor exp[βqIΔζ], where qI is the ion charge and β=(kBT)1, leading to the intrinsic binding constant:

Kint(I)=Kapp(I)exp[βqIΔζ].
(2)

Note that for vanishing surface charge, the zeta potential evaluates to zero, and therefore, Kapp(I)=Kint(I).

2.2. Finite-Size Effects

When lipid membranes are simulated under periodic boundary conditions, structural parameters, such as the area per lipid, are influenced by the periodic membrane copies along the bilayer normal. Due to their long-range nature, this finite-size artifact is predominantly mediated by electrostatic interactions between the surface charges of the bilayers. Denoting the area per lipid pertaining to infinite bilayer separation (i.e., infinite z-dimension of the simulated system) and zero surface tension with aL,o, and that pertaining to the simulated system of finite z-dimension with aL,o, the free energy associated with the change in area per lipid from aL,o to aL,o, is:

ΔGFS=κA(aL,o,aL,o)22aL,o.
(3)

The area compressibility κA in Equation (3) is given by:

κA=β1aL,o(NL/2)δaL,o2,
(4)

where δaL,o2 is the mean square fluctuation of the area per lipid. In Equation (3), it is assumed that κA does not change upon increasing the bilayer separation to infinity. The area per lipid at the infinite bilayer separation may be expressed in terms of the disjoining pressure [38]:

Pdis=κAdw1(1aL,o1aL,o,),
(5)

where dw is the water layer thickness. Setting Pdis equal to the disjoining pressure present in an MD simulation, i.e., [38]:

PMD=κAdw1(1aL,o1aL),
(6)

allows solving Equation (5) for aL,o,. Comparison of aL,o, with the actual area per lipid observed during the simulation gives an estimate of the magnitude of finite-size effects on the simulated membrane structure.

2.3. Inter-Bilayer Forces

The disjoining pressure given by Equation (6) may also be used to establish a connection between membrane surface charge and membrane structural parameters (thickness and area). It characterizes the repulsion between periodic bilayer copies due to electrostatic forces between the associated surface charges. Comparing it with the theoretically-expected repulsion, captured here by the theoretically-expected disjoining pressure [35]:

PTH=β1(ρNa+(zm)+ρCl(zm)2(ρNa+(zm)ρCl(zm))1/2),
(7)

where ρNa+(zm) and ρCl(zm) are sodium and chloride ion number densities in the center of the water layer, gives insight into the mechanism underlying area and thickness alteration of the membrane, i.e., whether this alteration is really related to the change in surface charge that occurs when the membrane is immersed in an electrolyte solution rather than in pure water.

2.4. Membrane Rigidity and Area Compressibility

Modeling a lipid bilayer membrane as an array of thin elastic shells and describing the shell deformations with elasticity theory [39], one can express the area compressibility κA (Equation (17)) and bending stiffness kc (bending rigidity; Equation (18)) in terms of the two-dimensional Lamé constants λ and μ. If one assumes a vanishing shear modulus, one can then derive a relation between bending stiffness, area compressibility and membrane thickness, captured by a coefficient B [40]:

B=kcκA1DHH2,
(8)

where DHH is the bilayer thickness. It was seen that the magnitude of B is a measure for the extent of interdigitation present between the two bilayer leaflets [41]. A value B=1/48 is characteristic of membranes showing little or no interdigitation, i.e., leaflets that can slide past each other with little hindrance, whereas a larger value, B=1/12, indicates coupled leaflet motion due to interdigitation. Note that the quadratic dependence of kc on DHH in Equation (8) derives from a very simple membrane model. A more involved description of the membrane accounting for conformational entropy of the fatty acid chains suggests kcDHH5/2 [42].

3. Results and Discussion

Based on the simulations described in Section 4.1, the original and modified force-field descriptions for ion-lipid interactions were compared in terms of their influence on structural properties (membrane area per lipid, area compressibility, bilayer thickness, water layer thickness, lipid order parameters), the thermodynamics of ion-membrane binding (apparent binding constants, coordination numbers), lipid diffusion coefficients, electrostatic potential variation along the bilayer normal and membrane rigidity.

In simulations with the modified force-field version, the area per lipid increases in comparison to the original force field by about 16 (SNaCl(m1,wc)) or 14% (SNaCl(m2,wc); Table 1). It slightly exceeds the area per lipid of the membrane patch in pure water (0.607 nm2), namely by about 0.02 and 0.01 nm2 in simulations SNaCl(m1,wc) and SNaCl(m2,wc), respectively (Table 1). The results obtained with the modified force-field version entailing enhanced ion binding are thus in closer agreement to experimental data [4]. The area per lipid in simulations with the Parrinello–Rahman barostat is overall similar to that observed in simulations with the weak-coupling barostat (differences of 0–0.009 nm2).

Table 1
Area per lipid aL (Equation (9)), area compressibility κA (Equation (17)), bilayer thickness DHH (Section 4.3.1) and water layer thickness dw (Equation (10)) for systems containing membrane patches of 64 lipids per leaflet. The simulation acronyms ...

As a consequence of lateral membrane expansion (along with solvent incompressibility), the thickness of the solvent layer dw decreases in simulations SNaCl(m1,wc) and SNaCl(m2,wc) compared to the simulation with the original force field (Table 1). In addition, the bilayer thickness DHH decreases and basically adopts values observed for a membrane patch in pure water. For example, in simulations SNaCl(m1,wc) and SNaCl(m2,wc), DHH evaluates to 3.86 and 3.95 nm, respectively, i.e., it is very similar to the thickness found in simulation Swat(wc), 3.98 nm, and smaller than that found in simulation SNaCl(o,wc), 4.31 nm. Experimentally, the presence of aqueous NaCl leads to an increase in DHH above NaCl concentrations of 0.6 M [4]. This increase in bilayer thickness is in qualitative agreement with simulation SNaCl(o,wc), but not with the simulations involving the new force-field versions. In comparison to simulation SNaCl(o,wc), the binding of approximately equal amounts of cations and anions (SNaCl(m1,wc), SNaCl(m2,wc)) induces a contraction of the membrane along the bilayer normal and a lateral expansion. Intuitively, one can try to explain this in terms of repulsion between charged bilayer copies as present in simulation SNaCl(o,wc). Here, predominantly sodium ions bind below the membrane shear plane, whereas chloride ions are mainly located above the shear plane (Figure 1). Therefore, the membrane presents a net positive charge. An intuitive reasoning would suggest that the resulting electrostatic repulsion between periodic bilayer copies increases solvent layer thickness. Furthermore, due to solvent incompressibility, the xy-plane of the simulation box has to contract, which would be reflected in a decreased area per lipid (Table 1). However, an alternative mechanism is not based on electrostatic grounds, but ion-specific effects arising from direct ion-lipid interactions, i.e., the way in which ions bind to the lipid molecules may directly cause lateral compression or extension of the membrane. This is discussed in more detail below.

Figure 1
Profiles of water (ρOW(z)) and ion (ρI(z); I= Na+ or Cl) number densities normal to the membrane, along with integrated ion density profiles giving the cumulative number of ions (NI(z); I= Na+ or Cl), for simulation S ...

An increase of sodium-lipid and a decrease of chloride-lipid nonpolar repulsion, as achieved here by a scaling of corresponding Lennard–Jones C12 parameters (Section 4.2), brings both cation and anion binding in better agreement with the experimental findings. Apparent ion-membrane binding constants from electrophoresis and isothermal titration measurements [14] are 0.437±0.046 M1 for Na+ ions and 0.401±0.038 M1 for Cl ions. They may be compared to the corresponding values from the simulations performed here, which are reported in Table 2 and were calculated according to Equation (21) based on the simulated ion density profiles along the bilayer normal, shown in Figure 2 and Figure 3 for the modified force-field versions m1 and m2, respectively. Sodium ion binding is overestimated by around a factor two in the original force field, i.e., Kapp=0.87±0.16 M1, using the arithmetic mean of the binding constants pertaining to the two leaflets (0.71 and 1.02 M1; Table 2), with the indicated error denoting the standard deviation. However, it is essentially in line with experimental data when the modified force fields m1 or m2 are used. The binding constants are Kapp=0.35±0.01 and 0.39±0.11 M1, respectively, again resulting from averaging over the two leaflets and using the standard deviation as an error estimate (Table 2). Chloride ion binding is significantly underestimated by the original force field (Kapp=0.011±0.002 M1). By virtue of decreased chloride-ion lipid-tail nonpolar repulsion, the modified force-field versions allow those ions to cross the membrane shear plane (Figure 2 and Figure 3), thus leading to an increase in their apparent binding constants (Kapp=0.30±0.01 and 0.30±0.09 M1 in simulations SNaCl(m1,wc) and SNaCl(m2,wc), respectively; Table 2). Although chloride ion binding is still weaker than sodium binding, the structural features of the membrane basically correspond to those of a neutral membrane patch (see above and Table 1).

Figure 2
As Figure 1, but for simulation SNaCl(m1,wc), i.e., a simulation based on a modified force field (Table 7). The corresponding apparent binding constants are reported in Table 2.
Figure 3
As Figure 1, but for simulation SNaCl(m2,wc), i.e., a simulation based on a modified force field (Table 7). The corresponding apparent binding constants are reported in Table 2.
Table 2
Apparent binding constants Kapp (Equation (21)) for sodium and chloride ions, computed from ions binding to either leaflet (which is why two values are reported; Section 4.3.2), the corresponding numbers of bound ions N˜b(Na+) and N˜b ...

Analysis of coordination numbers of the ions reveals that the new force-field descriptions enhance the first-shell coordination of chloride ions with lipid tail atoms by about three orders of magnitude, while leaving first-shell coordination of sodium ions with lipid tail atoms essentially unaffected (Table 3). The first-shell coordination number of sodium ions with carbonyl- and phosphate-oxygen atoms is reduced by approximately two orders of magnitude in the new force-field description (Table 3). In particular, the binding to the carbonyl oxygen atoms is reduced in that it occurs almost exclusively in a solvent-separated fashion, as evidenced by almost negligible first-coordination shell peaks (heights of 0.21 and 0.43 for simulations SNaCl(m1,wc) and SNaCl(m2,wc), respectively) in comparison to the second-coordination shell peaks (heights of 3.07 and 3.90 for simulations SNaCl(m1,wc) and SNaCl(m2,wc), respectively). Corresponding peak heights for simulation SNaCl(o,wc) are 143.1 and 1.56, respectively, i.e., here, an enhanced direct binding in comparison to solvent-separated binding of sodium ions to the carbonyl oxygen atoms is observed. The force-field dependence of the relative balance of direct- and solvent-separated binding mechanisms of cations to lipid molecules has been noted before [43,44] and is of importance for the parameterization of ions in coarse-grained force-field descriptions (inclusion vs. omission of the first hydration shell in the beads representing ions). Although ion-water Lennard–Jones interaction parameters were not altered, the first-shell ion-water coordination numbers are slightly lower (chloride) or higher (sodium) when the new force-field description is employed (Table 3). Concerning second-shell coordination numbers, similar qualitative trends as for the first shell are observed, except for sodium-ion binding to the carbonyl oxygen atoms, where the new force-field description leads to an enhancement of the second-shell coordination number by about a factor two (Table 3).

Table 3
Average number of atoms j found in the first (Nij(1); Equation (13)) and second (Nij(2); Equation (14)) coordination shells of atom i, reported for simulations SNaCl(o,wc), SNaCl(m1,wc) and SNaCl(m2,wc). Atoms i are Na+ or Cl ions, and atoms ...

The alteration of ion-membrane binding is reflected in the electrostatic potential along the bilayer normal. In the absence of ions, the lipid headgroups (with the negatively-charged phosphate groups being closer to the membrane interior than the positively-charged choline groups) lead to a contribution of the electrostatic potential inside the membrane, which is negative with respect to the value in the bulk water (Figure 4c). Water not only dielectrically screens this electrostatic potential contribution, but overcompensates it, such that the electrostatic potential contribution due to water, as well as the total electrostatic potential are positive inside the membrane (Figure 4a,c). It was suggested before that this effect arises from the structure of water at the surface of hydrophobic media [45]; however, this explanation might be misleading because of the limitations of Equation (16) and the fact that preferential water orientation at hydrophobic surfaces is not per se decisive, but must be interpreted relative to a system with isotropically-oriented water molecules [46]. For the old force field, the asymmetric binding of the cations and the anions leads to a contribution of the ions to the electrostatic potential, which is positive inside the membrane (Figure 4b). Both water and lipids dielectrically screen (partially compensate) the effect of the ions on the electrostatic potential, which means that the contribution of both water and lipids to the electrostatic potential inside the membrane (with respect to the water) is decreased compared to the ion-free case. For the new force-field versions showing similar binding for both ion species, the effect of the ions on the electrostatic potential inside the membrane is much smaller. Hence, the change of the total electrostatic potential compared to the ion-free case, as well as the corresponding contributions from water and lipids are much smaller to negligible for the new force-field versions.

Figure 4
Symmetrized electrostatic potential ϕ(z) along the bilayer normal (Equation (16)), for simulations Swat(wc), SNaCl(o,wc), SNaCl(m1,wc) and SNaCl(m2,wc), shown as black, red, green and blue curves, respectively. The underlying atom-based charge ...

Membrane-associated net charge, as (considering apparent ion binding constants; Table 2) strongly present in simulation SNaCl(o,wc) and nearly vanishing in simulations SNaCl(m1,wc) and SNaCl(m2,wc), can cause a separation of membrane bilayer copies. This was already discussed above in terms of structural properties, e.g., an increase in the thickness of the water layer dw. However, changes in membrane area or water layer thickness can occur as well via an entirely different mechanism, e.g., because ions interacting with the membrane can push apart or pull closer neighboring lipid molecules. The question about what is the physical mechanism underlying the area-per-lipid increase in the present study can be answered through analyzing the disjoining pressure PTH, which arises from differential ion concentrations in the center of the water layer (Equation (7)). Here, it is found that PTH is indeed almost increased three-fold with the original force field in comparison to the modified force-field versions, the latter presenting more balanced association of cations and anions with the membrane (Table 4). The disjoining pressure was also calculated as a function of area compressibility, water layer thickness and area-per-lipid change upon ion binding (PMD; Equation (6)). The resulting pressures PMD obtained from the different simulations suggest that the forces influencing the distance between successive bilayer copies are indeed repulsive when the original force field is employed (positive disjoining pressure; Table 4), but attractive when the modified force-field versions are used (negative disjoining pressure; Table 4). These conclusions are independent of whether the water bilayer thickness and area-per-lipid change upon ion binding are taken from simulations with the weak-coupling or the Parrinello–Rahman barostat (Table 4).

Table 4
Disjoining pressure between successive bilayer copies, evaluated via Equation (6) (PMD) and Equation (7) (PTH). The simulation acronyms are explained in Table 7.

Note that PMD and PTH differ by about three orders of magnitude (Table 4). This indicates that the change in area per lipid induced by the new force-field description is actually not due to the alteration of membrane-associated net charge, but must rather be understood as a consequence of ion-specific effects on the arrangement of lipid molecules. Sodium ions can be envisioned as decreasing the area per lipid, because they may be coordinated by oxygen atoms pertaining to multiple lipid molecules at the same time, which decreases the distance between lipid molecules. Because of their more smeared out charge, chloride ions have a somewhat higher affinity to nonpolar membrane regions than sodium ions. The association of chloride ions with more hydrophobic membrane regions leads to a larger exposure of the fatty-acid chains to water, which causes the area per lipid to increase.

Membrane bending rigidity was calculated according to Equation (18). The results are reported in Table 5. The bending stiffness observed in the simulation involving pure water as the solvent (2.8× 1020 J) is similar to those found in comparable simulation studies, e.g., 4× 1020 J for the DPPC bilayer studied in [47]. According to experimental investigations, the bending fluctuations of a POPC bilayer increase in the presence of aqueous NaCl (above concentrations of about 0.5 M) [4], and concomitantly, the bending stiffness kc is expected to decrease. Here, the presence of salt in the solvent leads to contradictory results for the change in bending rigidity in comparison to pure water when old or new force-field descriptions are used, namely values of kc=1.7, 3.3 or 1.7 ×1020 J in simulations LNaCl(o,wc), LNaCl(m1,wc) and LNaCl(m2,wc), respectively (Table 5). Because the error of kc for the salt-free simulation is rather large, it is difficult to assess the effect of membrane-ion binding on kc, as well as whether the effect agrees with the experiment. Comparing kc in the presence of salt for the different force fields, the large difference between the similar force-field versions m1 and m2 is striking. This might be due to ill convergence of the simulations or a poor fit to the data points for the spectral intensity of undulation because of an insufficient number of data points in the considered low-wavenumber regime, i.e., a too large grid cell edge length (Equation (18) and Figure 5). Note that the grid employed for the analysis of membrane undulations was rather coarse (Section 4.3.1).

Figure 5
Spectral intensity of membrane undulations u2(q) as a function of the wavenumber q for simulations Lwat(wc), LNaCl(o,wc), LNaCl(m1,wc) and LNaCl(m2,wc), shown as black, red, green and blue circles, respectively. u2(q) ...
Table 5
Bending stiffness kc (Equation (18)) for simulations of large membrane patches in pure water or a NaCl solution described with the original or modified force-field versions. The average box-edge length along the (quadratic) membrane plane Lx ...

The coefficient B (Equation (8)) characterizing leaflet interdigitation is considerably smaller than the value B=1/48 commonly associated with the absence of interdigitation (Table 5), indicating that the two leaflets move in an uncoupled fashion. The underestimation of B in comparison to the theoretically-expected value of 1/48 may be due to an underestimation of kc or/and an overestimation of κA or/and of DHH in the present simulations. For the membrane in pure water, the latter quantities may be compared with readily-available experimental data. Experimentally, kc, κA and DHH are found to be 5×1020 J [48], 139 kJ·mol1·nm2 [49] and 3.83 nm [49], respectively, for a DPPC bilayer at ambient conditions. These values result in B=0.015 (Equation (8)), which is close to B=1/480.02. In comparison to the above experimental data, the simulated membrane shows a bending stiffness that is underestimated by 44% and an area compressibility that is overestimated by 45%, while the bilayer thickness is in fairly good agreement with the experimental value (Table 1 and Table 5). The membrane force field appears to render the membrane too soft (increased susceptibility to bending deformation and lateral compression), which results in the very small B-coefficient observed here. Considering the relatively large error in kc (Table 5), the calculated values of B are also affected by a rather large statistical uncertainty.

Comparing the effect of ions on lipid ordering and diffusion between the old and new force-field descriptions, the membrane properties appear to be closer to the pure-water situation when the latter is used. Deuterium order parameters of lipid tail CHn groups experience a minute decrease through sodium and chloride binding in simulations SNaCl(m1,wc) and SNaCl(m2,wc) in comparison to simulation Swat(wc), whereas they are elevated on average by 30.5% for the oleoyl chains and by 30.3% for the palmitoyl chains in simulation SNaCl(o,wc), i.e., lipid tail order is much enhanced through the unbalanced binding of sodium and chloride ions as occurring with the old force-field description (Figure 6). Similarly, the lipid center of mass diffusion coefficients (across all three examined time scales) are closer to the pure-water situation with the new force-field description in comparison to the old one (Table 6). With the latter, the long, medium and short time scale diffusion coefficients are reduced by 46.2, 31.7 and 13.0%, respectively, in comparison to the the diffusion coefficients found for the membrane in pure water. Thus, the excess of sodium-ion binding in simulation SNaCl(o,wc) appears to lead to a drastic slowing down of the long time scale lipid motion (Table 6).

Figure 6
Deuterium order parameter Si (Equation (11)) of the CHn groups i of the oleoyl (a) and palmitoyl (b) chains for simulations Swat(wc), SNaCl(o,wc), SNaCl(m1,wc) and SNaCl(m2,wc), shown as black, red, green and blue circles, respectively. The first ...
Table 6
Long- (2–6 ns), medium- (100–500 ps) and short-timescale (2–10 ps) lipid center of mass diffusion coefficients (Equation (15)) for simulations Swat(wc), SNaCl(o,wc), SNaCl(m1,wc) and SNaCl(m2,wc). The time intervals refer to the ...

Neutron scattering experiments deliver diffusion coefficients on a short timescale of 1–10 ps [50,51]; however, the influence of salts is unknown to date.

Experimental data based on fluorescence correlation spectroscopy (FCS) show that on millisecond timescales, lipid diffusion decreases with increasing NaCl concentration [19]. Such timescales are not routinely accessed by MD simulations. Simulations by Böckmann et al. [19] showed a decrease in lipid diffusion upon the addition of NaCl on nanosecond timescales. When inferring diffusion constants both from the simulations, as well as the experiments, analyses were performed based on a normal diffusion model. However, in fact, lipids undergo anomalous diffusion [52]. Therefore, diffusion constants inferred assuming normal diffusion are effective diffusion constants, which depend on the timescale considered. Table 6 shows that the effective diffusion constants decrease with increasing timescale. The smallest timescales considered in the table (2–10 ps) are similar to the timescales accessed by neutron scattering experiments (1–10 ps). Interestingly, the simulated diffusion coefficients on the 2–10-ps timescale are more similar to experimental values obtained at elevated temperatures ([50] reports a translational diffusion coefficient of 14×108 cm2·s1 for a system at 30 C and 400 × 108 cm2·s1 for a system at 63 C; [51] reports in-plane diffusion coefficients on the order of 15 × 108–600 × 108 cm2·s1 for a system at 60 C). As diffusion constants inferred from experiments and simulations typically depend on the timescale considered, also the effect of salt on the diffusion constants may depend on the timescales. Though the FCS experiments show a decrease in lipid diffusion on millisecond timescales [19], to our knowledge, no experimental data on the influence of ions on the effective diffusion constants on (sub)nanosecond timescales are available. On such timescales, the old force field also employed in [19] shows a decrease in lipid diffusion, whereas the force fields with the more realistic ion-membrane binding behavior introduced here do not show a substantial salt effect on lipid diffusion.

4. Computational Details

4.1. MD Simulations

All simulations were performed with the gromacs 4.6.5 engine [53,54]. Five different systems were simulated: (i) a 1 M aqueous NaCl solution consisting of 121 Na+ ions, 121 Cl ions and 6715 water molecules (NACL); (ii) a (small) hydrated bilayer consisting of 8×8 POPC molecules per leaflet and 5120 water molecules (Swat); (iii) a (small) bilayer consisting of 8×8 POPC molecules per leaflet solvated in a NaCl solution of 50 Na+ ions, 50 Cl ions and 5020 water molecules (SNaCl); (iv) a (large) hydrated bilayer consisting of 16×16 POPC molecules per leaflet and 20480 water molecules (Lwat); (v) a (large) bilayer consisting of 16×16 POPC molecules per leaflet solvated in a NaCl solution of 200 Na+ ions, 200 Cl ions and 20,080 water molecules (LNaCl). The concentration of the NaCl solution in the systems containing the POPC bilayer is about 0.5 M. Those systems were simulated using different force fields (for ion-ion and ion-lipid interactions), barostatting schemes and simulation times. An overview is provided in Table 7.

Table 7
Overview of the performed simulations. The first column provides an acronym for the different simulations, consisting of a label for the simulated system (NACL, Swat, SNaCl, Lwat or LNaCl; Section 4.1) and an indication of the employed force field (o, ...

The simulations were performed under periodic boundary conditions at constant temperature (300 K) and pressure (1 bar). If not specified differently, thermostatting and barostatting were performed using the weak-coupling scheme [55] with coupling times of 0.1 ps and 1 ps, respectively, as well as an isothermal compressibility appropriate for water (4.5 × 105 bar1) for the barostat. The membrane and solvent (water and, if present, ions) components were separately coupled to the temperature bath. The barostatting was applied semi-isotropically, except for the simulations of the aqueous NaCl solution, where the simulation box was scaled isotropically. The Linear Constraint Solver (LINCS) algorithm [57] was used to constrain all bond lengths to their reference values. The equations of motion were integrated according to the leap-frog algorithm [58], with a step size of 0.002 ps.

A buffered (Verlet) pair list was created based on a maximum allowed error for pair interactions per particle of 0.005 kJ·mol1·ps1. It was updated every five time steps. Van der Waals interactions were modeled with the Lennard–Jones potential, truncated at a distance of 1.2 nm, with the potential being shifted such that the energy vanishes at the truncation distance. Electrostatic interactions were computed using the particle-mesh-Ewald algorithm [59], based on a real-space cutoff of 1.2 nm, a grid spacing of 0.12 nm and a fourth order interpolation scheme.

Except for simulations involving the Parrinello–Rahman barostat [56] and for simulations of the NaCl solution, all simulations started from an energy-minimized conformation of a POPC bilayer immersed in solvent and were run for 210 ns, of which the first 90 ns were discarded for equilibration (Table 7). The simulation of the NaCl solution was equilibrated for 1 ns and was simulated for 9 ns for production. Simulations with the Parrinello–Rahman barostat [56] were only simulated for 120 ns, with the initial structure taken as the configuration obtained after 120 ns from the corresponding simulation with the weak-coupling barostat [55]. From these 120 ns, the first 30 ns were discarded for equilibration (Table 7).

For all simulations, except Lwat(wc), LNaCl(o,wc), LNaCl(m1,wc) and LNaCl(m2,wc), coordinates and energies were written to file every 2 ps. For the latter systems, coordinates and energies were written to file every 8 and 2 ps, respectively.

4.2. Force Field

Three different force-field parameterizations were employed that differ in the Lennard–Jones interactions between sodium and chloride ions, as well as between those ions and POPC molecules. First, an “unmodified” force-field parameter set that was also used in previous studies [5,14,15,19,36] was employed (simulations Swat(wc), Swat(pr), Lwat(wc), SNaCl(o,wc), SNaCl(o,pr) and LNaCl(o,wc); Table 7). This parameter set involves Berger lipid parameters [60], Groningen Molecular Simulation 87 (GROMOS87) parameters for sodium and chloride ions, which are equivalent to these parameters in the GROMOS 43a1 force field [61], and the Simple Point Charge (SPC) water model [62], along with application of the geometric-mean combination rule for the Lennard–Jones C6 and C12 parameters.

Two modified versions of this force field (indicated with “m1” and “m2” in Table 7) were created to investigate the effect of altered ion-lipid binding on membrane structural, mechanical, dynamic and electrostatic properties. The common features of these two parameter sets are: (i) an increase of the Lennard–Jones C12 coefficient for Na+ and Cl interactions by a factor sNaCl=1.33 to reproduce the same number of chloride ions (within the distance given by the first minimum in the sodium-chloride radial distribution function) around a sodium ion as observed by Hess and van der Vegt [63] for a 1 M aqueous NaCl solution involving the Kirkwood–Buff force field of Weerasinghe and Smith [64] (indicated as nCIP=0.1 in [63]); (ii) a decrease of the Lennard–Jones C12 coefficient for Cl interactions with the lipid tail atoms (atom types LP2, LP3 and LH1 in the Berger force field) by a factor sCl=0.057; (iii) an increase of the Lennard–Jones C12 coefficient for Na+ interactions with all lipid atoms by a factor sNa+. Versions m1 and m2 differ in the choice for sNa+, which was set to 3.2 or 3.0, respectively. Thus, the force-field modifications imply a circumvention of combination rules for the above-mentioned atom pairs. Note that these force-field modifications were not done with the aim to quantitatively reproduce the experimental apparent ion-membrane binding constants, which were determined as [14] 0.437±0.046 M1 for Na+ ions and 0.401±0.038 M1 for Cl ions for POPC vesicles in a 100 mM aqueous NaCl solution. Rather, it was attempted to achieve approximately similar Na+ and Cl binding and to investigate the effect on membrane properties.

4.3. Analysis of Simulations

4.3.1. Structural, Mechanical, Dynamic and Electrostatic Membrane Properties

Area per lipid, bilayer thickness, area compressibility and membrane bending rigidity were calculated. Note, in the following, that the POPC membrane lies parallel to the xy-plane of the coordinate system. Prior to any analysis, the membrane was centered around the origin of the computational box. The area per lipid aL is given as:

aL=LxLyNL/2,
(9)

where Lx and Ly are the box-edge lengths along the x and y dimensions and NL is the number of lipid molecules in the total bilayer (Table 7). The angular brackets denote arithmetic averaging over simulation frames. The bilayer thickness DHH is obtained by measuring the distance between maxima in the density profile of all phosphorus atoms along the z-axis. In combination with the average box edge length along the z-dimension, Lz, the water layer thickness dw is obtained as:

dw=LzDHH
(10)

and used in the calculation of the disjoining pressure PMD (Equation (6)). Here, the quantity aL,o was set to the area per lipid observed in a simulation without ions (Swat(wc)). The calculation of the disjoining pressure according to Equation (7) involves concentrations ρNa+(zm) and ρCl(zm) in the center of the solution moiety of the system, zm. These concentrations were read from the corresponding ion density profiles along the z-axis of the computational box. The density profiles were constructed using 500 bins along the z-axis, which amounts to bin widths of approximately 0.016–0.018 nm.

Deuterium order parameters Si for CHn groups i in the fatty acid tails were calculated. They are defined as:

Si=123cos2θi1,
(11)

where θi is the angle between the bilayer normal and the C–D bond vector associated with the carbon atom in CHn group i. The angular brackets imply averaging over both lipid molecules and trajectory frames. As the Berger lipid force field [60] uses a united-atom description of the fatty acid tails, the positions of the deuterium atoms were reconstructed based on the appropriate bond hybridization geometries.

Radial distribution functions gij(r) between ions i and other atoms j in the system were calculated as:

gij(r)=Nj(r)4πr2Δrρj,
(12)

where Nj(r) is the number of particles j in a shell of width Δr located at distance r from particle i and ρj is the number density of particles j. Here, sodium or chloride ions were considered as particles i, and water oxygen atoms, lipid tail carbon atoms, lipid ester carbonyl oxygen atoms, or lipid ester phosphate oxygen atoms were considered as particles j. The number of particles j within the first coordination shell of particles i is:

Nij(1)=4πρj0R1drr2gij(r),
(13)

where R1 indicates the position of the first minimum of gij(r). The number of particles j within the second coordination shell of particles i is:

Nij(2)=Nij(1)+4πρj0R2drr2gij(r),
(14)

where R2 indicates the position of the second minimum of gij(r). The bin width Δr was set to 0.002 nm.

Diffusion coefficients D for the center of mass of lipid molecules were calculated based on the Einstein relation:

D=12nddt[r(t)r(to)]2,
(15)

where n=2 for two-dimensional diffusion of the lipid molecules in the bilayer plane and [r(t)r(to)]2 is the mean square displacement of the center of mass of lipid molecules at time tto. Here, the angular brackets imply averaging over lipid molecules and time origins to. In practice, diffusion coefficients were calculated separately for the two monolayers and successively averaged. An error estimate was obtained from the propagation of uncertainty in ddt[r(t)r(to)]2 for each monolayer. Furthermore, because lipid diffusion in bilayers is anomalous [52], three time scales (2–6 ns, 100–500 ps and 2–10 ps) in which the lipid molecules show different diffusion behavior (i.e., a plot of mean square displacement vs. time has different slopes) were analyzed.

The (lattice-sum) electrostatic potential along the bilayer normal, ϕ(z), was calculated by double integration of the atom-based charge density ρa(z) (see, e.g., [45]):

ϕ(z)=Cϵo10Lz/2dzzLz/2dzρa(z),
(16)

where ϵo is the permittivity of vacuum and C an arbitrary integration constant. In practice, the integration was done numerically based on a histogram of the atom-based charge density, with a bin width of about 0.093 nm. Different subsets of the system were considered as contributing to ρa(z). The integration constant C was set such that ϕ(z) is zero in the bulk when ρa(z) involved all atoms, only ions or only lipid atoms and is zero at the bilayer center when ρa(z) involved only water atoms. The electrostatic potential was symmetrized with respect to the bilayer center. Note that the summation over atoms in Equation (16) may be considered inappropriate [65,66], but it was used here in the absence of an unambiguous definition for “molecule-based” summation involving lipid molecules.

Membrane area compressibility κA is calculated according to Equation (4), which, using aL in Equation (9), may be rewritten:

κA=β1aL(NL/2)δaL2.
(17)

Note that this compressibility is affected by membrane undulations. Hence, it is an effective compressibility that has contributions from the true (intrinsic) area compressibility describing a change in aL along the membrane surface that may actually be curved due to undulations and a second contribution, because here, an area-projection into the xy-plane is analyzed, considering that aL is given by Equation (9) [67]. Because the fluctuations of the area per lipid δaL2 are spurious in simulations involving the weak-coupling barostat, the Parrinello–Rahman barostat was used for simulations from which κA was determined (simulations Swat(pr), SNaCl(o,pr), SNaCl(m1,pr) and SNaCl(m2,pr); Table 7).

Simulations of a membrane patch four-times the size of the original membrane (simulations Lwat(wc), LNaCl(o,wc), LNaCl(m1,wc) and LNaCl(m2,wc); Table 7) were performed to calculate the bending stiffness. Specifically, the bending stiffness kc was determined from a fit of [47]:

u2(q)=β1Lx2kcq4+γ˜q21
(18)

to the spectral intensity of membrane undulations u2(q) as a function of the wavenumber q, where the microscopic surface tension coefficient γ˜ also results from the fit. The fit was carried out in the low-wavenumber regime [47], here up to q=2.0 nm1. The correlation coefficients of the fits were 0.99, 0.91, 0.93 and 0.97 for simulations Lwat(wc), LNaCl(o,wc), LNaCl(m1,wc) and LNaCl(m2,wc), respectively. The undulation intensity u2(q) was obtained as described by Lindahl and Edholm [47]. In detail, here, the carbon atoms connecting the phosphoglycerol headgroup and the fatty-acid tails were used to analyze membrane undulation. A shifted set of z-coordinates ({z˜i}) of those atoms was created by subtracting the arithmetic mean of the z-coordinates of the subset of these atoms located in the same bilayer, i.e., z˜i=ziNa1j=1Nazj, where zi is the z-coordinate of a considered carbon atom i and the sum runs over all Na of the considered carbon atoms j that are in the same bilayer leaflet as atom i (including atom i). The considered carbon atoms were assigned to two-dimensional histograms according to their x- and y-coordinates, i.e., the histogram was constructed parallel to the xy-plane. Twelve grid cells per dimension were chosen, corresponding to grid cell edges of about 1.0–1.1 nm, i.e., grid cell areas of 1.0–1.2 nm2 which, considering the observed area per lipids, cover 1–2 lipid molecules. Per frame, for each grid cell xi,yj, an average z˜¯(xi,yj) of z˜i over all atoms located in the respective grid cell was calculated. If no atoms were situated in a given cell in a given simulation frame, a value z˜¯(xi,yj) for the respective cell was determined as the average of z˜¯(xk,yl) values pertaining to the eight neighboring grid cells (under the consideration of periodic boundary conditions along the x- and y-dimensions of the computational box). The resulting discrete undulation function z˜¯(xi,yj) was transformed to reciprocal space by means of the Fastest Fourier Transform in the West 3 (FFTW3) library [68], and the resulting average spectral intensity u2(q) per wavenumber was obtained from averaging over all simulation frames.

Whenever time series of the quantity of interest were available, statistical errors were assessed via block averaging [69]. Otherwise, statistical errors were estimated as the standard error pertaining to two samples from dividing a total simulation into two halves (i.e., the standard deviation divided by the square root of two), or from the propagation of uncertainties.

4.3.2. Thermodynamics of Ion-Membrane Binding

Apparent binding constants of sodium and chloride ions to the POPC bilayer were calculated according to Equation (1). This requires knowledge of the number of membrane-bound ion species and of the mean salt concentration in the bulk. The former was obtained as:

N˜b(Na+)=0zsdzρNa+(z)aLNL/2
(19)

for sodium ions and, equivalently, as:

N˜b(Cl)=0zsdzρCl(z)aLNL/2
(20)

for chloride ions, where ρNa+(z) and ρCl(z) are the number densities of sodium and chloride ions, respectively, along the z-axis of the system and zs is the location of the shear plane, defined as that distance along the z-axis where the water oxygen number density is equal to half of its bulk value. This choice is suggested from the profile of the average velocity of matter parallel to a charged bilayer exposed to an external electric field parallel to the bilayer [45]. N˜b(Na+) and N˜b(Cl) differ from the corresponding quantities in Equation (1) in that they refer to only one of the two bilayer leaflets. Therefore, the integration in Equations (19) and (20) can in principle be performed for the two leaflets of the bilayer, giving rise to two values for each of N˜b(Na+) and N˜b(Cl), which, due to symmetry, are expected to be equal within the statistical error. Note that if N˜b(Na+) and N˜b(Cl) are used in Equation (1), α(I) has to be calculated as α˜(I)=N˜b(I)NL/2, and:

Kapp(I)=(c¯blk)1α˜(I)1α˜(I).
(21)

The mean molar salt concentration in the bulk c¯blk corresponds to the geometric mean of ρNa+(zm) and ρCl(zm), where zm is the center of the water layer along the z-axis, in units of moles per liter, i.e., assuming ρNa+(z) and ρCl(z) to be given in units of nm3, c¯blk reads:

c¯blk=1024NA1ρNa+(zm)ρCl(zm)1/2,
(22)

where NA is Avogadro’s constant.

5. Conclusions

The strength of ion-membrane binding is governed by the equilibrium between solvation and the membrane-interaction of the ions and was not accurately reproduced by the lipid force field from Berger et al. [60] in combination with the SPC water model [62] and GROMOS87 Lennard–Jones parameters for sodium and chloride ions [61]. This might be due to force field shortcomings associated with: (i) the use of ad hoc combination rules for the calculation of heteroatomic Lennard–Jones interaction parameters based on homoatomic ones, which was seen to have a spurious impact on ion-ion interactions; (ii) the use of an ad hoc 12th power of interatomic separation to describe Pauli repulsion in the Lennard–Jones potential, which might be inappropriate for very small ions exerting a large electric field on the surroundings (see also [70]); (iii) the lack of explicit polarization (implicit polarization is, to a certain extent, only accounted for through the permanent molecular dipole moment of the SPC water model), which might be especially important for the nonpolar lipid tail carbon atoms [71] owing to the very high polarizability of alkanes [72].

In the present study, the first problem was addressed through the introduction of scaling factors for the Lennard–Jones C12 parameter of ion-lipid and ion-ion interactions. The resulting modified force-field description reproduces better the experimental binding constants for the association of sodium and chloride ions with POPC bilayers, as well as the structural properties of a POPC bilayer in an aqueous sodium-chloride solution. The latter are comparable to salt-free conditions. However, ion aggregation at the bilayer-solution interface still occurs, indicating that points (ii) and (iii) above provide opportunities for future force-field improvement.

Acknowledgments

The authors thank the Innovationsfonds Forschung of the University of Freiburg for financial support and Chris Oostenbrink for critical reading of the manuscript, discussions, support and generosity in terms of computational resources.

Author Contributions

Author Contributions

M.M.R. and V.K. conceived of the MD simulations. M.M.R. performed the MD simulations. M.M.R., C.K. and V.K. analyzed the data. M.M.R., C.K. and V.K. wrote the paper.

Conflicts of Interest

Conflicts of Interest

The authors declare no conflict of interest.

References

1. Lösche M., Möhwald H. Electrostatic interactions in phospholipid membranes. II. Influence of divalent ions on monolayer structure. J. Colloid Interface Sci. 1989;131:56–67. doi: 10.1016/0021-9797(89)90145-8. [Cross Ref]
2. Clarke R.J., Lüpfert C. Influence of anions and cations on the dipole potential of phosphatidylcholine vesicles: A basis for the Hofmeister effect. Biophys. J. 1999;76:2614–2624. doi: 10.1016/S0006-3495(99)77414-X. [PubMed] [Cross Ref]
3. Renoncourt A., Vlachy N., Bauduin P., Drechsler M., Touraud D., Verbavatz J.M., Dubois M., Kunz W., Ninham B.W. Specific alkali cation effects in the transition from micelles to vesicles through salt addition. Langmuir. 2007;23:2376–2381. doi: 10.1021/la062837z. [PubMed] [Cross Ref]
4. Pabst G., Hodzic A., Štrancar J., Danner S., Rappolt M., Laggner P. Rigidification of neutral lipid bilayers in the presence of salts. Biophys. J. 2007;93:2688–2696. doi: 10.1529/biophysj.107.112615. [PubMed] [Cross Ref]
5. Klasczyk B., Knecht V., Lipowsky R., Dimova R. Interactions of alkali metal chlorides with phosphatidylcholine vesicles. Langmuir. 2010;26:18951–18958. doi: 10.1021/la103631y. [PubMed] [Cross Ref]
6. Ferber U.M., Kaggwa G., Jarvis S.P. Direct imaging of salt effects on lipid bilayer ordering at sub-molecular resolution. Eur. Biophys. J. 2011;40:329–338. doi: 10.1007/s00249-010-0650-7. [PubMed] [Cross Ref]
7. Zimmermann R., Küttner D., Renner L., Kaufmann M., Werner C. Fluidity modulation of phospholipid bilayers by electrolyte ions: Insights from fluorescence microscopy and microslit electrokinetic experiments. J. Phys. Chem. A. 2012;116:6519–6525. doi: 10.1021/jp212364q. [PubMed] [Cross Ref]
8. Rios G.M., Belleville M.-P., Paolucci-Jeanjean D. Membrane engineering in biotechnology: Quo vamus? Trends Biotechnol. 2007;25:242–246. doi: 10.1016/j.tibtech.2007.04.003. [PubMed] [Cross Ref]
9. Saeui C.T., Mathew M.P., Liu L., Urias E., Yarema K.J. Cell surface and membrane engineering: Emerging technologies and applications. J. Funct. Biomater. 2015;6:454–485. doi: 10.3390/jfb6020454. [PMC free article] [PubMed] [Cross Ref]
10. Satoh K. Determination of binding constants of Ca2+, Na+ and Cl ions to liposomal membranes of dipalmitoylphosphatidylcholine at gel phase by particle electrophoresis. Biochim. Biophys. Acta. 1995;1239:239–248. doi: 10.1016/0005-2736(95)00154-U. [PubMed] [Cross Ref]
11. Eisenberg M., Gresalfi T., Riccio T., McLaughlin S. Adsorbtion of monovalent cations to bilayer membranes containing negative phospholipids. Biochemistry. 1979;18:5213–5223. doi: 10.1021/bi00590a028. [PubMed] [Cross Ref]
12. Evans D.F., Wennerstrom H. The Colloidal Domain. 2nd ed. VCH; Weinheim, Germany: 1999.
13. Seimiya T., Ohki S. Ionic structure of phospholipid membranes, and binding of calcium ions. Biochim. Biophys. Acta. 1973;298:546–561. doi: 10.1016/0005-2736(73)90073-4. [PubMed] [Cross Ref]
14. Knecht V., Klasczyk B. Specific binding of chloride ions to lipid vesicles and implications at molecular scale. Biophys. J. 2013;104:818–824. doi: 10.1016/j.bpj.2012.12.056. [PubMed] [Cross Ref]
15. Klasczyk B., Knecht V. Validating affinities for ion-lipid association from simulation against experiment. J. Phys. Chem. A. 2011;115:10587–10595. doi: 10.1021/jp202928u. [PubMed] [Cross Ref]
16. Garcia-Celma J.J., Hatahet L., Kunz W., Fendler K. Specific anion and cation binding to lipid membranes investigated on a solid supported membrane. Langmuir. 2007;23:10074–10080. doi: 10.1021/la701188f. [PubMed] [Cross Ref]
17. Collins K.D. Ions from the Hofmeister series and osmolytes: Effects on proteins in solution and in the crystallization process. Methods. 2004;34:300–311. doi: 10.1016/j.ymeth.2004.03.021. [PubMed] [Cross Ref]
18. Leontidis E., Aroti A. Liquid expanded monolayers of lipids as model systems to understand the anionic Hofmeister series. 2. Ion partitioning is mostly a matter of size. J. Phys. Chem. B. 2009;113:1460–1467. doi: 10.1021/jp809444n. [PubMed] [Cross Ref]
19. Böckmann R.A., Hac A., Heimburg T., Grubmüller H. Effect of sodium chloride on a lipid bilayer. Biophys. J. 2003;85:1647–1655. doi: 10.1016/S0006-3495(03)74594-9. [PubMed] [Cross Ref]
20. Böckmann R.A., Grubmüller H. Multistep binding of divalent cations to phospholipid bilayers: A molecular dynamics study. Angew. Chem. Int. Ed. 2004;43:1021–1024. doi: 10.1002/anie.200352784. [PubMed] [Cross Ref]
21. Lin Y.-W., Liao L.-F. Probing interactions between uranyl ions and lipid membrane by molecular dynamics simulation. Comput. Theor. Chem. 2011;976:130–134. doi: 10.1016/j.comptc.2011.08.016. [Cross Ref]
22. Valley C.C., Perlmutter J.D., Braun A.R., Sachs J.N. NaCl interactions with phosphatidylcholine bilayers do not alter membrane structure but induce long-range ordering of ions and water. J. Membr. Biol. 2011;244:35–42. doi: 10.1007/s00232-011-9395-1. [PubMed] [Cross Ref]
23. Jarerattanachat V., Karttunen M., Wong-ekkabut J. Molecular dynamics study of oxidized lipid bilayers in NaCl solution. J. Phys. Chem. B. 2013;117:8490–8501. doi: 10.1021/jp4040612. [PubMed] [Cross Ref]
24. Sachs J.N., Woolf T.B. Understanding the Hofmeister effect in interactions between chaotropic anions and lipid bilayers: Molecular dynamics simulations. J. Am. Chem. Soc. 2003;125:8742–8743. doi: 10.1021/ja0355729. [PubMed] [Cross Ref]
25. Sachs J.N., Nanda H., Petrache H.I., Woolf T.B. Changes in phosphatidylcholine headgroup tilt and water order induced by monovalent salts: Molecular dynamics simulations. Biophys. J. 2004;86:3772–3782. doi: 10.1529/biophysj.103.035816. [PubMed] [Cross Ref]
26. Vácha R., Siu S.W.I., Petrov M., Böckmann R.A., Barucha-Kraszewska J., Jurkiewicz P., Hof M., Berkowitz M.L., Jungwirth P. Effect of alkali cations and halide anions on the DOPC lipid membrane. J. Phys. Chem. A. 2009;113:7235–7243. doi: 10.1021/jp809974e. [PubMed] [Cross Ref]
27. Vácha R., Jurkiewicz P., Petrov M., Berkowitz M.L., Böckmann R.A., Barucha-Kraszewska J., Hof M., Jungwirth P. Mechanism of interaction of monovalent ions with phosphatidyl lipid membranes. J. Phys. Chem. B. 2010;114:9504–9509. doi: 10.1021/jp102389k. [PubMed] [Cross Ref]
28. Phillips R.B., Kondev J., Theriot J. Physical Biology of the Cell. Garland Science; New York, NY, USA: 2009.
29. Kagawa R., Hirano Y., Taiji M., Yasuoka K., Yasui M. Dynamic interactions of cations, water and lipids and influence on membrane fluidity. J. Membr. Sci. 2013;435:130–136. doi: 10.1016/j.memsci.2013.02.006. [Cross Ref]
30. Zhou Y., Raphael R.M. Solution pH alters mechanical and electrical properties of phosphatidylcholine membranes: Relation between interfacial electrostatics, intramembrane potential, and bending elasticity. Biophys. J. 2007;92:2451–2462. doi: 10.1529/biophysj.106.096362. [PubMed] [Cross Ref]
31. Babu Boggara M., Faraone A., Krishnamoorti R. Effect of pH and ibuprofen on the phospholipid bilayer bending modulus. J. Phys. Chem. B. 2010;114:8061–8066. doi: 10.1021/jp100494n. [PubMed] [Cross Ref]
32. Mitkova D., Marukovich N., Ermakov Y.A., Vitkova V. Bending rigidity of phosphatidylserine-containing lipid bilayers in acidic aqueous solutions. Colloids Surf. A Physicochem. Eng. Asp. 2014;460:71–78. doi: 10.1016/j.colsurfa.2013.12.059. [Cross Ref]
33. Winterhalter M., Helfrich W. Effect of surface charge on the curvature elasticity membranes. J. Phys. Chem. 1988;92:6865–6867. doi: 10.1021/j100335a004. [Cross Ref]
34. Winterhalter M., Helfrich W. Bending elasticity of electrical charged bilayers: Coupled monolayers, neutral surfaces, and balancing stresses. J. Phys. Chem. 1992;96:327–330. doi: 10.1021/j100180a060. [Cross Ref]
35. Berkowitz M.L., Raghavan K. Interaction forces between membrane surfaces. Role of electrostatic concepts. In: Blank M., Vodyanoy I., editors. Biomembrane Electrochemistry. Volume 235. American Chemical Society; Washington, DC, USA: 1994. pp. 3–25.
36. Pandit S.A., Bostick D., Berkowitz M.L. Molecular dynamics simulation of a dipalmitoylphosphatidylcholine bilayer with NaCl. Biophys. J. 2003;84:3743–3750. doi: 10.1016/S0006-3495(03)75102-9. [PubMed] [Cross Ref]
37. De Jong D.H., Schäfer L., de Vries A.H., Marrink S.J., Berendsen H.J.C., Grubmüller H. Determining equilibrium constants for dimerization reactions from molecular dynamics simulations. J. Comput. Chem. 2011;32:1919–1928. doi: 10.1002/jcc.21776. [PubMed] [Cross Ref]
38. Smirnova Y.G., Aeffner S., Risselada H.J., Marrink S.J., Müller M., Knecht V. Interbilayer repulsion forces between tension-free lipid bilayers from simulation. Soft Matter. 2013;9:10705–10718. doi: 10.1039/c3sm51771c. [Cross Ref]
39. Landau L.D., Lifshitz E.M. Theory of Elasticity. 3rd ed. Elsevier; Oxford, UK: 1986.
40. Helfrich W. Out-of-plane fluctuations of lipid bilayers. Z. Naturforsch. 1975;30:841–842. [PubMed]
41. Bermúdez H., Hammer D.A., Discher D.E. Effect of bilayer thickness on membrane bending rigidity. Langmuir. 2004;20:540–543. doi: 10.1021/la035497f. [PubMed] [Cross Ref]
42. Szleifer I., Kramer D., Ben-Shaul A., Roux D., Gelbart W.M. Curvature elasticity of pure and mixed surfactant films. Phys. Rev. Lett. 1988;60:1966–1969. doi: 10.1103/PhysRevLett.60.1966. [PubMed] [Cross Ref]
43. Zhang W., Róg T., Gurtovenko A.A., Vattulainen I., Karttunen M. Atomic-scale structure and electrostatics of anionic palmitoyloleoylphosphatidylglycerol lipid bilayers with Na+ counterions. Biophys. J. 2007;92:1114–1124. [PubMed]
44. Von Deuster C.I.E., Knecht V. Competing interactions for antimicrobial selectivity based on charge complementarity. Biochim. Biophys. Acta. 2011;1808:2867–2876. doi: 10.1016/j.bbamem.2011.08.005. [PubMed] [Cross Ref]
45. Knecht V., Klasczyk B., Dimova R. Macro-versus microscopic view on the electrokinetics of a water-membrane interface. Langmuir. 2013;29:7939–7948. doi: 10.1021/la400342m. [PubMed] [Cross Ref]
46. Reif M.M., Hünenberger P.H. Origin of asymmetric solvation effects for ions in water and organic solvents investigated using molecular dynamics simulations: The Swain acity-basity scale revisited. J. Phys. Chem. B. 2016;120:8485–8517. doi: 10.1021/acs.jpcb.6b02156. [PubMed] [Cross Ref]
47. Lindahl E., Edholm O. Mesoscopic undulations and thickness fluctuations in lipid bilayers from molecular dynamics simulations. Biophys. J. 2000;79:426–433. doi: 10.1016/S0006-3495(00)76304-1. [PubMed] [Cross Ref]
48. Evans E., Rawicz W. Entropy-driven tension and bending elasticity in condensed-fluid membranes. Phys. Rev. Lett. 1990;64:2094–2097. doi: 10.1103/PhysRevLett.64.2094. [PubMed] [Cross Ref]
49. Nagle J.F., Tristram-Nagle S. Structure of lipid bilayers. Biochim. Biophys. Acta. 2000;1469:159–195. doi: 10.1016/S0304-4157(00)00016-2. [PMC free article] [PubMed] [Cross Ref]
50. Tabony J., Perly B. Quasielastic neutron scattering measurements of fast local translational diffusion of lipid molecules in phospholipid bilayers. Biochim. Biophys. Acta. 1991;1063:67–72. doi: 10.1016/0005-2736(91)90354-B. [PubMed] [Cross Ref]
51. König S., Pfeiffer W., Bayerl T., Richter D., Sackmann E. Molecular dynamics of lipid bilayers studied by incoherent quasi-elastic neutron scattering. J. Phys. II Fr. 1992;2:1589–1615. doi: 10.1051/jp2:1992100. [Cross Ref]
52. Jeon J.-H., Martinez-Seara Monne H., Javanainen M., Metzler R. Anomalous diffusion of phospholipids and cholesterols in a lipid bilayer and its origins. Phys. Rev. Lett. 2012;109 doi: 10.1103/PhysRevLett.109.188103. [PubMed] [Cross Ref]
53. Van der Spoel D., Lindahl E., Hess B., Groenhof G., Mark A.E., Berendsen H.J.C. GROMACS: Fast, flexible, and free. J. Comput. Chem. 2005;26:1701–1718. doi: 10.1002/jcc.20291. [PubMed] [Cross Ref]
54. Hess B., Kutzner C., van der Spoel D., Lindahl E. GROMACS 4: Algorithms for highly efficient, load-balanced, and scalable molecular simulation. J. Chem. Theory Comput. 2008;4:435–447. doi: 10.1021/ct700301q. [PubMed] [Cross Ref]
55. Berendsen H.J.C., Postma J.P.M., van Gunsteren W.F., di Nola A., Haak J.R. Molecular dynamics with coupling to an external bath. J. Chem. Phys. 1984;81:3684–3690. doi: 10.1063/1.448118. [Cross Ref]
56. Parrinello M., Rahman A. Crystal structure and pair potentials: A molecular-dynamics study. Phys. Rev. Lett. 1980;45:1196–1199. doi: 10.1103/PhysRevLett.45.1196. [Cross Ref]
57. Hess B., Bekker H., Berendsen H.J.C., Fraaije J.G.E.M. LINCS: A linear constraint solver for molecular simulations. J. Comput. Chem. 1997;18:1463–1472. doi: 10.1002/(SICI)1096-987X(199709)18:12<1463::AID-JCC4>3.0.CO;2-H. [Cross Ref]
58. Hockney R.W. The potential calculation and some applications. Methods Comput. Phys. 1970;9:135–211.
59. Darden T., York D., Pedersen L. Particle mesh Ewald: An Nlog(N) method for Ewald sums in large systems. J. Chem. Phys. 1993;98:10089–10092. doi: 10.1063/1.464397. [Cross Ref]
60. Berger O., Edholm O., Jähnig F. Molecular dynamics simulations of a fluid bilayer of dipalmitoylphosphatidylcholine at full hydration, constant pressure, and constant temperature. Biophys. J. 1997;72:2002–2013. doi: 10.1016/S0006-3495(97)78845-3. [PubMed] [Cross Ref]
61. Van Gunsteren W.F., Billeter S.R., Eising A.A., Hünenberger P.H., Krüger P., Mark A.E., Scott W.R.P., Tironi I.G. Biomolecular Simulation: The GROMOS96 Manual and User Guide. Verlag der Fachvereine; Zürich, Switzerland: 1996.
62. Berendsen H.J.C., Postma J.P.M., van Gunsteren W.F., Hermans J. Interaction models for water in relation to protein hydration. In: Pullman B., editor. Intermolecular Forces. Reidel; Dordrecht, The Netherlands: 1981. pp. 331–342.
63. Hess B., van der Vegt N.F.A. Cation specific binding with protein surface charges. Proc. Natl. Acad. Sci. USA. 2009;106:13296–13300. doi: 10.1073/pnas.0902904106. [PubMed] [Cross Ref]
64. Weerasinghe S., Smith P.E. A Kirkwood-Buff derived force field for sodium chloride in water. J. Chem. Phys. 2003;119:11342–11349. doi: 10.1063/1.1622372. [Cross Ref]
65. Vorobjev Y.N., Hermans J. A critical analysis of methods of calculation of a potential in simulated polar liquids : Strong arguments in favor of “Molecule-based” summation and of vacuum boundary conditions in Ewald summations. J. Phys. Chem. B. 1999;103:10234–10242. doi: 10.1021/jp984211q. [Cross Ref]
66. Kastenholz M.A., Hünenberger P.H. Computation of methodology-independent ionic solvation free energies from molecular simulations: I. The electrostatic potential in molecular liquids. J. Chem. Phys. 2006;124 doi: 10.1063/1.2172593. [PubMed] [Cross Ref]
67. Waheed Q., Edholm O. Undulation contributions to the area compressibility in lipid bilayer simulations. Biophys. J. 2009;97:2754–2760. doi: 10.1016/j.bpj.2009.08.048. [PubMed] [Cross Ref]
68. Frigo M., Johnson S.G. The design and implementation of FFTW3. Proc. IEEE. 2005;93:216–231. doi: 10.1109/JPROC.2004.840301. [Cross Ref]
69. Allen M.P., Tildesley D.J. Computer Simulation of Liquids. Oxford University Press; New York, NY, USA: 1987.
70. Reif M.M., Hünenberger P.H. Computation of methodology-independent single-ion solvation properties from molecular simulations. IV. Optimized Lennard–Jones parameter sets for the alkali and halide ions in water. J. Chem. Phys. 2011;134 doi: 10.1063/1.3567022. [PubMed] [Cross Ref]
71. Szklarczyk O.M., Bachmann S.J., van Gunsteren W.F. A polarizable empirical force field for molecular dynamics simulation of liquid hydrocarbons. J. Comput. Chem. 2014;35:789–801. doi: 10.1002/jcc.23551. [PubMed] [Cross Ref]
72. Anslyn E.V., Dougherty D.A. Modern Physical Organic Chemistry. University Science Books; Philadelphia, PA, USA: 2006.

Articles from Membranes are provided here courtesy of Multidisciplinary Digital Publishing Institute (MDPI)