A Poisson-Boltzmann approach for a lipid membrane in an electric field

The behavior of a non-conductive quasi-planar lipid membrane in an electrolyte and in a static (DC) electric field is investigated theoretically in the nonlinear (Poisson-Boltzmann) regime. Electrostatic effects due to charges in the membrane lipids and in the double layers lead to corrections to the membrane elastic moduli which are analyzed here. We show that, especially in the low salt limit, i) the electrostatic contribution to the membrane's surface tension due to the Debye layers crosses over from a quadratic behavior in the externally applied voltage to a linear voltage regime. ii) the contribution to the membrane's bending modulus due to the Debye layers saturates for high voltages. Nevertheless, the membrane undulation instability due to an effectively negative surface tension as predicted by linear Debye-H\"uckel theory is shown to persist in the nonlinear, high voltage regime.


Introduction
Bilayers formed from lipid molecules are an essential component of the membranes of biological cells. The mechanical properties of membranes at equilibrium are characterized by two elastic moduli, the surface tension and the bending modulus [1]. These moduli typically depend on electrostatic properties, and their modifications in case of charged membranes/surfaces in an electrolyte have been examined theoretically in the period 1980-90's as reviewed e.g. in Ref. [2]: they have been first derived by Winterhalter and Helfrich [3] within linearized Debye-Hückel (DH) approximation, then by Lekkerkerker [4] in the nonlinear Poisson-Boltzmann (PB) regime for charged monolayers and by Ninham et al. [5] for charged symmetric bilayers. Later on Helfrich et al. revisited the question of the electrostatic corrections to the bending modulus of charged symmetric bilayers [6].
Nowadays, the study of deformations of membranes or vesicles in electric fields is an active field of research linked to many biotechnological applications. For instance, the application of electric fields is used to produce artificial vesicles from lipid films (electroformation), or to create pores in vesicles (electroporation), which is an important route for drug delivery. Both processes are widely used experimentally although they are not well understood theoretically. The effects of electric fields on giant unilamellar vesicles have been reviewed recently in Ref. [7]. This system shows a rich panel of possible behaviors and morphological transitions depending on experimental conditions -electric field frequency, conductivities of the medium and of the membrane, salt concentration, etc. These observations are supported to a large extend by theoretical modeling [8].
Originally it was Helfrich et al. who pointed out that the deformation of lipid vesicles in electric fields can serve as a means to determine the electrostatic corrections to the membrane elastic moduli [9]. Besides this observation on the macroscopic scale other techniques can provide valuable insights into the moduli corrections such as AFM, impedance spectroscopy [10], neutron reflectivity [11] and X-ray scattering [12]. Recently [12], X-ray scattering experiments have been carried out on a system of two superposed lipid membranes in an AC electric field. In trying to analyze this data, we noted that these experiments have been carried out at relatively high voltages, in a regime where the linearized DH approach may not be applicable. In order to describe such a situation theoretically, we extend previous work [13,14,15] based on the DH approach to the nonlinear PB regime, which is more suitable for realistic situations in which the induced surface charges on the membrane are large.
In this paper, we present a simple approach to calculate electrostatic corrections in the elastic moduli of a quasi-planar lipid membrane. The membrane is assumed to be non-conductive to ions, non-permeable to water and electrically neutral, it is subjected to a normal DC electric field and embedded in an electrolyte described by the Poisson-Boltzmann equation. In this situation, the electric field leads to an accumulation of charges on both sides of the membrane, which affect the mechanical properties of the membrane. The electrostatic corrections to the elastic moduli can be used to analyze the instability of a lipid membrane in an applied DC electric field. In contrast to previous work, based on a free energy approach [4,6,16], our method is based on a calculation of the stress (or force) balance at the membrane surface using electrokinetic equations [17]. This method is related to the work of Kumaran [18] who used a similar approach in the context of equilibrium charged membranes. Two points are worth mentioning: first, our approach is able to describe the capacitive effects of the membrane and of the Debye layers while keeping the simplicity of the zero thickness approximation on which most of the literature on lipid membranes is based. Second, our approach can include non-equilibrium effects which can not be described within the free energy approach. For instance, in Refs. [13,14,15] we investigated the effects of ionic currents flowing through the membrane, which in turn affect the fluid flow near the membrane. Other types of non-equilibrium effects that could be included in that framework are those arising from the stochasticity of ion channels.

