# Dielectric permittivity of water confined in stacks of charged lipid layers: Extracting profiles from molecular dynamics simulations using a modified Poisson-Boltzmann equation

Ludovic Gardré, Swen Helstroffer, Pierre Muller, Fabrice Thalmann, Thierry

Charitat, Laurent Joly, Claire Loison

► To cite this version:

Ludovic Gardré, Swen Helstroffer, Pierre Muller, Fabrice Thalmann, Thierry Charitat, et al. Dielectric permittivity of water confined in stacks of charged lipid layers: Extracting profiles from molecular dynamics simulations using a modified Poisson-Boltzmann equation. The Journal of Chemical Physics, 2025, 163 (19), 10.1063/5.0274252. hal-05371265

HAL Id: hal-05371265

https://hal.science/hal-05371265v1

Submitted on 18 Nov 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>
# Dielectric permittivity of water confined in stacks of charged lipid layers: extracting profiles from molecular dynamics simulations using a modified Poisson-Boltzmann equation

Ludovic Gardré,¹ Swen Helstroffer,² Pierre Muller,² Fabrice Thalmann,² Thierry Charitat,² Laurent Joly,¹ and Claire Loison¹, a)

¹) Université Claude Bernard Lyon 1, CNRS, Institut Lumière Matière, UMR5306, Villeurbanne, France

²) UPR 22/CNRS, Institut Charles Sadron, Université de Strasbourg, 23 rue du Loess, BP 84047 67034 Strasbourg Cedex 2, France

**Abstract** Most organic and inorganic surfaces (e.g., glass, or lipid membranes) become charged in aqueous solutions. The resulting ionic distribution induces effective interactions between the charged surfaces, which depend on the permittivity of the confined aqueous solution. To explore this phenomenon in very confined systems, we perform all-atom molecular dynamics (MD) simulations of charged lipid membranes separated by a salted water layer of varying thickness. To extract an effective permittivity from the atomistic model, we compare the ion distributions of these simulations with the ones of a continuous, mesoscopic model: a modified Poisson-Boltzmann (mPB) equation taking into account a spatially varying dielectric permittivity and an explicit Born solvation energy for ions. Such mPB/MD comparisons, applied to lipid membranes at various hydration levels, reveal a sharp decrease of the permittivity upon dehydration, converging to a plateau value that we attribute to lipid headgroups. We discuss the limitations of the mPB model in the dehydrated lipid membranes, in particular through the appearance of ion/ion correlations, and compare our results with alternative computational methods. In our tension-free simulations of the fluid membranes, an increase of the area per lipid indicates that the permittivity decrease is accompanied by an intermembrane attraction. Our approach could be generalized to estimate a mesoscopic permittivity of liquids confined by other interfaces, provided ions follow Boltzmann statistics.

**Keywords** Confined water, dielectric permittivity, lipid layers, molecular dynamics

## I. INTRODUCTION

Water is a crucial component of the human body, constituting over half of its composition¹. In this complex biological system, water exists in confined spaces such as blood capillaries, within cells, or trapped between lipid membranes where it has been extensively studied². Confined water exhibits properties significantly different from those in the bulk³. The dielectric permittivity in particular is one of these. This macroscopic property of a continuous medium describes how it becomes polarized under an external electric field. The polarization arises from the complex microscopic polarization mechanisms within the medium and generally exhibits non-local behaviour. In bulk water, its value tends to be high due to the polar nature of the molecule. However, in a confined geometry, the permittivity becomes inhomogeneous and anisotropic. Relatively to a planar confining wall, a decrease of the normal component of the water dielectric permittivity has been observed experimentally⁴, and in models at various length scales⁵⁻⁷. Additionally, the observed attraction between like-charged lipid layers in aqueous monovalent electrolytes has been interpreted as the signature of a pronounced dielectric decrement of the confined water⁸⁻¹⁰.

Though it is experimentally very challenging⁴,¹¹,¹² to probe the dielectric permittivity in these regions, numerical simulations can give insights. Seminal works by Stern and Feller on zwitterionic lipid membranes¹³, or by Balenegger and Hansen on water models¹⁴, have proposed linear-response formalism to obtain the permittivity tensor from the polarization density, taking into account the dielectric boundary conditions in the considered geometries. In high-dielectric regions, obtaining a quantitative estimate of the component normal to the interface turns out to be difficult because of the numerical statistical uncertainties on the polarization density fluctuations are dramatically amplified on the permittivity measurement. More recently, Schaaf and Gekle¹⁵, and Netz and co-workers⁵,⁸,¹⁶ have successfully computed both the parallel and perpendicular components of the dielectric permittivity in water confined by various interfaces. Unfortunately, the fluctuation-dissipation formula cannot be directly applied to study confined water in systems containing ions, such as charged lipid membranes, which are very abundant in biology.

In this paper, we present an alternative approach to extract an effective dielectric permittivity profile of water confined between stacks of charged lipid membranes, i.e. in the presence of free charges. Our method uses molecular dynamics (MD) simulations and is based on a mean-field approximation with the use of a modified Poisson-Boltzmann (mPB) equation. These modifications com-

<footnote>a) Corresponding author: claire.loison@univ-lyon1.fr</footnote>
bine an effective local z-dependant dielectric permittivity¹⁷ with the addition of a solvation energy term in the Boltzmann distribution of the ions¹⁸. Indeed, because of the differences between lipid-ions and water-ions solvation, the ions are expelled from the lipid region. This has to be taken into account in the description of the potential energy of the ions.⁸

In a first part, we present the systems we studied, the details of the mPB model and the methods we used to calibrate and solve it. We then apply the model to investigate the value of the dielectric permittivity of water confined in stacks of charged lipid layers, in the presence of a monovalent salt. We present its results for a large range of hydration level of the lipid membranes, which allow to extract a dielectric permittivity at the centre of the confined water, at each hydration value. We then compare our approach to alternative ways of estimating the dielectric permittivity in such systems.

## II. METHODS

### A. Studied systems and molecular dynamics simulations

We have followed the standard setup for planar bilayers: a tensionless membrane patch with 3-dimensional periodic boundary conditions and the z-axis along the bilayer normal. Figs. 1(a) and (b) show snapshots of example systems at high hydration and at two different temperatures.

**Composition** — All systems were built from the same base elements: a lipid bilayer made of 1,2-dipalmitoyl-sn-glycero-3-phosphoserine (DPPS, negatively charged in the head) solvated in water, see Fig. 1. We neutralized the charge of the lipids with $Na^+$ counterions. In addition, to test the model in presence of both cations and anions, we included NaCl salt in the solvent, at a concentration of approximately 0.25 M. The amount of water plays an important role in these systems. It was controlled at the creation step by setting a hydration number (HN), corresponding to the number of water molecules per lipid (e.g., HN = 20 means there are 20 water molecules per lipid). For HN ≥ 11, the bilayer was made of two leaflets of 100 lipids. At lower HN, to keep the salt concentration constant while still having a statistically relevant number of salt ions in the system, we quadrupled the systems size so each leaflet of the bilayer contained 400 lipids. The compositions of the different systems are detailed in Tab. S3 in Supplementary Material.

**Construction of initial conformations** — For HN ≥ 11, the systems were built using charmm-gui¹⁹. This tool provided pre-equilibrated systems with bilayers in the fluid phase. We followed the charmm-gui default pre-equilibration process that consists of six runs: two 125-ps runs in the (N,V,T) ensemble, followed by one 125-ps run and three 500-ps runs in the (N,P,T) ensemble. Restraints are applied on the lipid atoms positions and

on the dihedrals inside the lipid molecules. These restraints are gradually removed run after run throughout the whole equilibration process. For HN ≤ 10, systems were built iteratively: for each HN, we started by replicating the system at the HN just above (e.g., we replicated the HN = 11 system to create the HN = 10 system, and so on down to HN = 02). We then randomly removed water molecules and/or ions to match the desired HN and salt concentration.

**Force field** — We chose CHARMM36²⁰ for the lipid atoms and the modified TIP3P (mTIP3P)²¹ for water, since CHARMM36 was originally developed with this water force field. It is known that TIP3P permittivity does not match the experimental value for bulk water, but the typical decrease of permittivity when temperature is increased is reproduced²². We performed independent simulations on bulk mTIP3P water, either pure or with 0.25 M NaCl that confirm this result for the mTIP3P model (see Sect. S-IV in Supplementary Material). Thus, in this work, we focus on the relative variations of the permittivity, not the absolute values. We did not use explicit polarisable force fields²³, nor scaled-charge force fields²⁴. These may be considered if one would extend this study²⁵. Moreover, one has to remember that the impact of dehydration on the lipid bilayer structure is not yet perfectly reproduced by CHARMM36²⁶.

Short-ranged van der Waals interactions were modelled with a Lennard-Jones potential. A cutoff distance was set at 12 Å, and we smoothly switched forces to 0 between 10 and 12 Å. We used Lorentz-Berthelot combination rules to compute the cross-species parameters. Long ranged Coulombic interactions were computed with the Particle-Mesh Ewald method²⁷,²⁸. All simulations were carried out with a time step of 2 fs, using Gromacs 2021.2 with CUDA support²⁹.

