Monovalent Ions and Water Dipoles in Contact with Dipolar Zwitterionic Lipid Headgroups-Theory and MD Simulations

The lipid bilayer is a basic building block of biological membranes and can be pictured as a barrier separating two compartments filled with electrolyte solution. Artificial planar lipid bilayers are therefore commonly used as model systems to study the physical and electrical properties of the cell membranes in contact with electrolyte solution. Among them the glycerol-based polar phospholipids which have dipolar, but electrically neutral head groups, are most frequently used in formation of artificial lipid bilayers. In this work the electrical properties of the lipid layer composed of zwitterionic lipids with non-zero dipole moments are studied theoretically. In the model, the zwitterionic lipid bilayer is assumed to be in contact with aqueous solution of monovalent salt ions. The orientational ordering of water, resulting in spatial variation of permittivity, is explicitly taken into account. It is shown that due to saturation effect in orientational ordering of water dipoles the relative permittivity in the zwitterionic headgroup region is decreased, while the corresponding electric potential becomes strongly negative. Some of the predictions of the presented mean-field theoretical consideration are critically evaluated using the results of molecular dynamics (MD) simulation.

permittivity in the zwitterionic headgroup region is decreased, while the corresponding electric potential becomes strongly negative. Some of the predictions of the presented mean-field theoretical consideration are critically evaluated using the results of molecular dynamics (MD) simulation.