Model equations
We consider a steady (DC) voltage V between two electrodes at a fixed distance z = ±L/2, applied to an initially flat membrane located at z = 0. The membrane is embedded in an electrolyte of monovalent ions with densities n + and n − . The membrane is treated as non-conductive for both ion species and is (effectively) uncharged; thus we focus solely on capacitive effects. A point on the membrane is characterized within the Monge representation by a height function h(r ⊥ ), where r ⊥ is a two-dimensional in-plane vector in the membrane.
In the electrolyte, the electric potential φ obeys Poisson's equation Here e is the elementary charge, ǫ is the dielectric constant of the electrolyte and we have introduced half of the charge density, For the sake of simplicity, we assume a symmetric 1 : 1 electrolyte, so that far away from the membrane n + = n − = n * , and the total system is electrically neutral. The densities of the ion species obey the Poisson-Nernst-Planck equations with ionic current densities where k B T is the thermal energy. For simplicity we consider the case where both ion types have the same diffusion coefficient D and neglected various corrections for concentrated solutions [19].
As boundary conditions (BC), the potential at the electrodes is externally held at This BC is oversimplified for real electrodes, but captures the main effects of the electric field, see the discussion in Ref. [15]. The distance between the electrodes is assumed to be much larger than the Debye length, L ≫ λ D = κ −1 . In that case, the bulk electrolyte is quasi-neutral, n + = n − = n * , with negligible charge density (compared to the total salt concentration), so that far from the membrane The BC at the membrane is crucial to recover the correct physical behavior. Let n be the unit vector normal to the membrane. We use the Robin-type BC where is a length scale that contains the membrane thickness d and the ratio of the dielectric constant of the electrolyte, ǫ, and of the membrane, ǫ m . This BC was originally developed for electrodes sustaining Faradaic current [20,21] or charging capacitively [22]. In Refs. [23,14,15], this BC was derived, and was shown to properly account for the jump in the charge distribution which occurs near the membrane as a result of the dielectric mismatch between the membrane and the surrounding electrolyte.