**Simulation Ensemble** — The equilibration and production runs were performed in the (N, $P_n$, $\gamma$, T) ensemble, were $P_n$ is the normal pressure and $\gamma$ is the surface tension. The temperature was maintained with a Nosé-Hoover thermostat³⁰. Systems were kept in the $P_n = 1$ bar and $\gamma = 0$ ensemble with the semi-isotropic Parrinello-Rahman barostat³¹, i.e. with the constraints $L_x = L_y$, and $P_z = (P_x + P_y)/2$.

**Equilibration and Production** — Every system provided by charmm-gui was first run at T = 333.15 K during 1 ts. Then, we increased the temperature up to T = 353.15 K with a step of 5 K. For every temperature, we ran a 1 ts simulation before increasing the temperature.

Even if all initial configurations were initially prepared in the fluid phase, depending on the HN and temperature, we observed spontaneous phase transitions from fluid to gel (or vice versa). When it occurred, we extended the simulation length and discarded everything before the transition. Therefore, all observables shown in this work were computed over 1 ts phase-transition-free production runs. The error bars in Figs. 3, 6 and 7 were obtained
This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.
PLEASE CITE THIS ARTICLE AS DOI: 10.1063/5.0274252

<figure><img src="image_2.png" /><figcaption>FIG. 1. (a, b) Snapshots of systems at HN = 35, respectively at T = 353.15 K (bilayer in the fluid phase) and T = 333.15 K (bilayer in the gel phase). The black boxes represent the simulation boxes, replicated in all directions of space. For clarity, the lipids' hydrogen atoms are not drawn in the snapshots, nor are water molecules in the simulation boxes. The z coordinate is set at 0 at the centre of the left bilayer. It reaches zmax = Lz/2 at the middle of the interlayer water. (c) Detail of the DPPS lipid molecule. The gold and silver colours respectively correspond to the heads and the tails of the lipids as drawn in the snapshots. (d) Interlayer species.</figcaption></figure>

by dividing the production run into five trajectories of 200 ns, which were analysed separately and considered as independent. The standard error of the mean was multiplied by the appropriate Student's t-factor for a 95% confidence interval (i.e. ≈ 2.78 for five values).

*Lipid ordering* — Depending on its hydration and temperature, the equilibrated bilayer can be in the fluid or in the gel phase. At the highest temperature, the bilayer is in the fluid phase and consequently the lipids are disordered, see Fig. 1(a). At the lowest temperature, the bilayer is in the gel phase: the lipid tails are ordered, but not their heads, see Fig. 1(b).

Due to the structural difference between these two phases, the membrane has a different thickness which influences the box length in the *z* direction, meaning that *L*<sub>*z*</sub> is not the same in all systems and varies depending on the physico-chemical conditions (e.g., *T*, HN, etc.).

Even if the orientation of the lipid tilts may be differ-
ent in the two monolayers in the gel phase, see Fig. 1(b),

all the systems are considered symmetrical with respect
to the plane defined by z = zmax. In practice, when com-
puting charge densities from MD simulations, the centre
of mass of the bilayer is placed at the centre of the sim-
ulation box (in Fig. 1, this is equivalent to shifting the
simulation box by Lz/2 in the z direction). The densities
are then integrated along the z axis, then symmetrized
and centred with respect to the bilayer centre of mass.

B. Modified Poisson-Boltzmann model

The purpose of this work is to establish a method to
obtain, from MD simulations, an effective profile of the
dielectric permittivity of the ion-carrying water confined
between charged lipid membranes. To do so, we propose
to modify the Poisson-Boltzmann equation, a mean-field
approximation that describes how electrolyte solutions
behave in the vicinity of a charged interface³².

First, let us recall the standard Poisson-Boltzmann model, starting from Maxwell-Gauss equation in matter:

$$ \nabla \cdot D = \rho_{\text{free}}, \qquad (1) $$

where $D$ is the electric displacement field and $\rho_{\text{free}}$ is the free-charge density.

Considering a planar geometry (with the $z$ axis perpendicular to the charged plane), the problem is reduced to one dimension as a consequence of the systems symmetries and invariances:

$$ -\varepsilon_0 \varepsilon_{\perp} \frac{d^2 V}{dz^2} = \rho_{\text{free}}(z) = \rho_{\text{ions}}(z), \qquad (2) $$

where $V$ is the electric potential, $\varepsilon_0$ the vacuum permittivity, $\varepsilon_{\perp}$ is the perpendicular component of the relative permittivity (assumed homogeneous) given by the relation $D_{\perp} = \varepsilon_0 \varepsilon_{\perp} E_{\perp}$, and considering the total free-charge density is equal to the charge density of the ions $\rho_{\text{ions}} = \rho_{+} + \rho_{-}$, where $\rho_{+}$ and $\rho_{-}$ are the charge densities of the cations and anions, respectively. Considering a monovalent salt, and using a Boltzmann distribution of cations and anions:

$$ \rho_{\pm}(z) = \pm eC_{\text{salt}} \exp[\mp\beta eV(z)], \qquad (3) $$

with $\beta = (k_B T)^{-1}$ ($k_B$ the Boltzmann constant and $T$ the temperature) and $C_{\text{salt}}$ the ions concentration in a reservoir where $V$ vanishes, one can express $\rho_{\text{ions}}$ as follows:

$$ \begin{aligned} \rho_{\text{ions}}(z) &= \rho_{+}(z) + \rho_{-}(z) \\ &= eC_{\text{salt}} \left\{ \exp[-\beta eV(z)] - \exp[\beta eV(z)] \right\}. \end{aligned} \qquad (4) $$

Combining Eq. (2) and Eq. (4), one derives the standard 1D Poisson-Boltzmann equation:

$$ \frac{d^2 V}{dz^2} = \frac{2eC_{\text{salt}}}{\varepsilon_0 \varepsilon_{\perp}} \sinh[\beta eV(z)]. \qquad (5) $$

To simplify Eq. 5, one can introduce the dimensionless reduced electric potential $\phi$, defined as $\phi(z) = \beta eV(z)$, and two characteristic lengths: first, the Bjerrum length $l_B$ (the characteristic length at which the thermal energy equals the electrostatic interaction energy between two monovalent ions), given by³³:

$$ l_B = \frac{\beta e^2}{4\pi\varepsilon_0\varepsilon_{\perp}}; \qquad (6) $$

second, the Debye length $\lambda_D$ (the characteristic length at which the electrolyte solution screens the electric field)³³:

$$ \lambda_D = \frac{1}{\sqrt{8\pi l_B C_{\text{salt}}}}. \qquad (7) $$

One can then rewrite Eq. (5) as:

$$ \frac{d^2 \phi}{dz^2} = \frac{1}{\lambda_D^2} \sinh[\phi(z)]. \qquad (8) $$

In our modified Poisson-Boltzmann (mPB) model, we allowed the relative permittivity $\varepsilon_{\perp}$ to vary along the $z$ axis¹⁷. In addition, we included the lipid heads charge density $\rho_{\text{lipids}}$ as a fixed external contribution to the total free-charge density. We then rewrote Eq. (2) as follows:

$$ \begin{aligned} -\varepsilon_0 \frac{d}{dz} \left[ \varepsilon_{\perp}(z) \frac{dV}{dz} \right] &= \rho_{\text{free}}(z) \\ &= \rho_{\text{lipids}}(z) + \rho_{\text{ions}}(z). \end{aligned} \qquad (9) $$

Since $\varepsilon_{\perp}$ varies with $z$, so does the ions solvation energy. We then introduced the solvation energy $W_{\perp}$ in the Boltzmann factor, using the Born expression¹⁸:

$$ W_{\pm}(z) = \frac{e^2}{8\pi\varepsilon_0 r_{\pm}} \left( \frac{1}{\varepsilon_{\perp}(z)} - \frac{1}{\varepsilon_{\perp}^{\text{ref}}} \right), \qquad (10) $$

where $r_{\pm}$ is the radius of the ion and $\varepsilon_{\perp}^{\text{ref}}$ is a reference value for the relative permittivity of water (e.g., bulk water value). This potential describes the energetic cost to create an electric field in the surrounding dielectric medium due to the presence of the ion. It represents the energy for the dipoles in the solvent to reorient because of the presence of the charge. The increase of this energy when the ion enters the lipid region describes that the lipid heads and tails are more difficult to polarize than the water dipoles. This Born term therefore represents how the lipid-ion solvation differs from the water-ion solvation. Noticeably, it approximates solely dipole/ion electrostatic interactions. The cost to create a cavity for the ion is not included in this Born term. This model therefore considers that the free energy associated to the Lennard-Jones interactions of the ions remain approximately constant when the ions move relative to the lipidic membrane. The values of the parameters of the Born model will be discussed in section IID, devoted to the model calibration. In the end, the ions charge density became:

$$ \begin{aligned} \rho_{\text{ions}}(z) &= \rho_{+}(z) + \rho_{-}(z) \\ &= eC_{\text{salt}} \\ & \times \left\{ \exp[-\beta W_{+}(z) - \beta eV(z)] - \exp[-\beta W_{-}(z) + \beta eV(z)] \right\}. \end{aligned} \qquad (11) $$

### C. Solving the mPB model's equation