Introduction
The lipid bilayer is a basic building block of cell membranes. Although the cell membrane is a highly heterogeneous structure, composed of lipids, proteins, carbohydrates and other components [1][2][3], a large number of physical properties of biological (cell) membrane which are important in physiological processes [4,5] can be studied also in pure planar lipid bilayer systems [6,7]. Due to its double-sided accessability in experiments the planar lipid bilayer is suitable for experimental manipulation [8,9] even if some other pure lipids systems, like lipid vesicles, may mimic also three-dimensional geometry of the cell membrane.
Among others a planar lipid bilayer can be constructed on a small hole positioned on a barrier separating two compartments filled with electrolyte solution. Experimental setup usually consists of two teflon pieces with separate compartments and 100 µm-1 mm diameter hole in between. The diameter of the hole mostly depends on the procedure used in planar lipid bilayer formation [10]. In experiments, the electrodes are dipped into electrolyte solution approximately 1 cm away form the barrier. Electrolyte solution is basically an aqueous solution of salt composed of water molecules ( Figure 1) and monovalent positively and negatively charged ions. Such experimental setup is ordinarily used to determine various mechanical and electrical properties of planar lipid bilayers [9,11,12]. Figure 1. A single water molecule is considered as a sphere with permittivity n 2 and point-like rigid (permanent) dipole with dipole moment p at the center of the sphere [13], where n is the optical refractive index of water.
In general, some of the lipid molecules bear net negative charge, while other lipids like glycerophospholipids are electrically neutral [14]. Glycerophospholipids (see Figure 2) are composed of dipolar (zwitterionic) headgroup and two nonpolar tails [1,2,6,7,14]. The tails are hydrophobic fatty acids. Due to its hydrophilic nature [14] the headgroups of lipids are in contact with aqueous solution. The negative charges of the dipolar lipid headgroup ( Figure 2) are in contact with nonpolar tails on the one side (left side in Figure 2) and with electrolyte solution on the other (right side in Figure 2), thus its charge distribution is described in the model as negatively charged surface at x = 0 ( Figure 2). The positive charges of the headgroups in planar bilayer of dipolar lipids protrude further in the electrolyte solution as depicted in Figure 2. In the headgroup region of planar dipolar lipid layer the cations are attracted towards the negatively charged surface at x = 0, while anions are depleted from this region. Using the MD simulations, the accumulation of sodium cations near the phosphate groups and accumulation of anions near the choline groups in the DOPC bilayer were predicted by [19]. These theoretical predictions were supported also by the results of fluorescence spectroscopy experiments on zwitterionic phospholipid bilayers which show that the addition of salt into the solution being in contact with the lipid bilayer restricts the mobility of the hydrated lipid headgroups and their lateral diffusion when compared to pure water solution [19]. In the high electric field of dipolar lipid headgroups (see [20] and the references therein) the water dipoles are expected to be oriented towards the charged plane at x = 0 (see for example [13,21]). Due to accumulation of cations and saturation effect in polarization the electric potential decreases towards the x = 0 plane [22].
Most of the electrostatic models of electrolyte solution in contact with lipid surfaces [20,[23][24][25][26] assume that the dielectric permittivity in the electrolyte solution is constant. In the absence of an explicit consideration of orientational ordering of water molecules the assumption of constant permittivity is the consequence of the assumption of the constant number density of water molecules in the system [13]. But actually, close to the membrane surface the orientation and depletion of water molecules may result in strong spatial variation of permittivity [13,[27][28][29][30][31].
In this paper the effect of the nonhomogeneous volume charge distribution in the headgroup region of the planar dipolar lipid layer on the space dependent electric potential and permittivity is presented. An analytical mean-field model, based on the previously developed Langevin-Poisson-Boltzman (LPB) model [21], is introduced. The relative (dielectric) permittivity, related to the electric field strength was analysed in the headgroup region and its close vicinity. To test the predictions of the model, the realistic values of the input model parameters, previously determined in measurements or simulations for DPPC lipid molecules, were used. The results of an analytical model are compared with the results of molecular dynamic simulation (MD) of DPPC planar lipid bilayer.

Modified Langevin-Poisson-Boltzmann (MLPB) Model
The Langevin-Poisson-Boltzmann (LPB) model [21] is generalized to take into account the cavity field [13] in the saturation regime. In addition, the electronic polarization of the water is taken into account by assuming that the point-like rigid (permanent) dipole embedded in the center of the sphere with a volume equal to the average volume of a water molecule in the electrolyte solution ( Figure 1) [13,32]. The permittivity of the sphere is taken to be n 2 , where n = 1.33 is the optical refractive index of water. The relative (effective) permittivity of the electrolyte solution ε r can be then expressed as [13,32] : where P = |P| is the magnitude of the polarization vector due to a net orientation of permanent point-like water dipoles having dipole moment p, ε 0 is the permittivity of the free space, while E = |E| is the magnitude of the electric field strength. The absolute value of polarization P (x) is given by [13] : where n w (x) is the space dependency of the number density of water molecules, p 0 is the magnitude of the external dipole moment p e = (3/(2 + n 2 )) p (see also Figure 1) [13,32], In the following, for the sake of simplicity the finite volume of ions and water molecules is not taken into account as in [13] and consequently the number density of water molecules is considered to be constant all over the solution and equal to its bulk value n w0 , i.e., n w (x) = n 0w from where it follows : Combining Equations 1 and 4 yields space dependency of permittivity within modified Langevin-Poisson-Boltzmann model (MLPB model) in the form: In the limit of vanishing electric field strength (E(x) −→ 0 everywhere in the solution) the above equation for relative permittivity ε r (x) gives the classical Onsager expression: At room temperature (298 K) the above Equation 7 gives ε r = 78.5 for bulk solution. The parameters p 0 and n 0w /N A are 3.1 Debye and 55 mol/l, respectively.

Poisson Equation
In the model the dipolar lipid headgroup is described by two charges at fixed distance D, i.e., it is assumed that the headgroup has non-zero dipole moment. The negative charges of the phosphate groups of dipolar (zwitterionic) lipids are described by negative surface charge density σ = −e 0 /a 0 at x = 0 (see Figure 2), where a 0 is the area per lipid. The orientational ordering of water is taken into account assuming the spatial dependence of permittivity ε r (x) as described by Equation 6.
The corresponding Poisson equation in planar geometry can thus be written in the form (see e.g., [22]) : where ϕ(x) is the electric potential, ρ Zw (x) is the macroscopic (net) volume charge density of positive charges of dipolar (zwitterionic) headgroups, ρ ions (x) is the macroscopic (net) volume charge density of co-ions (n − ) and counter-ions (n + ) of the electrolyte solution (see Figure 3). Since we neglect the finite volumes of the salt ions and water molecules [21,22] the co-ions and counter-ions are assumed to be distributed according to Boltzmann distribution functions [20,[22][23][24][25][26]) : where e 0 is the unit charge and n 0 bulk number density of salt co-ions and counter-ions. The lipid headgroups can be oriented at various angles ω relative to the normal vector to the planar lipid layer (Figure 2), hence the volume charge density due to the positive charges of the lipid dipolar headgroups can be written in the form: where a 0 is the area per lipid, while P(x) is the probability density function indicating the probability that the positive charge of the dipolar lipid headgroup is located at the distance x from the negatively charged surface at x = 0 : where x ≤ D. Equation 13 neglects the finite volume of lipid headgroups. The normalization condition yields : Using Equations 11-13 and 15 it follows from Equation 8 : The boundary conditions are (see for example [22]) : where in Equation 17 the surface charge density σ = −e 0 /a 0 . Note that the area per lipid a 0 is different in gel and liquid phase. E = |ϕ ′ |. In numerical calculations the distance from the negatively charged surface x was limited to 12 nm, where the boundary condition stated in Equation 18 was applied. The modified LPB equation (Equation 16) was solved numerically as described in the Appendix A. All the results were obtained using a 0 = 0.48 nm 2 and a 0 = 0.60 nm 2 corresponding to DPPC lipid in gel-crystalline (below 314 K) and liquid-crystalline phase (above 314 K), respectively [33]. Other values of model parameters were: the dipole moment of water p 0 = 3.1 Debye, bulk concentration of salt n 0 /N A = 0.1 mol/l and concentration of water n 0w /N A = 55 mol/l. N A is Avogadro number.

Molecular Dynamics Simulations (MD)
The molecular dynamics (MD) model of DPPC planar lipid bilayer was constructed in NAMD program using all molecule performance CHARMM 36 force field. The model consists of 256 lipid units and 20174 water molecules. The solvent was 450 mM KCl modeled by 153 K + and 153 Cl − ions [34,35]. Chemical bonds between hydrogen and heavy atoms were constrained to their equilibrium value. Long-range electrostatic forces were taken into account using a fast implementation of the particle mesh Ewald (PME) method [36,37]. The model was examined at constant pressure (1.013 × 10 5 Pa) and constant temperature (232 K) employing Langevin dynamics and the Langevin piston method. The equations of motion were integrated using the multiple time-step algorithm. A time step of 2.0 fs was employed. Short-and long-range forces were calculated every one and two time steps, respectively.
The model was equilibrated and followed 30 ns. The last 15 ns of the simulation were used for extraction of dipole orientation angle. From the P and N atoms positions the dipole was determined for all 256 lipids in each of 1500 simulation frames and exported to Matlab2012b. The distribution of vector amplitude corresponding to distance D between charges was extracted as well as distribution of the angle ω between the dipole and normal vector to the planar lipid bilayer plane (Figure 2). To obtain the probability density P(x), projection of each headgroup dipole on the normal vector to the planar lipid bilayer plane was calculated. The average distance between P and N atoms (0.42 nm) was used as a parameter D in MLPB model (Equation 16).

Results
Electric potential ϕ and relative permittivity ε r as a function of the distance from the charged planar surface (x = 0) is presented in Figure 4. The results are presented for two values of the temperature: T = 310 K (a 0 = 0.48 nm 2 ) and T = 323 K (a 0 = 0.60 nm 2 ). It can be seen in Figure 4 that the relative permittivity ε r (x) is considerably decreased in the vicinity of charged planar surface. At the charged planar surface (x = 0) the value of ε r (x) drops to 44 at T = 310 K and to 55 at T = 323 K. The effect of the charged planar surface at x = 0 is very weak already at the distance x = D. Far away from the surface (x = 0) the values of ε r (x) is 75.5 at temperature T = 310 K and 72.6 for temperature T = 323 K. The electric potential in the vicinity of the charged planar surface is considerably negative. It is − 60 mV for temperature T = 310 K and − 54 mV for temperature T = 323 K.
Charge density profile of co-ions (n − ) and counter-ions (n + ) of the electrolyte solution for two temperatures 310 K and 323 K can be seen in Figure 3. Near the negatively charged planar surface at x = 0 one can observe strong accumulation of positively charged counter-ions and depletion of the negatively charged co-ions. With increasing distance from the headgroup region (0 < x ≤ D), i.e., for x larger than D, the number density of co-ions (n − ) decreases and the number density of counter-ions (n + ) increases. Far away from the charged planar surface, the concentration of counter-ions (n + ) equals the concentration of co-ions (n − ) corresponding to electroneutrality condition in bulk solution. At higher temperature (323 K) the DPPC has increased area per lipid a 0 = 0.60 nm 2 resulting in lower area density of the lipid molecules and hence less negative surface charge density at x = 0 plane as in the case of lower temperature (310 K) . Consequently, also the calculated ion number density profiles are lower for x > D (Figure 3).
The average headgroup orientation angle < ω > as a function of the temperature T can be seen on Figure 5. At DPPC phase transition temperature (314 K) the value of < ω > can not be calculated as phase transition effect is not included in MLPB model. The average dipole orientation angle < ω >, calculated from P(x) as described in Appendix B, is not temperature dependent. At DPPC liquid-crystalline phase average headgroup orientation angle < ω > = 69.36 degrees, which agrees with median angle < ω > between N and P dipole and normal vector obtained in MD simulation (< ω > = 68 degrees). The difference of < ω > between liquid-crystalline and gel-crystalline phase is a consequence of a different values of area per lipid a 0 in gel-crystalline phase and liquid-crystalline phase.   Figure 2) as a function of the temperature T . Area per lipid a 0 = 0.48 nm 2 below 314 K, corresponding to DPPC lipid gel-crystalline phase (full line) and a 0 = 0.60 nm 2 above 314 K, corresponding to DPPC lipid liquid-crystalline phase (dashed line). A gap near the phase transition temperature (314 K) is present, because phase transition effect is not included in MLPB model. Average dipole orientation angle < ω > was calculated from P(x) as described in Appendix B. The values of other model parameters are the same as in Figure 4.

Discussion and Conclusions
The comparison between the probability density P(x) calculated within MLPB model (Equation 13) and P(x) obtained in MD simulations, can be seen on Figure 6. In MLPB model the function P(x) is steeply increasing in the vicinity of x = 0 plane ( Figure 6, plot A) which is a consequence of the exponential Boltzmann factor in the function P(x) (Equation 13). This result is clearly not in accordance with MD simulation (Figure 6, plot B), where the function P(x) is saturated. The discrepancy between the predictions of MLPB model and MD simulations arises due to the fact that the finite volumes of lipid headgroups and finite volumes of ions and water molecules are not considered in MLPB model. Taking into account the finite volume of lipid headgroups within the lattice statistics approach (see also [13,24]) yields for probability density that the positive charge of the dipolar lipid headgroup is located at certain distance x from the x = 0 surface in the form : where Λ is determined from normalization Equation 14. The parameter α is equal to the ratio between the average volume of the positively charged parts of dipolar (zwitterionic) headgroups and the average volume of the salt solution in the headgroup region. Equation 21 predicts the saturation of the probability density function P(x), corresponding to the close packing of the lipid headgroups in accordance with the results of MD simulations. Figure 6 thus shows the dependence of P(x) calculated by using Equation 21 for the values of the parameter α. It can be seen that taking into account the finite volume of lipid headgroups leads to better agreement between the predicted P(x) dependencies within MLPB model and MD simulations. In the limit of α → ∞ (when all lattice sites are occupied by the headgroups) the probability density function P(x) becomes constant as expected. On the other hand, in the limit of small values of α (i.e., negligible volume of the headgroups) the probability density P(x), calculated by using Equation 21, approaches to the probability denisty P(x) determined by Equation 13 (curve A in Figure 6).  Figure 4.
Although the relative permittivity ε r is usually considered as a constant, it is shown in present paper that ε r can considerably change within the dipolar (zwitterionic) lipid headgroup region. Consequently, the electric potential in this region is substantially decreased.
To conclude, our model shows that the effect of decreasing potential and permittivity has an impact only in the headgroup region of dipolar zwitterionic lipids and its close vicinity. The average orientation angle of the zwitterionic lipid headgroup dipole (< ω >) predicted within the presented MLPB model is comparable with the results obtained in MD simulation.

A.1. Variant A
Equation 16 can be rewritten as : where ψ(x) is reduced electric potential : and The boundary conditions (Equations 17-20) can be rewritten as : Equation 22 was solved in Matlab using the standard function for the boundary value problems (bvp4c). The second part of the right hand side of Equation 22 is defined only within headgroup region, therefore it is equal to zero for x > D. The value of ε r (x) (Equation 6) and the integral ∫ exp(−ψ(x))dx and dε r (x)/dx in Equation 22 were calculated outside of bvp4c function. The dε r (x)/dx was derived from Equation 6, where the Langevin function L(x) was expanded for small values of electric field into Taylor series up to the cubic term L(x) ≈ x/3 − x 3 /45. ConsideringẼ(x) = dψ(x) dx , can be written : where the constants are : ) 2 n 0w p 2 0 β ε 0 (32) For high values of electric field the ε r (x) and the dε r (x)/dx are derived from Equation 6 with non-expanded Langevin function : where constants are the same as above (Equations 31-33).

A.2. Variant B
By inserting the Equation 6 and its derivative into Equation 8, the Equation 8 can be rewritten as : where the constants K1, K2, K3 are defined as above (Equations 31-33). Equation 36 was solved in Matlab using standard function for boundary value problem (bvp4c) with multipoint boundary values (Equations 25-28). The value of the integral ∫ exp(−ψ(x))dx in Equation 36 was calculated in iteration process outside of bvp4c function.

B. Average Orientation of Lipid Head-Groups
From Equations (13 and 15) it follows: Average orientation angle of the headgroup dipole (< ω >) can be then written as normalized distribution function: where the dipole orientation angle is (see also Figure 2) :