Poisson-Boltzmann approach for a membrane in an external potential
Here we show how the well-known solution of the Poisson-Boltzmann (PB) equation for a single charged plate in an electrolyte can be used to describe the present situation of a capacitive membrane with induced Debye layers in an external potential. In a steady state situation and when there is no electric current through the membrane, one obtains from Eqs. (3,4) After a direct integration using the BCs from above, one obtains and insertion in Poisson's equation then yields the Poisson-Boltzmann (PB) equation Linearization (for φ − V 2 ≪ 1) leads to the well-known Debye-Hückel (DH) equation, where and κ −1 = λ D is the Debye length that defines the characteristic length scale for charge relaxation in the electrolyte. The nonlinear PB equation (11) for the planar case can be integrated analytically [2], leading to for z > 0. The expressions for z < 0 can be obtained using the symmetry of the system: Just as in the classical PB solution for a single charged plate in contact with an electrolyte, the non-dimensional parameter c is determined by the BC for the potential at the membrane. For a flat charged surface surrounded by an electrolyte [2], c is given by a simple quadratic equation and can be expressed in terms of the ratio of the two characteristic length scales: the Gouy-Chapman length b = 2ǫk B T e|σ| (with |σ| the charge density of the surface) and the Debye length λ D = κ −1 . In contrast, in the case of a membrane in an electric field, we obtain from Eq. (7) the following nonlinear equation Note that two dimensionless parameters enter this equation: the ratio of electrostatic to thermal energy, and the dimensionless parameter [21] which quantifies the electrical coupling between the membrane and the Debye layers. More precisely, depending on whether this parameter is large or small with respect to one, the capacitance of the diffuse part of the double layer, C D = ǫκ, or the capacitance of the membrane, C m = ǫ m /d, dominates the overall voltage drop. The non-dimensional parameter c = c(λ m ,V ) which is determined by Eq. (16) is related to the potential at the membrane, φ(0 + ), and to the charge density at the membrane, ρ(0 + ), by the following relations From Eq. (19), one can see that the values of c must be restricted to the interval [0, 1].
In the limit of small voltages,V ≪ 1, there is a linear relation between c and the charge 8en * (2+κλm) , in accordance with the calculations based on the DH approximation of Ref. [15]. Fig. 1a) shows the solution c of Eq. (16) as a function of external voltage for various values of κ, which corresponds to varying the amount of salt since κ ∝ √ n * . Clearly, for low salt, the linear approximation remains valid only for rather small voltages (for instance it holds only for V 0.1V in pure water), while for high salt V = 5V is still in the linear regime. Fig. 1b) shows a comparison of the charge distribution for the nonlinear PB (blue) and the linear DH solution (green). As expected, the figure shows that the DH approximation underestimates the surface charge on the membrane layers as compared to the PB calculation. The figure also shows the distribution of the positive and negative ions, which both tend to n * far from the interface as a result of electroneutrality.
Although the present situation differs from the case of a single charged plate in an PB electrolyte, the structure of the solutions Eqs. (14,15) is very similar in both problems. Because of this, there is an equivalent to the Contact theorem [2], which relates in the single charged plate problem the surface charge density to the limiting value of the potential/ionic density at the plate: namely, one can give the effective surface charge σ eff for the charged plate problem that creates the same voltage/charge distributions as the capacitive membrane in the external field. This effective surface charge reads

Corrections to membrane elastic moduli
Surface tension. The electrostatic corrections to the surface tension can be calculated directly from the stresses acting on the membrane in the flat configuration (also called the base state) as explained in Refs. [24,15]. The total stress tensor reads which contains the pressure, the viscous stresses in the surrounding fluid (the electrolyte) and the Maxwell stress due to the electrostatic field. η is the viscosity of the electrolyte and v its velocity field. The electric field is given by E = −∇φ.
In the base state, the electric field is oriented in z-direction, and the condition Using Eqs. (14,15) and imposing p(z → ∞) = 0, this is readily solved by and similarly with z → −z for z < 0. Let us call S a closed surface englobing the membrane with the normal vector field n. The force acting on the surface S in the x-direction (chosen to be the direction of the lateral stress) can be calculated from the stress tensor as F x = S x · τ · n dS. Since τ is divergence free, the surface S can be deformed, for convenience to a cube of size L, and it is easy to see that the integral is non-zero only on the faces of the cube with the normal along ±x. With dS = Ldz and ∆Σ = F x /L for n = +x, we arrive at ∆Σ = From Eq. (22), one has τ xx (z) = −p(z) − ǫ 2 (∂ z φ) 2 , where φ(z) and p(z) are the potential and the pressure in the base state given above. Upon integration (using Lκ ≫ 1), one obtains the corrections to the surface tension, as a sum of two terms.
First, there is the external contribution due to the Debye layers, Second, there is the internal contribution due to the electric field inside the membrane (cf. Ref. [15]), which is given by Note that both corrections to the surface tension are negative, which means that these corrections can lead to an instability as soon as the total surface tension Σ tot = Σ 0 + ∆Σ κ + ∆Σ m (which is the sum of the bare tension Σ 0 plus the above corrections) becomes negative [25]. Bending modulus. To obtain the correction to the bending modulus we perform, as detailed in Ref. [15] for the DH case, a calculation of the potential at first order in the membrane height. Then, by solving the hydrodynamics problem of the electrolyte around the membrane (in Stokes approximation), one determines the pressure and obtains the total stress tensor. The growth rate of membrane fluctuations, s(k ⊥ ), where k ⊥ is the wave vector in the membrane plane, is then determined by imposing that the discontinuity of the normal-normal component of the total stress tensor at the membrane has to equal the membrane restoring force: Here F H is the standard Helfrich free energy, and Σ 0 and K 0 are the bare surface tension and the bare bending modulus of the membrane, respectively. Expanding the left hand side of Eq. (27) in powers of k ⊥ yields the growth rate of the form (29) Details of the calculations can be found in Appendix A. The surface tension corrections calculated above in Eqs. (25,26) can be recovered by this method, which provides a self-consistency check. For the bending modulus one obtains for the external contribution due to the Debye layers and the internal contribution due to the voltage drop at the membrane, respectively. The field inside the membrane is given by