In this section, we introduce the method we followed to solve the mPB model shown in Eq. (9). In this equation, the ions charge density $\rho_{\text{ions}}$ is given by the analytical expression from Eq. (11), while the lipid charge
density $\rho_{lipids}$ corresponds to values extracted from MD simulations beforehand. These values were centred and symmetrized around the centre of mass of the bilayer, allowing to focus on one single half of the system when solving the equation. Therefore, we set the origin of the z axis at the centre of the bilayer and the maximum z value then corresponded to the centre of the interlayer water channel.

There were two unknown profiles in the equation: the electric potential $V(z)$ and the relative permittivity $\varepsilon_{\perp}(z)$. We chose to model the latter with a parametrized curve. This led to having only one fully unknown profile, the electric potential, so we could numerically solve the equation.

To decide on the shape of the dielectric permittivity profile, we analysed a system made of a zwitterionic lipid bilayer and no ions, i.e. without any free charges. In this situation, the electric displacement field is constant in the solvent, and the dielectric response can be obtained as the linear response of the polarization density upon the perturbation of an electric field. Fluctuation-dissipation approaches can then be used to extract a dielectric permittivity profile from the local dipole moment fluctuations in equilibrium simulations⁸. We used MAICoS (https://www.maicos-analysis.org), which provided us a profile of $\varepsilon_{\perp}^{-1}(z)$ that was close to a sigmoid (see Fig. S1). This shape was similar to the inverse permittivity profiles obtained by Schlaich et al. for neutral lipid bilayers such as Digalactosyldiacylglycerol (DGDG)³⁴. Supposing that the shape of the dielectric permittivity profile would be similar in a lipidic system with free charges, we chose to model its inverse (to stay consistent with the reference provided by MAICoS) using one sigmoid and one Gaussian, to allow for more degrees of freedom:

$$ \begin{aligned} & \varepsilon_{\perp}^{-1} (z, z_1, K, \sigma, \mu, C, \varepsilon_{\perp}^{\text{tails}}, \varepsilon_{\perp}^{\text{plateau}}) \\ & = \frac{1}{\varepsilon_{\perp}^{\text{plateau}}} + \left( \frac{1}{\varepsilon_{\perp}^{\text{tails}}} - \frac{1}{\varepsilon_{\perp}^{\text{plateau}}} \right) \\ & \qquad \times \frac{1}{1 + \exp[K(z - z_1)]} \\ & \qquad + \frac{C}{\sigma\sqrt{2\pi}} \exp\left[-\frac{(z - \mu)^2}{2\sigma^2}\right]. \end{aligned} \quad (12) $$

$z_1$ and $K$ are parameters for the sigmoid component: the first one controls the coordinate of its inflection point while the other one controls its stiffness. $\sigma$, $\mu$ and $C$ are parameters for the Gaussian component. The first two respectively correspond to its standard deviation and its expected value, while $C$ controls its height. Finally, $\varepsilon_{\perp}^{\text{tails}}$ and $\varepsilon_{\perp}^{\text{plateau}}$ respectively control the left and the right horizontal asymptotes. An example of this curve is shown in Fig. S2 in the supplementary material. The effect of the addition of a second Gaussian function in Eq. 12 is

also displayed in Fig. S3. We considered the results with a single Gaussian more robust.

Equation 12 reflects our choice for a profile without large or numerous oscillations. It cannot describe complex profiles with several oscillations or negative values, as the ones that have been calculated for flat interfaces such as graphene³⁴. Indeed, we expect that the roughness of the lipidic interface diminishes strong density oscillations and associated permittivity oscillations. The Born energy also describes a solvation cost that is averaged over slices parallel to the interface. Therefore, our definition for $\varepsilon_{\perp}^{-1}$ is thought as a "coarse-grained dielectric constant". It may be linked to the one proposed by Borgis et al.⁷. This quantity, noted $\tilde{\varepsilon}_{\perp}(z)$ in their article, is obtained by averaging the polarization density over slices of the thickness of the same size as a water molecule. Borgis et al. show that the z-averaging leads to smooth profiles in $\tilde{\varepsilon}_{\perp}(z)$, that remains strictly positive.

Since we parametrized the inverse of the dielectric permittivity, we modified Eq. (9) to obtain an equation depending on $\varepsilon_{\perp}^{-1}$ and its derivative with respect to $z$. The complete system of equations numerically solved is given in section S-I in the supplementary material.

In this solving process, we added an optimization procedure, to obtain the best values for the parameters of $\varepsilon_{\perp}^{-1}$ (Eq. (12)). To do so, we fitted the ionic charge densities computed with Eq. (11) on the ionic charge densities extracted from the MD simulation.

### D. Calibrating the mPB model

The calibration process served the purpose of fixing the value of two parameters: $C_{\text{salt}}$ and $\varepsilon_{\perp}^{\text{ref}}$.

Three parameters were set prior to the calibration: the ionic radii of Na⁺ and Cl⁻ and the asymptotic value for the relative permittivity in the lipid tails $\varepsilon_{\perp}^{\text{tails}} = 2$. This later choice should be discussed. Values obtained from MD simulations are usually reported around 1¹³,³⁵, in agreement with the value we obtain in the analysis of the DPPC bilayer with MAICoS, as visible in Fig. S1. However, theoretical calculation³⁶ and experimental measurements¹¹ report a value closer to 2. Fig. S5 shows that this choice has practically no impact on the permittivity profile in the region where ions are present.

The Born radii are in fact effective parameters, and may differ from the commonly measured ionic crystal radius³⁷⁻³⁹. We have used the Born radii of these ions calculated from their free energies of hydration by Schmid and al.⁴⁰: $r_+ = 0.187$ nm and $r_- = 0.186$ nm, after checking that the values chosen had no significant impact on the results.

We calibrated the model for each of the five temperatures. Every calibration was realised in two stages, using the HN = 45 system which contained enough water to have a solvent reservoir between the membranes. The
first stage was to determine the ions concentration in the reservoir $C_{\text{salt}}$. Using the ionic Boltzmann distributions as in Eq. (11), one can get rid of the electrostatic energy by multiplying the cations and anions concentrations (respectively $C_+$ and $C_-$):

$$C_+(z) \times C_-(z) = C_{\text{salt}} \exp[-\beta eV(z) - \beta W_+(z)] \times C_{\text{salt}} \exp[\beta eV(z) - \beta W_-(z)] = C_{\text{salt}}^2 \exp[-\beta(W_+(z) + W_-(z))]. \quad (13)$$

In the reservoir, the ions concentration $C_{\text{salt}}$ is constant. Moreover, one can also set $W_{\pm} = 0$ in the reservoir by defining the relative permittivity of the reservoir water as the reference value $\epsilon_{\perp}^{\text{ref}}$ (Eq. (10)), which will be obtained in the second stage of the calibration. In the reservoir, Eq. (13) then gives:

$$C_{\text{salt}} = \sqrt{C_{+} \times C_{-}}. \quad (14)$$

Fig. S4 in the supplementary material shows the ionic concentrations $C_+$ and $C_-$ in the most hydrated system at $T = 353.15$ K, and the square root of their product that indeed reaches a constant value in the reservoir. The value of $C_{\text{salt}}$ was then obtained using Eq. 14 in the reservoir, i.e. at large $z$ values. The values $C_{\text{salt}}$ are reported Tab. S2. They are all close to $0.20 \pm 0.01$ M.

Once this first step was done, we solved the mPB model's equation by following the method described in the precedent section IIC, except we added the following constraint during the optimization process: $\epsilon_{\perp}^{\text{plateau}} = \epsilon_{\perp}^{\text{ref}}$. The dielectric permittivity indeed reached a plateau in the water region (where $z$ values are close to $z_{\text{max}}$, as shown in Fig. 1). We identified this plateau value, which corresponds to the horizontal asymptote defined by the parameter $\epsilon_{\perp}^{\text{plateau}}$ in Eq. (12), as the dielectric permittivity of the reservoir water. Since this parameter is optimized, the constraint allowed to maintain $W_{\pm} = 0$ in the reservoir, ensuring self-consistency with the first stage of the calibration procedure.

Finally, we kept the optimized value of $\epsilon_{\perp}^{\text{plateau}}$ as the value of $\epsilon_{\perp}^{\text{ref}}$ in all the systems at the same temperature. The values $\epsilon_{\perp}^{\text{ref}}(T)$ for all temperatures are reported in Supplementary Material, Tab. S2. One expects these values of $\epsilon_{\perp}^{\text{ref}}(T)$ to approach the values of bulk mTIP3 at the relevant temperature and $C_{\text{salt}}$ salt concentration. The permittivity of bulk TIP3P water is known to decrease with increasing temperature$^{22}$. We observe the same tendency for our values of $\epsilon_{\perp}^{\text{ref}}(T)$. The permittivity of bulk TIP3P water is also known to decrease with the addition of salt$^{41}$. We also observed the same tendency for our values of $\epsilon_{\perp}^{\text{ref}}$, that are about 10 to 15% lower than the values for pure mTIP3P at the given temperature, measured in independent simulations (see

Sect. S-IV). We conclude that our reference permittivities for the reservoirs are reasonable.

### III. RESULTS

#### A. Solutions of the mPB model and corresponding profiles

After the calibration was done, we were able to solve the mPB model's equation (Eq. (S4) in the supplementary material). For each system, we obtained three profiles: the electric potential $V(z)$, the dielectric permittivity $\epsilon_{\perp}^{\text{mPB}}(z)$ and the ionic charge densities $\rho_{\pm}(z)$ predicted by the model (computed using the two latter and Eq. (11)). Examples of these results for a DPPS bilayer at 353.15 K are shown in Fig. 2, both at high hydration (HN = 35, left side) and low hydration (HN = 10, right side). As illustrated in Fig. 1, thanks to the symmetry of the system, all plots are simplified to show only one half of the system. The lipid tails are located in the left side, the lipid heads are in the region where $\rho_{\text{lipids}}$ is not null and the water is on the right side (with an overlap in the hydrophilic heads).

Fig. 2 (a, d) show the estimated relative dielectric permittivity profiles $\epsilon_{\perp}^{\text{mPB}}$ in the systems, obtained by inverting Eq. (12). The general shape of the profiles naturally depends on the parametrized curve we chose to describe its inverse: one sigmoid and one Gaussian. In particular, thanks to the Gaussian component, the dielectric permittivity can display a small bump. This can notably be the case in the lipid heads region, where the structure is more complex because of the presence of free charges (ions) and the complex distribution of charges and dipoles in the lipid heads. An example solution of the mPB model without the Gaussian component is given in Fig. S2 in the supplementary material to demonstrate its importance.

At high hydration number, the dielectric permittivity reaches a plateau in the confined interlayer water. This confirms that the bilayer is highly hydrated, therefore at the centre of the channel the water molecules are not constrained by the interactions with the lipids. With further dehydration however, this plateau is never reached, as seen in Fig. 2 (d). In this hydration state, all the remaining water molecules are located in the lipid heads region, where their orientation is somewhat constrained by electrostatic interactions, which is why a plateau is never reached.

Since the system is considered symmetrical with respect to the plane defined by $z = 0$ (other half of the bilayer, as visible in Fig. 1) and periodic through the MD periodic boundary conditions, we imposed when solving the model's equation that the derivative of the reduced electric potential is null at the borders to ensure its continuity, as visible in Fig. 2 (b, e).

For highly hydrated systems, we found a magnitude of about 100 mV for the electric potential difference $\Delta V$ be-
<figure><img src="image_3.png" /><figcaption>FIG. 2. Results from the modified Poisson Boltzmann model at 353.15 K, at HN = 35 (left) and HN = 10 (right). (a, d) Relative dielectric permittivity profile. The lipid in the background represents the approximative z position of the tails and the heads in the lipid layer. (b, e) Reduced electric potential φ(z) and electric potential V(z). (c, f) Charge densities. Markers show charge densities extracted from the MD simulation (for clarity, only one data point out of seven is plotted). Lines show charge densities obtained with the mPB model, i.e. with Eq. (11). According to Eq. (9), both blue plots include the lipid charge density extracted from the MD simulation.</figcaption></figure>

tween the lipid heads and the water, which corresponds to Δφ ≈ 3.3. Though we did not manage to find an other value for a DPPS membrane in the literature to compare with the value given by the model, we found some values for a 1,2-dipalmitoyl-sn-glycero-3-phosphocholine (DPPC) membrane. This lipid has a different head than that of DPPS and is zwitterionic instead of being negatively charged. Stern and Feller carried out simulations of DPPC membranes¹³ and found a magnitude of around 1 V (Δφ ≈ 33), while experimental values⁴² gave a smaller magnitude of 227 mV (Δφ ≈ 7.5). One of the reasons the authors gave to explain such difference was the inaccuracy of the force field they used, which was CHARMM27⁴³. Despite being a different lipid, we expected the order of magnitude of DPPS membrane's electric potential difference to be somewhat close from that of DPPC membrane's. To strengthen this hypothesis, one can use Grahame equation³², which derives from the standard Poisson-Boltzmann model and links the surface charge to the surface reduced potential difference, in the case of a planar solid wall. For a monovalent salt, it

writes:

$$ \sinh\left(\frac{\Delta\phi}{2}\right) = 2\pi\ell_B\lambda_D\frac{\sigma_s}{e}, \quad (15) $$

where $\ell_B$ is the Bjerrum length defined by Eq. 6 and $\lambda_D$ is the Debye length defined by Eq. 7. Applied in our systems, where $\sigma_s$ is on the order of 1 e per 0.5 nm², taking $\varepsilon_⊥ \simeq 60$ and $T = 353$ K, Eq. 15 gives a value $\Delta\phi \approx 5$, which is consistent with the value we measured. Since we used a much more recent version of the CHARMM force field than Stern and al., our results tend to validate their statement against the accuracy of CHARMM27.

With the profiles obtained in Fig. 2 (a) and (b), the model computes the ionic charge densities using Eq. (11). The results are shown in Fig. 2 (c) and compared with the charge densities obtained from the MD simulations. For both anions and cations, the charge densities computed by the model are in good agreement with MD data. In addition, since the solution for the electric potential appears to be consistent within the PB model, these results suggest that the profile of the dielectric permittivity $\varepsilon_⊥^{mPB}$ is a good approximation of the permittivity in the system.
In comparison to the classical PB model, the ions are more repelled by the interface. In our mPB, this repulsion is due to the solvation energy, via the Born energy term. Fig. S6 shows that the Born energy dominates over the electrostatic potential when entering the lipid region. In other mPB formulations, this repulsion was also taken into account by other forms of potentials$^{8,16}$, but they also dominate over the electrostatic energy nearby the interface.

Similarly, with the profiles obtained in Fig. 2 (d) and (e), the model computes the ionic charge densities shown in Fig. 2 (f). This result shows that even in a low hydration state, the charge densities computed by the model stay in good agreement with MD data. Moreover, the model is also able to reproduce more complex shapes of charge densities as they appear in systems where the bilayer is in the gel phase. Fig. S2 (a) in the supplementary material shows results from the mPB model on a system at HN = 35 and T = 333.15 K, which corresponds to a bilayer in the gel phase.

## B. Dielectric permittivity at the centre of the interlayer water channel

Once the mPB model has been solved in all systems, at all temperatures, one can extract the relative dielectric permittivities at the centre of the water channels by evaluating $\epsilon_{\perp}^{\text{mPB}}(z_{\text{max}})$, where $z_{\text{max}} = L_z/2$ as shown in Fig. 1. The resulting values are plotted in Fig. 3 for bilayers in the fluid phase and in the gel phase.

The value at the centre of the water layer corresponds to the maximum of the dielectric permittivity profile. If one wanted to define an effective permittivity of the whole water slit, one would typically take the average of its inverse$^7$ $\epsilon_{\perp}^{-1}$. Averaging over $\epsilon_{\perp}^{-1}$ results in an effective permittivity that is dominated by the lower values of $\epsilon_{\perp}$. Therefore, a typical effective permittivity of the whole water slit would be lower than $\epsilon_{\perp}^{\text{mPB}}(z_{\text{max}})$. The value of such an effective average of the dielectric permittivity strongly depends on the definition of its thickness and position, which are not easily defined in the case of lipid bilayers. We therefore resorted to investigate the maximum of the permittivity profiles.

In the gel phase, the dielectric permittivity is rather constant upon dehydration while staying at high HN. However, it starts to raise at some point, which is counter-intuitive. The same phenomenon appears for the fluid phase, however in this case the raise seems to be continuous upon all the dehydration process. At the time we write this paper, we do not have a clear explanation for this increase. On the one hand, it could be because the real dielectric permittivity is actually a non-local property, but it is described as a local one in the model. In other words, $\epsilon_{\perp}^{\text{mPB}}$ is an effective dielectric permittivity profile, optimized to describe the ionic charge densities, but not the real profile itself. It depends on the

parametrized curve we chose to model it. Therefore, it is possible that this curve does not allow to describe the dielectric permittivity profile with enough accuracy and that the model compensates with odd values of the dielectric permittivity. Noticeably, in the study by Borgis et al (Stockmayer fluid confined by Lennard-Jones planar walls$^7$), this quantity also reaches a maximum well above the bulk value of the dielectric constant nearby the confining interface.

On the other hand, this increase could come from the mPB model itself. We brought two modifications to the standard PB model: we assumed a varying dielectric permittivity along the *z* axis and added the Born energy in the Boltzmann distribution of the ions. The model may lack some other interactions in this complex environment.

While the model gives high values of the dielectric permittivity at large hydration, we observe a strong decrease of the dielectric permittivity, down to a plateau value at HN = 11 which is around half its value when the attraction arises. It is generally assumed that HN = 10 is a threshold value, corresponding to the water molecules bound to the lipid heads$^{44}$. For this reason, the value of $\epsilon_{\perp}^{\text{mPB}}$ that remains close to 50 for the dehydrated systems (HN < 12) can be interpreted as the dielectric permittivity in lipid head region, that becomes almost independent

<figure><img src="image_4.png" /><figcaption>FIG. 3. Maximum dielectric permittivity versus hydration number at different temperatures, for bilayers in the fluid phase (top) and bilayers in the gel phase (bottom). Filled markers represent bilayers in the fluid phase and open markers represent bilayers in the gel phase. The purple region shows where we identify an attractive regime between the membranes (see Sect. IV D). The mPB model cannot be used in the hatched region because the ions no longer follow a Boltzmann distribution.</figcaption></figure>
of the water content.

The plot of the maximum dielectric permittivity as a function of the interlayer water thickness is shown in Fig. S10 in the supplementary material.

Under HN = 8, the mPB model fails due to a lack of statistics. Because of the low number of anions and since the dynamics in the water channel between the membranes slows down a lot when reaching this dehydration state (even at high T), the charge density profile of the anions no longer follows a Boltzmann distribution.

## IV. DISCUSSION

### A. Importance of ion/ion correlations in confined electrolytes between charged membranes

Since this work is based on a Boltzmann description of the ion distribution, with energies that depend solely on the z component, the ion/ion electrostatic interactions are taken into account at a mean-field level. But, depending on the nature of the ions, their number or the permittivity, position correlations between ions may appear. To determine whether these correlations must be taken into account, one can compute the coupling constant Ξ of the ions confined between membranes with a given surface charge density σs as follows$^{45}$:

$$ \Xi = 2\pi q^3 \ell_B^2 \frac{\sigma_s}{e}, \quad (16) $$

where e is the elementary charge, q is the valence of the ions, and $\ell_B$ is the Bjerrum length (Eq. 6).

Low Ξ values correspond to a weak coupling regime where mean-field theories, like the Poisson-Boltzmann model, are adequate. However, these models tend to fail in a strong coupling (SC) regime, i.e. for high Ξ values. Several experimental observations have highlighted the appearance of ion/ion correlations in confined electrolytes between charged membranes. In particular, SC regime can predict the attraction between stacks of identically charged lipid layers for Ξ values higher than about 10. For example, in the system of Ref. 46, which consists of charged lipid vesicles with Ca$^{2+}$ counterions, Ξ was estimated to be on the order of 20 due to the large valence of the ions (assuming a water dielectric constant of 80). A similar attraction has been also experimentally observed by Mukhina and al.$^9$ and Helstroffer et al.$^{10}$, this time in the presence of monovalent counterions. The authors interpreted this attractive behaviour as being due to a reduced dielectric constant in the confined water; they found that values between 10 and 30 gave a coupling constant Ξ high enough to justify an attraction between the membrane in their system, within the SC theory. Therefore, it is important to investigate the appearance of ion/ion correlation effects in our simulations. Unfortunately, the calculation of the coupling constant Ξ is not obvious, since we do not have a homogeneous and isotropic permittivity $\epsilon_r$.

### B. Validity of the mPB model

To estimate when ion/ion couplings may appear locally, we have computed the plasma parameter $\Gamma$, which compares the electrostatic interaction energy between two ions and the thermal energy. It is given by$^{33}$:

$$ \Gamma(z) = \frac{\ell_B(z)}{d_{\text{ion}}(z)}, \quad (17) $$

where $\ell_B(z)$ is the local Bjerrum length defined by Eq. 6 using the local $\epsilon_\perp(z)$, and $d_{\text{ion}}(z)$ is the typical distance between two ions and can be computed as follows: $d_{\text{ion}}(z) = \sqrt[3]{C_+(z) + C_-(z)}$. The value $\Gamma = 1$ is the threshold at which correlations between the ions start to become non negligible$^{33,47}$. Since these are not considered in the PB model, its accuracy decreases when $\Gamma$ exceeds this value.

Fig. 4 shows the evolution of the plasma parameter in systems at T = 353.15 K, at different HN. At high hydration, the plasma parameter exceeds 1 in the lipid heads, but quickly decays to a lower value in the water. At HN = 14 and below, $\Gamma$ stays over 1 in the water. Because of the correlations, the results given by the mPB model in this hydration range will be less accurate, though they can still be analysed qualitatively. In addition, the plasma parameter could be underestimated since it is computed using the dielectric permittivity given by the mPB model that neglects correlations and some anisotropy effects.

### C. Comparison with alternative ways of estimating the relative dielectric permittivity

In this section, we compare the results obtained from the mPB model described in this work with alternative ways of estimating the dielectric permittivity.

**Born dielectric permittivity**— By getting rid of the electrostatic energy term the same way we did in Eq. (13),

<figure><img src="image_5.png" /><figcaption>FIG. 4. Evolution of the plasma parameter at different HN. Above the dashed line at $\Gamma = 1$, ion correlations start to become non negligible.</figcaption></figure>
i.e. by multiplying the ionic concentrations, one can compute a dielectric permittivity profile $\epsilon_{\perp}^{\text{Born}}$ which would be the effective permittivity in the Born model. This "Born dielectric permittivity" derives from Eq. (10) and is given by:

$$ \frac{1}{\epsilon_{\perp}^{\text{Born}}(z)} = \frac{1}{\epsilon_{\perp}^{\text{ref}}} + \frac{16\pi\epsilon_0}{\beta e^2} \frac{r_+r_-}{r_+ + r_-} \times \ln\left(\frac{C_{\text{salt}}}{\sqrt{C_+(z) \times C_-(z)}}\right). \quad (18) $$

Since the mPB model presented in this work highly relies on the Born model to compute the effective permittivity profile, comparing our results with the profile of $\epsilon_{\perp}^{\text{Born}}$ can give an insight on the self-consistency of the mPB model. Fig. 5 shows a comparison between $\epsilon_{\perp}^{\text{mPB}}$ and $\epsilon_{\perp}^{\text{Born}}$. Eq. (18) is not valid when the ionic concentrations are null, which is why $\epsilon_{\perp}^{\text{Born}}$ has no value in the tails region. Elsewhere however, it appears that the two profiles are quite similar, and even more at low HN. Noticeably, in the mPB model, we have neglected the change in steric interactions between the water and lipid region, i.e. the variation of the free energy to form a cavity for the ions in the different locations of the system. The consistency between the profiles of $\epsilon_{\perp}^{\text{Born}}(z)$ and $\epsilon_{\perp}^{\text{mPB}}(z)$ illustrated in Fig. 5 supports this assumption.

*Fluctuation-dissipation Approach*— The dielectric permittivity is commonly calculated using the fluctuation-dissipation theorem that links the fluctuations of the total dipole moment $M$ to the dielectric permittivity$^{48}$. In

<figure><img src="image_6.png" /><figcaption>FIG. 5. Comparison between $\epsilon_{\perp}^{\text{mPB}}$, $\epsilon_{\perp}^{\text{Born}}$ and $\epsilon_{\perp}^{\text{water density}}$ for two systems at T = 353.15 K. (a) HN = 35. (b) HN = 10.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>FIG. 6. Dielectric response of the water molecules versus hydration number at different temperatures. Black dots correspond to $\epsilon_{\perp}^{\text{mPB}}(z_{\text{max}})$ at 353.15 K. Filled markers represent bilayers in the fluid phase and open markers represent bilayers in the gel phase. The purple region shows where we identify an attractive regime between the membranes (see Sect. IV D).</figcaption></figure>

a homogeneous medium containing solely water dipoles, the relation is:

$$ \epsilon^{\text{dip. fluct.}} = 1 + \frac{\langle M^2 \rangle - \langle M \rangle^2}{3\epsilon_0 n_w V_w k_B T}, \quad (19) $$

where $n_w$ is the number of water molecules, and $V_w = 29.7 \AA^3$ is the volume of a water molecule. Here we propose to estimate the dielectric permittivity through this formula, beyond its validity domain, since the system is inhomogeneous, and contains other dipoles than the water molecules. In practice, this approach should be considered as a way to assess the total water dipole fluctuations in a normalized representation.

Results of this method are shown in Fig. 6. The plot of $\epsilon^{\text{dip. fluct.}}$ as a function of the interlayer water thickness is shown in Fig. S11 in the supplementary material. At high hydration (HN ≥ 35), both mPB model and this method give similar results. However, they differ more and more upon dehydration, because of two main reasons. The first one is that $\epsilon^{\text{dip. fluct.}}$ only takes into account the water molecules, while $\epsilon_{\perp}^{\text{mPB}}$ acts as an effective permittivity in the system, and does not strictly corresponds to the dielectric permittivity of the sole water molecules. The second reason is that by essence, $\epsilon^{\text{dip. fluct.}}$ is monotonically decreasing upon dehydration, as is linked to the decrease of the rotational degree of freedom of water molecules. At very low HN, the only remaining water molecules are surrounded by charged components (being either the ions or the lipid heads). The dipole moment of the water molecules is then constrained by these interactions and cannot really fluctuate, explaining the extremely low values of the permittivity obtained with this method in such regime.

*Correlation between permittivity and water density*—
Finally, the dielectric permittivity profile shapes shown in Figs. 2(a) and 2(d) suggest that the permittivity profile follows the density profiles of water, $\rho_w(z)$. Such an ad-hoc scaling was already proposed by M. Miettinen$^{49}$ for cationic bilayers$^{50,51}$. We model here this correlation with a function that interpolates between the two extreme values $\epsilon_{tails}$ to $\epsilon_{\perp}^{\text{ref}}$ with a variation proportional to $\rho_w/\rho_w^{\text{bulk}}$:

$$\epsilon_{\perp}^{\text{density}}(z) = \epsilon_{\text{tails}} + \frac{\rho_w(z)}{\rho_w^{\text{bulk}}} [\epsilon_{\perp}^{\text{ref}} - \epsilon_{\text{tails}}]. \quad (20)$$

The values for $\epsilon_{\perp}^{\text{ref}}$ are reported in Tab. S2. The $\rho_w^{\text{bulk}}$ were chosen as the water densities at $z_{max}$ for HN = 45.

The comparison of $\epsilon_{\perp}^{\text{density}}$ profile with the one of $\epsilon_{\perp}^{\text{mPB}}$ for two different HN at $T = 353.15$ K are reported in Fig. 5. It turns out that the profiles given by Eq. 20 are qualitatively similar. A first difference is that $\epsilon_{\perp}^{\text{density}}(z_{max})$ cannot describe the increase of $\epsilon_{\perp}^{\text{mPB}}(z_{max})$ upon dehydration since its maximum remains $\epsilon_{\perp}^{\text{ref}}$ whatever the hydration level. Moreover, Eq. 20 underestimates both $\epsilon_{\perp}^{\text{mPB}}$ and $\epsilon_{\perp}^{\text{Born}}$. Especially in the region where the water density is very low, the contribution of the lipid dipoles to the permittivity $\epsilon_{\perp}^{\text{mPB}}$ is not negligible.

### D. Resulting force between the bilayers

In the case of charged lipid bilayers with monovalent salt, the force between the bilayers is generally expected to be repulsive$^{52}$, but recent observations of DPPS membranes with NaCl salt suggest that this interaction could also be attractive$^{9,10}$. This attractive behaviour was associated to a decrease of permittivity in the water layer.

In this context, we also aimed at determining if the membranes in our simulations were attracting or repelling each others. The force between lipid bilayers can be extracted from MD simulations, for example by calculating the water chemical potential$^{53,54}$. It would be very interesting to perform such calculations in the present case, but it is beyond the scope of the present article. Given the relatively high computational cost, we resorted to a different approach, valid for membrane simulations without surface tension$^{55}$. It involves looking at the Helmholtz free energy $F$, which depends on the volume $V$, the temperature $T$ and the total number of particles $N$ in the system. Smirnova and al. expressed, in the (NVT) ensemble, the excess of free energy $\Delta F$, with respect to a single fully hydrated (in bulk water) lipid bilayer with no surface tension, as follows:

$$\Delta F = \Delta F_{\text{bilayer}} + \Delta F_{\text{water}} + g(d_w)A, \quad (21)$$

where $\Delta F_{\text{bilayer}}$ is the excess free energy contribution coming from the elastic deformation of the bilayer, $\Delta F_{\text{water}}$ is the one coming from the compression of the bulk of water surrounding the membrane, $A$ is the surface

of the interface and $g$ is the interface potential. In stacked bilayers systems, the resulting force between lipid bilayers is generally considered as the sum of three main contributions. The first one comes from the van der Waals interactions, which are attractive. The second one, the hydration repulsion, prevents the membranes from sticking to each others: it is a short-ranged interaction and decays exponentially with the thickness of the interlayer confined water$^{56}$. The third one comes from the electrostatic interactions. We suppose here that these interactions between membranes can be described within the last term $g(d_w)A$, i.e. proportional to the membrane area, and as a function of the water thickness $d_w$. This assumption turned out to be relevant for the neutral lipid membranes$^{55}$, and we assume it to hold for charged membranes as well. The potential $g$ depends on the thickness of the water thickness $d_w$, which can be estimated from the HN as:

$$d_w(\text{HN}) = \frac{2 \times \text{HN} \times V_w}{A_{\text{lipid}}}, \quad (22)$$

where $V_w$ is the water molecular volume and $A_{\text{lipid}}$ is the area per lipid. We did not consider the dependence of this molecular volume on the temperature here as we were not interested in quantitative measurements but rather in the qualitative behaviour in the simulation.

Smirnova and al. stated that since the variation of the simulation box volume is negligible, the variations of Gibbs free enthalpy $G$ are similar to the Helmholtz free energy variations: $\Delta G \approx \Delta F$. They were then able to express the repulsion pressure $P(d_w) = -g'(d_w)$, the derivative of the interface potential $g$ can be expressed as a function of the area variation $\Delta A = A - A_0$:

$$P(d_w) = -\frac{K_A}{A_0} \left[ \frac{\Delta A(d_w)}{d_w} + \int_{d_w}^{\infty} \frac{\Delta A(d)}{d^2} dd \right], \quad (23)$$

where $K_A$ is the area compressibility modulus, $A_0$ is the reference area at the largest HN, and $A(d_w)$ is the area for a given water thickness, i.e. a given HN. The attractive vs repulsive behaviour between the membranes depends on the sign of this pressure. It will be repulsive if $P(d_w) > 0$ and attractive otherwise. Noticeably, in the work by Smirnova et al., Eq. 23 was used under several conditions$^{55}$: (1) there is no fluid/gel phase change in the bilayer. (2) $K_A$ does not vary upon dehydration. (3) the integral in Eq. 23 is supposed to be negligible and the sign of the area variation $\Delta A(d_w)$ directly determines the sign of the pressure $P(d_w)$. Since $K_A$, $d_w$ and $A_0$ are positive, a larger area denotes an attractive behaviour while a smaller area stands for a repulsion.

To constrain our analysis within the condition (1), we shall solely discuss here simulations in the fluid phase. Concerning the condition (2), Fig. S7 shows that for the simulations in the fluid phase, the compressibility modulus is practically constant in the range $11 \le \text{HN} \le 45$. Concerning the condition (3), Fig. S8 comparing the two terms of Eq. 23 shows that the second term is negligible
<figure><img src="image_8.png" /><figcaption>FIG. 7. Area per lipid versus hydration number at different temperatures. Filled markers represent bilayers in the fluid phase and open markers represent bilayers in the gel phase. The black dotted line shows the reference area per lipid $A_0$ at $T = 353.15$ K. At this temperature, following Eq. (23), we identify the purple region as the hydration range at which the attractive behaviour appears.</figcaption></figure>

at larger HN, but it can become non negligible at the lowest HNs. Nevertheless, in the range of hydration where $K_A$ is constant, at larger distances $d$, $\Delta A(d)$ does not change sign. In this hydration range, the two terms of Eq. 23 are of the same sign. Therefore, we analyse in the following the sign of $\Delta A(d_w)$ to infer the sign of $P(d_w)$.

Fig. 7 shows the area per lipid as a function of the HN, at the five studied temperatures. The same plot as a function of the interlayer water thickness $d_w$ is shown in Fig. S9 in the Supplementary Material. Since the number of lipids per layer remains constant, it is equivalent to study the total area variation or the area per lipid variation. It appears that at high HN (typically above 20), there is no significant change in the area per lipid upon dehydration. This indicates that in this HN range the bilayers are fully hydrated and they do not interact significantly. However, if we keep dehydrating the system, two thresholds appear. The first one is between HN = 20 and HN = 15. After it is reached, we see an increase in the area per lipid. The second threshold appears around HN = 10-11, below which the opposite phenomenon occurs: the area per lipid starts to decrease. These changes in the area per lipid can even induce phase transitions in the bilayers. In the coldest systems between HN = 20 and HN = 15, the bilayer undergoes a transition from the gel phase to the fluid phase. At lower HN, it undergoes a transition from the fluid phase to the gel phase in all the simulations.

One has to take care of these phase changes when interpreting the differences in area per lipid using Eq. (23). In this work, we safely focus on the system at 353.15 K, for which the bilayers stay in a fluid phase down to HN = 6. At high hydration, the system is in an equilibrium state and there is no significant force acting on the bilayers.

However, after the first threshold is met, the area per lipid is higher than the reference value. This corresponds to a negative value for the repulsion pressure, i.e. an attractive interaction between the bilayers. If the dehydration continues after the second threshold, the conditions necessary for Eq. 23 do not hold any longer in our simulations. But the strong decrease in area might be the hint that the short-range hydration repulsion well known in the literature for any kind of lipid bilayers⁵⁶ becomes more important.

To conclude, even if we do not calculate explicitly the interface potential, the area per lipid variations indicate some inter-membrane attraction at $T = 353.25$ K in the range of $11 \le HN \le 16$, which is the same HN-range where the permittivity of the water layer drops.

V. CONCLUSION

We have presented a 1D modified Poisson-Boltzmann model, which allows the dielectric permittivity to vary along the $z$ axis and takes into account the solvation energy of the ions (Born energy) in their Boltzmann distribution.

In the model's equation, the two free charge densities come from the ions and from the lipids. The first one has an analytical description (a Boltzmann distribution) while the second one is added as an external fixed contribution obtained from MD simulations. Two unknowns remain in the equation: the electric potential and the dielectric permittivity. We left the electric potential as the unknown of the equation and we made an educated guess on the shape of the dielectric permittivity profile. With the help of a previous work that could only compute this profile in systems without free charges⁵, we modelled the dielectric permittivity profile as a parametrized curve and optimized the parameters so that the analytical description of the ions in the model fitted at best the charge densities obtained from MD simulations.

Since this model can estimate the relative dielectric permittivity in a system containing free charges, we were able to use it to investigate water confined between charged lipid bilayers. The results show a self-consistence of the model, giving a dielectric permittivity profile close from the one that directly derives from the Born model. Moreover, when focusing on the dielectric permittivity at the centre of the interlayer channel, the mPB model shows a strong decrease of the dielectric constant when entering a hydration regime where the membranes attract each others. This observation is in agreement with the common interpretation of the phenomenon in the literature. At intermediate hydration levels, the mPB model results in a local maximum of the dielectric constant that is higher than the bulk value. This may be due to the fact that our profiles are by construction smoothed, coarse-grained, and treat the dielectric permittivity as a local property. At high hydration, it is in good agree-
ment with a classic way of computing the dielectric constant from the fluctuations of the water molecules dipole moment. However upon dehydration, it diverges rather quickly from this method, which only considers the water molecules.

In the end, the model presented in this work is generic and thus can be applied in a large variety of systems where an electrolyte solution is confined between charged walls. We expect the mPB model to be more adapted to the simpler case of homogeneous planar walls, as the thermal fluctuations of the membranes are not taken into account in the model, since the lipid charge density profile is averaged over the simulation length. The curve to adjust the inverse permittivity profile would have to be adapted to these new systems. This approach could also be useful to estimate the dielectric permittivity profile in biological systems, where the biological fluid of interest is not too much confined so the distributions of the salt it contains can be modelled with a Boltzmann distribution.

## SUPPLEMENTARY MATERIAL

The supplementary material includes (I) details on the numerical solving of the mPB model, (II) the inverse dielectric permittivity profile for a DPPC bilayer, (III) inverse dielectric permittivity profiles for the DPPS bilayers, and their variation when the fitting function includes 0, 1 or 2 Gaussians, (IV) details on the calibration of the mPB model and the independent permittivity measurements on bulk water, (V) a comparison of Born and electrostatic energy profiles, (VI) the area compressibility modulus as a function of hydration level at different temperatures, (VII) the dependence of the area per lipid, and relevant permittivities as a function of $d_w$, and (VII) the list of simulated systems and their compositions.

## ACKNOWLEDGMENTS

We thank Hélène Berthoumieux, Ananya Debnath, Marcus Müller and Alexander Schlaich for fruitful discussions. This study was supported by the French Agence Nationale de la Recherche through the Banana_Slip project (ANR-21-CE30-0026). We used computer and storage resources by GENCI at IDRIS thanks to the grants A0110807662 and A0150807662 on the supercomputer Jean Zay's V100 partition. The authors thank the reviewers for their valuable feedback and insightful suggestions.

## AUTHOR DECLARATIONS

*Conflict of interest*— The authors have no conflict to disclose.

## DATA AVAILABILITY

The data that support the findings of this study are available from the corresponding author upon reasonable request. Some complete trajectories, some input and output files have been published on Zenodo repository⁵⁷, to permit a reference in MDverse⁵⁸ and NMRLipids⁵⁹ catalogs. Other trajectories are available upon reasonable request to the authors.

## REFERENCES

¹R. M. Forbes, A. R. Cooper, and H. H. Mitchell. The composition of the adult human body as determined by chemical analysis. *The Journal of biological chemistry*, 203:359–366, 1953.

²H. Binder. Water near lipid membranes as seen by infrared spectroscopy. *European Biophysics Journal*, 36:265–279, 2007.

³N. E. Levinger. Water in confinement. *Science*, 298:1722–1723, 2002.

⁴L. Fumagalli, A. Esfandiar, R. Fabregas, S. Hu, P. Ares, A. Jardananan, Q. Yang, B. Radha, T. Taniguchi, K. Watanabe, G. Gomila, K. S. Novoselov, and A. K. Geim. Anomalously low dielectric constant of confined water. *Science*, 360:1339–1342, 2018.

⁵A. Schlaich, E. W. Knapp, and R. R. Netz. Water dielectric effects in planar confinement. *Physical Review Letters*, 117, 2016.

⁶Jean-François Olivierier, James T. Hynes, and Damien Laage. Confined waters dielectric constant reduction is due to the surrounding low dielectric media and not to interfacial molecular ordering. *The Journal of Physical Chemistry Letters*, 12(17):4319–4326, 2021.

⁷D. Borgis, D. Laage, L. Belloni, and G. Jeanmairet. Dielectric response of confined water films from a classical density functional theory perspective. *Chemical Science*, 14:11141–11150, 2023.

⁸A. Schlaich, A. P. Dos Santos, and R. R. Netz. Simulations of nanoseparated charged surfaces reveal charge-induced water reorientation and nonadditivity of hydration and mean-field electrostatic repulsion. *Langmuir*, 35:551–560, 2019.

⁹T. Mukhina, A. Hemmerle, V. Rondelli, Y. Gerelli, G. Fragneto, J. Daillant, and T. Charitat. Attractive interaction between fully charged lipid bilayers in a strongly confined geometry. *Journal of Physical Chemistry Letters*, 10:7195–7199, 2019.

¹⁰Sven Helstroffer, Ludovic Gardré, Giovanna Fragneto, Arnaud Hemmerle, Léo Henry, Laurent Joly, Fabrice Thalmann, Claire Loison, Pierre Muller, and Thierry Charitat. The role of confined water in the emergence of electrostatic strong coupling as revealed by nanoseparated charged lipid layers. *The Journal of Physical Chemistry Letters*, 16(33):8369–8376, 2025.

¹¹G. Gramse, A. Dols-Perez, M. A. Edwards, L. Fumagalli, and G. Gomila. Nanoscale measurement of the dielectric constant of supported lipid bilayers in aqueous solutions with electrostatic force microscopy. *Biophysical Journal*, 104:1257–1262, 2013.

¹²R. Wang, M. Souilamas, A. Esfandiar, R. Fabregas, S. Benaglia, H. Nevison-Andrews, Q. Yang, J. Normansell, P. Ares, G. Ferrari, A. Principi, A. K. Geim, and L. Fumagalli. In-plane dielectric constant and conductivity of confined water. 7 2024.

¹³Harry A. Stern and Scott E. Feller. Calculation of the dielectric permittivity profile for a nonuniform system: Application to a lipid bilayer simulation. *The Journal of Chemical Physics*, 118(7):3401–3412, 02 2003.

¹⁴V. Ballenegger and J. P. Hansen. Dielectric permittivity profiles of confined polar fluids. *Journal of Chemical Physics*, 122, 2005.

¹⁵Christian Schaaf and Stephan Gekle. Spatially resolved dielectric constant of confined water and its connection to the nonlocal nature of bulk water. *The Journal of Chemical Physics*, 145(8):084901, 08 2016.
16 Douwe Jan Bonthuis, Stephan Gekle, and Roland R. Netz. Dielectric profile of interfacial water and its effect on double-layer capacitance. *Phys. Rev. Lett.*, 107:166102, Oct 2011.

17 D. M. Huang, C. Cottin-Bizonne, C. Ybert, and L. Bocquet. Aqueous electrolytes near hydrophobic surfaces: Dynamic effects of ion specificity and hydrodynamic slip. *Langmuir*, 24:1442–1450, 2008.

18 A. Poddar, D. Maity, A. Bandopadhyay, and S. Chakraborty. Electrokinetics in polyelectrolyte grafted nanofluidic channels modulated by the ion partitioning effect. *Soft Matter*, 12:5968–5978, 2016.

19 S. Jo, T. Kim, V. G. Iyer, and W. Im. Charmm-gui: A web-based graphical user interface for charmm. *Journal of Computational Chemistry*, 29:1859–1865, 2008.

20 J. B. Klauda, R. M. Venable, J. A. Freites, J. W. O'Connor, D. J. Tobias, C. Mondragon-Ramirez, I. Vorobyov, A. D. MacKerell, and R. W. Pastor. Update of the charmm all-atom additive force field for lipids: Validation on six lipid types. *Journal of Physical Chemistry B*, 114:7830–7843, 2010.

21 A. D. Jr. MacKerell, D. Bashford, M. Bellott, R. L. Jr. Dunbrack, J. D. Evansck, M. J. Field, S. Fischer, J. Gao, H. Guo, S. Ha, D. Joseph-McCarthy, L. Kuchnir, K. Kuczera, F. T. K. Lau, C. Mattos, S. Michnick, T. Ngo, D. T. Nguyen, B. Prodhom, W. E. Reiher, B. Roux, M. Schlenkrich, J. C. Smith, R. Stote, J. Straub, M. Watanabe, J. Wiórkiewicz-Kuczera, D. Yin, and M. Karplus. All-atom empirical potential for molecular modeling and dynamics studies of proteins. *The Journal of Physical Chemistry B*, 102(18):3586–3616, 1998.

22 S. Izadi, R. Anandakrishnan, and A. V. Onufriev. Building water models: A different approach. *Journal of Physical Chemistry Letters*, 5:3863–3871, 2014.

23 K. Vanommeslaeghe and A. D. Mackerell Jr. Charmm additive and polarizable force fields for biophysics and computer-aided drug design. *Biochimica et Biophysica Acta - General Subjects*, 1850:861–871, 2015.

24 E. Duboué-Dijon, M. Javanainen, P. Delcroix, P. Jungwirth, and H. Martinez-Seara. A practical guide to biologically relevant molecular simulations with charge scaling for electronic polarization. *Journal of Chemical Physics*, 153, 2020.

25 J. Melcr, T. M. Ferreira, P. Jungwirth, and O. H. S. Ollila. Improved cation binding to lipid bilayers with negatively charged pops by effective inclusion of electronic polarization. *Journal of Chemical Theory and Computation*, 16:738–748, 2020.

26 Milla Kurki, Alexey M. Nesterenko, Nicolai E. Alsaker, Tiago M. Ferreira, Sami Kyllönen, Antti Poso, Piia Bartos, and Markus S. Miettinen. Solid-state nmr validation of opls4: Structure of pc-lipid bilayers and its modulation by dehydration. *The Journal of Physical Chemistry B*, 128(50):12483–12492, 2024.

27 Ulrich Essmann, Lalith Perera, Max L. Berkowitz, Tom Darden, Hsing Lee, and Lee G. Pedersen. A smooth particle mesh ewald method. *The Journal of Chemical Physics*, 103(19):8577–8593, 11 1995.

28 Tom Darden, Darrin York, and Lee Pedersen. Particle mesh ewald: An nlog(n) method for ewald sums in large systems. *The Journal of Chemical Physics*, 98(12):10089–10092, 06 1993.

29 M. J. Abraham, T. Murtola, R. Schulz, S. Páll, J. C. Smith, B. Hess, and E. Lindahl. Gromacs: High performance molecular simulations through multi-level parallelism from laptops to supercomputers. *SoftwareX*, 1-2:19–25, 2015.

30 Glenn J. Martyna, Michael L. Klein, and Mark Tuckerman. Noséhoover chains: The canonical ensemble via continuous dynamics. *The Journal of Chemical Physics*, 97(4):2635–2643, 08 1992.

31 M. Parrinello and A. Rahman. Polymorphic transitions in single crystals: A new molecular dynamics method. *Journal of Applied Physics*, 52(12):7182–7190, 12 1981.

32 J. N. Israelachvili. *Intermolecular and Surface Forces (Third Edition)*. Academic Press, 2011.

33 C. Herrero and L. Joly. The poisson-boltzmann equation in micro- and nanofluidics: A formulary. *Physics of Fluids*, 36, 2024.

34 Philip Loche, Cihan Ayaz, Emanuel Wolde-Kidan, Alexander Schlaich, and Roland R. Netz. Universal and nonuniversal aspects of electrostatics in aqueous nanoconfinement. *The Journal of Physical Chemistry B*, 124(21):4365–4371, 2020. PMID: 32364728.

35 Hugh Nymeyer and Huan Xiang Zhou. A method to determine dielectric constants in nonhomogeneous systems: Application to biological membranes. *Biophysical Journal*, 94:1185–1193, 2008.

36 Wei-Tze Huang and David G. Levitt. Theoretical calculation of the dielectric constant of a bilayer membrane. *Biophysical Journal*, 17:111–128, 1977.

37 Alexander A. Rashin and Barry Honig. Reevaluation of the born model of ion hydration. *The Journal of Physical Chemistry*, 89:5588–5593, 1985.

38 Spiros Kournopoulos, Mirella Simões Santos, Srikanth Ravipati, Andrew J. Haslam, George Jackson, Ioannis G. Economou, and Amparo Galindo. The contribution of the ion-ion and ion-solvent interactions in a molecular thermodynamic treatment of electrolyte solutions. *Journal of Physical Chemistry B*, 126:9821–9839, 2022.

39 Gabriel M. Silva, Bjørn Maribo-Mogensen, Xiaodong Liang, and Georgios M. Kontogeorgis. Improving the born equation: Origin of the born radius and introducing dielectric saturation effects. *Fluid Phase Equilibria*, 576, 2024.

40 R. Schmid, A. M. Miah, and V. N. Sapunov. A new table of the thermodynamic quantities of ionic hydration: values and some applications (enthalpyentropy compensation and born radii). *Physical Chemistry Chemical Physics*, 2:97–102, 2000.

41 Chunyi Zhang, Shuwen Yue, Athanassios Z. Panagiotopoulos, Michael L. Klein, and Xifan Wu. Why dissolving salt in water decreases its dielectric permittivity. *Phys. Rev. Lett.*, 131:076801, 08 2023.

42 K. Gawrisch, D. Ruston, J. Zimmerberg, V. A. Parsegian, R. P. Rand, and N. Fuller. Membrane dipole potentials, hydration forces, and the ordering of water at membrane surfaces. *Biophysical Journal*, 61:1213–1223, 1992.

43 S. E. Feller and A. D. MacKerell. An improved empirical potential energy function for molecular simulations of phospholipids. *Journal of Physical Chemistry B*, 104:7510–7515, 2000.

44 S. Malik and A. Debnath. Dehydration induced dynamical heterogeneity and ordering mechanism of lipid bilayers. *Journal of Chemical Physics*, 154, 2021.

45 A. G. Moreira and R. R. Netz. Strong-coupling theory for counter-ion distributions. *Europhysics Letters*, 52:705–711, 2000.

46 K. Komorowski, A. Salditt, Y. Xu, H. Yavuz, M. Brennich, R. Jahn, and T. Salditt. Vesicle adhesion and fusion studied by small-angle x-ray scattering. *Biophysical Journal*, 114:1908–1920, 2018.

47 Y. Levin. Electrostatic correlations: from plasma to biology. *Reports on Progress in Physics*, 65:1577–1632, 2002.

48 J. G. Kirkwood. The dielectric polarization of polar liquids. *The Journal of Chemical Physics*, 7:911–919, 1939.

49 Markus Miettinen. *Computational Modeling of Cationic Lipid Bilayers in Saline Solutions*. Aalto University, 2010.

50 Andrey A. Gurtovenko, Markus Miettinen, Mikko Karttunen, and Ilpo Vattulainen. Effect of monovalent salt on cationic lipid membranes as revealed by molecular dynamics simulations. *The Journal of Physical Chemistry B*, 109(44):21126–21134, 2005. PMID: 16853736.

51 Markus S. Miettinen, Andrey A. Gurtovenko, Ilpo Vattulainen, and Mikko Karttunen. Ion dynamics in cationic lipid bilayer systems in saline solutions. *The Journal of Physical Chemistry B*, 113(27):9226–9234, 2009. PMID: 19534449.

52 J. Marra. Direct measurement of the interaction between phosphatidylglycerol bilayers in aqueous electrolyte solutions. *Biophysical Journal*, 50:815–825, 1986.

53 E. Schneck, F. Sedlmeier, and R. R. Netz. Hydration repulsion between biomembranes results from an interplay of dehydration and depolarization. *Proceedings of the National Academy of Sci-
<figure><img src="image_9.png" /></figure>

ences of the United States of America, 109:14405–14409, 2012.

<footnote n="54">54 A. Schlaich, J. O. Daldrop, B. Kowalik, M. Kandu, E. Schneck, and R. R. Netz. Water structuring induces nonuniversal hydration repulsion between polar surfaces: Quantitative comparison between molecular simulations, theory, and experiments. *Langmuir*, 40:7896–7906, 2024.</footnote>

<footnote n="55">55 Y. G. Smirnova, S. Aeffner, H. J. Risselada, T. Salditt, S. J. Marrink, M. Müller, and V. Knecht. Interbilayer repulsion forces between tension-free lipid bilayers from simulation. *Soft Matter*, 9:10705–10718, 2013.</footnote>

<footnote n="56">56 W. Xu and F. Pincet. Quantification of phase transitions of lipid mixtures from bilayer to non-bilayer structures: Model, experimental validation and implication on membrane fusion. *Chemistry and Physics of Lipids*, 163:280–285, 2010.</footnote>

<footnote n="57">57 Ludovic Gardré, 10.5281/zenodo.16318887 2025.</footnote>

<footnote n="58">58 Johanna KS Tiemann, Magdalena Szczuka, Lisa Bouarroudj, Mohamed Oussaren, Steven Garcia, Rebecca J Howard, Lucie Delemotte, Erik Lindahl, Marc Baaden, Kresten Lindorff-Larsen, Matthieu Chavent, and Pierre Poulain. Mdverse, shedding light on the dark matter of molecular dynamics simulations. *eLife*, 12:RP90061, 08 2024.</footnote>

<footnote n="59">59 Anne M. Kiirikki, Hanne S. Antila, Lara S. Bort, Pavel Buslaev, Fernando Favela-Rosales, Tiago Mendes Ferreira, Patrick F.J. Fuchs, Rebeca Garcia-Fandino, Ivan Gushchin, Batuhan Kav, Norbert Kuerka, Patrik Kula, Milla Kurki, Alexander Kuzmin, Anusha Lalitha, Fabio Lolicato, Jesper J. Madsen, Markus S. Miettinen, Cedric Mingham, Luca Monticelli, Ricky Nencini, Alexey M. Nesterenko, Thomas J. Piggot, Ángel Piñeiro, Nathalie Reuter, Suman Samantray, Fabián Suárez-Lestón, Reza Talandashti, and O. H.Samuli Ollila. Overlay databank unlocks data-driven analyses of biomolecules for all. *Nature Communications*, 15:1–15, 2024.</footnote>