Discussion
Let us now discuss the nonlinear electrostatic effects on the membrane elastic moduli in the limits of low and high applied voltages. Low voltage regime. In the low voltage limit,V ≪ 1, a solution of Eq. (16) to linear order in c, yields Here ρ(0 + ) is half the charge density at the membrane, corresponding to the quantity called ρ m in Ref. [15] for the DH case. Expanding φ(z) and ρ(z), as well as the corrections to the moduli ∆Σ κ , ∆Σ m , ∆K κ and ∆K m for small c, one exactly recovers all of the results given in Ref. [15]. Specifically, all corrections to the moduli scale quadratically with the external voltage, ∝ V 2 . This is due to the fact that both the potential and the induced charge are proportional to the applied voltage. High voltage regime. In the opposite limit,V ≫ 1 implies c → 1. Introducing α = 1−c and expanding Eq. (16) for small α, one gets − 2λm α +4 ln α+V +λ m −4 ln 2 = 0. For high values ofV one can neglect the last two (constant) terms and gets where W (x) is Lambert's function, i.e. y = W (x) is the solution of ye y = x. As a result, the nonlinear electrostatics strongly effects the corrections to the moduli from the Debye layers. In fact, the external contribution to the surface tension scales as (2+λm) 2 : in the high voltage regime, the external surface tension correction scales linearly with the voltage instead of quadratically. The crossover fromV 2 toV is salt dependent, i.e. depends on the value ofλ m = λ m κ, cf. Fig. 1. In contrast, the external surface tension correction ∆Σ m remains a quadratic function of the applied voltage. Fig. 2 displays both corrections to the surface tension as a function of voltage for low salt (pure water). Fig. 2a) shows that due to the crossover to a linear voltage dependence for the external surface tension correction, rapidly ∆Σ m (the blue curve) wins and ∆Σ κ (red curve) becomes negligible for V 1.5V . Fig. 2b)  that |∆Σ κ | > |∆Σ m | for all voltages (and both proportional to V 2 ), while already above V 0.4V the internal contribution exceeds the external contribution. For the bending modulus corrections, the differences between the DH and PB models are even more striking: As a function of the applied voltage, the external contribution levels off at a constant value In contrast, the internal contribution ∆K m continues to grow quadratically in the high voltage limit. Fig. 3a) and b) show comparisons of the nonlinear PB solutions (solid) and the linear DH solutions (dashed) for both contributions to the bending modulus. Note, however, the different scales: Although the external contribution saturates, see Fig. 3a), this value is larger by two orders of magnitude than the still growing value from the internal contribution, cf. Fig. 3b). In conclusion, for pure water and voltages ≃ 0.2−5V , the total correction to the bending modulus will appear constant within this voltage interval. Only for still higher voltages the quadratic growth due to the internal contribution ∆K m will dominate. Membrane instability. Let us briefly discuss the consequences for the membrane undulation instability. As mentioned above, an instability develops in this system when the total surface tension, the sum of the bare value and the electrostatic corrections, becomes negative. The threshold value for the voltage, V c , has been calculated in Ref. [15] within linearized DH theory and reads (in case of a non-conductive membrane) This curve is shown in Fig. 4 as the blue line. The red line, in contrast, shows the nonlinear PB result given by numerical solution of 0 = Σ tot = Σ 0 + ∆Σ κ + ∆Σ m with Eqs. (25,26) for the surface tension corrections and c(λ m ,V ) calculated from Eq. (16). In the high salt limit, the Debye layers shrink to zero and as a result the external contribution vanishes. For the internal contribution, the inside field is exactly calculated from the DH approach, as can be seen from the behavior of the parameter c in Figure 1a. Thus, both curves merge and attain the limiting value given by Eq. (36), In contrast, for low salt, the linear result overestimates the threshold value, since it underestimates the induced charges at the membrane, as shown in Fig. 1b).
In addition to the threshold voltage for the instability, the most unstable wave number associated to the instability will be affected by the electrostatic nonlinearities of the PB equation. This wave vector corresponds to the maximum of the growth rate given by Eq. (29) with respect to k ⊥ . In the low salt regime, since the system is more unstable in the presence of electrostatic nonlinearities, cf. Fig. 4, in general the wave vector will be increased as compared to the predictions from the DH approach.

Conclusions
We have investigated the nonlinear electrostatic effects of an external DC electric field on a purely capacitive membrane, which is non-conductive for the ions and bears no fixed charges, in an electrolyte. We have calculated in the nonlinear Poisson-Boltzmann regime the corrections to the membrane elastic moduli -both the external ones due to the Debye layers surrounding the membrane and the internal corrections due to the electric field inside the membrane. Strong deviations from the linear Debye-Hückel behavior have been found in the low salt regime at already rather moderate voltages. In particular we have shown that the external contribution to the surface tension crosses over from a quadratic dependence on the externally applied voltage as predicted by the linearized theory to a linear voltage dependence. In contrast, the internal contribution remains quadratic and becomes dominating at high voltages. The external contribution to the bending modulus even saturates for high voltages, while the internal contribution remains quadratic in voltage.
In addition, our work confirms that the surface tension still grows in absolute value with the voltage, which means that the membrane undulation instability present in the DH theory (due to an effectively negative surface tension) persists in the nonlinear PB regime. The nonlinearities affect the threshold in voltage and the characteristic wavelength of the instability.
The method presented here can serve as a starting point for further extensions. One example would be to include fixed charges in addition to induced charges, similarly as studied in [26]. Another possible direction would be to improve the description of the membrane, cf. the recent work in Ref. [27] where thickness fluctuations of the membrane and fluctuations of the lipid dipole orientations within the membrane are accounted for in a comprehensive continuum model for a membrane in a normal DC electric field. Other possible extensions could include other types of non-linear effects, for example due to membrane elasticity, due to inclusions of proteins such as ion channels or pumps in the membrane [28], and also various relevant non-equilibrium effects, coupling electrostatics and hydrodynamics as in induced charge electro-osmosis [29]. Finally, in the biological context, the heterogeneity of the bilayer composition is another feature which is beyond the present model and which is likely to be important.
A little care has to be taken for the correction due to the internal field, cf. Ref. [15]. Note that we accounted for the finite membrane thickness d in Eq. (A.8) above. As the internal field is constant, and due to symmetry, one can write The potential inside the membrane can be written (using again the symmetry, as well as ρ = 0 inside; cf. also Ref. [14] for details) φ m 1 k ⊥ , d 2 can be calculated approximately by using the outside solution, Eq. (A.2), and imposing the BC at the membrane, leading to On the right hand side, now all quantities are known. Finally, one expands Eq. (A.11) and all the other quantities entering the total stress discontinuity, Eq. (A.7), up to fourth order in k ⊥ . From Eq. (A.9) one then obtains the growth rate s(k ⊥ ) of membrane fluctuations.